2.1 Governing equation for the non-axisymmetric flow

We seek to model the system of equations (2.2)–(2.4) for small Ekman numbers. We choose the QG approximation to model the evolution of the velocity field (Or & Busse Reference Or and Busse1987; Brummell & Hart Reference Brummell and Hart1993; Cardin & Olson Reference Cardin and Olson1994; Plaut & Busse Reference Plaut and Busse2002; Aubert et al. Reference Aubert, Gillet and Cardin2003; Jones, Rotvig & Abdulrahman Reference Jones, Rotvig and Abdulrahman2003; Morin & Dormy Reference Morin and Dormy2004; Schaeffer & Cardin Reference Schaeffer and Cardin2005; Gillet & Jones Reference Gillet and Jones2006; Calkins et al. Reference Calkins, Aurnou, Eldredge and Julien2012; Teed, Jones & Hollerbach Reference Teed, Jones and Hollerbach2012). The QG approximation reduces the 3-D system to a 2-D system by taking advantage of the small variation of the flow along $z$ due to the rapid rotation (Gillet & Jones Reference Gillet and Jones2006). This approximation is only justified in the case of small slope of the boundaries, such as the Busse (Reference Busse1970) annulus. In the case of a sphere, the approximation is therefore not rigorously justified in any asymptotic limit. Consequently, the QG model is intended as a simplified model of convection in a rapidly rotating sphere that allows us to investigate unexplored regions of the parameter space. When possible, comparisons with theoretical, experimental and 3-D numerical models show that the QG model correctly reproduces key properties of the full system (Aubert et al. Reference Aubert, Gillet and Cardin2003; Morin & Dormy Reference Morin and Dormy2004; Gillet & Jones Reference Gillet and Jones2006; Gillet et al. Reference Gillet, Brito, Jault and Nataf2007; Plaut et al. Reference Plaut, Lebranchu, Simitev and Busse2008).

The QG model assumes that the fluid dynamics is dominated by the geostrophic balance, i.e. the Coriolis force balances the pressure gradient. The geostrophic velocity $\boldsymbol{u}^{g}$ is invariant along $z$ and $\boldsymbol{u}^{g}=(u_{s}^{g},u_{\unicode[STIX]{x1D719}}^{g},0)$ in cylindrical polar coordinates. By taking the $z$ -component of the curl of the momentum equation (2.2) and averaging it along $z$ , we obtain the equation for the axial vorticity, $\unicode[STIX]{x1D701}^{g}=(\unicode[STIX]{x1D735}\times \boldsymbol{u}^{g})\boldsymbol{\cdot }\boldsymbol{e}_{z}$ ,

(2.9) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}^{g}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}^{g}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D701}^{g}-\left(\frac{2}{Ek}+\unicode[STIX]{x1D701}^{g}\right)\left\langle \frac{\unicode[STIX]{x2202}u_{z}}{\unicode[STIX]{x2202}z}\right\rangle =\unicode[STIX]{x1D6FB}_{e}^{2}\unicode[STIX]{x1D701}^{g}-Ra\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\right\rangle , & & \displaystyle\end{eqnarray}$$

with

(2.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}_{e}^{2}A\equiv \frac{1}{s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(s\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}s}\right)+\frac{1}{s^{2}}\frac{\unicode[STIX]{x2202}^{2}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{2}}, & & \displaystyle\end{eqnarray}$$

and

(2.11) $$\begin{eqnarray}\displaystyle \langle A\rangle \equiv \frac{1}{2H}\int _{-H}^{+H}A\,\text{d}z, & & \displaystyle\end{eqnarray}$$

where $H=\sqrt{1-s^{2}}$ is the axial distance from the spherical boundary to the equatorial plane.

The axial velocity $u_{z}$ is assumed to be linear in $z$ . Three-dimensional linear convection models with internal heating show that the variation in $z$ of $u_{z}$ is not exactly linear, but it is a reasonable assumption because the variations in $z$ with respect to the linear profile are small compared with the variations in $s$ and $\unicode[STIX]{x1D719}$ (Gillet & Jones Reference Gillet and Jones2006).

The velocity $\boldsymbol{u}^{g}$ can be described by a streamfunction $\unicode[STIX]{x1D713}$ that describes the non-axisymmetric (i.e.  $\unicode[STIX]{x1D719}$ -dependent) geostrophic flow with the addition of an axisymmetric azimuthal flow,

(2.12) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{g}=\frac{1}{H}\unicode[STIX]{x1D735}\times (H\unicode[STIX]{x1D713}\boldsymbol{e}_{z})+\overline{u_{\unicode[STIX]{x1D719}}^{g}}\boldsymbol{e}_{\unicode[STIX]{x1D719}}, & & \displaystyle\end{eqnarray}$$

where

(2.13) $$\begin{eqnarray}\displaystyle \overline{A}\equiv \frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}A\,\text{d}\unicode[STIX]{x1D719}. & & \displaystyle\end{eqnarray}$$

This choice for the streamfunction implies that the divergence of $\boldsymbol{u}^{g}$ in the equatorial plane is non-zero because of the return axial flow due to the slope of the boundaries,

(2.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{e}\boldsymbol{\cdot }\boldsymbol{u}^{g}=-\unicode[STIX]{x1D6FD}u_{s}^{g}, & & \displaystyle\end{eqnarray}$$

where

(2.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}_{e}\boldsymbol{\cdot }\boldsymbol{A}\equiv \frac{1}{s}\frac{\unicode[STIX]{x2202}sA_{s}}{\unicode[STIX]{x2202}s}+\frac{1}{s}\frac{\unicode[STIX]{x2202}A_{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}, & & \displaystyle\end{eqnarray}$$

and

(2.16) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}=\frac{1}{H}\frac{\text{d}H}{\text{d}s}=-\frac{s}{H^{2}}. & & \displaystyle\end{eqnarray}$$

The third term on the left-hand side of (2.9) requires us to determine $u_{z}$ at the boundary $z=\pm H$ :

(2.17) $$\begin{eqnarray}\displaystyle u_{z}|_{\pm H}=\pm \frac{1}{H}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}|_{\pm H}\pm \unicode[STIX]{x1D6FD}Hu_{s}^{g}, & & \displaystyle\end{eqnarray}$$

where the normal vector at the boundary is $\boldsymbol{n}=\boldsymbol{e}_{r}$ . The normal component, $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}|_{\pm H}$ , is the Ekman pumping induced by the viscous boundary layer and is determined by asymptotic methods (Greenspan Reference Greenspan1968),

(2.18) $$\begin{eqnarray}\displaystyle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}|_{z=\pm H}=-\frac{Ek^{1/2}}{2}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \left(\frac{\boldsymbol{n}\times \boldsymbol{u}\pm \boldsymbol{u}}{\sqrt{|\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{ e}_{\boldsymbol{z}}|}}\right)=Ek^{1/2}P(s,u_{s}^{g},u_{\unicode[STIX]{x1D719}}^{g}). & & \displaystyle\end{eqnarray}$$

This analytical formulation is valid for a linear Ekman layer and for variations of the velocity that are slower than the rotation period and on larger length scales than $\mathit{O}(Ek^{1/2})$ . The function $P$ is derived for a spherical boundary in Schaeffer & Cardin (Reference Schaeffer and Cardin2005).

Our numerical code solves the equation for the evolution of the non-axisymmetric streamfunction

(2.19) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}{\mathcal{L}}\unicode[STIX]{x1D713}+(\boldsymbol{u}^{g}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{e}){\mathcal{L}}\unicode[STIX]{x1D713}-\left(\frac{2}{Ek}+{\mathcal{L}}\unicode[STIX]{x1D713}\right)\left(\frac{Ek^{1/2}P}{2H^{2}}+\unicode[STIX]{x1D6FD}u_{s}^{g}\right)=\unicode[STIX]{x1D6FB}_{e}^{2}{\mathcal{L}}\unicode[STIX]{x1D713}-Ra\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\right\rangle ,\qquad & & \displaystyle\end{eqnarray}$$

where

(2.20) $$\begin{eqnarray}\displaystyle {\mathcal{L}}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D701}^{g}=-\unicode[STIX]{x1D6FB}_{e}^{2}\unicode[STIX]{x1D713}+\frac{1}{s}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(\frac{s^{2}\unicode[STIX]{x1D713}}{H^{2}}\right). & & \displaystyle\end{eqnarray}$$

The no-slip and impenetrable boundary conditions imply that

(2.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}s}=0\quad \text{at }s=1. & & \displaystyle\end{eqnarray}$$

At the centre of the sphere, the solution must be regular. Using a decomposition in Fourier modes,

(2.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}(s,\unicode[STIX]{x1D719},t)=\mathop{\sum }_{m=1}^{+\infty }\hat{\unicode[STIX]{x1D713}}^{m}(s,t)\text{e}^{\text{i}m\unicode[STIX]{x1D719}}, & & \displaystyle\end{eqnarray}$$

the regularity condition is

(2.23) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D713}}^{m}=O(s^{m})\quad \text{at }s=0. & & \displaystyle\end{eqnarray}$$

2.2 Governing equation for the zonal flow

In our model, the streamfunction $\unicode[STIX]{x1D713}$ only describes the non-axisymmetric motions, so the axisymmetric azimuthal flows, or zonal flows, are treated separately. We take the $\unicode[STIX]{x1D719}$ - and $z$ -averages of the $\unicode[STIX]{x1D719}$ -component of the momentum equation to obtain

(2.24) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\overline{u_{\unicode[STIX]{x1D719}}^{g}}}{\unicode[STIX]{x2202}t}+\overline{u_{s}^{g}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D719}}^{g}}{\unicode[STIX]{x2202}s}}+\overline{\frac{u_{s}^{g}u_{\unicode[STIX]{x1D719}}^{g}}{s}}+\frac{2}{Ek}\langle \overline{u_{s}}\rangle =\unicode[STIX]{x1D6FB}^{2}\overline{u_{\unicode[STIX]{x1D719}}^{g}}-\frac{\overline{u_{\unicode[STIX]{x1D719}}^{g}}}{s^{2}}. & & \displaystyle\end{eqnarray}$$

Note that the geostrophic balance imposes that $\overline{u_{s}^{g}}=0$ . The fourth term on the left-hand side of (2.24) involves the $z$ -dependent radial velocity, which corresponds to the Ekman pumping term. Using the incompressibility of the fluid, it can be shown (Aubert et al. Reference Aubert, Gillet and Cardin2003) that

(2.25) $$\begin{eqnarray}\displaystyle \langle \overline{u_{s}}\rangle =\frac{Ek^{1/2}}{2H^{3/2}}\overline{u_{\unicode[STIX]{x1D719}}^{g}}. & & \displaystyle\end{eqnarray}$$

The no-slip boundary condition at the outer sphere implies

(2.26) $$\begin{eqnarray}\displaystyle \overline{u_{\unicode[STIX]{x1D719}}^{g}}=0\quad \text{at }s=1. & & \displaystyle\end{eqnarray}$$

By symmetry at the centre,

(2.27) $$\begin{eqnarray}\displaystyle \overline{u_{\unicode[STIX]{x1D719}}^{g}}=0\quad \text{at }s=0. & & \displaystyle\end{eqnarray}$$

2.3 Governing equation for the temperature

The dimensionless equation for the evolution of the temperature perturbation in three dimensions is

(2.28) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}^{3d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E9}=\frac{1}{Pr}(2ru_{r}^{3d}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6E9}), & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}^{3d}$ is the velocity in three dimensions. In cylindrical polar coordinates,

(2.29) $$\begin{eqnarray}\displaystyle \boldsymbol{u}^{3d}=(u_{s}^{g},u_{\unicode[STIX]{x1D719}}^{g},Ek^{1/2}zP+\unicode[STIX]{x1D6FD}zu_{s}^{g}). & & \displaystyle\end{eqnarray}$$

The temperature is fixed at the outer boundary so

(2.30) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E9}=0\quad \text{at }r=1. & & \displaystyle\end{eqnarray}$$

At the centre, the heat flux of the spherically averaged temperature perturbation is zero,

(2.31) $$\begin{eqnarray}\displaystyle \left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}\right\rangle _{{\mathcal{S}}}=\int _{{\mathcal{S}}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}r}\,\text{d}{\mathcal{S}}=0\quad \text{at }r=0, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{S}}$ is a spherical surface. By symmetry, the non-spherically symmetric components of $\unicode[STIX]{x1D6E9}$ are zero,

(2.32) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E9}-\langle \unicode[STIX]{x1D6E9}\rangle _{{\mathcal{S}}}=0\quad \text{at }r=0. & & \displaystyle\end{eqnarray}$$

2.4 Numerical method

In summary, the governing equations of our hybrid QG–3-D model are

(2.33) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}{\mathcal{L}}\unicode[STIX]{x1D713}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{e}){\mathcal{L}}\unicode[STIX]{x1D713}-\left(\frac{2}{Ek}+{\mathcal{L}}\unicode[STIX]{x1D713}\right)\left(\frac{Ek^{1/2}P}{2H^{2}}+\unicode[STIX]{x1D6FD}u_{s}\right)=\unicode[STIX]{x1D6FB}_{e}^{2}{\mathcal{L}}\unicode[STIX]{x1D713}-Ra\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\right\rangle ,\qquad & \displaystyle\end{eqnarray}$$

(2.34) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\overline{u_{\unicode[STIX]{x1D719}}}}{\unicode[STIX]{x2202}t}+\overline{u_{s}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}s}}+\overline{\frac{u_{s}u_{\unicode[STIX]{x1D719}}}{s}}+\frac{\overline{u_{\unicode[STIX]{x1D719}}}}{Ek^{1/2}H^{3/2}}=\unicode[STIX]{x1D6FB}^{2}\overline{u_{\unicode[STIX]{x1D719}}}-\frac{\overline{u_{\unicode[STIX]{x1D719}}}}{s^{2}}, & \displaystyle\end{eqnarray}$$

(2.35) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}^{3d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E9}=\frac{1}{Pr}(2ru_{r}^{3d}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6E9}), & \displaystyle\end{eqnarray}$$

$g$

Equations (2.33) and (2.34) are solved on a 2-D grid in the equatorial plane. A second-order finite difference scheme is implemented in radius with irregular spacing (finer near the outer boundary). In the azimuthal direction, the streamfunction and zonal velocity are expanded in Fourier modes. Equation (2.35) is solved on a 3-D grid. Similarly to the 2-D grid, a finite difference scheme is used in radius. The temperature is expanded in spherical harmonics $Y_{l}^{m}$ in the angular coordinates with $l$ representing the latitudinal degree and $m$ the azimuthal mode. This allows for a simple implementation of the boundary conditions at the centre. A Crank–Nicolson scheme is implemented for the time integration of the diffusion terms and an Adams–Bashforth procedure is used for the other terms. Our hybrid approach is numerically efficient compared with a fully 3-D approach when considering small Prandtl numbers, i.e. when the diffusion of the temperature happens on shorter time scales than the viscous dissipation. In this case, the 3-D grid of the temperature requires less radial points and azimuthal modes than the 2-D grid of the velocity, so smaller grid steps and time step can be used for the velocity. The most demanding nonlinear simulations were run at $Ek=10^{-8}$ . The numerical resolution required to resolve the simulation run at $Ek=10^{-8}$ and $Pr=10^{-1}$ for the largest Rayleigh number ( $Ra\approx 6Ra_{c}$ ) is $1600$ radial points, $300$ Fourier modes $m$ for the 2-D velocity grid and $800$ radial points, $200$ degrees $l$ and $200$ modes $m$ of the spherical harmonics for the 3-D temperature grid.

The buoyancy term in (2.33) requires the integration of the temperature along $z$ for each point of coordinates $(s,\unicode[STIX]{x1D719})$ on the 2-D grid at each time step. To do so, we use a linear interpolation in $r$ and $\unicode[STIX]{x1D703}$ from the 3-D spherical grid onto a 3-D cylindrical projection of the velocity grid. The spectral decomposition in $\unicode[STIX]{x1D719}$ for the velocity and temperature allows a straightforward interpolation in spectral space. To calculate the 3-D velocity used for the advection of temperature, we use a linear interpolation in $s$ . Again, the interpolation in $\unicode[STIX]{x1D719}$ is straightforward in spectral space.

In order to benchmark the numerical model, we performed calculations with the linearised version of the code to determine the critical parameters at the linear onset of convection. The linear results are given in appendix A and compared with results from the asymptotic theory and from published numerical studies using either a 3-D or a QG model. The comparison shows that our numerical results tend to the asymptotic values as the Ekman number decreases. In the following, the critical Rayleigh number at the linear onset of convection calculated with our numerical code is denoted $Ra_{c}$ .

2.5 Definition of the output parameters

To quantify the global properties of convection, we use the following quantities calculated from the output of the simulations. In our dimensionless units, the Reynolds number corresponds to the root mean square (r.m.s.) value of the velocity,

(2.36) $$\begin{eqnarray}\displaystyle Re^{\ast }(t)=\left(\frac{3}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{1}(u_{s}^{2}(s,\unicode[STIX]{x1D719},t)+u_{\unicode[STIX]{x1D719}}^{2}(s,\unicode[STIX]{x1D719},t))2H(s)s\,\text{d}s\,\text{d}\unicode[STIX]{x1D719}\right)^{1/2}, & & \displaystyle\end{eqnarray}$$

and

(2.37) $$\begin{eqnarray}\displaystyle Re=\frac{1}{\unicode[STIX]{x0394}t}\int _{\unicode[STIX]{x0394}t}Re^{\ast }(t)\,\text{d}t, & & \displaystyle\end{eqnarray}$$

where $u_{\unicode[STIX]{x1D719}}$ includes the zonal velocity. We measure the convective Reynolds number, $Re_{c}^{\ast }$ , as in (2.36) but including only the non-axisymmetric velocity. Its time average is denoted $Re_{c}$ . Similarly, we measure the zonal Reynolds number, $Re_{0}^{\ast }$ , including only the axisymmetric velocity, and its time average is denoted $Re_{0}$ . The time averages are calculated over a long period of time after the system reaches a statistical equilibrium.

We define the radial profile of the time-averaged zonal velocity as

(2.38) $$\begin{eqnarray}\displaystyle U_{0}(s)=\frac{1}{\unicode[STIX]{x0394}t}\int _{\unicode[STIX]{x0394}t}\overline{u_{\unicode[STIX]{x1D719}}}(s,t)\,\text{d}t, & & \displaystyle\end{eqnarray}$$

and the radial profile of the r.m.s. radial velocity,

(2.39) $$\begin{eqnarray}\displaystyle U_{s}(s)=\frac{1}{\unicode[STIX]{x0394}t}\int _{\unicode[STIX]{x0394}t}\left(\frac{1}{2\unicode[STIX]{x03C0}s}\int _{0}^{2\unicode[STIX]{x03C0}}u_{s}^{2}(s,\unicode[STIX]{x1D719},t)s\,\text{d}\unicode[STIX]{x1D719}\right)^{1/2}\,\text{d}t. & & \displaystyle\end{eqnarray}$$

The Nusselt number, which is the ratio of the total heat flux, $Q_{t}$ , to the conductive heat flux $Q_{\unicode[STIX]{x1D705}}$ , is often used to measure the efficiency of the convective heat transport. In a convective system with internal heating and fixed temperature boundaries, the mean total outward heat flux must balance heat production, so $Q_{t}=r_{o}S/3$ . To estimate the heat flux that is conducted along the spherically averaged temperature gradient, we define the equivalent internal heating, $S_{eq}$ , that would correspond to the purely conductive state with a temperature at the centre $T_{c}=T_{s}(r=0)+\langle \unicode[STIX]{x1D6E9}(r=0)\rangle _{{\mathcal{S}}}$ , where $\langle \unicode[STIX]{x1D6E9}(r=0)\rangle _{{\mathcal{S}}}$ is an output of the simulation. Assuming that $Q_{\unicode[STIX]{x1D705}}=r_{o}S_{eq}/3$ , we can define the Nusselt number, $Nu$ , as

(2.40) $$\begin{eqnarray}\displaystyle Nu=\frac{Q_{t}}{Q_{\unicode[STIX]{x1D705}}}=\frac{S}{S_{eq}}=\frac{1}{PrT_{c}}, & & \displaystyle\end{eqnarray}$$

where $S_{eq}=PrT_{c}S$ and $T_{c}=Pr^{-1}+\langle \unicode[STIX]{x1D6E9}(r=0)\rangle _{{\mathcal{S}}}$ . The decrease of the temperature at the centre due to convection is thus a good measure of the efficiency of the convective heat transport (Goluskin Reference Goluskin2016). We find that longer time integrations are required to ensure the saturation of $Nu$ than for the saturation of the kinetic energy.

We also measure the r.m.s. value of the axisymmetric temperature perturbation as

(2.41) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E9}_{0}^{\ast }(t)=\left(\frac{3}{2}\int _{0}^{1}\int _{0}^{\unicode[STIX]{x03C0}}\overline{\unicode[STIX]{x1D6E9}}^{2}r^{2}\sin \unicode[STIX]{x1D703}\,\text{d}r\,\text{d}\unicode[STIX]{x1D703}\right)^{1/2}. & & \displaystyle\end{eqnarray}$$

3.1 Supercritical convection for $Ek=10^{-5}$

We begin our study with simulations performed at $Ek=10^{-5}$ and $Pr=10^{-1}$ in order to link the results obtained with our hybrid QG–3-D model to previous studies using either 3-D or QG models.

Figure 2.

Snapshots of the streamfunction for $Ek=10^{-5}$ and $Pr=10^{-1}$ for (a) $Ra=1.13Ra_{c}$ , (b) $Ra=1.19Ra_{c}$ and (c) $Ra=7.06Ra_{c}$ . The time-averaged radial profiles of the zonal velocity, $U_{0}$ (solid line), and the r.m.s. radial velocity, $U_{s}$ (dashed line), are plotted according to the axis on the right-hand side.

For this Ekman number, the bifurcation from the basic conduction state is supercritical at $Ra=Ra_{c}$ . For values of $Ra$ just above $Ra_{c}$ , the convective flow consists of a thermal Rossby wave, a well-known flow pattern in rotating convection (e.g. Busse Reference Busse1970). Figure 2(a) shows a snapshot of the streamfunction (colour) for $Ra=1.13Ra_{c}$ . The linear thermal Rossby wave develops around the radius $s\approx 0.5$ and propagates eastward. For small Prandtl numbers, the thermal Rossby wave has a large tilt in the equatorial cross-section (Zhang Reference Zhang1992a ).

Figure 3.

Evolution of the global properties of convection for $Ek=10^{-5}$ and $Pr=10^{-1}$ . The vertical and horizontal error bars represent the standard deviation in the time series.

Figure 3(a) shows the evolution of $Re$ as a function of the normalised Rayleigh number, $Ra/Ra_{c}-1$ . For $Ra/Ra_{c}<1.2$ , the evolution of $Re$ with $Ra/Ra_{c}-1$ can be approximated by a power law $Re\sim (Ra/Ra_{c}-1)^{0.42}$ . A sudden transition occurs at $Ra/Ra_{c}=1.2$ , after which the slope of the curve becomes steeper. Around this transition, or kink, the convection changes from quasi-steady to time dependent. The sequence of bifurcations that follows the kink has been previously documented in both 3-D models (Zhang Reference Zhang1992a ; Sun, Schubert & Glatzmaier Reference Sun, Schubert and Glatzmaier1993; Tilgner & Busse Reference Tilgner and Busse1997; Grote & Busse Reference Grote and Busse2001; Simitev & Busse Reference Simitev and Busse2003) and QG models (Schnaubelt & Busse Reference Schnaubelt and Busse1992; Morin & Dormy Reference Morin and Dormy2004). Our hybrid QG model reproduces this sequence: periodic vacillations ( $Ra=1.18Ra_{c}$ ), chaotic fluctuations with a localisation of the convection ( $Ra=1.19Ra_{c}$ , figure 2 b) and bursts of convection ( $Ra=1.21Ra_{c}$ ). These various states are attributed to the interactions of the zonal flow and the axisymmetric temperature with the convective columns (Teed et al. Reference Teed, Jones and Hollerbach2012). Interestingly, Morin & Dormy (Reference Morin and Dormy2004) find that the kink occurs closer to the onset of convection as the Ekman number decreases, implying that quasi-steady convection in the form of thermal Rossby waves is expected to occur within a vanishingly small range of Rayleigh numbers for $Ek\rightarrow 0$ . For higher Rayleigh number, the time dependence of the convection becomes irregular. Figure 2(c) shows the streamfunction for $Ra=7.06Ra_{c}$ . The thermal Rossby wave is no longer visible and the velocity comprises a wide range of length scales.

The kink at $Ra=1.2Ra_{c}$ is also clearly visible in the plot of $Nu$ as a function of $Ra/Ra_{c}-1$ shown in figure 3(b). Near the onset, Busse & Or (Reference Busse and Or1986) predicted that the evolution of the convective heat transport with the Rayleigh number follows a power law $Nu\sim Ra/Ra_{c}-1$ , which has later been confirmed by Gillet & Jones (Reference Gillet and Jones2006) with numerical results for $Pr>1$ . For $Ra<1.2Ra_{c}$ , our numerical results can be approximated by a power law of weaker exponent, $0.88$ . The slow increase of $Nu$ with the Rayleigh number for $Pr<1$ was previously observed by Gillet & Jones (Reference Gillet and Jones2006). After the kink, we have too few numerical points to determine whether the points still follow a power law.

In figure 2(a), we plot the time-averaged radial profile of the zonal velocity, $U_{0}$ , for $Ra=1.13Ra_{c}$ . Zonal flows are produced from the nonlinear interactions of the non-axisymmetric velocity components (2.34). The production of zonal flow by thermal Rossby waves has been extensively studied in the literature (e.g. Busse & Hood Reference Busse and Hood1982; Plaut et al. Reference Plaut, Lebranchu, Simitev and Busse2008). In a spherical geometry, the divergence of the Reynolds stresses produces a retrograde zonal jet in the central part and a prograde jet in the outer part. The correlation of the velocity components is caused by the tilt of the thermal Rossby waves. Figure 2(c) shows $U_{0}$ for $Ra=7.06Ra_{c}$ . The zonal flow still has a double jet structure with a prograde jet on the outer part. The peak velocity of the zonal flow is three times larger than the maximum of the radial profile of the r.m.s. radial velocity, $U_{s}$ , also shown in the figure.

Weakly nonlinear analysis predicts that the amplitude of the zonal flow scales as the square of the amplitude of the convective velocity (Busse & Hood Reference Busse and Hood1982). Figure 3(c) shows the ratio $Re_{0}/Re_{c}$ as a function of the convective velocity $Re_{c}$ . This ratio varies linearly with $Re_{c}$ just above the onset, in agreement with the weakly nonlinear analysis. After the transition from quasi-steady convection to time-dependent convection, the evolution of $Re_{0}/Re_{c}$ with $Re_{c}$ departs from the weakly nonlinear prediction and the slope becomes shallower. The transition occurs when $Re_{0}$ becomes comparable to $Re_{c}$ , in agreement with the idea that the bifurcation to time-dependent convection is due to the presence of zonal flows.

3.2 Subcritical convection for $E=10^{-8}$

We now present results obtained at the small Ekman numbers that are out of reach of earlier and recent 3-D models (e.g. Yadav et al. Reference Yadav, Gastine, Christensen, Duarte and Reiners2016).

Figure 4.

(a) Time series of $Re^{\ast }$ (thick black line) and $Re_{0}^{\ast }$ (red) for $Ek=10^{-8}$ , $Pr=10^{-2}$ and $Ra=1.06Ra_{c}$ . (b,c) Streamfunction (colour, snapshot) and radial profiles of the zonal velocity (snapshot, solid line, plotted according to the right-hand side axis) during (b) the growing phase at $t=1.36\times 10^{-3}$ and (c) the saturated phase at $t=1.52\times 10^{-3}$ .

Figure 4(a) shows the time series of the Reynolds number, $Re^{\ast }$ , and the zonal Reynolds number, $Re_{0}^{\ast }$ , for $Ek=10^{-8}$ , $Pr=10^{-2}$ and $Ra=1.06Ra_{c}$ . After a small temperature perturbation is added at $t=0$ , the kinetic energy grows exponentially and eventually saturates around a mean value. Despite the proximity to the onset of convection, the values of the Reynolds numbers are large in the saturated phase, on average $Re\approx 15200$ and $Re_{0}\approx 9600$ . Figure 4(b,c) show snapshots of the streamfunction during the growing phase and the saturated phase respectively. During the growth, the flow has the distinctive pattern of a thermal Rossby wave. The azimuthal modulation of the non-axisymmetric pattern is caused by the interaction of Rossby waves of different azimuthal wavenumbers. A snapshot of the zonal flow is also shown as a solid black line. The zonal flow has the familiar double jet structure seen at the onset of convection at $Ek=10^{-5}$ (figure 2 a). During the saturated phase, the flow is markedly different. The flow is now organised into two regions in the equatorial plane. In the inner region ( $s<0.6$ ), the pattern of the thermal Rossby waves is lost and the convective flow is vigorous with a wide range of length scales. In the outer part ( $s>0.6$ ), the velocity has a smaller amplitude and the characteristic elongated and tilted pattern of the Rossby waves is visible. This type of zonation of the convection was previously described by Sumita & Olson (Reference Sumita and Olson2000), Aubert et al. (Reference Aubert, Gillet and Cardin2003), Miyagoshi, Kageyama & Sato (Reference Miyagoshi, Kageyama and Sato2010) for Rayleigh numbers several times above critical. In the outer part, the vortex stretching term due to the slope of the boundaries is predominant over the other forces, so Rossby waves are easily excited there. The zonal flow is also modified: its radial structure is more complex with three main jets of similar radial widths, where a jet is defined by a zonal velocity of same sign. The innermost and outermost jets remain prograde while the middle jet is retrograde. The outermost prograde zonal jet is driven by the Reynolds stresses from the tilted thermal Rossby waves in the outer region.

Figure 5.

(a) $Re$ as a function of $Ra/Ra_{c}$ for $Ek=10^{-8}$ and $Pr=10^{-2}$ . The limit between the initial value of $Re^{\ast }$ leading to either the decay or the growth of the solution is plotted as triangles. (b) Same as (a) for all the cases of subcritical convection ( $Ek=10^{-8}$ , $Pr\in [10^{-2},10^{-1}]$ and $Ek=10^{-7}$ , $Pr=10^{-2}$ ).

Figure 5(a) shows the values of $Re$ when the Rayleigh number is varied around the onset of convection for $Ek=10^{-8}$ and $Pr=10^{-2}$ . As previously noted in figure 4(a), $Re$ becomes large (of the order of $10^{4}$ ) when $Ra$ is only just slightly increased above $Ra_{c}$ . Remarkably, when the Rayleigh number is decreased during a simulation from $Ra\geqslant Ra_{c}$ to $Ra<Ra_{c}$ , the convection does not decay and remains vigorous up to a value of the Rayleigh number, $Ra_{s}$ , where the convection abruptly shuts down. At $Ra=Ra_{s}$ , the value of $Re$ is discontinuous: for Rayleigh numbers above $Ra_{s}$ , we have $Re>5000$ , whereas for all $Ra$ below $Ra_{s}$ , $Re=0$ . This behaviour ( $Ra_{s}<Ra_{c}$ ) is typical of a subcritical onset of nonlinear convection. For $Ek=10^{-8}$ and $Pr=10^{-2}$ , we find that $Ra_{s}/Ra_{c}=0.60$ . Figure 6(a,b) show snapshots of the streamfunction for $Pr=10^{-2}$ and two Rayleigh numbers on the subcritical branch, $Ra/Ra_{c}=0.60$ and $0.83$ . For $Ra=Ra_{s}=0.60Ra_{c}$ , the convective flow is similar to the saturated case described previously in figure 4(c). The main difference is that the amplitude of the flow in the inner convective region is concentrated in a localised structure, that drifts in the azimuthal direction with time. The localisation of the convection might indicate that the system struggles to maintain convection for this Rayleigh number. However the kinetic energy of the convective and zonal flows do not decay during the entire time integration of the simulation, which corresponds to $200$ convective turnover time scales (where one turnover time scale is calculated as $1/Re_{c}$ ), or equivalently, $2$ global thermal diffusion time scales. For larger Rayleigh numbers, the localisation in the inner convective region disappears and the convection fills all longitudes, as observed for $Ra=0.83Ra_{c}$ in figure 6(b). The radial profile $U_{s}$ is plotted in figure 6(a,b). The r.m.s. radial velocity peaks in the inner convective region, and then, monotonically decreases in the outer Rossby wave region. The boundary between the convective inner region and the outer Rossby wave region moves outwards when $Ra$ increases. The radial profile of the zonal flow, $U_{0}$ , is also plotted in figure 6(a,b). As described in the saturated case of figure 4(c), the zonal flow displays three jets in both cases. The peak of the middle retrograde jet roughly coincides with the maximum of $U_{s}$ . The zonal flows represent a significant portion of the total flow with the amplitude of $Re_{0}$ increasing from $0.5Re_{c}$ to $1.3Re_{c}$ when $Ra/Ra_{c}$ varies from $0.60$ to $3$ . No significant change in the flow is observed when the Rayleigh number is varied across the linear critical value $Ra_{c}$ , and all of the output parameters evolve continuously, similarly to $Re_{c}$ in figure 5(a). No continuous branch of convection is found at $Ra=Ra_{c}$ when the convection is started from a small perturbation as described for $Ek=10^{-5}$ .

Figure 6.

Streamfunction (colour, snapshot) and radial profiles of $U_{s}$ (dashed line) and $U_{0}$ (solid line) (time averages) for $Ek=10^{-8}$ , $Pr=10^{-2}$ : (a) $Ra=Ra_{s}=0.60Ra_{c}$ and (b) $Ra=0.83Ra_{c}$ .

Within the interval $Ra_{s}\leqslant Ra<Ra_{c}$ , the system has two stable solutions: one is located on the subcritical branch shown in red in figure 5(a) and is found by using an initial condition with a large enough amplitude; the other solution corresponds to $Re=0$ and is found by starting the simulation from a small perturbation. In order to quantify how large the finite amplitude of the initial condition must be to access the subcritical branch, we perform a series of simulations where the initial condition is taken from a snapshot of the saturated solution on the subcritical branch for given parameters and we initially divide the amplitude of the velocity by a given factor. By a trial and error procedure, we determine the factor that separates a solution where the kinetic energy decays to zero from a solution where the kinetic energy grows back to the value on the subcritical branch. The limit between the values of $Re^{\ast }$ of the initial condition leading to either the decaying solution or the growing solution are shown in figure 5(a). As expected, we find that for values of $Ra$ just above $Ra_{s}$ , the finite amplitude of the initial condition must be much larger than for values of $Ra$ just below  $Ra_{c}$ .

By varying the Prandtl number between $[10^{-2},10^{-1}]$ and the Ekman number between $[10^{-8},10^{-5}]$ , we find that subcritical convection occurs for $Ek=10^{-8}$ at $Pr\leqslant 10^{-1}$ and $Ek=10^{-7}$ at $Pr\leqslant 10^{-2}$ . Table 1 gives the values of $Ra_{s}/Ra_{c}$ for the subcritical cases and figure 5(b) shows $Re$ as a function of $Ra/Ra_{c}$ for these cases. Note that long time integrations are required to determine $Ra_{s}$ . The values of $Ra_{s}$ are obtained from simulations where the convection does not decay for the entire duration of the calculation, which is at least $200$ convective turnover time scales. We find that the value of $Re$ at $Ra=Ra_{s}$ is of the order of or larger than $10^{3}$ . The ratio $Ra_{s}/Ra_{c}$ decreases when $Ek$ decreases, and also when $Pr$ decreases. Consequently, our results suggest that the subcritical behaviour is amplified by low Ekman and Prandtl numbers.

Table 1.

Values of the ratio of the Rayleigh number at the nonlinear onset of convection to the critical Rayleigh number at the linear onset in the subcritical cases.

The subcritical behaviour of convection in internal heating models in the limit of small $Ek$ was anticipated from the work of Soward (Reference Soward1977), who found weakly nonlinear solutions close to the critical Rayleigh number at the linear onset predicted by the local theory of Busse (Reference Busse1970), which is significantly smaller than the true global value (Jones et al. Reference Jones, Soward and Mussa2000). These nonlinear solutions are maintained in the subcritical domain because the zonal shear produced by nonlinear effects opposes the process of phase mixing, which is responsible for the decay of the linear solution. Subcritical convection for small Ekman numbers was also predicted by Plaut et al. (Reference Plaut, Lebranchu, Simitev and Busse2008) from a weakly nonlinear analysis of a quasi-geostrophic model with internal heating. Plaut et al. found that the bifurcation to the thermal Rossby waves changes from supercritical to subcritical for $Ek\leqslant 6.7\times 10^{-7}$ (in our scaling) at $Pr=1$ . Their analysis demonstrates that the nonlinearities in the temperature equation are responsible for the transition. These nonlinearities are produced by the interactions of the zonal flow with the temperature perturbation of the wave and by the interactions of the axisymmetric temperature with the velocity of the wave. The difference between our numerical results and the weakly nonlinear analysis of Plaut et al. (Reference Plaut, Lebranchu, Simitev and Busse2008), which predicts subcriticality for larger $Ek$ at $Pr=1$ , is possibly due to their treatment of the temperature in two dimensions. As far as we know, the numerical results of the present paper are the first to demonstrate the existence of subcritical rotating convection in a fully nonlinear model in spherical geometry. The requirements of small Ekman numbers and small Prandtl numbers (at least for $Ek\in [10^{-8},10^{-7}]$ ) can explain why subcriticality has not been observed in previous numerical models (e.g. by Morin & Dormy (Reference Morin and Dormy2004) with QG simulations at $Ek=2\times 10^{-7}$ and $Pr=1$ ). Although subcriticality in the fully 3-D system has been predicted from the work of Soward (Reference Soward1977) for small $Ek$ , its presence in our results might yet be a consequence of our quasi-geostrophic approximation. Future 3-D simulations will be needed to determine whether this is the case.

3.3 A hysteresis loop for $E=10^{-6}$ and $Pr=10^{-2}$

For $Ek=10^{-6}$ , the bifurcation at the onset of convection is supercritical for all the Prandtl numbers studied here ( $Pr\in [10^{-2},10^{-1}]$ ). For $Pr=10^{-2}$ , we observe an interesting behaviour near the onset: two stable branches of convection co-exist within a limited range of Rayleigh numbers. This behaviour can be seen in the evolution of $Re$ and $Nu-1$ with $Ra/Ra_{c}-1$ in figure 7. On the upper branch, $Re$ is one order of magnitude larger than on the lower branch. We thus refer to the upper branch as strong and the lower branch as weak. The two branches are stable and co-exist for Rayleigh numbers between $1.09\leqslant Ra/Ra_{c}\leqslant 1.27$ . The co-existence of two stable branches for a range of $Ra$ leads to a hysteresis loop. The convection selects the weak or the strong branch depending on the dynamical history of the system. For instance, the system remains on the strong branch if the initial condition is taken from a simulation with $Ra>1.27Ra_{c}$ and $Ra$ is then decreased down to $1.09Ra_{c}$ . For $Ra>1.27Ra_{c}$ , only the strong branch exists according to our calculations. When the system is located on the weak branch, an increase of the Rayleigh number above $1.27Ra_{c}$ therefore leads the system to jump to the strong branch.

Figure 7.

Evolution of (a) $Re$ and (b) $Nu-1$ as a function of the Rayleigh number for $E=10^{-6}$ and $Pr=10^{-2}$ .

In order to quantify the amplitude of the velocity required to jump from one branch to the other, we follow the same numerical procedure as described in § 3.2. As an initial condition, we use a snapshot of a solution on the weak branch that we multiply by a given factor and determine the initial amplitude of the flow leading to either a return to the weak branch or a growth to the strong branch. The triangles in figure 7(a) indicate the limit between the initial value of the r.m.s. velocity leading to the growth or the decay of the energy. As expected a larger amplitude of the initial velocity is required to jump to the strong branch for $Ra=1.09Ra_{c}$ than $Ra=1.27Ra_{c}$ . In fact, for $Ra=1.27Ra_{c}$ , we observe that, when the system is on the weak branch, a fluctuation of the kinetic energy can lead the system to jump to the strong branch regime. Figure 8 shows this jump in a time series of the Reynolds number, $Re^{\ast }$ . The two branches are well separated as there is no slow variation on the weak branch towards the strong branch, but only a sudden jump that lasts less than $3$ convective turnover time scales. The ratio of zonal to total Reynolds numbers (plotted according to the right axis) is approximately $0.45$ on the weak branch and increases to approximately $0.8$ after the jump to the strong branch.

Figure 8.

Time series of $Re^{\ast }$ (thick black line, left axis) and $Re_{0}^{\ast }/Re^{\ast }$ (red, right axis) for a simulation near the edge of the hysteresis loop $E=10^{-6}$ , $Pr=10^{-2}$ and $Ra=1.27Ra_{c}$ .

Figure 9.

Streamfunction (colour), $U_{0}$ (solid line) and $U_{s}$ (dashed line) for $Ek=10^{-6}$ , $Pr=10^{-2}$ and $Ra=1.19Ra_{c}$ on (a) the weak branch and (b) the strong branch.

Typical convective flows on the weak and strong branches are shown in figure 9(a,b) for the same Rayleigh number $Ra=1.19Ra_{c}$ . On the weak branch, the flow consists in thermal Rossby waves. The zonal flow has a double jet structure and its amplitude is comparable to the amplitude of the radial flow. The characteristics of this flow are similar to the case described near the onset of convection for $Ek=10^{-5}$ (§ 3.1). On the strong branch, the thermal Rossby waves are not visible anymore and a wide range of length scales are observed. The zonal flow has a larger amplitude than the radial velocity. Not only do the solutions look different on the two branches, the global properties of the convection also follow different scaling laws with the input parameters. The evolution of $Nu$ with the Rayleigh number (figure 7 b) is less steep on the strong branch compared with the weak branch. On the weak branch we find that the points approximately follow a power law $Nu\sim (Ra/Ra_{c}-1)^{\unicode[STIX]{x1D6FC}}$ with $\unicode[STIX]{x1D6FC}=1$ , whereas on the strong branch $\unicode[STIX]{x1D6FC}=0.6$ .

3.4 Nonlinear oscillations for $E=10^{-7}$ and $Pr=10^{-1}$

For $E=10^{-7}$ and $Pr=10^{-1}$ , we observe an interesting behaviour between the linear onset of convection and the saturated branch of convection, for values of $Ra$ between $Ra_{c}$ and $1.025Ra_{c}$ . This behaviour is best observed in the time series of $Re^{\ast }$ shown in figure 10 for Rayleigh numbers in this interval. $Re^{\ast }$ first grows at a rate that depends on $Ra$ and then reaches a value of approximately $400$ . After a few fluctuations around this mean value, $Re^{\ast }$ decreases rapidly. The system subsequently enters a new growing phase and the oscillation repeats itself. The duration of the oscillation depends on the growth rate, and so, larger Rayleigh numbers have shorter period. For $Ra=1.025Ra_{c}$ , the system eventually stays on the saturated branch after a few oscillations, after which it never goes back to the oscillating phase.

Figure 10.

Time series of the Reynolds number, $Re^{\ast }$ , for $E=10^{-7}$ and $Pr=10^{-1}$ and Rayleigh numbers near the onset.

Figure 11.

Phase diagram of (a) $Re_{0}^{\ast }/Re_{c}^{\ast }$ versus $Re_{c}^{\ast }$ and (b) $\unicode[STIX]{x1D6E9}_{0}^{\ast }$ versus $Re_{c}^{\ast }$ for the same time series as figure 10.

This behaviour is similar to nonlinear oscillations such as relaxation oscillations found in many dynamical systems (Fauve Reference Fauve, Godrèche and Manneville1998). These oscillations are quasi-periodic and may be better understood using phase diagrams. Figure 11(a) shows the phase diagram of the ratio of the zonal to convective Reynolds numbers, $Re_{0}^{\ast }/Re_{c}^{\ast }$ , versus the convective Reynolds number, $Re_{c}^{\ast }$ , for the same time series as figure 10. Figure 11(b) shows the phase diagram for the r.m.s. axisymmetric temperature, $\unicode[STIX]{x1D6E9}_{0}^{\ast }$ , versus $Re_{c}^{\ast }$ . During the exponential growing phase $Re_{c}^{\ast }$ (labelled (A)), all trajectories follow the same path and $Re_{0}^{\ast }$ increases as the square of $Re_{c}^{\ast }$ , in agreement with weakly nonlinear analysis near the onset of convection (see § 3.1). The system reaches a maximum (B), where $Re_{0}^{\ast }$ has a similar amplitude to $Re_{c}^{\ast }$ . A dissipative phase (C) then follows during which $Re_{0}^{\ast }$ and $Re_{c}^{\ast }$ initially decay at the same rate. The three quantities, $Re_{c}^{\ast }$ , $Re_{0}^{\ast }$ and $\unicode[STIX]{x1D6E9}_{0}^{\ast }$ , all start to decay at the same time. While $Re_{c}^{\ast }$ and $Re_{0}^{\ast }$ enter together the growing phase (A), the decay of $\unicode[STIX]{x1D6E9}_{0}^{\ast }$ persists for longer. The three phases (A, B, C) are fairly similar for all the Rayleigh numbers, although the system can follow different trajectories during the dissipative phase. For the largest Rayleigh number ( $Ra=1.025Ra_{c}$ ), the system exits the oscillating regime to reach a stagnation point (D), which is clearly distinct from the phase (B) according to the measurements of $\unicode[STIX]{x1D6E9}_{0}^{\ast }$ . During the phase (D), the characteristics of the flow are similar than the ones described in § 3.2 for $Ek=10^{-8}$ , so this saturated regime corresponds to the strong branch of convection.

Relaxation oscillations are a well-known phenomenon in rotating spherical convection for Rayleigh numbers a few times above onset (e.g. Grote & Busse Reference Grote and Busse2001; Morin & Dormy Reference Morin and Dormy2004; Teed et al. Reference Teed, Jones and Hollerbach2012). However, unlike in these higher $Ra$ relaxation oscillations, which are caused by the disruption of the convection by the zonal flow and the axisymmetric temperature gradient (Teed et al. Reference Teed, Jones and Hollerbach2012), we do not observe here a time lag between the maximum of $Re_{0}^{\ast }$ and $Re_{c}^{\ast }$ , so the oscillations of figure 10 are produced by a different mechanism.

Figure 12.

Snapshots of the streamfunction and the zonal flow for $E=10^{-7}$ , $Pr=0.1$ and $Ra=1.004Ra_{c}$ . Between the first and last snapshots, a time corresponding to about $10^{4}$ rotation periods has elapsed.

Figure 12 shows a series of snapshots of the streamfunction (colour) and the zonal flow (solid black line) during the different phases of the nonlinear oscillation for $Ra=1.004Ra_{c}$ . During the growing phase (A) ( $t=0.076$ ), the flow consists of a thermal Rossby wave and a double zonal jet, both developing around the radius $s=0.5$ , reminiscent of the solution on the weak branch. When the system reaches the maximum (B) ( $t=0.078$ ), the thermal Rossby wave and the zonal flow drift inwards. A narrow prograde zonal jet starts to develop around $s=0.35$ ( $t=0.0781$ ) and rapidly gains a very large amplitude ( $t=0.0783$ ). The non-axisymmetric flow is maximum inside the narrow newly formed prograde jet and keeps the same azimuthal wavenumber as the thermal Rossby wave. The amplitudes of both the zonal and non-axisymmetric flows then decay rapidly (phase (C)). A new oscillation begins with the growth of a thermal Rossby wave and its associated double zonal jet around the radius $s=0.5$ , while the zonal flow of the past oscillation finishes to be dissipated ( $t=0.082$ ). The inward migration of the non-axisymmetric and zonal flows during the phase (B) is probably the cause for the disruption of the convection and the oscillatory behaviour. The mechanism responsible for the inward migration remains unclear as the competitive effects that either inhibit ( $\unicode[STIX]{x1D6FD}$ effect) or promote (gravity along $s$ ) convection depend on radius.

Figure 13.

Space–time diagram of (a) $\overline{u_{\unicode[STIX]{x1D719}}}$ and (b) $\unicode[STIX]{x1D6E5}$ for $E=10^{-7}$ , $Pr=0.1$ and $Ra=1.004Ra_{c}$ .

We note that a secondary effect might arise from the inward migration of the zonal flow in the form of a shear instability. Indeed, the Rossby wave that is confined within the intense narrow prograde jet ( $t=0.0783$ in figure 12) could have been excited by such an instability. A sufficient and necessary condition for the barotropic instability of a shear flow in a inviscid Boussinesq fluid is that the quantity $\unicode[STIX]{x1D6E5}=2\unicode[STIX]{x1D6FD}Ek^{-1}-\text{d}\overline{\unicode[STIX]{x1D701}}/\text{d}s$ , where $\overline{\unicode[STIX]{x1D701}}=\text{d}\overline{u_{\unicode[STIX]{x1D719}}}/\text{d}s+\overline{u_{\unicode[STIX]{x1D719}}}/s$ , changes sign at some radius (Kuo Reference Kuo1949; Ingersoll & Pollard Reference Ingersoll and Pollard1982). Guervilly, Cardin & Schaeffer (Reference Guervilly, Cardin and Schaeffer2012) showed that the threshold of the instability obtained with numerical simulations at finite $Ek$ approaches asymptotically the inviscid theoretical prediction as the Ekman number is decreased and a good quantitative agreement is obtained at $Ek=10^{-7}$ . Figure 13 shows the space–time diagram of $\overline{u_{\unicode[STIX]{x1D719}}}$ and $\unicode[STIX]{x1D6E5}$ between the time $t=0.0778$ and $t=0.08$ for the same simulation as figure 12. $\unicode[STIX]{x1D6E5}$ is always negative, although it comes close to zero when the narrow prograde jet is at its maximum. It is nonetheless possible that the linear stability criteria is violated, but this cannot be strictly identified in the nonlinear simulation because the zonal flow would be rapidly modified by the growth of the shear instability. We emphasise that if the zonal flow is shear unstable, it is a secondary effect and not the cause of the inward migration.

We find similar nonlinear oscillations in the vicinity of the linear onset of convection for $Ek=10^{-6}$ and $Pr=10^{-1}$ , but in a narrower window of Rayleigh numbers (between $Ra=Ra_{c}$ and $Ra=1.0015Ra_{c}$ ). For larger Rayleigh numbers, the convection is similar to the solution found on the weak branch with thermal Rossby waves. This therefore indicates that the nonlinear oscillations are not a consequence of the existence of a strong branch of convection.