The onset of thermo-compositional convection in rotating spherical shells

TL;DR

This paper investigates the linear onset of thermo-compositional convection in rotating spherical shells, revealing how differences in thermal and compositional diffusivities influence the transition between various convective regimes.

Contribution

It introduces an effective Rayleigh number and a mixing parameter to analyze the linear stability of double-diffusive convection, providing new insights into the dynamics across different regimes.

Findings

Identification of critical curves in the Ra-α space

Analysis of asymptotic behaviors and regime transitions

Mapping of convection modes across diffusivities and rotation rates

Abstract

Double-diffusive convection driven by both thermal and compositional buoyancy in a rotating spherical shell can exhibit a rather large number of behaviours often distinct from that of the single diffusive system. In order to understand how the differences in thermal and compositional molecular diffusivities determine the dynamics of thermo-compositional convection we investigate numerically the linear onset of convective instability in a double-diffusive setup. We construct an alternative equivalent formulation of the non-dimensional equations where the linearised double-diffusive problem is described by an effective Rayleigh number, $\text{Ra}$, measuring the amplitude of the combined buoyancy driving, and a second parameter, $\alpha$, measuring the mixing of the thermal and compositional contributions. This formulation is useful in that it allows for the analysis of several limiting cases and reveals dynamical similarities in the parameters space which are not obvious otherwise. We analyse the structure of the critical curves in this $\text{Ra}-\alpha$ space, explaining asymptotic behaviours in $\alpha$, transitions between inertial and diffusive regimes, and transitions between large scale (fast drift) and small scale (slow drift) convection. We perform this analysis for a variety of diffusivities, rotation rates and shell aspect ratios showing where and when new modes of convection take place.

Tables

Table 1. Table 1: Non-dimensional model parameters. In the text we occasionally refer to both the Prandtl and the Schmidt numbers as “Prandtl numbers” for brevity.

Parameter	Definition
Coriolis Number	$\tau =2\Omega d^2/\nu$
Thermal Rayleigh Number	$R_t={\alpha }_T\gamma d^4{T}^{\ast }/{\nu }^2$
Compositional Rayleigh Number	$R_c={\alpha }_C\gamma d^4{C}^{\ast }/{\nu }^2$
Thermal Prandtl Number	$\mathrm{Pr}=\nu /\kappa$
Compositional Prandtl Number (Schmidt Number)	$\mathrm{Sc}=\nu /D$
Radius Ratio	$\eta =r_i/r_o$

Equations

$${\text{upartial}}_t\mathbf{u}=-\tau \hat{\mathbf{k}}\times \mathbf{u}-\nabla \pi +({\mathrm{R}}_{\mathrm{t}}\Theta +{\mathrm{R}}_{\mathrm{c}}\chi )\mathbf{r}+{\nabla }^2\mathbf{u},$$

$$\nabla \cdot \mathbf{u}=0,$$

$${\text{upartial}}_t\Theta ={\mathrm{Pr}}^{-1}{\nabla }^2\Theta -\mathbf{u}\cdot \nabla T,$$

$${\text{upartial}}_t\chi ={\mathrm{Sc}}^{-1}{\nabla }^2\chi -\mathbf{u}\cdot \nabla C,$$

$$\mathbf{u}={\mathbf{u}}_P+{\mathbf{u}}_T=\nabla \times \nabla \times \mathcal{S}(\mathbf{r},t)\mathbf{r}+\nabla \times \mathcal{T}(\mathbf{r},t)\mathbf{r},$$

$$\mathcal{X}(\mathbf{r},t)=\tilde{\mathcal{X}}(\mathbf{r})\exp (\mathrm{i}t(\omega -\mathrm{i}\Gamma )),$$

$$(\mathrm{i}\omega +\Gamma )({\nabla }^2-2r{\text{upartial}}_r)({\nabla }_H^2\tilde{\mathcal{S}})$$

$$=-\tau \mathbf{r}\cdot \nabla \times \nabla \times (\hat{\mathbf{k}}\times \tilde{\mathbf{u}})+{\nabla }_H^2({\mathrm{R}}_{\mathrm{t}}\tilde{\Theta }+{\mathrm{R}}_{\mathrm{c}}\tilde{\chi })+\mathbf{r}\cdot \nabla \times (\nabla \times ({\nabla }^2\tilde{\mathbf{u}})),$$

$$(\mathrm{i}\omega +\Gamma )({\nabla }_H^2\tilde{\mathcal{T}})=-\tau \mathbf{r}\cdot \nabla \times (\hat{\mathbf{k}}\times \tilde{\mathbf{u}})+\mathbf{r}\cdot \nabla \times ({\nabla }^2\tilde{\mathbf{u}}),$$

$$(\mathrm{i}\omega +\Gamma )\tilde{\Theta }={\mathrm{Pr}}^{-1}{\nabla }^2\tilde{\Theta }-\tilde{\mathbf{u}}\cdot \nabla T,$$

$$(\mathrm{i}\omega +\Gamma )\tilde{\chi }={\mathrm{Sc}}^{-1}{\nabla }^2\tilde{\chi }-\tilde{\mathbf{u}}\cdot \nabla C,$$

$${\nabla }_H^{2}f={\nabla }^{2}f-\frac{1}{r^2}\frac{\text{upartial}}{\text{upartial}r}\left(r^2\frac{\text{upartial}f}{\text{upartial}r}\right),$$

$$\tilde{\mathcal{S}}={\text{upartial}}_{rr}\tilde{\mathcal{S}}={\text{upartial}}_r(\tilde{\mathcal{T}}/r)=\tilde{\Theta }=\tilde{\chi }=0\text{at}r=r_i,r_o.$$

$$(\mathrm{i}\omega +\Gamma ){\left[A\right]}_{n,l}\left[\begin{array}{c}{\tilde{s}}_{n,l,m}\\ {\tilde{t}}_{n,l,m}\\ {\tilde{\theta }}_{n,l,m}\\ {\tilde{\chi }}_{n,l,m}\end{array}\right]={\left[B\right]}_{n,l}\left[\begin{array}{c}{\tilde{s}}_{n,l,m}\\ {\tilde{t}}_{n,l,m}\\ {\tilde{\theta }}_{n,l,m}\\ {\tilde{\chi }}_{n,l,m}\end{array}\right],$$

$${\left[A\right]}_{n,l}=\left[\begin{array}{cccc}\square & \varnothing & \varnothing & \varnothing \\ \varnothing & \square & \varnothing & \varnothing \\ \varnothing & \varnothing & \square & \varnothing \\ \varnothing & \varnothing & \varnothing & \square \end{array}\right],{\left[B\right]}_{n,l}=\left[\begin{array}{cccc}\square & \square & \square & \square \\ \square & \square & \varnothing & \varnothing \\ \square & \varnothing & \square & \varnothing \\ \square & \varnothing & \varnothing & \square \end{array}\right],$$

$$2n+l-m+2\le 3+2N,$$

$$\left(\mathrm{R_{t}}\ {\varTheta}+\mathrm{R_{c}}\ \chi\right){\bm{r}}=\mathrm{Ra}\,\big{(}\cos(\alpha){\varTheta}+\sin(\alpha)\chi\big{)}{\bm{r}},$$

$$\mathrm{Ra}=\sqrt{{{\mathrm{R}}_{\mathrm{t}}}^2+{{\mathrm{R}}_{\mathrm{c}}}^2}.$$

$$\alpha =\text{atan2}({\mathrm{R}}_{\mathrm{c}},{\mathrm{R}}_{\mathrm{t}}),$$

$$\Psi =\Theta \cos \alpha +\chi \sin \alpha ,{\Psi }^{\prime }=\Theta \cos \alpha -\chi \sin \alpha ,$$

$$\Xi =T\cos \alpha +C\sin \alpha ,{\Xi }^{\prime }=T\cos \alpha -C\sin \alpha ,$$

$$\Theta \cos \alpha =\frac{\Psi +{\Psi }^{\prime }}{2},\chi \sin \alpha =\frac{\Psi -{\Psi }^{\prime }}{2}.$$

$${\text{upartial}}_t\mathbf{u}=-\tau \hat{\mathbf{k}}\times \mathbf{u}-\nabla \pi +\mathrm{Ra}\Psi \mathbf{r}+{\nabla }^2\mathbf{u},$$

$${\text{upartial}}_t\Psi ={{\mathrm{P}}_{+}}^{-1}{\nabla }^2\Psi +{{\mathrm{P}}_{-}}^{-1}{\nabla }^2{\Psi }^{\prime }-\mathbf{u}\cdot \nabla \Xi ,$$

$${\text{upartial}}_t{\Psi }^{\prime }={{\mathrm{P}}_{+}}^{-1}{\nabla }^2{\Psi }^{\prime }+{{\mathrm{P}}_{-}}^{-1}{\nabla }^2\Psi -\mathbf{u}\cdot \nabla {\Xi }^{\prime },$$

$${{\mathrm{P}}_{+}}^{-1}=\frac{{\mathrm{Pr}}^{-1}+{\mathrm{Sc}}^{-1}}{2},{{\mathrm{P}}_{-}}^{-1}=\frac{{\mathrm{Pr}}^{-1}-{\mathrm{Sc}}^{-1}}{2}.$$

$${{\mathrm{P}}_{+}}^{-1}\approx \frac{1}{2}\max ({\mathrm{Pr}}^{-1},{\mathrm{Sc}}^{-1}),{{\mathrm{P}}_{-}}^{-1}\approx \pm \frac{1}{2}\max ({\mathrm{Pr}}^{-1},{\mathrm{Sc}}^{-1}),$$

$${\text{upartial}}_t\Psi \approx \frac{\max ({\mathrm{Pr}}^{-1},{\mathrm{Sc}}^{-1})}{2}{\nabla }^2(\Psi \pm {\Psi }^{\prime })-\mathbf{u}\cdot \nabla \Xi ,$$

$${\text{upartial}}_t{\Psi }^{\prime }\approx \frac{\max ({\mathrm{Pr}}^{-1},{\mathrm{Sc}}^{-1})}{2}{\nabla }^2({\Psi }^{\prime }\pm \Psi )-\mathbf{u}\cdot \nabla {\Xi }^{\prime },$$

$$\begin{array}{rl}{R}_{\text{crit}}=7.252{\left(\frac{P\tau }{1+P}\right)}^{4/3}& (1-\eta {)}^{7/3},m_{\text{crit}}=0.328{\left(\frac{P\tau }{1+P}\right)}^{1/3}(1-\eta {)}^{-2/3},\\ {\omega }_{\text{crit}}& =-0.762{\left(\frac{{\tau }^2}{P(1+P{)}^2}\right)}^{1/3}(1-\eta {)}^{2/3},\end{array}$$

$$\begin{array}{rl}\underset{P\to \infty }{\lim }{R}_{\text{crit}}=7.252{\tau }^{4/3}(1& -\eta {)}^{7/3},\underset{P\to \infty }{\lim }{m}_{\text{crit}}=0.328{\tau }^{1/3}(1-\eta {)}^{-2/3}\\ & \underset{P\to \infty }{\lim }{\omega }_{\text{crit}}\to 0.\end{array}$$

$$\frac{\text{upartial}{R}_{\text{crit}}}{\text{upartial}P}>0,\frac{\text{upartial}{m}_{\text{crit}}}{\text{upartial}P}>0,\frac{\text{upartial}{\omega }_{\text{crit}}}{\text{upartial}P}>0,$$

$${\mathrm{Ra}}_c=\frac{{\mathrm{Ra}}_0}{\cos \alpha +\sin \alpha }.$$

$$\mathrm{Sc}=\mathrm{Pr}(1+\delta ).$$

$${{\mathrm{P}}_{+}}^{-1}\approx \frac{{\mathrm{Pr}}^{-1}}{2}\left(2-\delta \right),{{\mathrm{P}}_{-}}^{-1}\approx \frac{{\mathrm{Pr}}^{-1}}{2}\delta ,$$

$${\text{upartial}}_t\Psi$$

$${\text{upartial}}_t{\Psi }^{\prime }$$

Full text

1 Introduction

Convection in a rotating spherical fluid shell provides one of the fundamental models for understanding the large-scale motions and the magnetic fields observed in many geophysical, planetary and astrophysical systems, see (Jones 2011, Glatzmaier 2013, Busse and Simitev 2015). Thermal and compositional convection occur, for instance, in Earth’s outer core which is composed mostly of iron and nickel, alloyed to lighter elements, supposedly silicon, sulphur and oxygen (Jeanloz 1990). Heat is continually lost to outer space establishing a secular-cooling thermal gradient that can drive thermal convection in the outer core. In addition, secular cooling leads to freezing of iron onto the inner core, a process in which both latent heat and light material are released and give rise to additional thermal and chemical buoyancy (Jacobs 1953, Braginsky 1963). Both buoyancy components are important when modelling the geodynamo. Indeed, recent estimates of thermal conductivity for iron at core conditions (Pozzo et al. 2012, Davies et al. 2015) confirm that outer core convection cannot be driven by thermal buoyancy alone. It is estimated that the compositional contribution to the buoyancy flux in the Earth’s core is around 80% (Lister and Buffett 1995). To model thermo-compositional convection Braginsky and Roberts (1995) suggested that temperature and concentration can be combined into a single “co-density” field. The co-density formulation has since been widely used in numerical simulations of the geodynamo and planetary dynamos, see reviews (Jones 2011, Christensen and Wicht 2015). However, the co-density formulation requires that a single set of effective boundary conditions, a single effective distribution of sources and a single effective value of the co-density diffusivity must be used. While the co-density formalism has been adapted to mimic thermo-compositional boundary conditions in the Core (Hori et al. 2012), stratification (Olson et al. 2017, Takehiro and Sasaki 2018), the Ganymede snow zone (Christensen 2015) and other situations, it is nor clear how well the it captures the dynamics of the distinct heat and composition fields. Only few modelling studies of core dynamics have employed a genuine ’double-diffusive’ formulation where both temperature and composition are included as separate fields. Breuer et al. (2010) found that the convective flow patterns differ significantly depending on the dominant driving component with abrupt changes in differential rotation, mean kinetic energy and transport efficiencies indicating a regime change when the relative contribution of compositional driving exceeds 20%. Trümper et al. (2012) extended the latter work and concluded that the main reason for this flow behaviour is the pronounced difference in diffusivities as measured by the Prandtl numbers of the thermal component and the compositional component. Trümper et al. (2012) also considered the effect of distinct boundary conditions for temperature and concentration and obtained some preliminary results on the onset of convection using an initial value code. Some comparison with (Trümper et al. 2012) is made in our article. In a double-diffusive model of Mercury’s dynamo Manglik et al. (2010) observed that when thermal and compositional buoyancy are of equal intensity, finger convection penetrates the upper layer enhancing the poloidal magnetic field, a significant difference compared to co-density cases. Exploring a geodynamo model Takahashi (2014) confirmed that the morphology of the poloidal field is determined by the balance of thermal and compositional driving with a dipolar magnetic field maintained when thermal buoyancy comprises less than 60% of the total driving and non-dipolar fields otherwise prevailing due to helicity reduction and concluded that the fraction of power injection by thermal convection in the present geodynamo is below this threshold.

A systematic study of the linear onset of double-diffusive convection in rotating spherical shells is required for understanding the behaviour of turbulent dynamos driven by thermal and compositional buoyancy. Indeed, it is often found that the properties of convection at onset provide much insight to finite-amplitude convection (Simitev and Busse 2003, Busse and Simitev 2006). Much is known about the onset of purely thermal convection – we refer to the recent monograph of Zhang and Liao (2017) which provides an extensive list of references. Interest in double-diffusive convection has been motivated mostly by oceanic applications. The most fascinating feature of double-diffusive convection is that instability can occur in a fluid where the vertical density gradient is statically stable. Surprisingly, this instability is due to diffusion, a typically stabilising effect. When the instability is driven by the small-diffusivity component it is called fingering; when the instability driven by the large-diffusivity component it is known as diffusive convection; see (Radko 2013). The linear onset in a two-component fluid layer constrained between two horizontal boundaries was studied by (Stern 1960, Nield 1967, Veronis 1968, Baines and Gill 1969) who naturally focussed on the fingering and diffusive regimes. The effect of rotation with respect to a vertical axis was considered by Pearlstein (1981) who found that a non-rotating layer can be destabilized by rotation, and that a rotating layer can be destabilized by the addition of a bottom-heavy solute gradient. These results were summarized in the context of core convection by Fearn et al. (1988). The effect of rotation perpendicular to gravity and the effect of inclined top and bottom boundaries is crucial in our case. These effects were considered by Busse (2002) and by Simitev (2011) in the setting of a rapidly rotating cylindrical annulus with conical caps. Busse (2002) established that when rotation and gravity are mutually perpendicular the second buoyancy component can balance the Coriolis force and so significantly facilitate convection. Simitev (2011) found that due to additional “double-diffusive” eigenmodes, the neutral curves for the onset of instability in a rotating annulus are typically multi-valued and form regions of instability in the parameter space that may be entirely disconnected from each other. It was also observed that while known asymptotic expressions for the critical Rayleigh number and frequency derived by Busse (2002) describe the onset of convection over an extended range of non-asymptotic parameter values but do not capture the full complexity of the critical curves. Net et al. (2012) studied numerically the influence of an externally enforced compositional gradient on the onset of convection of a mixture of two components in a rotating fluid spherical shell. Both positive and negative compositional gradients were considered in the latter study and it was found that the influence of the mixture is significant in both cases, in particular the critical values of the thermal Rayleigh number, of the frequency and of the wave number depends strongly on the direction of the compositional gradient. The aim of the present article is to extend and complement the latter studies and to contribute to a detailed linear analysis of double-diffusive convection in rotating spherical shells including fingering, diffusive as well as direct top-heavy regimes. Indeed, there is much to learn since both the parameter space and the space of valid modelling assumptions that can be made is very large. For instance, our model is mathematically similar to those of Busse (2002) and Simitev (2011) but is set in a spherical geometry. In contrast to Net et al. (2012), we consider the stress-free case for the velocity boundary conditions and internal rather than differential heating. Similarly to Net et al. (2012) we allow for both a stabilizing compositional gradient which may occur for instance, in lower main-sequence stars with heated helium-rich core surrounded by lighter hydrogen layers (Kippenhahn et al. 2012), as well as destabilising compositional gradients relevant, for instance, in the in the case of the Earth’s core where solidification with the release of light components takes place. In extension to Net et al. (2012) we consider negative values of the thermal expansion coefficient a situation relevant in some cases (Squyres et al. 1983, Röttger et al. 1994). In comparison to the linear analysis of (Trümper et al. 2012) we use an eigenvalue solver that offers many additional details on the onset of convective instability.

The article is organised as follows. We start by reviewing the mathematical set-up used to estimate the parameters of the system at onset in section 2. In section 3 we introduce a formalism that allows us to measure the combined effect of the thermal and compositional buoyancy and to better understand the convective processes at play. Section 4 proceeds to describe the competition of eigenmodes that leads to the formation of the global critical curves for onset. The remaining sections 5, 6 and 7 are devoted to describing the onset of convection depending on the Prandtl and Coriolis numbers and aspect ratio, respectively. A summary of the results and conclusions is presented in section 8.

2 Mathematical formulation and numerical solution

We investigate the onset of thermo-compositional convection in a rotating spherical shell. The shell has a thickness $d=r_{o}-r_{i}$, where $r_{o}$ and $r_{i}$ are the inner and the outer radii, respectively, and rotates about an axis aligned with the $z$-direction at a constant rate ${\varOmega}$. The unit vectors pointing in the $z$-direction and in the radial direction are denoted by ${\hat{\bm{k}}}$ and ${\hat{\bm{r}}}$, respectively, and ${\bm{r}}$ is the position vector with magnitude $r$. We assume that the spherical shell is full of incompressible fluid solution with constant kinematic viscosity $\nu$, thermal diffusivity $\kappa$, and chemical diffusivity $D$. We take the density $\rho$ of the fluid to depend linearly on changes in composition and temperature with first order expansion coefficients $\alpha_{C}$ and $\alpha_{T}$, respectively. We employ the Boussinesq approximation in that variations in density are assumed important only when affecting the gravitational force $-\rho\gamma{\bm{r}}$, with $\gamma$ a constant. In order to isolate the effects induced by differences in thermal and chemical diffusivities we follow (Busse 2002, Net et al. 2012) and disregard any differences in source-sink distribution and in boundary conditions for the temperature $T$ and the concentration $C$. Static profiles $T(r)$ and $C(r)$ with radial gradients $\upartial_{r}T=-\beta_{T}r$ and $\upartial_{r}C=-\beta_{C}r$ then exist assuming the temperature and concentration are fixed at the boundaries and have uniformly distributed sources with constant densities $\beta_{T}$ and $\beta_{C}$, respectively.

Using $d^{2}/\nu$ as the unit of time, $d$ as the unit of length, $T^{*}=\beta_{T}d^{2}\nu/\kappa$ as the unit of temperature and $C^{*}=\beta_{C}d^{2}\nu/D$ as the unit of concentration we arrive at the following linearised equations in adimensional units,

[TABLE]

where $\bm{u}$ is the flow velocity, $\pi$ is an effective pressure including all terms that can be written in gradient form, and ${\varTheta}$ and $\chi$ are the temperature and the compositional anomalies from the static reference states $T$ and $C$, respectively. The non-dimensional parameters are defined in Table 1. Note that here the compositional scale, $C^{*}$, is inversely proportional to $D$ in contrast with the work of Simitev (2011) who scales $C^{*}$ with $\kappa$. As a consequence, the compositional Rayleigh numbers reported bellow must be divided by the Lewis number, $\mathrm{Le}=\kappa/D$, when compared to the latter work.

Except when otherwise mentioned, the equations are solved for a spherical shell with an inner to outer radius ratio of $\eta=0.35$.

Exploiting the solenoidality of the velocity field we use the poloidal-toroidal decomposition

[TABLE]

and represent $\bm{u}$ in terms of a poloidal and a toroidal scalar functions ${\cal S}({\bm{r}},t)$ and ${\cal T}({\bm{r}},t)$, respectively. All unknown scalar quantities, ${\cal X}({\bm{r}},t)\equiv[{\cal S},{\cal T},{\varTheta},\chi]^{\top}({\bm{r}},t)$ are then assumed to obey a linear Fourier mode ansatz in time

[TABLE]

where $\omega$ is the frequency of oscillation (or drift rate), ${\varGamma}$ is the growth rate, and ⊤ denotes transpose. When ${\varGamma}$ is negative perturbations decay and the system is considered stable, otherwise instability occurs. Operating on equation (1a) by ${{\bm{r}}}{\bm{\cdot}}{\bm{\nabla}}\times$ and by ${{\bm{r}}}{\bm{\cdot}}{\bm{\nabla}}\times{\bm{\nabla}}\times$ four scalar equations are obtained

[TABLE]

In order to find the locus of marginal convective stability where ${\varGamma}=0$ the scalar quantities $\tilde{\cal X}({\bm{r}})$ are further expanded in terms of spherical harmonics for the angular part. Due to the linearity of the equations and of the orthogonality properties of spherical harmonics, individual azimuthal wave numbers $m$ decouple and can be investigated one by one. To complete the spatial discretisation of the problem in the radial direction we follow (Zhang and Busse 1987, Ardes et al. 1997) and expand $\tilde{\cal X}({\bm{r}})$ in trigonometric functions obeying the boundary conditions. After computing the appropriate Galerkin projection integrals numerically equations (3) take the matrix form

[TABLE]

for fixed $m$, with summation implied over the degree of the associated Legendre polynomials $l$ and the index of the radial functions $n$. Matrices $A$ and $B$ are of the form

[TABLE]

where squares represent non-null blocks. A triangular truncation of the sums is chosen such that there are the same number of radial functions as associated Legendre polynomials (Zhang and Busse 1987, Ardes et al. 1997),

[TABLE]

with $N$ an integer bigger than 2. $N$ represents the required resolution for the calculation and has been set to values higher than 10 and as close as feasibly possible to $m$. Equation (4) is then solved for the complex eigenvalues $({\mathrm{i}}\omega+{\varGamma})$ using standard numerical eigenvalue methods implemented in the LAPACK library111The LAPACK Library, Linear Algebra PACKage (LAPACK) http://www.netlib.org/lapack..

Once the eigenvalue problem (4) is solved for a set of fixed parameter values non-trivial numerical extremization and continuation problems must be tackled in order to follow the marginal stability curves in the parameter space. The numerical code for the solution of the problem can be obtained via (Silva and Simitev 2018).

3 An effective Rayleigh-number formalism

Multivalued critical curves for the onset of thermo-compositional convection in related problems were found in (Simitev 2011, Net et al. 2012) and occur in the present setting as well as illustrated in figure 1. It is difficult to find the extrema of multivalued curves and to this end we introduce below a new adimensional parameter $\mathrm{Ra}$ in terms of which the critical curves become single valued. This new parameter can be used to measure criticality with respect to the onset of convection similarly to the purely thermal case. An added advantage is that the associated transformation of the linearised governing equations makes it easy to interpret two important limiting cases as discussed within this section. To introduce $\mathrm{Ra}$ we note that the momentum equation (1a) can be transformed by representing the buoyancy force in the form

[TABLE]

with $\alpha$ varying between $-\pi/2$ and $3\pi/4$ and $\mathrm{Ra}$ being positive. We refer to $\mathrm{Ra}$ as the effective Rayleigh number. It is related to the thermal and the compositional Rayleigh numbers by

[TABLE]

We refer to $\alpha$ as the Rayleigh mixing angle since it corresponds to an angle in the $\mathrm{R_{t}}$-$\mathrm{R_{c}}$  plane. It is also defined in terms of the thermal and the compositional Rayleigh numbers,

[TABLE]

where the function $\text{atan2}(y,x)$ is the four-quadrant inverse tangent defined for $x\in\mathbb{R}$, $y\in\mathbb{R}$ as the principal argument $\text{Arg}(z)$ of the complex number $z=x+iy$, a notation used in many programming languages. The case of purely thermal convection corresponds to $\alpha=0$, the case of purely compositional convection is obtained at $\alpha=\pi/2$ and the typical co-density approach is corresponds to $\alpha=\pi/4$ with $\mathrm{Pr}=\mathrm{Sc}$.

We also remark that our parametrisation has several advantages over previous parametrisations such as the once proposed by Breuer et al. (2010) and Trümper et al. (2012). These authors considered a total Rayleigh number constructed simply as the sum of the thermal and compositional Rayleigh numbers. This is a very good approximation for the cases when the Prandtl and the Schmidt numbers are equal $\mathrm{Pr}=\mathrm{Sc}$ but fails when they are even only marginally different. In addition to converting critical curves to single valued functions, and providing a single positive measure of criticality, our reparametrisation, maps transitions between convective regimes and large Rayleigh number regions to asymptotes in the $\alpha$. In more practical terms, plotting on logarithmic scales can be employed due to the introduction of the strictly positive effective Rayleigh number.

The Rayleigh angle $\alpha$ introduced above appears in some respects similar to the Turner angle Tu defined by Ruddick (1983) but the two are not directly related as they arise in different models of double-diffusive convection. Tu is used in the oceanographic literature to determine the local stability of an inviscid saline water column (McDougall et al. 1988) and arises as a parameter in the so called “unbounded gradient layer” model where density increases downwards and instabilities are due to local gradients of temperature and concentration (Radko 2013). In this case double-diffusive convection is characterised by an intrinsic length scale which makes it possible to non-dimensionalise this model so that Tu (or equivalently the so called stability density ratio $R_{\rho}$) is the only parameter needed to control instability and identify “salt-fingering” and “diffusive” regimes. In contrast, we are motivated by planetary core applications where “bottom-light” convection with concurrently destabilizing temperature and composition is arguably more important. The natural length scale of vigorous bottom-light convection is not a function only of the ratio of local gradients of temperature and concentration but is also determined by other factors including geometry and driving. To account for these we employ a “vertically-bounded layer” model of the type first used by Veronis (1968). Here buoyancy forces are driven by the global variation or temperature and concentration across the volume of the spherical shell and two Rayleigh numbers are now needed to describe the convective instability. The effective Rayleigh number $\mathrm{Ra}$ and the Rayleigh angle $\alpha$ are a convenient reparametrisation of those as discussed above.

To take advantage of our reparametrisation, we now define a new set of dynamical variables ${\varPsi}$ and ${\varPsi}^{\prime}$ and static reference profiles ${\varXi}$ and ${\varXi}^{\prime}$ using the transformations

[TABLE]

with the inverses being

[TABLE]

By substituting these expressions into equations (1a–1d) and adding and subtracting equations (1c) and (1d), we arrive at the transformed problem

[TABLE]

with effective Prandtl numbers defined as

[TABLE]

Note that, while the number $\mathrm{P_{+}}$ is strictly positive, the number $\mathrm{P_{-}}$ can be negative thus providing means of concentrating rather than diffusing the fields ${\varPsi}$ or ${\varPsi}^{\prime}$. This mechanism can drive convection at lower than expected Rayleigh numbers and at very large scales as we shall see.

Formulation (11) allows us to consider two important limiting cases. Firstly, we note that the transformed momentum equation (11a) only depends on the field ${\varPsi}$. A consequence of this decoupling is that, in the case of $\mathrm{P_{-}}^{-1}$ approaching 0, the system will behave as if driven only by one buoyancy generating field. This situation arises when $\mathrm{Pr}$ is close or equal to $\mathrm{Sc}$ and thus, we obtain the usual co-density approximation. We consider this case in more depth in sections 5.1 and 5.2.

A second limiting case occurs when the Prandtl numbers are significantly different in magnitude. Then the new parameters can be approximated as

[TABLE]

and equations (11b,c) can be rewritten approximately as

[TABLE]

with the symbol $\pm$ taking the positive sign when $\mathrm{Pr}^{-1}\gg\mathrm{Sc}^{-1}$, and the negative sign otherwise. With all other parameters fixed, the changes to the critical value of $\mathrm{Ra}$ can be due only to different behaviours of the system with respect to $\alpha$. We consider this case in more detail in section 5.3 where we will also explore the effects of $\mathrm{Pr}$ on the curves $\mathrm{Ra}_{c}(\alpha)$.

4 Construction of global critical curves

The introduction of the equation of concentration in system (1) leads to the occurrence of additional eigenmodes not present in the purely thermal case. This is due to the increased dimension of the of the eigenvalue problem (4) compared to the corresponding purely thermal problem, as is well illustrated in (Simitev 2011). As an example, figure 2 shows a comparison between selected eigenmodes of a thermo-compositional case and the eigenmodes of the corresponding purely thermal case for two particular azimuthal wave numbers $m=4$ and $m=7$ with $\mathrm{R_{c}}=4\times 10^{4}$, $\tau=1.2\times 10^{3}$ and $\mathrm{Sc}=10$. Two types of eigenmodes are observed – (a) modes that occur in both the thermo-compositional and the purely thermal case and appear to correspond to each other, and (b) additional modes that have no counterparts in the purely thermal case. These latter modes appear to represent the “fingering” and “diffusive” instability modes known from classical studies of double-diffusive convection, see (Radko 2013), but modified by rotation and spherical geometry, see further below for discussion. The purely double-diffusive modes of type (b) can become unstable in parameter regions disconnected from the main region of convective instability as indicated by their positive growth rate shown in the left panel of figure 2. This was also observed in the case of the cylindrical annulus (Simitev 2011). Whether such a disconnected region will remain present in the global critical curve depends on the overlap with the instability regions of the eigenmodes of other wave numbers. This leads us to a discussion of the construction of global critical curves below.

Instability arises whenever an eigenmode attains a positive growth rate, as illustrated by the shaded regions in figure 2. Thus, for each azimuthal wave-number $m$ a composite critical curve is obtained which consists of pieces corresponding to the critical curves of the structural eigenmodes that become unstable first. Figures 1 and 3 present examples of critical curves for fixed values of the wave number $m$. The location of the transition between instability due to modes of type (a) and instability due to modes of type (b) is easy to recognise in these figures. For instance, one such transition occurs in figure 1(a) for $m=4$ at $\mathrm{R_{c}}=-4.8\times 10^{4}$ where the curve is not smooth. We identify the arising purely double-diffusive mode as being analogous to a diffusive instability in a non-rotating layer of cold fresh water over warm salty water because to the left of this transition $\mathrm{R_{c}}$ is negative and $\mathrm{R_{t}}$ is positive which in our model corresponds to cooler and lighter fluid over warmer and heavier fluid. A second such transition occurs for the same $m$ at $\mathrm{R_{c}}=5.8\times 10^{4}$ where the curve folds back on itself and becomes three-valued. We identify the arising purely double-diffusive mode as being analogous to a fingering instability in a non-rotating layer of warm salty water over cold fresh water because the substantial part this branch occupies a region where $\mathrm{R_{c}}$ is positive and $\mathrm{R_{t}}$ is negative which in our model corresponds to warmer and heavier fluid over cooler and lighter fluid. A similar fingering branch where the critical curve folds back on itself is visible in figures 2(a) and 5(a) of (Net et al. 2012) in the case of a shell with rigid boundaries and with differential rather than internal driving, as well as in figure 4(b) of (Simitev 2011) in the rotating annulus geometry. In the figure of Simitev (2011) the diffusive branch is also reported. While fingering instability occurs via exchange of stabilities in the classical doubly-diffusive plane layer model irrespective of rotation, in the spherical geometry any mode including fingering is overstable in order to overcome the Taylor-Proundman constraint. Overstability takes the form of travelling waves similarly to purely thermal motions (Busse 1970). Figure 1 also illustrates the appearance of critical curves when plotted in the $\alpha-\mathrm{Ra}$ plane.

Finally, the global critical curve can be conveniently constructed in the $\alpha-\mathrm{Ra}$ plane as the lower envelope of the composite curves for each wave number $m$. In other words, at each value of $\alpha$ instability occurs at the smallest value of $\mathrm{Ra}$ among the critical curves for the separate values of $m$. This is illustrated in the example shown in figure 3. There, we can also see that the critical curves for different values of $m$ intersect and overlap allowing convection for only certain $m$-values and not for others. Some azimuthal wave-numbers will contribute to large portions of the global critical curve (e.g. $m=9$ in figure 3); others will contribute only point-wisely (e.g. $m=15$) and others not at all (e.g. $m=1$). Distinct branches of the global critical curve can be easily identified and these correspond to the various convective regimes discussed above as well as further in the text.

In addition to the purely double-diffusive modes of type (b), here we are also interested in the top-heavy convection modes of type (a) as these are most likely the modes responsible for the generation of magnetic field in the core. These modes appear analogous to purely thermal convection modes as discussed in relation to figure 2. To confirm this and to help interpret parameter dependences, in subsequent sections we compare certain branches of our numerical solutions to known closed-form approximations of the critical Rayleigh number, azimuthal wave number and drift-rate for the onset of columnar convection in rapidly rotating axisymmetric systems. The main difficulties for analytical study of convection in rotating spheres and shells are the significant variation in the Coriolis force with the angle between gravity and angular velocity and the inclination of the spherical boundaries with respect to the axis of rotation. In the asymptotic limit $\tau\to\infty$ a so called “local” theory of viscous convection in an internally heated sphere was developed by Roberts (1968) and Busse (1970). Soward (1977) showed that the critical mode predicted by the local theory is not physically realized as it decays due to phase mixing. In particular, in the case of internal heating convection onsets deep inside the shell and the true value of the critical Rayleigh number may exceed the local theory estimate by about 25% (Jones 2015). Following Yano (1992), a “global” asymptotic theory of the onset of viscous convection was constructed for rapidly-rotating spheres (Jones et al. 2000) and for spherical shells (Dormy et al. 2004) and results were shown to agree with accurate numerical solutions of the corresponding linear problems. Because the global theory does not in general yield expressions in convenient closed analytical form, we use here the following less accurate local theory approximations

[TABLE]

where $P$ refers to either the thermal Prandtl number $\mathrm{Pr}$ or the Schmidt number $\mathrm{Sc}$. These expressions are derived from equations (2.7), (4.6) and (4.7) of Yano (1992) by re-scaling in terms of the dimensionless parameters, length and time scales used in the present paper. In particular, the factor $(1-\eta)=(r_{o}-r_{i})/r_{o}=d/r_{o}$ enters due to length scale conversion: the radius of the outer surface $r_{o}$ is used as length scale by Yano (1992) while the shell thickness $d=r_{o}-r_{i}$ is used in our work. In turn, the expressions of Yano (1992) provide a first-order correction of the local asymptotic results of Busse (1970) and are thus expected to better account for the effects of finite inclination of the outer spherical boundary. While expressions (15) are not valid asymptotic results for the configuration studied here we consider them adequate for the purpose of identifying trends in the numerical results presented in our paper. They provide an approximate idea of the dependence of the critical values for the onset of convection on the parameters of the problem and similar expressions have been used in this sense in other studies e.g. (Simitev and Busse 2003).

In the limit of $P\to\infty$ (i.e. either $\mathrm{Pr}\to\infty$ or $\mathrm{Sc}\to\infty$) we obtain

[TABLE]

The derivatives with respect to $P$ of equation (15) are all positive

[TABLE]

with all derivatives vanishing as $P\to\infty$. Therefore, we expect that in the region $0\leq\alpha\leq\pi/2$ where $\mathrm{R_{t}}$ and $\mathrm{R_{c}}$ are positive and specifically at $\alpha=0$ (purely thermal case) and $\alpha=\pi/2$ (purely compositional case) the values of $\mathrm{Ra}$ and $m$ will increase towards limiting values while $\omega$ will decrease towards zero as $\mathrm{Pr}$ is increased. Equations (15) also show that $R_{\text{crit}}\propto\tau^{4/3}$, $m_{\text{crit}}\propto\tau^{1/3}$ and $\omega_{\text{crit}}\propto\tau^{2/3}$. Therefore, we expect that the critical Rayleigh and wave numbers and drift-rate will increase as $\tau$ is increased. Finally, the $\eta$-dependence in equations (15) is $R_{\text{crit}}\propto(1-\eta)^{{7/3}}$, $m_{\text{crit}}\propto(1-\eta)^{-2/3}$ and $\omega_{\text{crit}}\propto(1-\eta)^{2/3}$. Therefore, we expect that the critical Rayleigh number and drift-rate will decrease and the critical azimuthal wave number will increase as $\eta$ is increased. Equations (15) are appropriate only for sufficiently large $P\tau\gg 1$. At small values of $P\tau$ convection takes the form of eddies attached to the outer equatorial boundary that can be understood as inertial oscillations modified by viscous friction and thermal buoyancy effects. Asymptotic expressions for this case can be found in Zhang (1994), Busse and Simitev (2004) but will not be used here since $P\tau$ is kept relatively large.

5 Dependence on the Prandtl number

5.1 The case of equal Prandtl numbers

The simplest case to consider is the case of equal Prandtl and Schmidt numbers as it is essentially equivalent to the case of purely thermal convection. When temperature and composition have the same diffusivities the thermal and the compositional Prandtl numbers are identical, i.e. $\mathrm{Pr}=\mathrm{Sc}$, and composition and temperature will evolve in a similar manner. In terms of equations (11), the parameter $\mathrm{P_{-}}$ vanishes and as a result equations (11b) and (11c) decouple from each other, the momentum equation (11a) depends only on ${\varPsi}$ while ${\varPsi}^{\prime}$ becomes a passive tracer for the flow. System (11) reduces to the familiar co-density model. The critical Rayleigh number dependence takes the form of a straight line $\mathrm{R_{t}}=-\mathrm{R_{c}}+\mathrm{Ra}_{0}$ in the $\mathrm{R_{t}}-\mathrm{R_{c}}$ plane, as proposed by Breuer et al. (2010), and in the $\mathrm{Ra}-\alpha$ plane we obtain

[TABLE]

Figure 4(a) shows the critical curves $\mathrm{Ra}_{c}$ as a function of $\alpha$ for $\mathrm{Pr}$ varying between $\mathrm{Pr}=10^{-5}$ (lower curves) to $\mathrm{Pr}=10^{4}$ (higher curves) and for values of the Coriolis number and the shell aspect ratio fixed to $\tau=10^{4}$ and $\eta=0.35$ (the dependence on $\tau$ and $\eta$ is discussed in later sections). All curves are symmetric with respect to $\alpha=\pi/4$ (equal Rayleigh numbers) and approach infinity with asymptotes at $\alpha=-\pi/4$ and $\alpha=3\pi/4$. Equations (11a) and (11b) show that the larger $\mathrm{Pr}$ becomes, the smaller its influence is. The coefficient $\mathrm{Ra}_{0}$ in (18) should then become independent of $\mathrm{Pr}$ at large values of this parameter. We will refer to this regime as the advective regime. At Prandtl numbers much smaller than one the process of diffusion takes over and prevents small scale convection from growing. Then much larger scales dominate and the onset of convection occurs at lower values of the Rayleigh numbers. We will refer to this regime as the diffusive regime. It is easy to recognise that the advective and the diffusive regimes are, in fact, the familiar columnar convection regime (Busse 1970) and inertial convection (Zhang and Busse 1987, Zhang 1994) regime, respectively. The transition between these two regimes occurs at about $\mathrm{Pr}=0.1$ for $\tau=10^{4}$ (Ardes et al. 1997, Busse and Simitev 2004). To confirm this we plot the critical values of $\mathrm{Ra}_{0}$, of the azimuthal wave number $m$ and of the drift rate $\omega$ as a function of $\mathrm{Pr}$ in figure 4(b). The values of $\mathrm{Ra}_{0}$ and the $m$ saturate in the limits $\Pr\to 0$ and $\Pr\to\infty$, with small values of $m$ ($m=2$) for $\Pr\to 0$ and large values of $m$ ($m\geq 14$) for $\Pr\to\infty$. The transition between the diffusive and the advective regime is easy to identify in figure 4(b). We have also plotted for comparison the approximations (15) for the critical parameter of columnar convection. As noted previously, approximations for the critical parameters of inertial convection are also available (Zhang 1994, Busse and Simitev 2004) but have not be included in figure 4(b).

The critical modes exhibit the typical spacial features of inertial convection for smaller values of $\Pr$ and of columnar convection for larger values of $\Pr$. For lower values of the $\Pr$, convection takes the form of large non-spiralling convection eddies attached to the outer boundary of the shell near the equator and exhibit a strong clockwise drift. Above the abrupt transition at about $\mathrm{Pr}=0.1$ convection columns with much smaller azimuthal scale appear. Near the transition the structure of columns resembles the double-humped mode described by Ardes et al. (1997). By $\mathrm{Pr}=0.2$ columnar convection is fully established and show a typical spiralling shape. Convection columns are detached from outher spherical surface and are arrange in a cartridge belt immediately outside of the tangent cylinder. Inertial and columnar convection regimes with properties similar to these described here for $\tau=10^{4}$ exist also for more general values of $\tau$ (Zhang 1994, Simitev and Busse 2003).

5.2 Small departures from the $\mathrm{Pr}=\mathrm{Sc}$ regime

We now depart from the co-density case of a “single-diffusive” fluid characterised by equal Prandtl numbers. In this section we develop a first order approximation to equations (11) that will allow for the analysis of small departures from the regime of $\mathrm{Pr}=\mathrm{Sc}$. We consider a small number $|\delta|\ll 1$ such that the Prandtl and Schmidt numbers are similar in value

[TABLE]

To first order in $\delta$ the modified Prandtl numbers introduced in equations (12) become

[TABLE]

with $\mathrm{P_{+}}^{-1}$ being always much larger than $\mathrm{P_{-}}^{-1}$ but the latter still finite. With these approximations equations (11b) and (11c) become to $O(\delta)$,

[TABLE]

Using definitions (9) and the reference temperature and composition profiles equations (21) can be rewritten in the form

[TABLE]

where $\beta=\beta_{\mathrm{T}}=\beta_{\mathrm{C}}$ was assumed. Equations (22) are only coupled to each other by the last term in  (22a) which is multiplied by the small quantity $\delta$ and so the diffusion of $\chi$ has only a small contribution to the evolution of the buoyancy profile ${\varPsi}$. It is now easy to understand the physical effects of deviating from the case of equal Prandtl numbers as discussed bellow. These are illustrated in figures 5 and 6 for values of the Schmidt number $\mathrm{Sc}$ just below and just above $\mathrm{Pr}=1$, respectively, and for fixed values of the Coriolis number and the shell aspect ratio $\tau=10^{4}$ and $\eta=0.35$.

When $\mathrm{Sc}<\mathrm{Pr}$, so that $\delta<0$, and $0<\alpha<\pi/2$ (positive compositional and thermal Rayleigh numbers) equation (22a) shows that the growth in time of the compositional component $\chi$ has a small but destabilising effect on the buoyancy profile ${\varPsi}$. Thus convection can now occur for lower values of $\mathrm{Ra}$ as $|\delta|$ becomes larger. This is evident in the bottom panel of figure 5(a) where, upon decreasing delta from 0 to $-0.5$, a reduction in $\mathrm{Ra}$ is indeed observed. The magnitude of this reduction increases as $\alpha$ is increased from 0 to $\pi/2$. The middle panel of figure 5(b) shows this effect in the $\mathrm{R_{t}}-\mathrm{R_{c}}$ plane. More negative values of $\delta$ correspond to smaller values of $\mathrm{Sc}$ so that less energy is required for the onset of convection and both critical Rayleigh numbers (thermal and compositional) decrease. This confirms the results of Simitev (2011) where for positive Rayleigh numbers, an increase in $\mathrm{R_{c}}$ leads to a reduction in the critical value of $\mathrm{R_{t}}$ required for the onset of convection.

When $\mathrm{Pr}<\mathrm{Sc}$, so that $\delta>0$, and $0<\alpha<\pi/2$ the diffusion of $\chi$ has a stabilising effect on the buoyancy profile, ${\varPsi}$, with convection occurring at higher critical values of $\mathrm{Ra}$ as $\delta$ increases. This effect is evident in the bottom panel of figure 6(a), where an increase in $\mathrm{Ra}$ is seen as $\delta$ is increased. Similarly, this effect becomes larger as $\alpha$ is increased from 0 to $\pi/2$ as illustrated in the lower panel of figure 6(a) and in the middle panel of figure 6(b).

For sufficiently small negative values of alpha, $\alpha<0$, $|\alpha|\ll 1$ and $\mathrm{Sc}<\mathrm{Pr}$ the sign of $\sin\alpha$ is negative and the diffusion of $\chi$ has a stabilising effect as $\delta\sin\alpha$ is positive. This is illustrated in the top panel of figure 5(b) where as $\alpha$ changes from positive to negative and we see that decreasing $\delta$ increases the critical value of the thermal Rayleigh number $\mathrm{R_{t}}$ whereas when $\alpha>0$ decreasing $\delta$ decreases the critical $\mathrm{R_{t}}$. Similarly, when $\alpha>0$, $|\alpha|\ll 1$ and $\mathrm{Sc}>\mathrm{Pr}$ the diffusion of $\chi$ leads to a reduction of the buoyancy profile. This is seen most clearly in the upper panel of figure 6(b). Increasing $\delta$ leads to a reduction of the critical value of $\mathrm{R_{t}}$. This phenomenon can be explained by the fact that, when $\alpha<0$, the compositional buoyancy is now acting against the thermal component of buoyancy. Therefore, an increase in $\mathrm{Sc}$ inhibits the effect of the compositional component and thus lowers the critical value of $\mathrm{R_{t}}$.

In summary, the departure from equal Prandtl numbers has a small but noticeable affect on the onset of convection, manifested by a deviation from the “co-density” values of $\mathrm{Ra}$ obtained in the case $\mathrm{Pr}=\mathrm{Sc}$. However, the regions where $\mathrm{R_{t}}$ and $\mathrm{R_{c}}$ have opposite signs, namely $-\pi/2<\alpha<0$ and $\pi/2<\alpha<\pi$, will be most strongly affected as a delicate balance exists there between stabilising and destabilising forces.

The lower panels of figures 5(a) and 6(a) show that the asymptotic behaviour exhibited by the critical curves $\mathrm{Pr}=\mathrm{Sc}$ as $\alpha\rightarrow-\pi/4$ and $\alpha\rightarrow 3\pi/4$ is followed on only one side by the critical curves with $\mathrm{Pr}\neq\mathrm{Sc}$. When $\mathrm{Sc}>\mathrm{Pr}$ the asymptote at $\alpha\rightarrow-\pi/4$ is approached, when $\mathrm{Sc}<\mathrm{Pr}$ the asymptote at $\alpha\rightarrow 3\pi/4$ asymptote is approached. Since the shift of the asymptote is due to coupling between equations (22), we consider small departures from the values $\alpha=-\pi/4$ and $3\pi/4$ by setting $\alpha=-\pi/4+\epsilon$ and $3\pi/4+\epsilon$, where $|\epsilon|\ll 1$. Expanding $\cos\alpha$ and $\sin\alpha$ in Taylor series around these values equation (22a) becomes

[TABLE]

where, the upper and lower signs of the $\pm$ term refers to $\alpha=-\pi/4$ and $3\pi/4$, respectively. We now consider the two cases $\mathrm{Pr}>\mathrm{Sc}$ ($\delta<0$) and $\mathrm{Pr}<\mathrm{Sc}$ ($\delta>0$) in turn. Firstly, when $\mathrm{Pr}>\mathrm{Sc}$ and $\alpha=3\pi/4+\epsilon$ with $\epsilon>0$ (the negative of $\pm$ is taken) the diffusion of $\chi$ is destabilising and the advective term $\epsilon\sqrt{2}\beta u_{r}$ is stabilising. Clearly when $\delta=\epsilon=0$ the buoyancy profile simply diffuses away. As $\delta$ is decreased the gradient of the buoyancy profile is increased and thus $\epsilon$ can take greater values until the advective term becomes greater than the diffusion of $\chi$. This is supported by the bottom panel of figure 5(a) where the asymptote of $\mathrm{Ra}$ is shifted further to the right as $|\delta|$ increases. This behaviour is also observed in the bottom panel of figure 5(b) where a shift is seen away from the line $\alpha=3\pi/4$ towards $\alpha=\pi$ as $\delta$ decreases from [math] to $-0.5$ and thus lower critical values of $\mathrm{R_{c}}$ are required. Secondly, when $\mathrm{Pr}<\mathrm{Sc}$ the asymptote shift is mirror reflected at $\alpha=\pi/4$ and occurs to the left near $\alpha=-\pi/4$. When $\delta>0$ and $\alpha=-\pi/4+\epsilon$ with $\epsilon<0$, the positive sign of equation (23) shows again that the diffusion of $\chi$ and the advective term are destabilising and stabilising, respectively. Clearly, when the diffusive term can overcome the advective term, convection can occur. Therefore, as $\delta$ is increased the asymptote will shift, with negative values of $\epsilon$ of higher magnitude able to still support convection. This continues with increasing $\epsilon$ until the magnitude of the advective term becomes too large and convection is unable to occur. The bottom panels of figures 6(a) and figures 6(b) support this and show a shift of the asymptote from $\alpha=-\pi/4$ to $\alpha=-\pi/2$ and thus lower critical values of $\mathrm{R_{t}}$ are required.

An interesting feature is the jump from small scale to large scale convection immediately beyond the co-density asymptotes ($\alpha=-\pi/4,3\pi/4$) of the case $\mathrm{Pr}=\mathrm{Sc}$. The azimuthal wave-number decreases significantly beyond the co-density asymptote for non-zero $\delta$ as seen in the middle panels of figures 5a and 6a and a large drop in magnitude of the drift-rate $\omega$ is associated with this (see the top panels of figures 5a and 6a). However the wave-number and drift-rate quickly rise towards the new asymptotes. Another interesting effect occurs on the opposite side to the asymptote shifts. Here the co-density asymptotes are still in effect ($\delta=0$) but as $|\delta|$ is increased there is a reduction in the critical Rayleigh number. Again we the switch from small scale to large scale convection around these asymptotes.

Finally, we mention that the calculations shown in figures 5 and 6 are for the specific value $\mathrm{Pr}=1$ but our simulations show that similar results are obtained for other values of $\mathrm{Pr}$.

5.3 Large differences between Prandtl numbers

In this section we investigate the case when the values of the Prandtl and the Schmidt numbers, $\mathrm{Pr}$ and $\mathrm{Sc}$, are significantly different from each other. First, recall that equations (13) and (14) model the evolution of the buoyancy field in this case. For brevity, we will only analyse the case of $\mathrm{Sc}\gg\mathrm{Pr}$ as the case of $\mathrm{Sc}\ll\mathrm{Pr}$ is symmetric with respect to $\alpha=\pi/4$. In particular, we fix the values of the Schmidt number to $\mathrm{Sc}=25\mathrm{Pr}$, the Corriolis number to $\tau=10^{4}$ and the shell aspect ratio to $\eta=0.35$ in this section. Here equation (14a) becomes

[TABLE]

The evolution of ${\varPsi}$ is now essentially only affected by the diffusion of the physical field which is characterised by the smaller of the values of $\mathrm{Pr}$ and $\mathrm{Sc}$. Since convection regimes are rather different for values of the Prandtl numbers above and below 0.1 we present our analysis in these two cases. At Prandtl numbers smaller than 0.1 convection is generally large scale and three distinct regions of $\mathrm{Ra}_{c}(\alpha)$ can be identified. At Prandtl numbers greater than 0.1 instead of three, four distinct regions emerge all which are functions of $\alpha$. Figure 7 shows the results of our calculations.

5.3.1 Small Prandtl number case

At the lowest value of $\mathrm{Pr}$ in figure 7(a), $\mathrm{Pr}=10^{-5}$, the critical curve $\mathrm{Ra}$ is approximately symmetric with respect to $\alpha=\pi/4$ for values of the mixing angle from $\alpha=-\pi/4$ to $\alpha=3\pi/4$. In this range the null curve is straight line of negative gradient when plotted in the $\mathrm{R_{c}}-\mathrm{R_{t}}$ plane. At values $\mathrm{Pr}=10^{-5}$ and $\mathrm{Sc}=2.5\times 10^{-4}$ there is little difference in how thermal and compositional components affect convection which is large scale ($m=2$) with a prograde drift ($\omega<0$) and is of equatorially attached type. Figure 8 shows the stream function ($1/r\upartial{\cal S}/\upartial\theta$) in the equatorial plane for the case of $\mathrm{Pr}=10^{-5}$ (top left). This situation is easy to understand in terms of equation (11b) as any small scale buoyancy anomalies are rapidly diffused away due to small values of the Prandtl numbers; the advection of the background profiles cannot overcome such strong diffusion.

The behaviour in the region between mixing angles $-\pi/2<\alpha$ and $\alpha<-\pi/4$ is also easy to interpret – here the buoyancy profile ${\varXi}$ has a negative sign ($\cos\alpha+\sin\alpha<0$) and thus acts to stabilise the system. In fact, convection is only possible because the advection profile ${\varXi}^{\prime}$ is now approaching a maximum and therefore generates enough ${\varPsi}^{\prime}$ to fuel the production of ${\varPsi}$ in equation (11b) and consequently the velocity ${\bm{u}}$ in equation (11a). This, however, comes at the cost of requiring a higher value of $\mathrm{Ra}$ for the onset of convection to take place.

5.3.2 Intermediate Prandtl number case

Increasing $\mathrm{Pr}$ from $10^{-5}$ to $10^{-2}$ in figure 7(a) results in a general shift of the critical $\mathrm{Ra}$ curve upwards and to the left, as most clearly seen by comparing the curves $\mathrm{Pr}=10^{-5}$ and $\mathrm{Pr}=10^{-2}$. Wave number values increase gently from $m=3$ to $m=7$ so that convection is now on a slightly smaller scale. At $\mathrm{Pr}=10^{-1}$ an overall increase in the critical value of $\mathrm{Ra}$ occurs due to the transition from inertial convection to columnar convection. An interesting effect is that, as we switch from thermally controlled to compositionally controlled convection at about $\alpha=\pi/3$, there is a transition from large scale ($m=6$) to a significantly smaller scale mixed convection ($m=15$) which has already been observed by Trümper et al. (2012) with $\mathrm{Le}=30$. As $\alpha$ increases towards $\pi/2$ (purely compositional convection), the value of $m$ decreases again to $m=10$ in analogy to purely thermal convection with internally distributed heat sources and a Prandtl number of 2.5.

Figure 8 illustrates the corresponding patterns of convection with the case $\mathrm{Pr}=10^{-3}$ plotted in the top right panel and the case $\mathrm{Pr}=10^{-1}$ plotted in the lower left panel.

5.3.3 Large Prandtl number case

The onset of convection for large values of the Prandtl number, $\mathrm{Pr}>0.1$, is illustrated in figure 7(b) where the drift rate, the most unstable wave number and the critical effective Rayleigh number are plotted as a function of the mixing angle $\alpha$. As $\mathrm{Pr}$ is increased $\mathrm{Ra}$ and $m$ approach limiting values and $\omega$ appears to tend to zero as expected from the approximation (15). Four different regimes of convection can be identified. These are described in turn below starting from negative values of $\alpha$. We focus on the curve $P=1$ in order to be able to quote specific parameter values.

For values of $\alpha$ smaller than approximately $-\pi/4$ motions are due to a diffusive instability where the most unstable mode is an additional double-diffusive mode as described in relation to figures 1 and 2. The gradient of the buoyancy profile ${\varXi}$ changes sign with respect to the gradients of $T(r)$ and $C(r)$ and is shallower. A similar situation was analysed in section 5.2 where an asymptote located at $\alpha=-\pi/4$ has now shifted to $\alpha=-\pi/2$. The diffusion contribution to the buoyancy field ${\varPsi}$ is still dominated by the temperature but it is strongly reduced since $\cos\alpha<1/2$. The buoyancy force, on the other hand, is dominated by the concentration but acts in reversed direction pushing regions of positive ${\varPsi}$ down and regions of negative ${\varPsi}$ up. This gives rise to large scale ($m=3$), slow convection drifting counter-clockwise as illustrated in figure 7 (bottom right plot, panel (a)). The onset of convection occurs at lower values of the effective Rayleigh number than for purely thermal convection.

For values of $\alpha$ approximately between $-\pi/4$ and $\pi/4$ motions are due to single-diffusive eigenmodes of purely thermal nature. Indeed, when values of $\alpha$ are close to zero $\cos\alpha\approx 1-\alpha^{2}/2$ and $\sin\alpha\approx\alpha$ so that the temperature remains the only quantity being diffused in equation (14a) and it is also the main source of buoyancy. The system behaves like a single-diffusive purely thermal convection at $\mathrm{Pr}=10$ and exhibits values of the azimuthal wave number, drift rate and effective Rayleigh number as expected from equations (15), with the effective Rayleigh number modified according to equation (18). The spatial pattern of the flow is illustrated in figure 7 (bottom right plot, panel (b)). Abrupt transitions to neighbouring regimes occur at both ends of this region, namely an abrupt transition to diffusive-instability at $\alpha\approx-\pi/4$ and an abrupt transition to chemical convection at $\alpha\approx\pi/12$.

For values of $\alpha$ approximately between $\pi/4$ and $\pi/2$ motions are due to motions are due to single-diffusive eigenmodes of purely chemical nature. These are modes otherwise similar to the thermal convective modes of the regime to the left but characterised by a different value of the Prandtl number – the value of the Schmidt number is $\mathrm{Sc}=250$. This relatively large change in the values of the Prandtl numbers results in a large abrupt jump in wave number to much smaller scales, reduction of the drift rate by an order of magnitude and a discrete change in slope of the $\mathrm{Ra}$ curve. Because these are effectively single-diffusive modes the critical values are well approximated by equations (15) with $\mathrm{Ra}$ modified by equation (18). The spatial pattern of the flow is illustrated in figure 7 (bottom right plot, panel (c)).

For values of $\alpha$ greater than approximately $\pi/2$ motions are due to a fingering instability where the most unstable mode is an additional double-diffusive mode as described in relation to figures 1 and 2. For $\alpha>\pi/2$ the factor $\cos\alpha$ becomes negative and diffusion then acts to concentrate anomalies of ${\varPsi}$ around the anomalies of ${\varTheta}$ instead of dispersing them. However, because ${\varXi}$ is now smaller than either of $T$ or $C$, there will be less buoyancy produced by advection. This leads a relatively gradual transition between the purely chemical and the fingering regime as illustrated in figure 7(b) where a continuous decrease in the value of $m$ towards a constant value of $m=3$ and a continuous change in the slope of the $\mathrm{Ra}$ curve are observed. Only the drift rate shows a small discrete jump to smaller values. The spatial pattern of the flow is illustrated in figure 7 (bottom right plot, panel (d)). The asymptote on the far right side remains at $3\pi/4$ as was the case in discussed in relation to figure 6.

6 Dependence on the Coriolis number

So far we have analysed cases at fixed values of the Coriolis number, $\tau$. In this section we investigate the dependence of the onset of convection on this adimensional parameter. In the case of purely thermal convection the thermal Rayleigh number $\mathrm{R_{t}}$ of columnar convection is proportional to $\tau^{4/3}$ as approximated by equation (15). Here, we will show that the effective Rayleigh number for double-buoyant, double-diffusive convection obeys the same scaling in regions where the lowest diffusion quantity dominates in generating buoyancy including in the regions where both Rayleigh numbers are positive. We will also establish numerically an approximate scaling valid in the regions of diffusive and fingering regimes.

Figure 9(a) shows critical parameter values for the onset of convection as a function of $\alpha$ for selected values of $\tau$ and for fixed valued of the Prandtl numbers $\mathrm{Pr}=1$ and $\mathrm{Sc}=100$ and shell aspect ratio $\eta=0.35$. Since both Prandtl numbers are larger than unity the discussion presented in Section 5.3.3 is relevant. As the mixing angle $\alpha$ is varied for a fixed value of $\tau$ the four convective regimes identified in the later section are evident from the ‘boat-shaped’ plots which appear similar to those in figure 7(b). We observe a region symmetric around $\alpha=0$ where instability is due to purely thermal eigenmodes since buoyancy has a predominantly thermal component. The region where buoyancy has a predominantly chemical component occupies a range of mixing angle values approximately between $\alpha=\pi/4$ and $\alpha=\pi/2$. Both purely thermal and purely compositional modes take the form of columnar convection featuring a cartridge belt of $z$-aligned columns centred at mid shell as expected at these relatively large values of the Prandtl numbers. The transition between the two regimes is characterised by large abrupt jumps in the values of the drift rate $omega$ and the wave number $m$ due to the relative difference in the values of $\mathrm{Pr}$ and $Sc$. The single-diffusive regimes are flanked by regions of very large scale convection appear approximately in the octants $\alpha\in[-\pi/2,-\pi/4]$ and $\alpha\in[\pi/2,3\pi/4]$ corresponding to diffusive and fingering instabilities, respectively. These regions are characterised by Rayleigh numbers of opposite signs. This shape and convective structure were discussed in Section 5.3.

Figure 9(b) shows explicitly the dependence of the critical parameter values for the onset of convection on the Coriolis number $\tau$ for several selected values of $\alpha$ sampling the four distinct convective regimes. In order to compare trends with the approximations (15) of Yano (1992) the scaled quantities $\widehat{}\mathrm{Ra}=\mathrm{Ra}\tau^{-4/3}$ and $\widehat{\omega}=\omega\tau^{-2/3}$ have been actually plotted. The curves for all values of $\alpha$ sampled in figure 9(b) which belong to the regions of purely-thermal or purely-chemical instability follow the scaling of approximations (15) namely, $\mathrm{Ra}\propto\tau^{4/3}$ and $\omega\propto\tau^{2/3}$ as demonstrated by the curves becoming horizontal for the larger values of $\tau$. Although the wave number $m$ is not scaled in this figure it also also follows a thermal-type scaling $m\propto\tau^{1/3}$ that we expect from equations (15). The purely For instance, a purely thermal case occurs when $\alpha=0$ (filled triangles) so that $\mathrm{Ra}=\mathrm{R_{t}}$ and in this case the particular curves in figure 9(b) correspond to the curves for $\Pr=1$ in figure 3 of (Ardes et al. 1997), although there is a slight discrepancy as $\eta=0.4$ in the latter work. Similarly, when $\alpha=\pi/2$ (filled circles) convection is purely chemical and if the values of $\mathrm{R_{c}}$, $\mathrm{Sc}$ are used the asymptotic trends for columnar convection cn be reproduced. A different scaling emerges in the two essentially double-diffusive regions sampled by points $\alpha=-3\pi/8$ and $\alpha=5\pi/8$. These are the regions previously described as arising from large discrepancies in Prandtl numbers. The most unstable azimuthal wave number in these regions seems to be independent of the Coriolis number, $\tau$, and the critical effective Rayleigh number is now $\mathrm{Ra}\propto\tau$. Note that, because these new regions obey a shallower scaling of $\mathrm{Ra}$ with $\tau$, as this parameter grows, those large scale regions will have a lower critical value of $\mathrm{Ra}$ when compared to purely thermal or chemical convection.

Significant deviations from the scaling discussed of equation (15) seem to occur ony for small values of the Coriolis number. This is expected for two reasons. Firstly, at small $\tau$ the end of the range of validity of the approximation (15) is approached, and secondly, inertial effects start to play a role when $\tau\mathrm{Pr}=O(10^{3})$, (Simitev and Busse 2003, Busse and Simitev 2004).

7 Dependence on the shell thickness

In this last section the dependence on the shell thickness is explored. The Earth’s inner core grows due to secular cooling and different convective regimes may occur at various moments in geological time.

Similarly to past sections, we keep all but two of the parameters at fixed values, in particular $\mathrm{Pr}=1$, $\mathrm{Sc}=100$ and $\tau=10^{4}$, while varying the mixing angle $\alpha$ and the parameter of interest, in this case the aspect ratio $\eta$. Te critical values of the Rayleigh number, the most unstable wave number and the drift rate are plotted as functions of the two variable parameters in figure 10. As the value of the mixing angle $\alpha$ is varied in figure 10(a) the four convective regimes familiar from the past two sections are observed at all values of $\eta$. The boundaries between these regimes vary little with $\eta$ – the boundary between the diffusive and the thermal columnar regime is located near $-\pi/4$, the boundary between the thermal columnar regime and the chemical columnar regime is located near $\pi/4$, and the boundary between the chemical columnar regime and the fingering regime is located near $\pi/2$, all of them sifting to smaller values with increasing $\eta$. In order to investigate the dependence on $\eta$ directly and to compare trends with the approximation (15), we plot in figure 10(b) the scaled quantities $\overline{\omega}=\omega(1-\eta)^{-2/3}$, $m$ and $\overline{}\mathrm{Ra}=\mathrm{Ra}(1-\eta)^{-7/3}$ as function of $\eta$. The approximation appears to hold well for thick spherical shells $\eta\in[0.1,0.5]$ but fails for thin shells where the effects of the curved boundaries start to dominate and columnar structures cannot easily form. The approximation for $\mathrm{Ra}$, naturally, also fails for values of $\alpha$ that are not in the columnar regimes. In these regimes when the system is stabilised by one quantity but destabilised by the other the most unstable wave number seems to be proportional to the mean radius and the effective Rayleigh number, to the square of the radius at 44% of the shell.

8 Discussion and conclusions

In this paper we identify deviations of the onset of doubly buoyant, doubly diffusive spherical rotating convection from the purely thermal case. We introduce a new representation for the thermal and compositional Rayleigh numbers that makes it possible to define a strictly positive pre-factor for the buoyancy. This pre-factor is the effective Rayleigh number $\mathrm{Ra}$ that, in this representation, has a unique critical value at fixed values of all other parameters. In contrast, when either of $\mathrm{R_{t}}$ or $\mathrm{R_{c}}$ serves as a critical parameter it may take negative values and the curve of marginal stability may be single, double- or triple-valued. In our new representation, the trade-off between thermal and compositional buoyancy is determined by a mixing angle $\alpha$ which takes positive values in the right half of the $\mathrm{R_{t}}-\mathrm{R_{c}}$ plane.

Four different regimes of convection are identified depending on the values of the mixing angle $\alpha$ and the Prandtl and Schmidt numbers. In the case when the values of the Prandtl number and the Schmidt number are equal convection behaves like in a single-diffusive, one-buoyancy system for all values of the mixing angle. In the $\mathrm{R_{t}}-\mathrm{R_{c}}$ plane the critical curve is a straight line with slope $-1$. In the $\mathrm{Ra}-\alpha$ plane the marginal curve behaves like $1/(\cos\alpha+\sin\alpha)$ and approaches vertical asymptotes at $-\pi/4$ and $3\pi/4$. When the values of the Prandtl number and the Schmidt number are slightly different from each other a fingering and a diffusive convection regimes emerge due to the different rates of diffusion of the two components of the fluid mixture. A small-deviations approximation allows us to gain some insight into the nature of these regimes. When the values of the Prandtl and the Schmidt numbers are significantly different from each other an additional transition between a thermal columnar and a thermal chemical regime occurs if the values of $\mathrm{Pr}$ are large. For small values of $\mathrm{Pr}$ both the thermal and the chemical modes are in an inertial regime without a clearly identifiable transition. In more detail, when both Prandtl numbers are below 0.1 very large scale inertial convection dominates independently of $\alpha$. When both Prandtl numbers are above 0.1, rotating convection takes place with a variety of wave numbers depending on $\alpha$. In regions of parameter space where $\mathrm{R_{t}}$ or $\mathrm{R_{c}}$ are negative, only very large scale convection seems to be possible. In terms of mixing angle these regions correspond to $\alpha<0$ or $\alpha>\pi/2$ and the modes seem related to the diffusive and fingering regimes of in thermohaline systems. Convection is not possible when the quantity with the smaller Prandtl number and a negative Rayleigh number dominates buoyancy. When the quantity with the smaller Prandtl number and a negative Rayleigh number does not dominate the buoyancy budget, convection is very large scale and has a negative drift rate (prograde). This structure of convection appears to be qualitatively preserved for all values of the Coriolis number $\tau$ and the shell thickness $\eta$ investigated.

Most of this study was carried out for and relatively large values of Prandtl and Schmidt numbers and for moderate values of the Coriolis number $\tau$. This choice offers a reasonable trade-off between computational expense and the possibility to compare the results to known trends for columnar convection (Busse 1970, Yano 1992). Both thermal and chemical columnar modes which are observed for $\tau\mathrm{Pr}>O(10^{3})$ appear to satisfy the trends of the approximation given by equation (15), as long as the regime of large scale convection with either $\mathrm{R_{t}}<0$, $\mathrm{R_{c}}>0$ or $\mathrm{R_{t}}>0$, $\mathrm{R_{c}}<0$ is not considered. In those regions, $\mathrm{Ra}\propto\tau$, the preferred azimuthal wave number $m$ shows little variation with respect to $\tau$ and the drift rates scales as $|\omega|\propto\tau^{2/3}$. The columnar scaling is also observed when varying $\eta$ but only for relatively thick shells $\eta\in[0.1,0.5]$ beyond which the approximation given by equation (15) starts to break down. Again this thermal scaling is only correct as long as the regime of large scale convection with either $\mathrm{R_{t}}<0$, $\mathrm{R_{c}}>0$ or $\mathrm{R_{t}}>0$, $\mathrm{R_{c}}<0$ is not considered.

Acknowledgements

This work was supported by the Leverhulme Trust [grant number RPG-2012-600].