Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores

Guervilly, C.; Cardin, P.

doi:10.1093/gji/ggx315

Abstract

We study rapidly rotating Boussinesq convection driven by internal heating in a full sphere. We use a numerical model based on the quasi-geostrophic approximation for the velocity field, whereas the temperature field is 3-D. This approximation allows us to perform simulations for Ekman numbers down to 10⁻⁸, Prandtl numbers relevant for liquid metals (∼10⁻¹) and Reynolds numbers up to 3 × 10⁴. Persistent zonal flows composed of multiple jets form as a result of the mixing of potential vorticity. For the largest Rayleigh numbers computed, the zonal velocity is larger than the convective velocity despite the presence of boundary friction. The convective structures and the zonal jets widen when the thermal forcing increases. Prograde and retrograde zonal jets are dynamically different: in the prograde jets (which correspond to weak potential vorticity gradients) the convection transports heat efficiently and the mean temperature tends to be homogenized; by contrast, in the cores of the retrograde jets (which correspond to steep gradients of potential vorticity) the dynamics is dominated by the propagation of Rossby waves, resulting in the formation of steep mean temperature gradients and the dominance of conduction in the heat transfer process. Consequently, in quasi-geostrophic systems, the width of the retrograde zonal jets controls the efficiency of the heat transfer.

Core, Numerical modelling, Planetary interiors, Heat generation and transport

1 INTRODUCTION

Convection is the main heat transport process in the liquid cores of planets and is thought to be responsible for the generation of planetary magnetic fields. Convection is strongly affected by the rapid rotation of the planet via the action of the Coriolis force. Owing to the very low fluid viscosity, the convective flows are turbulent, although the nonlinear inertial effects are relatively weak compared with the Coriolis force. Under these conditions, and in the absence of magnetic fields, the primary dynamical balance is established between the Coriolis force and the pressure gradient and is called geostrophic balance. Geostrophic flows are invariant along the rotation axis, and so, in spherical geometry, they can only be axisymmetric and azimuthal (i.e. zonal). Convective flows, which are directed along the direction of gravity, cannot be exactly geostrophic, but nevertheless form tall columnar flows aligned with the rotation axis (Jones 2015); such flows are commonly referred to as ‘quasi-geostrophic’ (QG). These columnar convective flows produce coherent Reynolds stresses that drive geostrophic zonal flows (Gilman 1977; Busse & Hood 1982). Stress-free boundary conditions, where boundary friction is absent, favour the emergence of strong zonal flows (e.g. Aurnou & Olson 2001). In models with relatively small viscosity (which can be measured by the Ekman number, the ratio of the rotation period to the global viscous timescale), the zonal flows develop persistent multiple jets of alternating sign inside the tangent cylinder (e.g. Heimpel et al. 2005; Gastine et al. 2014). The observation of intense jets in geophysical and astrophysical objects (e.g. Schou et al. 1998; Porco et al. 2003; Livermore et al. 2017) has prompted much effort dedicated to their study, and in particular, their width and amplitude (e.g. Christensen 2002; Gillet et al. 2007; Read et al. 2015; Cabanes et al. 2017). Although zonal flows (and shear flows in general) do not transport heat outwards, they strongly affect the convection because they can deflect and shear the convective flows, thereby reducing the efficiency of the heat transfer (e.g. Aurnou et al. 2008; Goluskin et al. 2014; von Hardenberg et al. 2015; Yadav et al. 2016). In the present paper, we explore the effect of intense, multiple zonal jets on the convective heat transport in turbulent rotating convection for small Ekman numbers.

The numerical modelling of turbulent rotating flows is extremely challenging as it necessitates a wide range of dynamical length and time scales. Numerical models must therefore employ Ekman numbers that are several orders of magnitude larger than those found in natural objects. However, in the absence of magnetic fields, the lengthscale of the convective flows scales with the Ekman number, at least at the linear onset of convection. The coherence of the Reynolds stresses, and hence the width and amplitude of the zonal flows, might well be affected by the convective lengthscale, and thus by the Ekman number. In order to approach turbulent rotationally constrained convection at small Ekman numbers, we alleviate part of the computational limitations by using a QG approximation that was developed by Busse & Or (1986) for thermal convection in the annulus geometry of Busse (1970) with curved boundaries. The model neglects the variations of the axial vorticity of the flow along the rotation axis, which allows to compute the velocity in 2-D. This is an important limitation to the full dynamics of rotating convection (e.g. Calkins et al. 2013), but the rationale of using this QG model is that it allows the exploration of currently inaccessible regions of the parameter space, thereby informing future 3-D studies. Variations of the QG model have been successfully applied in numerous studies in spherical geometry (e.g. Cardin & Olson 1994; Morin & Dormy 2004; Calkins et al. 2012). Where possible, results from these studies have been successfully benchmarked against asymptotic theories (Gillet & Jones 2006; Labbé et al. 2015), 3-D numerical models (Aubert et al. 2003; Plaut et al. 2008), and laboratory experiments (Aubert et al. 2003; Schaeffer & Cardin 2005; Gillet et al. 2007).

Following the model constructed in Guervilly & Cardin (2016), we use a hybrid numerical model that couples the QG velocity to a 3-D implementation of the temperature in the whole sphere, in order to account for the spherical symmetry of the basic temperature background. The buoyancy driving is controlled by the temperature averaged along the direction of the rotation axis, which, contrary to QG models using a 2-D temperature field (Busse & Or 1986), is not assumed to be equal to the temperature in the equatorial plane. Solving the temperature in 3-D will allow us to assess the influence of the 3-D temperature on the QG dynamics. This implementation is particularly appropriate to model fluids with small Prandtl numbers (the ratio of the viscosity to the thermal diffusivity) that are typical of liquid metals (⁠|$\mathcal {O}(10^{-1})$|⁠) by permitting the use of a 3-D grid for the temperature that is coarser than the 2-D grid used for the velocity. For simplicity, we consider only thermal convection in a full sphere without a solid inner core. The thermal convection is driven by a homogenous internal heating, which is more relevant for the early history of the Earth's core.

The existence of a so-called strong branch of convection driven by internal heating, as first suggested by the weakly nonlinear analysis of Soward (1977), was recently found numerically by Guervilly & Cardin (2016) with the hybrid QG-3D model and by Kaplan et al. (2017) with a fully 3-D model for Ekman numbers smaller than |$\mathcal {O}(10^{-7})$| and Prandtl numbers smaller than unity. The bifurcation is subcritical at the onset of convection and the strong branch is characterized by Reynolds numbers greater than 1000 near the onset and strong zonal flows. In this paper, we focus on the production of zonal flows on this strong branch of convection for Ek ∈ [10⁻⁸, 10⁻⁷] and Pr ∈ [10⁻², 10⁻¹].

The layout of the paper is as follows. In Section 2, we detail the formulation of the hybrid QG-3D model. In Section 3, we describe the radial dependence of the convection and zonal flows and quantify the dependence of the jet width, convective lengthscale and zonal flow velocity on the model parameters. The drift and stability of the zonal flows is discussed in Section 4 and their mechanism of formation in Section 5. The effect of the zonal flows on the heat transport is presented in Section 6. Finally, a discussion of the results is given in Section 7.

2 MATHEMATICAL FORMULATION

We study Boussinesq thermal convection driven by internal heating in a rotating sphere. The rotation vector is |$\Omega \boldsymbol{e}_z$|⁠, where Ω is constant. The acceleration due to gravity is radial and linear, |$\boldsymbol{g}=g_0 r \boldsymbol{e}_r$|⁠. The radius of the sphere is r_o and no inner core is present. The fluid has kinematic viscosity ν, thermal diffusivity κ, density ρ, heat capacity at constant pressure C_p, and thermal expansion coefficient α, all of which are constant. We consider a homogeneous internal volumetric heating S. In the absence of convection, the static temperature profile T_s is calculated by solving the diffusive heat equation and can be written as

\begin{equation} T_s (r) = T_o + \frac{S}{6\kappa \rho C_p} (r_o^2-r^2), \end{equation}

(1)

where T_o is the imposed temperature at the boundary, r = r_o. The governing equations are solved in dimensionless form, obtained by scaling lengths with r_o, times with |$r_o^2/\nu$|⁠, and temperature with |$\nu S r_o^2/(6\rho C_p\kappa ^2)$|⁠. The system of dimensionless equations is:

\begin{eqnarray} \frac{\mbox{$\partial $} \mathbf {u}}{\mbox{$\partial $} t} + \left(\mathbf {u}\cdot \boldsymbol{\nabla } \right)\mathbf {u}+ \frac{2}{\mbox{{Ek}}}\mathbf {e}_z \times \mathbf {u}= - \boldsymbol{\nabla } p + \boldsymbol{\nabla }^2 \mathbf {u}+ \mbox{{Ra}}\Theta r\boldsymbol{e}_r , \end{eqnarray}

(2)

\begin{eqnarray} \nabla \cdot \mathbf {u}= 0, \end{eqnarray}

(3)

\begin{eqnarray} \frac{\mbox{$\partial $} \Theta }{\mbox{$\partial $} t} + \mathbf {u}\cdot \nabla \Theta -\frac{2}{\mbox{{Pr}}} r u_r = \frac{1}{\mbox{{Pr}}} \nabla ^2 \Theta , \end{eqnarray}

(4)

where u is the velocity field, p the pressure and Θ the temperature perturbation relative to the static temperature (1).

The dimensionless numbers are, the Ekman number,

\begin{equation} \mbox{{Ek}}=\frac{\nu }{\Omega r_o^2}, \end{equation}

(5)

the Rayleigh number,

\begin{equation} \mbox{{Ra}}=\frac{\alpha g_0 S r_o^6}{6 \rho C_p \nu \kappa ^2}, \end{equation}

(6)

and the Prandtl number,

\begin{equation} \mbox{{Pr}}=\frac{\nu }{\kappa }. \end{equation}

(7)

At r = r_o, the boundary condition for the velocity is no-slip and impenetrable and the temperature is fixed,

\begin{equation} \mathbf {u}= \boldsymbol{0}, \quad \Theta = 0\,\, {\rm at} \, \, r=r_o. \end{equation}

(8)

Throughout this paper, we use both spherical coordinates (r, θ, ϕ) and cylindrical polar coordinates (s, ϕ, z). The mathematical formulation and the numerical method are described in detail in Guervilly & Cardin (2016), where the linearized version of the code is benchmarked against theoretical and previous numerical results at the onset of convection. The governing equations and the assumptions of the model are briefly described below.

2.1 Governing equation for the non-axisymmetric flow

To model the system of eqs (2)–(4) for small Ekman and Rossby numbers, we use the QG approximation to model the evolution of the velocity field (e.g. Or & Busse 1987; Cardin & Olson 1994; Gillet & Jones 2006). The QG approximation reduces the 3-D system to a 2-D system by taking advantage of the small variations of the flow along z compared with variations in s and ϕ due to the rapid rotation. This approximation is only justified in the case of small slope of the boundaries, such as the Busse (1970) annulus. In the case of a sphere, the approximation is therefore not rigorously justified in any asymptotic limit. Consequently, our 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. 2003; Morin & Dormy 2004; Gillet & Jones 2006; Gillet et al. 2007; Plaut et al. 2008).

The QG model assumes that the fluid dynamics is dominated by the geostrophic balance, that is, the Coriolis force balances the pressure gradient at leading order. The leading-order velocity u^g is invariant along z and |$\mathbf {u}^g = (u^g_s,u^g_{\phi },0)$| in cylindrical polar coordinates. By taking the z-component of the curl of the momentum eq. (2) and averaging it along z, we obtain the equation for the axial vorticity, |$\zeta ^g = \left(\boldsymbol{\nabla }\times \mathbf {u}^g \right)\cdot \boldsymbol{e}_z$|⁠,

\begin{eqnarray} \frac{\mbox{$\partial $} \zeta ^g}{\mbox{$\partial $} t} + \left(\mathbf {u}^g \cdot \boldsymbol{\nabla }\right)\zeta ^g - \left(\frac{2}{\mbox{{Ek}}} + \zeta ^g \right)\left\langle \frac{\partial u_z}{\partial z} \right\rangle = \nabla _e^2 \zeta ^g - \mbox{{Ra}}\left\langle \frac{\partial \Theta }{\partial \phi } \right\rangle ,\nonumber\\ \end{eqnarray}

(9)

with

\begin{equation} \nabla ^2_e A \equiv \frac{1}{s} \frac{\partial }{\partial s}\left(s \frac{\partial A}{\partial s} \right)+ \frac{1}{s^2} \frac{\partial ^2 A}{\partial \phi ^2}, \end{equation}

(10)

and

\begin{equation} {\left\langle A \right\rangle \equiv \frac{1}{2H} \int ^{+H}_{-H} A \text{d}z, } \end{equation}

(11)

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

The velocity u^g can be described by a streamfunction ψ that models the non-axisymmetric (i.e. ϕ-dependent) components with the addition of an axisymmetric azimuthal flow,

\begin{equation} \mathbf {u}^g = \frac{1}{H} \boldsymbol{\nabla }\times \left(H \psi \boldsymbol{e}_z \right)+ \overline{u_{\phi }^g} \boldsymbol{e}_{\phi }, \end{equation}

(12)

where

\begin{equation} \skew8\bar{A} \equiv \frac{1}{2\pi } \int _{0}^{2\pi }A \text{d}\phi . \end{equation}

(13)

We choose this formulation for the streamfunction to account for the non-zero divergence of u^g in the equatorial plane due to the return axial flow at the sloping boundaries,

\begin{equation} \boldsymbol{\nabla }_e \cdot \mathbf {u}^g = -\beta u_s^g, \end{equation}

(14)

where

\begin{equation} \boldsymbol{\nabla }_e \cdot \boldsymbol{A} \equiv \frac{1}{s} \frac{\partial s A_s}{\partial s} + \frac{1}{s} \frac{\partial A_{\phi }}{\partial \phi }, \end{equation}

(15)

and

\begin{equation} \beta = \frac{1}{H}\frac{dH}{ds} = - \frac{s}{H^2}. \end{equation}

(16)

The axial velocity u_z is assumed to be linear in z. The third term on the left-hand side of eq. (9) requires us to determine u_z at the boundary z = ±H:

\begin{equation} {\left. u_z\right|_{\pm H} = \pm \frac{1}{H} \left. \mathbf {u}\cdot \mathbf {n}\right|_{\pm H} \pm \beta H u_s^g,} \end{equation}

(17)

where the normal vector at the boundary is |$\boldsymbol{n} = \boldsymbol{e}_r$|⁠. The normal component, .u · n|_{± H}, is the Ekman pumping induced by the viscous boundary layer and is determined by asymptotic methods for a linear Ekman layer, |$\left. \mathbf {u}\cdot \mathbf {n}\right|_{z=\pm H} = \mbox{{Ek}}^{1/2} P(s,u_s^g, u_{\phi }^g)$| (Greenspan 1968). The analytical function P is derived for a spherical boundary in Schaeffer & Cardin (2005).

The numerical code solves the evolution equation of the non-axisymmetric streamfunction ψ. The no-slip and impenetrable boundary conditions imply that ψ = ∂_sψ = 0 at s = 1. We use the regularity condition |$\hat{\psi }^m=\mathcal {O}(s^m)$| at s = 0, where |$\hat{\psi }^m(s,t)$| is the Fourier mode of azimuthal wavenumber m (see Guervilly & Cardin 2016, for more detail).

2.2 Governing equation for the zonal flow

In our model, the streamfunction ψ only describes the non-axisymmetric motions, so the axisymmetric azimuthal flows, or zonal flows, are treated separately. We take the ϕ- and z-averages of the ϕ-component of the momentum equation to obtain

\begin{equation} \frac{\mbox{$\partial $} \overline{u_{\phi }^g}}{\mbox{$\partial $} t} + \overline{u_{s}^g \frac{\partial u_{\phi }^g}{\partial s}} + \overline{\frac{u_{s}^g u_{\phi }^g}{s}} +\frac{2}{\mbox{{Ek}}} \left\langle \overline{u_s} \right\rangle = \nabla ^2 \overline{u_{\phi }^g} - \frac{\overline{u_{\phi }^g}}{s^2} . \end{equation}

(18)

Note that the geostrophic balance imposes that |$\overline{u_{s}^g}=0$|⁠. The fourth term on the left-hand side of (18) 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. 2003) that

\begin{equation} {\left\langle \overline{u_s} \right\rangle = \frac{\mbox{{Ek}}^{1/2}}{2H^{3/2}} \overline{u_{\phi }^g}.} \end{equation}

(19)

The no-slip boundary condition at the outer sphere and the symmetry at the centre imply that |$\overline{u_{\phi }^g} = 0$| at s = 0, 1.

2.3 Governing equation for the temperature

The dimensionless equation for the evolution of the temperature perturbation in 3-D is

\begin{equation} \frac{\mbox{$\partial $} \Theta }{\mbox{$\partial $} t}+ \mathbf {u}^{3d} \cdot \boldsymbol{\nabla } \Theta = \frac{1}{\mbox{{Pr}}} \left(2 r u^{3d}_r+\nabla ^2 \Theta \right)\!. \end{equation}

(20)

where u^3d is the velocity in 3-D. In cylindrical polar coordinates,

\begin{equation} \mathbf {u}^{3d} = \left(u_s^g, u_{\phi }^g,\mbox{{Ek}}^{1/2} z P + \beta z u_s^g\right) . \end{equation}

(21)

The temperature is fixed at the outer boundary so Θ = 0 at r = 1. At the centre of the sphere, the non-spherically symmetric components of Θ are zero by symmetry and the spherically symmetric component of ∂_rΘ is zero.

2.4 Numerical method

In the following, the superscripts g are removed for clarity. The evolution equations for ψ and |$\overline{u_{\phi }}$| 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 variables are expanded in Fourier modes. The evolution equation for the temperature 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. Further detail about the numerical interpolations between the 2-D and 3-D grids used to compute the buoyancy term and the advection of the temperature can be found in Guervilly & Cardin (2016) and Guervilly (2010).

Table 1 gives the list of the simulations presented in this paper with some output quantities and the numerical resolutions. To quantify some of the global properties of convection, we often use the Reynolds number, which is calculated from the output of the simulations and corresponds to the time-averaged root mean square (r.m.s.) value of the velocity in dimensionless unit,

\begin{equation} \mbox{{Re}}= \frac{1}{\Delta t} \int _{\Delta t} \left(\frac{3}{4\pi } \int _0^{2\pi } \int _0^{1} \left(u_s^2 + u_{\phi }^2\right) 2 H(s) s \text{d}s \text{d}\phi \right)^{1/2} \text{d}t, \end{equation}

(22)

where u_ϕ includes the zonal velocity. We measure the convective Reynolds number, Re_c, as in eq. (22) but including only the non-axisymmetric velocity. Similarly, we measure the zonal Reynolds number, Re₀, including only the axisymmetric velocity. The integration time over which the time averages are calculated is indicated in Table 1 for each simulation.

Table 1.

Open in new tab

List of input and output parameters for all the simulations presented in the paper. Re_c and Re₀ are the convective and zonal Reynolds numbers, respectively. The columns labelled |$(N^u_s,M^u_{\rm max})$| and |$(N^t_r,M^t_{\rm max},L^t_{\rm max})$| give the numerical resolutions on the 2-D and 3-D grids, respectively. The last column gives the integration time used to compute the time averages in units of 1/Ω and, in brackets, in units of a convective turnover timescale, l_c/U_c, where U_c is the r.m.s. convective velocity (equivalent to Re_c in our dimensionless units) and l_c is the convective lengthscale computed from eq. (23) and averaged between 0.1 ≤ s ≤ 0.8.

Ek	Pr	Ra	Ra/Ra_c	Re_c	Re₀	\|$(N^u_s,M^u_{\rm max})$\|	\|$(N^t_r,M^t_{\rm max},L^t_{\rm max})$\|	Integration time
10⁻⁷	10⁻¹	6 × 10⁹	1.19	1012	445	(1100, 200)	(500, 150, 150)	3 × 10⁵ (479)
10⁻⁷	10⁻¹	1 × 10¹⁰	1.99	2070	1110	(1100, 200)	(500, 150, 150)	3 × 10⁵ (794)
10⁻⁷	10⁻¹	2 × 10¹⁰	3.97	3663	2994	(1200, 200)	(500, 150, 150)	10⁵ (414)
10⁻⁷	10⁻¹	3 × 10¹⁰	5.96	5093	5064	(1200, 200)	(500, 150, 150)	10⁴ (54)
10⁻⁷	10⁻¹	4 × 10¹⁰	7.95	6243	6840	(1200, 200)	(500, 150, 150)	10⁴ (61)
10⁻⁷	10⁻¹	5 × 10¹⁰	9.94	7731	9515	(1500, 260)	(600, 180, 180)	10⁴ (73)
10⁻⁷	10⁻²	1.9 × 10⁹	1.05	5215	5710	(1000, 160)	(400, 96, 96)	2 × 10⁴ (90)
10⁻⁷	10⁻²	3 × 10⁹	1.65	9027	11708	(1000, 160)	(400, 96, 96)	10⁴ (67)
10⁻⁷	10⁻²	4.8 × 10⁹	2.64	11487	21834	(1200, 180)	(400, 96, 96)	10⁴ (82)
10⁻⁷	10⁻²	8.5 × 10⁹	4.68	17548	40009	(1400, 200)	(400, 96, 96)	10⁴ (105)
10⁻⁸	10⁻¹	7.45 × 10¹⁰	0.96	855	243	(1600, 256)	(700, 200, 200)	2 × 10⁶ (633)
10⁻⁸	10⁻¹	7.8 × 10¹⁰	1.01	1047	291	(1600, 256)	(700, 200, 200)	2 × 10⁶ (785)
10⁻⁸	10⁻¹	1.5 × 10¹¹	1.93	3393	1128	(1800, 280)	(700, 200, 200)	5 × 10⁵ (400)
10⁻⁸	10⁻¹	2 × 10¹¹	2.58	4876	1893	(1800, 280)	(700, 200, 200)	5 × 10⁵ (506)
10⁻⁸	10⁻¹	3 × 10¹¹	3.87	7108	4724	(1800, 280)	(700, 200, 200)	2 × 10⁵ (286)
10⁻⁸	10⁻¹	5 × 10¹¹	6.44	9920	9341	(1900, 300)	(700, 200, 200)	10⁵ (191)
10⁻⁸	10⁻¹	7 × 10¹¹	9.02	12023	13176	(2000, 320)	(750, 200, 200)	2 × 10⁶ (4453)
10⁻⁸	10⁻²	2 × 10¹⁰	0.68	7200	4753	(1400, 220)	(500, 128, 128)	2 × 10⁵ (216)
10⁻⁸	10⁻²	3 × 10¹⁰	1.01	13498	10598	(1400, 220)	(500, 128, 128)	10⁵ (177)
10⁻⁸	10⁻²	5 × 10¹⁰	1.69	23288	23309	(1500, 240)	(500, 128, 128)	4 × 10⁴ (108)
10⁻⁸	10⁻²	8 × 10¹⁰	2.70	33163	43609	(1600, 260)	(500, 128, 128)	5 × 10⁴ (191)

Ek	Pr	Ra	Ra/Ra_c	Re_c	Re₀	\|$(N^u_s,M^u_{\rm max})$\|	\|$(N^t_r,M^t_{\rm max},L^t_{\rm max})$\|	Integration time
10⁻⁷	10⁻¹	6 × 10⁹	1.19	1012	445	(1100, 200)	(500, 150, 150)	3 × 10⁵ (479)
10⁻⁷	10⁻¹	1 × 10¹⁰	1.99	2070	1110	(1100, 200)	(500, 150, 150)	3 × 10⁵ (794)
10⁻⁷	10⁻¹	2 × 10¹⁰	3.97	3663	2994	(1200, 200)	(500, 150, 150)	10⁵ (414)
10⁻⁷	10⁻¹	3 × 10¹⁰	5.96	5093	5064	(1200, 200)	(500, 150, 150)	10⁴ (54)
10⁻⁷	10⁻¹	4 × 10¹⁰	7.95	6243	6840	(1200, 200)	(500, 150, 150)	10⁴ (61)
10⁻⁷	10⁻¹	5 × 10¹⁰	9.94	7731	9515	(1500, 260)	(600, 180, 180)	10⁴ (73)
10⁻⁷	10⁻²	1.9 × 10⁹	1.05	5215	5710	(1000, 160)	(400, 96, 96)	2 × 10⁴ (90)
10⁻⁷	10⁻²	3 × 10⁹	1.65	9027	11708	(1000, 160)	(400, 96, 96)	10⁴ (67)
10⁻⁷	10⁻²	4.8 × 10⁹	2.64	11487	21834	(1200, 180)	(400, 96, 96)	10⁴ (82)
10⁻⁷	10⁻²	8.5 × 10⁹	4.68	17548	40009	(1400, 200)	(400, 96, 96)	10⁴ (105)
10⁻⁸	10⁻¹	7.45 × 10¹⁰	0.96	855	243	(1600, 256)	(700, 200, 200)	2 × 10⁶ (633)
10⁻⁸	10⁻¹	7.8 × 10¹⁰	1.01	1047	291	(1600, 256)	(700, 200, 200)	2 × 10⁶ (785)
10⁻⁸	10⁻¹	1.5 × 10¹¹	1.93	3393	1128	(1800, 280)	(700, 200, 200)	5 × 10⁵ (400)
10⁻⁸	10⁻¹	2 × 10¹¹	2.58	4876	1893	(1800, 280)	(700, 200, 200)	5 × 10⁵ (506)
10⁻⁸	10⁻¹	3 × 10¹¹	3.87	7108	4724	(1800, 280)	(700, 200, 200)	2 × 10⁵ (286)
10⁻⁸	10⁻¹	5 × 10¹¹	6.44	9920	9341	(1900, 300)	(700, 200, 200)	10⁵ (191)
10⁻⁸	10⁻¹	7 × 10¹¹	9.02	12023	13176	(2000, 320)	(750, 200, 200)	2 × 10⁶ (4453)
10⁻⁸	10⁻²	2 × 10¹⁰	0.68	7200	4753	(1400, 220)	(500, 128, 128)	2 × 10⁵ (216)
10⁻⁸	10⁻²	3 × 10¹⁰	1.01	13498	10598	(1400, 220)	(500, 128, 128)	10⁵ (177)
10⁻⁸	10⁻²	5 × 10¹⁰	1.69	23288	23309	(1500, 240)	(500, 128, 128)	4 × 10⁴ (108)
10⁻⁸	10⁻²	8 × 10¹⁰	2.70	33163	43609	(1600, 260)	(500, 128, 128)	5 × 10⁴ (191)

Table 1.

Open in new tab

List of input and output parameters for all the simulations presented in the paper. Re_c and Re₀ are the convective and zonal Reynolds numbers, respectively. The columns labelled |$(N^u_s,M^u_{\rm max})$| and |$(N^t_r,M^t_{\rm max},L^t_{\rm max})$| give the numerical resolutions on the 2-D and 3-D grids, respectively. The last column gives the integration time used to compute the time averages in units of 1/Ω and, in brackets, in units of a convective turnover timescale, l_c/U_c, where U_c is the r.m.s. convective velocity (equivalent to Re_c in our dimensionless units) and l_c is the convective lengthscale computed from eq. (23) and averaged between 0.1 ≤ s ≤ 0.8.

Ek	Pr	Ra	Ra/Ra_c	Re_c	Re₀	\|$(N^u_s,M^u_{\rm max})$\|	\|$(N^t_r,M^t_{\rm max},L^t_{\rm max})$\|	Integration time
10⁻⁷	10⁻¹	6 × 10⁹	1.19	1012	445	(1100, 200)	(500, 150, 150)	3 × 10⁵ (479)
10⁻⁷	10⁻¹	1 × 10¹⁰	1.99	2070	1110	(1100, 200)	(500, 150, 150)	3 × 10⁵ (794)
10⁻⁷	10⁻¹	2 × 10¹⁰	3.97	3663	2994	(1200, 200)	(500, 150, 150)	10⁵ (414)
10⁻⁷	10⁻¹	3 × 10¹⁰	5.96	5093	5064	(1200, 200)	(500, 150, 150)	10⁴ (54)
10⁻⁷	10⁻¹	4 × 10¹⁰	7.95	6243	6840	(1200, 200)	(500, 150, 150)	10⁴ (61)
10⁻⁷	10⁻¹	5 × 10¹⁰	9.94	7731	9515	(1500, 260)	(600, 180, 180)	10⁴ (73)
10⁻⁷	10⁻²	1.9 × 10⁹	1.05	5215	5710	(1000, 160)	(400, 96, 96)	2 × 10⁴ (90)
10⁻⁷	10⁻²	3 × 10⁹	1.65	9027	11708	(1000, 160)	(400, 96, 96)	10⁴ (67)
10⁻⁷	10⁻²	4.8 × 10⁹	2.64	11487	21834	(1200, 180)	(400, 96, 96)	10⁴ (82)
10⁻⁷	10⁻²	8.5 × 10⁹	4.68	17548	40009	(1400, 200)	(400, 96, 96)	10⁴ (105)
10⁻⁸	10⁻¹	7.45 × 10¹⁰	0.96	855	243	(1600, 256)	(700, 200, 200)	2 × 10⁶ (633)
10⁻⁸	10⁻¹	7.8 × 10¹⁰	1.01	1047	291	(1600, 256)	(700, 200, 200)	2 × 10⁶ (785)
10⁻⁸	10⁻¹	1.5 × 10¹¹	1.93	3393	1128	(1800, 280)	(700, 200, 200)	5 × 10⁵ (400)
10⁻⁸	10⁻¹	2 × 10¹¹	2.58	4876	1893	(1800, 280)	(700, 200, 200)	5 × 10⁵ (506)
10⁻⁸	10⁻¹	3 × 10¹¹	3.87	7108	4724	(1800, 280)	(700, 200, 200)	2 × 10⁵ (286)
10⁻⁸	10⁻¹	5 × 10¹¹	6.44	9920	9341	(1900, 300)	(700, 200, 200)	10⁵ (191)
10⁻⁸	10⁻¹	7 × 10¹¹	9.02	12023	13176	(2000, 320)	(750, 200, 200)	2 × 10⁶ (4453)
10⁻⁸	10⁻²	2 × 10¹⁰	0.68	7200	4753	(1400, 220)	(500, 128, 128)	2 × 10⁵ (216)
10⁻⁸	10⁻²	3 × 10¹⁰	1.01	13498	10598	(1400, 220)	(500, 128, 128)	10⁵ (177)
10⁻⁸	10⁻²	5 × 10¹⁰	1.69	23288	23309	(1500, 240)	(500, 128, 128)	4 × 10⁴ (108)
10⁻⁸	10⁻²	8 × 10¹⁰	2.70	33163	43609	(1600, 260)	(500, 128, 128)	5 × 10⁴ (191)

Ek	Pr	Ra	Ra/Ra_c	Re_c	Re₀	\|$(N^u_s,M^u_{\rm max})$\|	\|$(N^t_r,M^t_{\rm max},L^t_{\rm max})$\|	Integration time
10⁻⁷	10⁻¹	6 × 10⁹	1.19	1012	445	(1100, 200)	(500, 150, 150)	3 × 10⁵ (479)
10⁻⁷	10⁻¹	1 × 10¹⁰	1.99	2070	1110	(1100, 200)	(500, 150, 150)	3 × 10⁵ (794)
10⁻⁷	10⁻¹	2 × 10¹⁰	3.97	3663	2994	(1200, 200)	(500, 150, 150)	10⁵ (414)
10⁻⁷	10⁻¹	3 × 10¹⁰	5.96	5093	5064	(1200, 200)	(500, 150, 150)	10⁴ (54)
10⁻⁷	10⁻¹	4 × 10¹⁰	7.95	6243	6840	(1200, 200)	(500, 150, 150)	10⁴ (61)
10⁻⁷	10⁻¹	5 × 10¹⁰	9.94	7731	9515	(1500, 260)	(600, 180, 180)	10⁴ (73)
10⁻⁷	10⁻²	1.9 × 10⁹	1.05	5215	5710	(1000, 160)	(400, 96, 96)	2 × 10⁴ (90)
10⁻⁷	10⁻²	3 × 10⁹	1.65	9027	11708	(1000, 160)	(400, 96, 96)	10⁴ (67)
10⁻⁷	10⁻²	4.8 × 10⁹	2.64	11487	21834	(1200, 180)	(400, 96, 96)	10⁴ (82)
10⁻⁷	10⁻²	8.5 × 10⁹	4.68	17548	40009	(1400, 200)	(400, 96, 96)	10⁴ (105)
10⁻⁸	10⁻¹	7.45 × 10¹⁰	0.96	855	243	(1600, 256)	(700, 200, 200)	2 × 10⁶ (633)
10⁻⁸	10⁻¹	7.8 × 10¹⁰	1.01	1047	291	(1600, 256)	(700, 200, 200)	2 × 10⁶ (785)
10⁻⁸	10⁻¹	1.5 × 10¹¹	1.93	3393	1128	(1800, 280)	(700, 200, 200)	5 × 10⁵ (400)
10⁻⁸	10⁻¹	2 × 10¹¹	2.58	4876	1893	(1800, 280)	(700, 200, 200)	5 × 10⁵ (506)
10⁻⁸	10⁻¹	3 × 10¹¹	3.87	7108	4724	(1800, 280)	(700, 200, 200)	2 × 10⁵ (286)
10⁻⁸	10⁻¹	5 × 10¹¹	6.44	9920	9341	(1900, 300)	(700, 200, 200)	10⁵ (191)
10⁻⁸	10⁻¹	7 × 10¹¹	9.02	12023	13176	(2000, 320)	(750, 200, 200)	2 × 10⁶ (4453)
10⁻⁸	10⁻²	2 × 10¹⁰	0.68	7200	4753	(1400, 220)	(500, 128, 128)	2 × 10⁵ (216)
10⁻⁸	10⁻²	3 × 10¹⁰	1.01	13498	10598	(1400, 220)	(500, 128, 128)	10⁵ (177)
10⁻⁸	10⁻²	5 × 10¹⁰	1.69	23288	23309	(1500, 240)	(500, 128, 128)	4 × 10⁴ (108)
10⁻⁸	10⁻²	8 × 10¹⁰	2.70	33163	43609	(1600, 260)	(500, 128, 128)	5 × 10⁴ (191)

3 STRUCTURE OF THE CONVECTIVE AND ZONAL FLOWS AND SCALING OF THE VELOCITY

3.1 Radial dependence of the convection

In this section, the Ekman and Prandtl numbers are fixed to Ek = 10⁻⁸ and Pr = 10⁻¹. For these parameters, the stable solution is located on a strong branch of convection, which is discontinuous at the onset of convection (Guervilly & Cardin 2016; Kaplan et al. 2017). All cases presented in this paper are located on the strong branch. This branch is distinct from the weak branch of convection, which occurs for |$\mbox{{Ek}}\mbox{{Pr}}\gtrsim \mathcal {O}(10^{-8})$| and is continuous at the onset of convection. At the onset of convection, solutions on the weak branch take the form of propagating structures that are tilted in the prograde direction. These structures are known as thermal Rossby waves and have been extensively studied in the literature (e.g. Busse 1970; Zhang 1992). Near the onset of convection, the flows on the strong branch are starkly different and are described in detail below.

For Ek = 10⁻⁸ and Pr = 10⁻¹, the nonlinear convection is maintained below the linear onset of convection (quantified by the critical Rayleigh number Ra_c), down to a value Ra = 0.96Ra_c (Guervilly & Cardin 2016). We vary the Rayleigh number from this lower value to approximately 9Ra_c. Fig. 1 shows snapshots of the radial and azimuthal velocities in the equatorial plane for Ra/Ra_c = 0.96, Ra/Ra_c = 1.93 and Ra/Ra_c = 9.02. In all three cases, two dynamical regions can be distinguished: an inner region, where the convection is vigorous with values of the radial velocity up to 3000 for the lowest Ra and up to 30 000 for the largest Ra, and an outer region, where the radial flow has smaller amplitude and the flow has finer structures that are tilted in the prograde direction. This radial dependence of the convection, sometimes referred to as dual convection, was previously described in laboratory experiments (Sumita & Olson 2000) and QG (Aubert et al. 2003) and 3-D numerical models (Miyagoshi et al. 2010). The limit between the two regions is located around s = 0.5 for the lowest Ra and s = 0.8 for the largest Ra, so the limit moves outwards when the convection becomes more vigorous. In the inner region, the convective flows are strongly time dependent, especially for large Ra, and are subject to frequent nonlinear interactions. The contours of the radial velocity tend to be directed radially, contrary to the tilted contours of the outer region. For large Ra, the azimuthal lengthscales of the radial flow decreases with increasing radius. This is likely due to the increase of the slope of the boundary β with radius: the vortex stretching term in the axial vorticity equation, which depends on β, impedes the radial motion of wide vortices. The azimuthal extent of the convective flows clearly increases with Ra, which indicates the presence of an upscale energy transfer as expected in β-plane turbulence (e.g. Davidson 2013). For all Rayleigh numbers, the radial velocity is weak in the central region because the gravity goes to zero at the centre. In the outer region, β is large and the vortex stretching term is the dominant source of the axial vorticity, so this outer region is dominated by the propagation of Rossby waves. The nonlinear interactions are weaker in this region.

Snapshots of the radial velocity (left) and the azimuthal velocity (right) in the equatorial plane for Ek = 10−8 and Pr = 10−1 for (a) Ra/Rac = 0.96, (b) Ra/Rac = 1.93 and (c) Ra/Rac = 9.02.

Figure 1.

Snapshots of the radial velocity (left) and the azimuthal velocity (right) in the equatorial plane for Ek = 10⁻⁸ and Pr = 10⁻¹ for (a) Ra/Ra_c = 0.96, (b) Ra/Ra_c = 1.93 and (c) Ra/Ra_c = 9.02.

Open in new tab Download slide

For all Ra, the azimuthal flow has a visible axisymmetric (i.e. zonal) component. Fig. 2 shows the time-averaged profiles of the zonal flow for different Ra. The zonal flow is prograde in the outermost region for all Ra. In the outer region dominated by Rossby waves, the Reynolds stresses due to the correlation of the velocity along the tilted contours produce a prograde jet in the outer part and a neighbouring inner retrograde jet (e.g. Busse & Hood 1982). The behaviour of the zonal flow in the inner convective region is different depending on Ra. For the lowest Ra in Fig. 1, azimuthal flows appear to spiral inward from mid-radius. These flows have an axisymmetric average that is positive in the centre and negative near s = 0.4. The time-averaged profile of the zonal flow for this Ra has therefore three jets of alternating sign. For Ra/Ra_c = 1.93, the azimuthal flows consist in a multitude of meandering narrow jets. Around mid-radius, the azimuthal average of the narrow meandering jets is not well-defined. Their axisymmetric average is mostly negative because the retrograde jets have stronger amplitude. In the centre, the azimuthal flow has a clear prograde direction and is wider than the meandering jets at larger radius. For Ra/Ra_c = 9.02, the multiple azimuthal jets have stronger velocity than the radial flow and they do not meander so their net axisymmetric average is well-defined. The radial profile of the zonal flow shows persistent multiple jets. The central jet remains prograde and is wider than the jets located at larger radius. Overall, for Ra > 3Ra_c, the zonal flows develop persistent multiple jets and the innermost and outermost jets always remain prograde. The region occupied by the in-between jets becomes wider as Ra increases and the limit between inner convective region and outer Rossby wave region moves outwards. A thin viscous boundary layer can be observed on the profiles of the zonal flow of the largest Rayleigh number as the boundary condition is no-slip at s = 1. This thin boundary layer is well resolved in our model as the radial grid is refined near the boundary.

Figure 2.

Radial profile of the zonal flow (time average) for different Rayleigh numbers (indicated at the top of each subplot) and Ekman and Prandtl numbers (indicated at the bottom). The vertical axis is the radius. The range of the horizontal axis is different for each subplot: the maximum of the zonal flow increases 80-fold between the smallest and largest Rayleigh numbers for Ek = 10⁻⁸ and Pr = 10⁻¹.

Open in new tab Download slide

For large Ra, the zonal velocity is larger than the radial velocity. In this case, the radial shear exerted by the zonal flow can be faster than the vortex turnover timescale, so the zonal flow has a dominant role in the dynamics of the convective vortices. The radial flow has a smaller radial extent at the larger Ra due to the presence of the multiple zonal jets of strong amplitude. The convective structures then change from narrow in ϕ and extended in r to fatter in ϕ and shortened in r as Ra increases.

3.2 Zonal jet width and convective lengthscale

To study the influence of the zonal jets on the convective flow and vice versa, we compare the width of the zonal jet and the convective lengthscale. We expect these two lengthscales to be correlated: on the one hand, the radial shear exerted by the zonal flow on the radial velocity limits the size of the convective flow; on the other hand, the size of the most energetic convective eddies controls the width of the zonal jet by controlling the mixing length—this mechanism is discussed in detail in Section 5. To give an estimate of the typical lengthscale of the convection, we compute the integral lengthscale of the non-axisymmetric flow, which is defined as

\begin{equation} l_c(s) = \pi s\left(\frac{\sum \limits _{m=1}^{M^u_{\rm max}} E(m,s)}{\sum \limits _{m=1}^{M^u_{\rm max}} m E(m,s)}\right), \end{equation}

(23)

where E(m, s) is the time-averaged kinetic energy of a mode m of the Fourier decomposition. The axisymmetric mode (m = 0) is excluded from this definition. The convective lengthscale l_c and the width of the jets (denoted l_j) are plotted as a function of the radius in Fig. 3 for Ra/Ra_c = 9.02. The regions of retrograde zonal jets are indicated in grey. Both the jet width and the convective scale tend to decrease with radius in the inner convective region (s < 0.8). l_c takes local minimum values in the core of the jets and maximum values in the flanks of the jets where the radial shear is largest. l_c corresponds to an azimuthal lengthscale, so this indicates that the shear from the zonal flow elongates the convective flow in the azimuthal direction as expected. The values of l_c and l_j are comparable, although the convective scale is slightly smaller than the jet width. We have thus verified that a strong correlation exists between the convective scale and the jet width.

Radial profiles of the Rhines scale lR, the integral scale of the convective flow lc and the width of the jets lj for Ek = 10−8, Ra/Rac = 9.02 at Pr = 10−1. The grey bands correspond to the regions where the zonal flow is retrograde.

Figure 3.

Radial profiles of the Rhines scale l_R, the integral scale of the convective flow l_c and the width of the jets l_j for Ek = 10⁻⁸, Ra/Ra_c = 9.02 at Pr = 10⁻¹. The grey bands correspond to the regions where the zonal flow is retrograde.

Open in new tab Download slide

In the framework of β-plane turbulence, a fundamental length scale of the flow is the Rhines scale (Rhines 1975), which is given by

\begin{equation} l_{R}(s) = \pi \left(\frac{\mbox{{Ek}}U}{2 |\beta |} \right)^{1/2}, \end{equation}

(24)

where U is the r.m.s. velocity of the flow. In the literature, U is interpreted as either an eddy velocity, a jet velocity or a total velocity (Ingersoll & Pollard 1982; Dritschel & McIntyre 2008). Here we will first consider that U is the r.m.s convective velocity, and later discuss how our results differ when using the total velocity instead. Note that Ek appears in this formula due to our choice of dimensional units. The Rhines scale may be considered as the scale separating dynamical regimes dominated by either the turbulence (at smaller scales) or by Rossby waves (at larger scales). The width of the zonal jets is often thought to be closely related to the Rhines scales, although this notion has been called into question (e.g. Sukoriansky et al. 2007; Scott & Dritschel 2012). Here we are interested in the predicted dependence of the Rhines scale on β and U (rather than on the actual value of l_R given by (24) that arbitrarily includes a factor |$\pi /\sqrt{2}$|⁠) as this dependence can be compared with the lengthscales computed from our data set. We first compare the radial dependence of the Rhines scale to l_c and l_j. We plot l_R in Fig. 3, where we used the r.m.s value of the non-axisymmetric velocity (which varies in radius) to estimate U. In the inner convective region, the convective scale is approximately 2 times larger than the Rhines scale on average but the decreasing trends observed for each of the lengthscales is similar: between s = 0.3 and s = 0.7, l_j, l_c and l_R are all approximately divided by two. The ratio l_c/l_R is approximately constant in the inner convective region, which shows that the Rhines scale adequately predicts the radial dependence of the convective scale and the jet width.

As the Rayleigh number (and thus U) increases, the Rhines scale predicts that the most energetic convective eddies become larger. This is indeed what we observe qualitatively on the snapshots of Fig. 1. This increase of the convective scale should be accompanied by an increase of the jet width. Fig. 2 shows that the jets indeed tend to become wider when the Rayleigh number increases. However this increase might be due to the widening of the inner convective region as it pushes the outer region outwards. Larger Ra, currently out of reach of our computational capabilities, would be necessary to observe a sizeable increase in the size of the jets.

We can extend our study on the zonal and convective lengthscales and the predictive value of the Rhines scale from the results of simulations performed at different Ek and Pr. We cannot presently run simulations at Ek < 10⁻⁸, so this study is restricted to higher Ek, namely Ek = 10⁻⁷. Our QG-3D model allows us to explore small Pr, so we also use results from calculations at Pr = 10⁻². The two panels on the right of Fig. 2 shows the zonal velocity for (Ek, Pr) = (10⁻⁸, 10⁻²) and (Ek, Pr) = (10⁻⁷, 10⁻¹) at the largest Ra performed (see Table 1). For (Ek, Pr) = (10⁻⁷, 10⁻¹), the zonal flow has 5 jets of alternating sign. The zonal jets widen when the Ekman number increases if this leads to a larger r.m.s. convective velocity Ek U according to the Rhines scale (24).

For (Ek, Pr) = (10⁻⁸, 10⁻²), the time-averaged zonal flow also has five alternating jets. By visual inspection of the snapshots of the velocity, it appears that meandering azimuthal flows are present in the inner region, similarly to the case Ra = 1.93 for (Ek, Pr) = (10⁻⁸, 10⁻¹). This suggests that larger Rayleigh numbers would be necessary to get multiple persistent jets. However, we were not able to perform calculations at larger Ra to confirm this. Smaller values of Pr lead to larger values of the convective velocity (see Table 1), and hence, to wider jets.

To compare more systematically the convective lengthscale with the Rhines scale, Fig. 4 shows the values of l_c versus l_R that have been radially averaged in the inner convective region (between 0.1 ≤ s ≤ 0.8), for increasing values of Ra and different Ek and Pr. Both lengthscales vary significantly with radius so we also indicate the standard deviation with vertical and horizontal bars. All the simulations are located on the strong branch of convection, which is discontinuous at the onset of convection, and l_c is always larger than the wavelength of the linear instability. The best fit to all the points is |$l_c\sim l_R^{0.73(\pm 0.04)}$|⁠. Our simulations therefore indicate that the convective lengthscale increases with the convective flow speed, but it follows a power law of smaller exponent (namely l_c ∼ U^0.37) than predicted by the Rhines scale (namely l_c ∼ U^0.5). This result is in agreement with the work of Gastine et al. (2016) using 3-D numerical simulations of rotating convection in a spherical shell with Pr = 1: they find that the convective lengthscale approaches the power law given by the Rhines scale when the Ekman decreases, but the exponent remains smaller than a half (0.45 for Ek = 3 × 10⁻⁷). Here we fitted all the points from different sets of Ekman and Prandtl numbers with one power law because, according to the Rhines scale argument, the dependence of the convective lengthscale on the parameters can be explained by a power law dependence on the flow velocity only, irrespective of the values of Ek and Pr. However there are some visible variations of the power law exponent for the different data sets, which is further indication that our data do not entirely corroborate the Rhines scale argument.

Figure 4.

Convective lengthscale versus Rhines scale. Both quantities have been averaged radially in the inner convective region between 0.1 ≤ s ≤ 0.8. The horizontal and vertical bars indicate the standard deviation. The dotted line indicates the best fit to all the data points.

Open in new tab Download slide

We now return to the issue of the interpretation of the velocity scale U in the Rhines scale formula (24). In our simulations, we find that the convective lengthscale approximately follows a power law l_c ∼ U^0.28 when U is interpreted as the r.m.s. total velocity. The power law exponent is therefore further away from the Rhines scale prediction when using the total velocity rather than the convective velocity.

3.3 Scaling of the zonal flow velocity

To complete this section, we now discuss how the amplitude of the zonal velocity scales with the convective velocity. Fig. 5(a) shows the evolution of the zonal Rossby number, Ro₀ = Re₀Ek, as a function of the convective Rossby number, Ro_c = Re_cEk, for Ek ∈ [10⁻⁸, 10⁻⁷] and Pr = [10⁻², 10⁻¹]. In this section we use the Rossby numbers to clarify the dependence of the velocity on the viscosity. For Ek Pr ≤ 10⁻⁹, the nonlinear onset of convection is subcritical (Guervilly & Cardin 2016). In this case, the zonal and convective Reynolds numbers are discontinuous and Re_c ≳ 10³. For each fixed value of Ek, the data points approximately fall on a straight line, indicating a power law dependence on Ro_c. The dashed line represents Ro₀ = Ro_c. The data points cross this line for moderate values of the Rayleigh numbers that depend on the Prandtl number: at Ra ≈ 1.5Ra_c for Pr = 10⁻² and Ra ≈ 6Ra_c for Pr = 10⁻¹. The zonal flows have therefore large amplitude for relatively modest values of the Rayleigh numbers when Pr < 1, despite the presence of the Ekman boundary friction in our model. They reach an amplitude comparable to the amplitude of the convective flows for smaller Ra as Pr is decreased. In this sense, lower Pr is favourable for the zonal flows.

(a) Zonal Rossby number Ro0 as a function of the convective Rossby number Roc. The dotted line shows the power law whose exponent best fits all the data points. Ro0 compensated by the power laws of (b) eq. (26) and (c) eq. (27) as a function of Roc.

Figure 5.

(a) Zonal Rossby number Ro₀ as a function of the convective Rossby number Ro_c. The dotted line shows the power law whose exponent best fits all the data points. Ro₀ compensated by the power laws of (b) eq. (26) and (c) eq. (27) as a function of Ro_c.

Open in new tab Download slide

A scaling law for the zonal flow amplitude can be obtained by considering that the dominant force balance in the zonal velocity eq. (18) is established between the nonlinear interactions of the convective velocities and the friction in the Ekman layer. This is consistent with measurements of the contributions to the zonal energy budget in our simulations. For instance, in the case Ek = 10⁻⁸, Pr = 10⁻¹ and Ra/Ra_c = 9.02, we measured that the boundary friction accounts for 90 per cent of the total viscous dissipation of the axisymmetric flow. In terms of scaling arguments, this balance implies

\begin{equation} \mbox{{Ro}}_0\sim \mbox{{Ek}}^{-1/2} \frac{\mbox{{Ro}}_c^2}{l_c} . \end{equation}

(25)

Using the scaling |$l_c\sim \mbox{{Ro}}_c^{0.37}$| deduced from our data in Section3.2, we obtain the scaling |$\mbox{{Ro}}_0\sim \mbox{{Ro}}_c^{1.63}$|⁠. The best fit to all the data points in Fig. 5(a) is |$\mbox{{Ro}}_0\sim \mbox{{Ro}}_c^{1.64(\pm 0.04)}$|⁠. This agreement is not surprizing because all the global quantities are calculated from the same data set, but it shows that our measurements of the global values of the convective lengthscales and the velocities are consistent and that eq. (25) is suitable to scale the zonal flow velocity.

It is interesting to compare the observed scaling of the zonal velocity with predictive scaling laws based on physical arguments that have been derived in the literature (e.g. Aubert et al. 2001). To do so, we need to examine possible scaling arguments for l_c. Considering first that l_c scales with the Rhines scale (24), |$l_c\sim \mbox{{Ro}}_c^{1/2}$|⁠, gives

\begin{equation} \mbox{{Ro}}_0\sim \mbox{{Ek}}^{-1/2} \mbox{{Ro}}_c^{3/2}. \end{equation}

(26)

The exponent 3/2 is in good agreement with the results of the 3-D numerical simulations of Kaplan et al. (2017) for Ek ∈ [10⁻⁷, 10⁻⁶] and Pr ∈ [3 × 10⁻³, 10⁻¹].

Gillet et al. (2007) proposed an alternative scaling for l_c obtained by using the zonal velocity as a typical flow velocity in the Rhines scale. This gives |$l_c\sim \mbox{{Ro}}_0^{1/2}$|⁠, and so,

\begin{equation} \mbox{{Ro}}_0\sim \mbox{{Ek}}^{-1/3} \mbox{{Ro}}_c^{4/3}. \end{equation}

(27)

Gillet et al. (2007) found that this scaling law provides a good fit for their numerical data obtained with a QG model at |$\mbox{{Ek}}=\mathcal {O}(10^{-6})$| and Pr = 0.025.

Finally, we can consider l_c ∼ Ek^1/3, which corresponds to the scaling of the azimuthal lengthscale of the linear convective instability (Jones et al. 2000). This estimate is obviously not satisfactory because the convective lengthscale must increase with the Rayleigh number as shown in Figs 1 and 4. For this estimate, we deduce the scaling law

\begin{equation} \mbox{{Ro}}_0\sim \mbox{{Ek}}^{-5/6} \mbox{{Ro}}_c^2. \end{equation}

(28)

The power law (26) provides the closest exponent to our data best fit and is based on plausible physical arguments. By contrast, the power law (28) is too steep and the dependence of the pre-factor on Ek is much weaker in the data. The power law (27) requires a much stronger dependence of the convective lengthscale on Ro_c (namely |$l_c\sim \mbox{{Ro}}_c^{2/3}$|⁠) than observed. Figs 5(b) and (c) show the dependence of Ro₀ compensated by the power laws (26) and (27), respectively, on Ro_c. The data points compensated by the power law (26) align on a plateau, confirming that the exponent of this power law provides a reasonable agreement with our data.

The pre-factor of the power law (26) predicts that a decrease of Ek by a decade should lead to an increase of Ro₀ by approximately a factor 3. The dependence on the Ekman number cannot be estimated accurately from our data points as they only sample one decade of Ek. Fig. 5(b) tentatively suggests that the power law (26) does not entirely explain the dependence on the Ekman number: in our simulations, Ro₀ only increases by approximately a factor 2 when Ek decreases from 10⁻⁷ to 10⁻⁸.

4 DRIFT AND STABILITY OF THE ZONAL JETS

In the rest of the paper, we only consider simulations run at Ek = 10⁻⁸ and Pr = 10⁻¹, where we obtain the largest number of zonal jets. The zonal flows have a strong influence on the convective flows, and hence on the heat transport as we will discuss in Section 6, so it is important to discuss the persistence of the zonal jets. In this section, we examine the stability of the jets.

Using a QG model of thermal convection, Rotvig (2007) shows that zonal jets drift inwards provided that the slope β has a significant dependence on radius, while QG models with constant β produce multiple jets that do not drift (e.g. Jones et al. 2003). Rotvig showed that the drift also occurs in a 3-D spherical model, so this effect is not restricted to QG models. The drift rate is found to increase with β and with the Rayleigh number. The drift of zonal flows is also observed in the rotating turntable experiment of Smith et al. (2014). In the experiment, β is positive and the drift is observed outwards. These studies thus indicate that the direction of the zonal jet drift is related to the sign of β.

Fig. 6 shows the space–time diagram of the zonal velocity for Ra/Ra_c = 0.96 and Ra/Ra_c = 9.02. For Ra = 0.96, the middle retrograde zonal jet drifts inwards periodically. In our system β < 0 so this inwards migration of the zonal flows is in agreement with the work of Rotvig (2007) and Smith et al. (2014).

Space–time diagram of the zonal velocity for (a) Ra/Rac = 0.96 and (b) Ra/Rac = 9.02 (Ek = 10−8 and Pr = 10−1). The time interval corresponds to (a) 106 rotation periods and 317 convective turnover timescales (as defined in Table 1) and (b) 1.2 × 106 rotation periods and 2672 convective turnover timescales.

Figure 6.

Space–time diagram of the zonal velocity for (a) Ra/Ra_c = 0.96 and (b) Ra/Ra_c = 9.02 (Ek = 10⁻⁸ and Pr = 10⁻¹). The time interval corresponds to (a) 10⁶ rotation periods and 317 convective turnover timescales (as defined in Table 1) and (b) 1.2 × 10⁶ rotation periods and 2672 convective turnover timescales.

Open in new tab Download slide

For Ra/Ra_c = 9.02, the picture is completely different. For s > 0.4 the zonal jets do not drift and the standard deviation of their amplitude is of the order of 2000 (compared with a mean amplitude of approximately 30 000) over the course of the simulation. However at radius s < 0.4, the central prograde jet and its retrograde neighbour drift outwards. The retrograde jet initially forms around s = 0.2 and eventually merges with the retrograde jet located at s ≈ 0.35, closing off the prograde jet in the process. This sequence is not quite periodic and takes between 5 and 20 zonal turnover timescale (based on the time- and volume-averaged zonal velocity). The central region of the equatorial plane has the smallest values of β and dβ/ds, so the direction and the location of the drift indicates that this mechanism is different from the inward drift observed for smaller Ra and large β.

The drift might be related to instabilities of the prograde jets near the centre, where β goes to zero. The Rayleigh-Kuo criterion states that a necessary condition for the barotropic instability of a shear flow in a inviscid Boussinesq fluid is that the quantity |$\Delta =2\beta \mbox{{Ek}}^{-1}-d\overline{\zeta }/ds$|⁠, where |$\overline{\zeta }=d\overline{u_{\phi }}/ds+\overline{u_{\phi }}/s$|⁠, changes sign at some radius (Kuo 1949). This indicates that large slopes have a stabilizing effect on the zonal flow, so it is plausible that prograde zonal flows near the centre are unstable. The criterion is valid for an inviscid fluid, so it only provides an indication of the zonal flow stability for small Ekman numbers. Nevertheless, Guervilly et al. (2012) showed that the threshold of the instability obtained with numerical simulations is in good quantitative agreement with the Rayleigh–Kuo criterion for Ek < 10⁻⁷. As well as calculating the quantity Δ in the time series of Ra/Ra_c = 9.02, we also compute the product of |$\overline{u_{\phi }}$| with the zonal average of the eddy momentum flux convergence given in the zonal velocity eq. (18),

\begin{equation} R = - \left(\overline{u^{\prime }_{s} \frac{\partial u^{\prime }_{\phi }}{\partial s}} + \overline{\frac{u^{\prime }_{s} u^{\prime }_{\phi }}{s}} \right)\overline{u_{\phi }}, \end{equation}

(29)

where the prime denotes the non-axisymmetric velocity component. R is the only source of energy of the zonal flow, so we expect this term to be positive when the eddies feed energy into the zonal velocity. Fig. 7 shows the space–time diagram of R (colour), where the isocontours of |$\overline{u_{\phi }}=0$| (black line) and Δ = 0 (green line) have been superposed. In the central region, R is negative in the inward flank of the drifting prograde jet, meaning that the zonal flows lose energy to the non-axisymmetric flows there. The instability criteria (i.e. Δ > 0) is violated near the centre where β is small and in a few places in the drifting prograde jet. This indicates that this prograde jet might be marginally stable, leading to a transfer of energy from the zonal flow to non-axisymmetric flows. As the zonal flow loses energy in the inward flank of the prograde jet, we expect the fluid parcels located there to move outwards to conserve their angular momentum. This mechanism could explain why the prograde jet moves outwards, and by doing so, pushes the neighbouring retrograde jet outwards.

$Space–time diagram of R (colour) for Ra/Rac = 9.02 (Ek = 10−8, Pr = 10−1). The isocontours $\overline{u_{\phi }}=0$ and Δ = 0 have been represented in black and green, respectively. The time interval corresponds to 3.5 × 104 rotation periods and 78 convective turnover timescales (as defined in Table 1).$

Figure 7.

Space–time diagram of R (colour) for Ra/Ra_c = 9.02 (Ek = 10⁻⁸, Pr = 10⁻¹). The isocontours |$\overline{u_{\phi }}=0$| and Δ = 0 have been represented in black and green, respectively. The time interval corresponds to 3.5 × 10⁴ rotation periods and 78 convective turnover timescales (as defined in Table 1).

Open in new tab Download slide

5 DYNAMICAL DIFFERENCE BETWEEN RETROGRADE AND PROGRADE JETS

We now discuss the mechanism of formation of persistent zonal flows based on the extensive literature on the subject, particularly in the context of the ocean and atmosphere dynamics (e.g. Vallis 2006). In our simulations, retrograde zonal flows are faster and sharper than the rounded prograde zonal flows. This asymmetry is indicative of an important dynamical difference between the two types of zonal jets, related to their formation mechanism and to the sign of β. The asymmetry is also observed in 3-D models of rotating spherical convection (e.g. Heimpel & Aurnou 2007; Gastine et al. 2014) and in QG models with constant β (e.g. Teed et al. 2012). In the 3-D simulations, the multiple zonal jets emerge inside the tangent cylinder, where β is positive, so the prograde flows are faster and sharper than the retrograde flows.

In QG flows, the zonal flows are directly related to the distribution of potential vorticity (PV; e.g. McIntyre 2003). In our system, the equation for the evolution of the PV, q, is

\begin{equation} \frac{\mbox{$\partial $} q}{\mbox{$\partial $} t} + \boldsymbol{u} \cdot \nabla q= D + F , \end{equation}

(30)

where

\begin{equation} q=\frac{\zeta + 2\mbox{{Ek}}^{-1}}{H}, \end{equation}

(31)

D represents the dissipation terms from the boundary and the bulk and F the buoyancy source. ζ includes the axisymmetric and non-axisymmetric components of the vorticity. In the absence of buoyancy sources and dissipative effects, q is a materially invariant scalar so it tends to be locally homogenized by the turbulence. The zonal flows appear as a consequence of the mixing of PV with

\begin{equation} \overline{u_{\phi }}(s) = \frac{1}{s} \int _0^s s^{\prime } \overline{\zeta } \text{d}s^{\prime } = \frac{1}{s} \int _0^s s^{\prime } \left(H \overline{q} - 2 \mbox{{Ek}}^{-1} \right)\text{d}s^{\prime } . \end{equation}

(32)

In a sphere (where H decreases with s, i.e. β < 0), regions of weak PV gradients correspond to prograde zonal jets, while regions of strong PV gradient correspond to retrograde zonal jets. The PV distribution controls the QG dynamics and, in particular, the propagation of Rossby waves. The equation for linear Rossby waves is

\begin{equation} \frac{\partial \zeta ^{\prime }}{\partial t} + u_s H \frac{\text{d} \overline{q}}{\text{d}s} = 0, \end{equation}

(33)

where |$\zeta ^{\prime }=\zeta -\overline{\zeta }$| and

\begin{equation} \frac{\text{d}\overline{q}}{\text{d}s} = \frac{\text{d}\overline{\zeta }/H}{\text{d}s} - \frac{2\beta }{\mbox{{Ek}}H}. \end{equation}

(34)

The restoring mechanism of the Rossby waves is therefore stronger in the retrograde jets than in the prograde jets. A strong mixing leads to the formation of a staircase of PV, that is, a succession of regions of homogeneous PV separated by regions of steep PV gradients (e.g. Scott & Dritschel 2012). This occurs because perturbations to strong PV gradients are radiated as Rossby waves, an effect referred to as Rossby wave elasticity (McIntyre 2008). Perturbations are therefore inhibited in the regions of strong PV gradients, while they are intensified in regions of weak PV gradients. This leads to a feedback mechanism that further steepens the gradients of PV (Dritschel & McIntyre 2008). In a sphere, the relation between PV and zonal flows (32) implies that wide regions of weak PV gradients appear as large rounded prograde jets, whereas narrow regions of steep PV gradients appear as sharp retrograde jets.

The radial profile of the axisymmetric PV is plotted in Fig. 8 for Ra/Ra_c = 9.02. In the range s ∈ [0, 0.8], the succession of weak and steep PV gradient regions is visible and forms a relatively mild PV staircase for this value of Ra. Larger Ra are required to obtain well mixed PV regions and a better formed staircase profile. The asymmetry of the jets is clearly visible in our simulations but retrograde jets are not particularly narrower than the prograde jets. This is because the chosen delimitation of the jets (⁠|$\overline{u_{\phi }}=0$|⁠) is, somewhat arbitrarily, defined with respect to the planetary rotation. The profile of |$\overline{q}$| shows that the regions of steep PV gradients (which correspond to the cores of the retrograde jets) are in fact much narrower than the regions of weak PV gradients. Another characteristic of the zonal flow is the robustness of the prograde jet at the centre in all our simulations: this is well explained by the mixing in the central region that leads to a local increase of q. The PV mixing mechanism therefore provides a good explanation for the most notable features of the zonal flows observed in our simulations. Nevertheless, we note that alternative–although not mutually exclusive–mechanisms for the formation of persistent zonal flows have been put forward in the literature, such as resonant triad interactions (e.g. Pedlosky 1987).

Figure 8.

Radial profile of (a) the zonal velocity and (b) the axisymmetric potential vorticity (time averages) for Ra/Ra_c = 9.02 (Ek = 10⁻⁸, Pr = 10⁻¹). The grey bands correspond to the regions where the zonal flow is retrograde.

Open in new tab Download slide

The PV distribution indicates that the non-axisymmetric dynamics inside the retrograde jets might be dominated by Rossby waves. Assuming for simplicity |$\zeta ^{\prime }=-\nabla _e^2\psi$| (see equation 12), the dispersion relation of the Rossby waves is

\begin{equation} \omega = H \frac{\text{d}\overline{q}}{\text{d}s} \frac{k_{\phi }}{|\boldsymbol{k}|^2}, \end{equation}

(35)

where ω is the frequency of the wave and |$\boldsymbol{k}=(k_s,k_{\phi })$| is the wavenumber vector. Consequently, the azimuthal phase velocity, v_p, and azimuthal group velocity, v_g, of the Rossby waves are

\begin{equation} v_p = H \frac{\text{d}\overline{q}}{\text{d}s} \frac{1}{|\boldsymbol{k}|^2}, \quad v_g = H \frac{\text{d}\overline{q}}{\text{d}s} \frac{k_s^2-k_{\phi }^2}{|\boldsymbol{k}|^4}. \end{equation}

(36)

In our system, the gradient of |$\overline{q}$| is positive so the Rossby waves always have a positive azimuthal phase speed. To determine whether Rossby waves are present in the retrograde jets, we can track the direction of the azimuthal drift of the velocity patterns. Fig. 9 shows the Hovmöller map (longitude–time) of the radial velocity u_s at a fixed radius in a prograde jet (⁠|$\overline{u_{\phi }}=11600$|⁠) and in a retrograde jet (⁠|$\overline{u_{\phi }}=-24600$|⁠). The solid black line represents the azimuthal drift due to the advection by the zonal velocity at this radius, |$\Delta t= s\Delta \phi /\overline{u_{\phi }}$|⁠. In the prograde jet, the radial velocity structures move in the prograde direction at a rate that is consistent with the advection by the zonal velocity, and even faster for some structures. In the retrograde jet, the radial velocity patterns mainly move in the prograde direction. These patterns must therefore correspond to Rossby waves. The amplitude of the patterns of u_s is modulated as they move in the prograde direction. These modulations appear on neighbouring patterns and seem to travel in the retrograde direction. This observation is consistent with Rossby waves for which k_s < k_ϕ so their azimuthal group velocity is negative. In the frame of reference (i.e. rotating at the rotation rate Ω), the Rossby waves move at a velocity |$v_p+\overline{u_{\phi }}$|⁠. From Fig. 9(b), we estimate that this velocity is approximately 20 000, which gives v_p ≈ 45 000. By using the local value of the gradient of |$\overline{q}$| at this radius, we find that |$|\boldsymbol{k}|^2\approx 84$| from eq. (36). From Fig. 9(b) we estimate that the azimuthal wavenumber is k_ϕ = m/s ≈ 70, and so, we deduce that k_s ≈ 46. This result is consistent with a negative azimuthal group velocity.

$Hovmöller map (longitude–time) of the radial velocity us at a fixed radius: (a) s = 0.51 in a prograde jet ($\overline{u_{\phi }}=11600$) and (b) s = 0.58 in a retrograde jet ($\overline{u_{\phi }}=-24600$) for Ra/Rac = 6.44 (Ek = 10−8, Pr = 10−1). The solid black line shows the drift due to the advection by the zonal velocity. The time interval corresponds to 2000 rotation periods and 4 convective turnover timescales (as defined in Table 1).$

Figure 9.

Hovmöller map (longitude–time) of the radial velocity u_s at a fixed radius: (a) s = 0.51 in a prograde jet (⁠|$\overline{u_{\phi }}=11600$|⁠) and (b) s = 0.58 in a retrograde jet (⁠|$\overline{u_{\phi }}=-24600$|⁠) for Ra/Ra_c = 6.44 (Ek = 10⁻⁸, Pr = 10⁻¹). The solid black line shows the drift due to the advection by the zonal velocity. The time interval corresponds to 2000 rotation periods and 4 convective turnover timescales (as defined in Table 1).

Open in new tab Download slide

We might expect that the dynamical difference between prograde and retrograde jets affects the temperature distribution because the Rossby waves might modify the transport properties of the flow. We study this problem in the next section.

6 EFFECT OF THE ZONAL FLOWS ON THE HEAT TRANSPORT

First, we can assess the effect of the zonal flow on the amplitude of the radial flow by comparing the profile of the zonal velocity to the radial profile of the r.m.s radial velocity, |$u_s^{\ast }$|⁠, which is calculated as

\begin{equation} u_s^{\ast }(s) = \frac{1}{\Delta t} \int _{\Delta t} \left(\frac{1}{2\pi s} \int _0^{2\pi } u_{s}^2(s,\phi ,t) s \text{d}\phi \right)^{1/2} \text{d}t. \end{equation}

(37)

Fig. 10(a) shows these profiles for Ra/Ra_c = 9.02. The zonal flow (black line) is plotted according to the left axis and |$u_s^{\ast }$| (blue line) is plotted according to the right axis. Both profiles are time-averaged. In the inner convective region (s ≲ 0.8), the local maxima of |$u_s^{\ast }$| correlates well with the extrema of |$\overline{u_{\phi }}$| (the cores of the jets), that is, the zeros of the radial shear |$|\partial _s \overline{u_{\phi }}|$|⁠, while the local minima of |$u_s^{\ast }$| correlates with maxima of |$|\partial _s \overline{u_{\phi }}|$| (the flanks of the jets). The reduction of |$u_s^{\ast }$| in the flank of a jet can reach 30 per cent of its value in the core of the neighbouring jets. The radial velocity is therefore impeded by the radial shear and maximized in the cores of jets, irrespective of their sign. Thus the dynamical difference between prograde and retrograde jets cannot be directly diagnosed on this profile.

Figure 10.

Radial profile of (a) the zonal velocity (black line) and the r.m.s radial velocity (blue), (b) the axisymmetric temperature along the rotation axis (black) and in the equatorial plane (blue), (c) the static heat flux (black), the conductive heat flux in the equatorial plane (blue) and the convective heat flux in the equatorial plane (red). All the quantities have been time-averaged. The grey bands correspond to the regions where the zonal flow is retrograde. The parameters are Ra/Ra_c = 9.02, Ek = 10⁻⁸ and Pr = 10⁻¹.

Open in new tab Download slide

To determine how the temperature field is affected by the zonal flows, Fig. 11(a) shows a meridional slice of the axisymmetric temperature, |$\overline{T}=\overline{\Theta }+T_s$|⁠, averaged in time for the same simulation. The isotherms have an ellipsoidal shape which is elongated towards the equator. To complement this figure, the radial profiles of the axisymmetric temperature along the rotation axis, |$\overline{T}(\theta =0)$|⁠, and in the equatorial plane, |$\overline{T}(\theta =\pi /2)$|⁠, are plotted in Fig. 10(b). In the equatorial plane, the temperature has a flatter profile than along the rotation axis for s < 0.8, which explains the ellipsoidal shape of the isotherms. The thermal boundary layer is more pronounced in the equatorial plane. Nevertheless it remains thick (much thicker than the Ekman layer) because Pr < 1 and the Rayleigh number is moderate. The disparity between the temperature along the axial and equatorial directions is due to preferential direction of the convective flows in rapidly rotating convection. This effect is also observed in 3-D models (e.g. Zhang 1991; Yadav et al. 2016), but is amplified here by the use of a QG model.

$Meridional cross-section of the axisymmetric average of (a) the temperature $\overline{T}=\overline{\Theta }+T_s$, (b) the convective heat flux $\overline{F}_{\text{cv}}$, and (c) the conductive heat flux $\overline{F}_{\text{cd}}$. All the quantities have been time-averaged. Same parameters as Fig. 10.$

Figure 11.

Meridional cross-section of the axisymmetric average of (a) the temperature |$\overline{T}=\overline{\Theta }+T_s$|⁠, (b) the convective heat flux |$\overline{F}_{\text{cv}}$|⁠, and (c) the conductive heat flux |$\overline{F}_{\text{cd}}$|⁠. All the quantities have been time-averaged. Same parameters as Fig. 10.

Open in new tab Download slide

On top of their ellipsoidal shape, the isotherms also have undulations of small amplitude. These undulations are located at the same distance from the rotation axis for each isotherm so they are likely due to the presence of zonal flows. This causal link is visible in Fig. 10(b) where we use the radial profile of |$\overline{T}(\theta =\pi /2)$| as a proxy for the z-averaged axisymmetric temperature. The temperature profile is relatively flat in the prograde zonal jets. By contrast it is significantly steeper in the core of the retrograde jets. This indicates that the mean temperature is affected by the PV distribution: in the core of the retrograde jets, the inhibition of the turbulence and the dominance of the Rossby waves imply that the temperature cannot be efficiently homogenized.

This effect can be quantified by measuring the heat fluxes carried by convection and conduction through the system. In the steady state, the volume average of the heat equation implies that

\begin{equation} \frac{1}{\mathcal {S}} \int _{\mathcal {S}} \left(F_{\text{cv}} + F_{\text{cd}} \right)\text{d}\mathcal {S} = F_s, \end{equation}

(38)

where the convective heat flux is

\begin{equation} F_{\text{cv}} = \Theta u_r, \end{equation}

(39)

the conductive heat flux,

\begin{equation} F_{\text{cd}} = - \frac{1}{\mbox{{Pr}}} \frac{\partial T}{\partial r}, \end{equation}

(40)

the static heat flux,

\begin{equation} F_s = - \frac{1}{\mbox{{Pr}}} \frac{\text{d}T_s}{\text{d}r} = \frac{2}{\mbox{{Pr}}^2} r, \end{equation}

(41)

and |$\mathcal {S}$| is a spherical surface. Fig. 11 shows the meridional slices of the axisymmetric averages of the convective heat flux, |$\overline{F}_{\text{cv}}$|⁠, and of the conductive heat flux, |$\overline{F}_{\text{cd}}$|⁠. The fluxes are time-averaged. Both fluxes have a banded structure aligned with z. The convective heat flux is mainly concentrated around the equatorial plane and is weak near the axis. The conductive heat flux is maximum near the equator. For this simulation at Pr = 10⁻¹, the conduction carries a larger part of the heat than the convection in most of the domain, despite the large values of the radial velocity (of the order of 10⁴).

The thermal contrast between prograde and retrograde jets can be clearly observed in the profiles of the axisymmetric heat fluxes in the equatorial plane shown in Fig. 10(c). For comparison, the static heat flux F_s is also shown. Note that the profile of |$\overline{F}_{\text{cv}}$| and |$\overline{F}_{\text{cd}}$| at a given latitude are not representative of the spherical averages so their sum is not equal to F_s at each radius. The local maxima of the convective flux are located in the prograde zonal jets. The decrease of the convective flux matches the decrease of |$u_s^{\ast }$| in the shear layers. However the convective flux systematically reaches minima in the core of the retrograde jets while the radial velocity recovers there. The conductive flux is boosted in the shear layers, but even more so in the retrograde jets, where it is much larger than the convective flux. The balance between the thermal processes is very different in the prograde and retrograde jets: the prograde jets are regions where about half of the heat is carried by the convection and the temperature is fairly uniform, whereas most of the heat is carried by the conductive flux in the retrograde jets. This occurs despite high values of the r.m.s radial velocity in the retrograde jets, so the radial flow must not be well correlated with the temperature perturbation in these regions. Fig. 12 show the Hovmöller map for the non-axisymmetric temperature perturbation, |$\Theta ^{\prime }=\Theta -\overline{\Theta }$|⁠, for the same parameters and radius as Fig. 9. In both prograde and retrograde jets, the temperature perturbation drifts in the direction of the zonal flow with a drift rate consistent with the advection by the zonal velocity. In the core of the retrograde jets, which are dominated by the propagation of Rossby waves, the radial velocity and temperature perturbation are visibly not well correlated. The inefficiency of the Rossby waves at transporting heat outwards explains the weak convective heat flux there. In the prograde jet, Θ^΄ is well correlated with u_s, which is consistent with the high convective heat flux.

Figure 12.

Hovmöller map of the non-axisymmetric temperature perturbation, Θ^΄, at a fixed radius: (a) s = 0.51 in a prograde jet and (b) s = 0.58 in a retrograde jet for the same parameters as Fig. 9. The solid black line shows the drift due to the advection by the zonal velocity. The time interval corresponds to 2000 rotation periods and 4 convective turnover timescales (as defined in Table 1).

Open in new tab Download slide

In summary, our results show that the core of the retrograde jets—and not only the regions of intense shear—act as primary bottlenecks to the convective heat transport.

7 DISCUSSION

In this paper, we have studied the convective structures and zonal flows that form in rotating thermal convection for values of the Prandtl number relevant for liquid metals (⁠|$\mbox{{Pr}}=\mathcal {O}(10^{-1})$|⁠) and low Ekman (⁠|$\mbox{{Ek}}=\mathcal {O}(10^{-8})$|⁠) and Rossby numbers (Ro_c < 10⁻²). In order to reach low values of the Ekman number, we have used a hybrid numerical model that couples a QG approximation for the velocity to a 3-D temperature field. Convection is driven by internal heating in a full sphere geometry. The model includes Ekman pumping to mimic no-slip boundary conditions. We focus on the intense zonal flows that emerge on the strong branch of convection, which was described by Guervilly & Cardin (2016) using the same hybrid QG-3D model and by Kaplan et al. (2017) in a fully 3-D model. Persistent multiple zonal jets of alternating sign form due to the mixing of PV and exert a strong feedback on the convection. An upscale energy transfer takes place and the integral convective lengthscale increases with the vigour of the convection. The convective lengthscale and the zonal jet width are closely related: the radial shear exerted by the zonal flow on the radial velocity limits the size of the convective eddies, while the typical mixing length of the PV depends on the size of the most energetic convective eddies. The convective lengthscale varies radially in agreement with the Rhines scale (Rhines 1975), that is, as the square root of β, where β measures the slope of the boundaries. However the convective scale increases more slowly with the convective speed (following a power law of exponent 0.37) than predicted by the Rhines scale (power law of exponent 0.5).

In our QG model, convection carries heat mostly in the direction perpendicular to the rotation axis. We have shown that the principal barrier to this convective heat transport is located in the cores of the retrograde zonal jets. This is due to the formation of a staircase of PV: steep and weak gradients of the PV correspond to retrograde and prograde zonal jets, respectively. The steep PV gradients inhibit the eddies and favour the propagation of Rossby waves, which are inefficient at carrying heat outwards. The occurrence of eddy-transport barriers associated with strong zonal jets is well-documented in the atmospheric dynamics context (Dritschel & McIntyre 2008). In our simulations, these barriers lead to the steepening of the mean temperature gradient in the core of the retrograde jets, and there, the heat is largely carried by conduction. The unfavourable effect of the shear layer in the flanks of the zonal jets on the convective heat transport is secondary by comparison. To illustrate the thermal signature of the retrograde jets at the surface, we plot in Fig. 13 the latitudinal profile of the axisymmetric heat flux at r = 1. The most noticeable feature is that the heat flux is maximal at the equator, which is expected as the convective transport is largely perpendicular to the rotation axis, and is also observed in 3-D models (e.g. Zhang 1991; Yadav et al. 2016). This enhanced heat transport in the equatorial regions compared with the polar regions is used in a number of models of the Earth's inner core growth to explain the observed seismic anisotropy (e.g. Yoshida et al. 1996; Deguen & Cardin 2009). Of greater interest here are the more subtle variations of the surface heat flux at higher latitudes: the cores of the retrograde jets are characterized by local maxima of the (conductive) surface heat flux. This occurs because the axisymmetric temperature gradient is steepest in these regions. To highlight the small-scale anomalies of the surface heat flux and quantify their amplitude, we filter out the coefficients of the spherical harmonics of degree l smaller than 30. The filtered profile is plotted according to the right axis of Fig. 13 in blue. It shows that the heat flux anomalies due to the presence of the zonal flows are sharp and of a small amplitude, approximately 1 per cent of the mean surface heat flux. The thermal signal associated with zonal jets provides useful information in the context of the gas giant planets because the power emitted at the surface of the planet can be measured. For Jupiter and Saturn, the emitted power is approximately uniform in latitude with variations of small amplitude that resemble the structure of the zonal flows at the surface (Pirraglia 1984; Li et al. 2010).

Figure 13.

Latitudinal profile of the axisymmetric heat flux (black line) from the North pole to the equator at the outer boundary r = 1 for Ek = 10⁻⁸, Pr = 10⁻¹ and Ra/Ra_c = 9.02. The blue line shows the axisymmetric heat flux where the spherical harmonics coefficient of degrees l ≤ 30 have been filtered out and is plotted according to the right axis. The grey bands correspond to the regions where the zonal flow is retrograde.

Open in new tab Download slide

For increasing thermal driving (out of reach of our current computational resources), we expect the convective eddies to further enhance the mixing of PV in the prograde jets. The PV staircase would then become sharper with narrower regions of steep PV gradients, and hence narrower retrograde jets. This sharpening of the PV staircase in rapidly rotating convection is observed by Verhoeven & Stellmach (2014) at large Rayleigh numbers (20 times supercritical) in 2-D numerical simulations using the anelastic approximation. In their study, the β effect is due to the compressibility of the fluid (Ingersoll & Pollard 1982; Glatzmaier et al. 2009) and leads to the formation of multiple zonal jets, similarly to the topographic β effect as studied here. They show that the sharpening of the PV staircase results in the sharpening of the entropy staircase (analogous to the temperature staircase in the Boussinesq approximation). In QG systems with vigorous convection, we thus expect the heat transport process to be heterogeneous with wide convective regions separated by narrow conducting bands; the efficiency of the heat transfer throughout the whole system would largely be controlled by the efficiency of the conducting process across the retrograde jets, and thus, by the width of these conducting bands. This process is partly analogous to the occurrence of layering in double-diffusive convection where the heat transfer is controlled by the flux through the interface between overturning layers (Turner 1985).

The main features of the nonlinear QG dynamics discussed in this paper do not crucially rely on the temperature being 3-D. Consequently, we expect that QG models of rotating convection using a 2-D temperature field would be able to reproduce our observations qualitatively. QG-2D models could be used to pursue this study at lower Ekman numbers and larger Rayleigh numbers. The numerical framework used in the hybrid QG-3D model is however well suited to explore the possible existence of dynamos driven by QG flows at low magnetic Prandtl numbers (e.g. Gillet et al. 2011). The dynamo problem indeed requires to treat the magnetic field in 3-D. We shall investigate QG dynamos in a forthcoming study.

Acknowledgements

CG was supported by the Natural Environment Research Council under grant NE/M017893/1. PC acknowledges the Agence Nationale de la Recherche for supporting this project under grant ANR TuDy. This work was undertaken on TOPSY of the HPC facilities at Newcastle University and on the facilities of N8 HPC Centre of Excellence, provided and funded by the N8 consortium and EPSRC (Grant EP/K000225/1) and coordinated by the Universities of Leeds and Manchester. We are grateful to the referees for helpful comments that improved the manuscript.

REFERENCES

Aubert

J.

,

Brito

D.

,

Nataf

H.-C.

,

Cardin

P.

,

Masson

J.-P.

,

2001

.

A systematic experimental study of rapidly rotating spherical convection in water and liquid gallium

,

Phys. Earth planet. Inter.

128

51

–

74

.

Aubert

J.

Gillet

N.

Cardin

P.

,

2003

.

Quasigeostrophic models of convection in rotating spherical shells

,

Geochem. Geophys. Geosyst.

4

(

7

),

1052, doi:10.1029/2002GC000456

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Aurnou

J.M.

,

Olson

P.L.

,

2001

.

Strong zonal winds from thermal convection in a rotating spherical shell

,

Geophys. Res. Lett.

28

2557

–

2560

.

Aurnou

J.

,

Heimpel

M.

,

Allen

L.

,

King

E.

,

Wicht

J.

,

2008

.

Convective heat transfer and the pattern of thermal emission on the gas giants

,

Geophys. J. Int.

173

(

3

),

793

–

801

.

Busse

F.H.

,

1970

.

Thermal instabilities in rapidly rotating systems

,

J. Fluid Mech.

44

441

–

460

.

Busse

F.H.

,

Hood

L.L.

,

1982

.

Differential rotation driven by convection in a rapidly rotating annulus

,

Geophys. Astrophys. Fluid Dyn.

21

(

1–2

),

59

–

74

.

Busse

F.H.

Or

A.C.

,

1986

.

Convection in a rotating cylindrical annulus: thermal Rossby waves

,

J. Fluid Mech.

166

173

–

187

.

Cabanes

S.

Aurnou

J.

Favier

B.

Le Bars

M.

,

2017

.

A laboratory model for deep-seated jets on the gas giants

,

Nat. Phys.

13

387

–

390

.

Calkins

M.

,

Aurnou

J.

,

Eldredge

J.

,

Julien

K.

,

2012

.

The influence of fluid properties on the morphology of core turbulence and the geomagnetic field

,

Earth planet. Sci. Lett.

359

55

–

60

.

Calkins

M.A.

,

Julien

K.

,

Marti

P.

,

2013

.

Three-dimensional quasi-geostrophic convection in the rotating cylindrical annulus with steeply sloping endwalls

,

J. Fluid Mech.

732

214

–

244

.

Cardin

P.

,

Olson

P.

,

1994

.

Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core

,

Phys. Earth planet. Inter.

82

235

–

259

.

Christensen

U.R.

,

2002

.

Zonal flow driven by strongly supercritical convection in rotating spherical shells

,

J. Fluid Mech.

470

115

–

133

.

Davidson

P.A.

,

2013

.

Turbulence in Rotating, Stratified and Electrically Conducting Fluids

Cambridge Univ. Press

.

Deguen

R.

,

Cardin

P.

,

2009

.

Tectonic history of the Earth's inner core preserved in its seismic structure

,

Nat. Geosci.

2

(

6

),

419

–

422

.

Dritschel

D.G.

,

McIntyre

M.E.

,

2008

.

Multiple jets as PV staircases: The Phillips effect and the resilience of eddy-transport barriers

,

J. Atmos. Sci.

65

855

–

874

.

Gastine

T.

,

Heimpel

M.

,

Wicht

J.

,

2014

.

Zonal flow scaling in rapidly-rotating compressible convection

,

Phys. Earth planet. Inter.

232

36

–

50

.

Gastine

T.

,

Wicht

J.

,

Aubert

J.

,

2016

.

Scaling regimes in spherical shell rotating convection

,

J. Fluid Mech.

808

690

–

732

.

Gillet

N.

,

Jones

C.A.

,

2006

.

The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder

,

J. Fluid Mech.

554

343

–

369

.

Gillet

N.

,

Brito

D.

,

Jault

D.

,

Nataf

H.-C.

,

2007

.

Experimental and numerical study of convection in a rapidly rotating spherical shell

,

J. Fluid Mech.

580

83

–

121

.

Gillet

N.

,

Schaeffer

N.

,

Jault

D.

,

2011

.

Rationale and geophysical evidence for quasi-geostrophic rapid dynamics within the Earth's outer core

,

Phys. Earth planet. Inter.

187

380

–

390

.

Gilman

P.A.

,

1977

.

Nonlinear dynamics of Boussinesq convection in a deep rotating spherical shell – I.

,

Geophys. Astrophys. Fluid Dyn.

8

93

–

135

.

Glatzmaier

G.

,

Evonuk

M.

,

Rogers

T.

,

2009

.

Differential rotation in giant planets maintained by density-stratified turbulent convection

,

Geophys. Astrophys. Fluid Dyn.

103

31

–

51

.

Goluskin

D.

,

Johnston

H.

,

Flierl

G.R.

,

Spiegel

E.A.

,

2014

.

Convectively driven shear and decreased heat flux

,

J. Fluid Mech.

759

360

–

385

.

Greenspan

H.P.

,

1968

.

The Theory of Rotating Fluids

University Press

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Guervilly

C.

,

2010

.

Dynamos numériques planétaires générées par cisaillement en surface ou chauffage interne

,

PhD thesis

Université Joseph Fourier Grenoble

,

https://tel.archives-ouvertes.fr/tel-00576177

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Guervilly

C.

,

Cardin

P.

,

2016

.

Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model

,

J. Fluid Mech.

808

61

–

89

.

Guervilly

C.

,

Cardin

P.

,

Schaeffer

N.

,

2012

.

A dynamo driven by zonal jets at the upper surface: applications to giant planets

,

Icarus

218

100

–

114

.

Heimpel

M.

,

Aurnou

J.

,

2007

.

Turbulent convection in rapidly rotating spherical shells: a model for equatorial and high latitude jets on Jupiter and Saturn

,

Icarus

187

(

2

),

540

–

557

.

Heimpel

M.

,

Aurnou

J.

,

Wicht

J.

,

2005

.

Simulation of equatorial and high-latitude jets on Jupiter in a deep convection model

,

Nature

438

193

–

196

.

Ingersoll

A.P.

,

Pollard

D.

,

1982

.

Motion in the interiors and atmospheres of Jupiter and Saturn - Scale analysis, anelastic equations, barotropic stability criterion

,

Icarus

52

62

–

80

.

Jones

C.A.

,

2015

.

Thermal and compositional convection in the outer core

, in

Treatise on Geophysics

2nd edn, pp.

115

–

159

, ed.

Schubert

G.

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Jones

C.A.

,

Soward

A.M.

,

Mussa

A.I.

,

2000

.

The onset of thermal convection in a rapidly rotating sphere

,

J. Fluid Mech.

405

157

–

179

.

Jones

C.A.

,

Rotvig

J.

,

Abdulrahman

A.

,

2003

.

Multiple jets and zonal flow on Jupiter

,

Geophys. Res. Lett.

30

(

14

),

1731, doi:10.1029/2003GL016980

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Kaplan

E.J.

Schaeffer

N.

Vidal

J.

Cardin

P.

,

2017

.

Subcritical thermal convection of liquid metals in a rotating sphere

,

Phys. Rev. Lett.

in press

, http://arxiv:1701.04623.

Google Scholar

OpenURL Placeholder Text

WorldCat

Kuo

H.-L.

,

1949

.

Dynamic instability of two-dimensional nondivergent flow in a barotropic atmosphere.

,

J. Atmos. Sci.

6

105

–

122

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Labbé

F.

,

Jault

D.

,

Gillet

N.

,

2015

.

On magnetostrophic inertia-less waves in quasi-geostrophic models of planetary cores

,

Geophys. Astrophys. Fluid Dyn.

109

(

6

),

587

–

610

.

Li

L.

et al. ,

2010

.

Saturn's emitted power

,

J. geophys. Res.

115

E11002, doi:10.1029/2010JE003631

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Livermore

P.W.

,

Hollerbach

R.

,

Finlay

C.C.

,

2017

.

An accelerating high-latitude jet in Earth's core

,

Nat. Geosci.

10

(

1

),

62

–

68

.

McIntyre

M.E.

,

2003

.

Potential vorticity

,

Encyclopedia Atmos. Sci.

2

685

–

694

.

McIntyre

M.E.

,

2008

.

Potential-vorticity inversion and the wave-turbulence jigsaw: some recent clarifications

,

Adv. Geosci.

15

47

–

56

.

Miyagoshi

T.

,

Kageyama

A.

,

Sato

T.

,

2010

.

Zonal flow formation in the Earth's core

,

Nature

463

793

–

796

.

Morin

V.

Dormy

E.

,

2004

.

Time dependent beta-convection in rapidly rotating spherical shells

,

Phys. Fluids

16

(

5

),

1603

–

1609

.

Or

A.C.

,

Busse

F.H.

,

1987

.

Convection in a rotating cylindrical annulus. II. Transitions to asymmetric and vacillating flow.

,

J. Fluid Mech.

174

313

–

326

.

Pedlosky

J.

,

1987

.

Geophysical Fluid Dynamics

Springer

.

Pirraglia

J.A.

,

1984

.

Meridional energy balance of Jupiter

,

Icarus

59

(

2

),

169

–

176

.

Plaut

E.

,

Lebranchu

Y.

,

Simitev

R.

,

Busse

F.H.

,

2008

.

On the Reynolds stresses and mean fields generated by pure waves: applications to shear flows and convection in a rotating shell

,

J. Fluid Mech.

602

303

–

326

.

Porco

C.C.

et al. ,

2003

.

Cassini imaging of Jupiter's atmosphere, satellites, and rings

,

Science

299

1541

–

1547

.

Read

P.L.

et al. ,

2015

.

An experimental study of multiple zonal jet formation in rotating, thermally driven convective flows on a topographic beta-plane

,

Phys. Fluids

27

(

8

),

085111, doi:10.1063/1.4928697

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Rhines

P.

,

1975

.

Waves and turbulence on a beta-plane

,

J. Fluid Mech.

69

417

–

443

.

Rotvig

J.

,

2007

.

Multiple zonal jets and drifting: thermal convection in a rapidly rotating spherical shell compared to a quasigeostrophic model

,

Phys. Rev. E

76

(

4

),

046306, doi:10.1103/PhysRevE.76.046306

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Schaeffer

N.

,

Cardin

P.

,

2005

.

Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers

,

Phys. Fluids

17

(

10

),

104111, doi:10.1063/1.2073547

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Schou

J.

et al. ,

1998

.

Helioseismic studies of differential rotation in the solar envelope by the solar oscillations investigation using the Michelson Doppler Imager

,

Astrophys. J.

505

(

1

),

390

–

417

.

Scott

R.K.

,

Dritschel

D.G.

,

2012

.

The structure of zonal jets in geostrophic turbulence

,

J. Fluid Mech.

711

576

–

598

.

Smith

C.A.

,

Speer

K.G.

,

Griffiths

R.W.

,

2014

.

Multiple zonal jets in a differentially heated rotating annulus

,

J. Phys. Oceanogr.

44

(

9

),

2273

–

2291

.

Soward

A.M.

,

1977

.

On the finite amplitude thermal instability of a rapidly rotating fluid sphere

,

Geophys. Astrophys. Fluid Dyn.

9

(

1

),

19

–

74

.

Sukoriansky

S.

,

Dikovskaya

N.

,

Galperin

B.

,

2007

.

On the arrest of inverse energy cascade and the Rhines scale

,

J. Atmos. Sci.

64

3312

–

3327

.

Sumita

I.

,

Olson

P.

,

2000

.

Laboratory experiments on high Rayleigh number thermal convection in a rapidly rotating hemispherical shell

,

Phys. Earth planet. Inter.

117

153

–

170

.

Teed

R.

,

Jones

C.

,

Hollerbach

R.

,

2012

.

On the necessary conditions for bursts of convection within the rapidly rotating cylindrical annulus

,

Phys. Fluids

24

(

6

),

066604, doi:10.1063/1.4711398

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Turner

J.S.

,

1985

.

Multicomponent convection

,

Annu. Rev. Fluid Mech.

17

(

1

),

11

–

44

.

Vallis

G.K.

,

2006

.

Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-scale Circulation

Cambridge Univ. Press

.

Verhoeven

J.

,

Stellmach

S.

,

2014

.

The compressional beta effect: a source of zonal winds in planets?

,

Icarus

237

143

–

158

.

von Hardenberg

J.

,

Goluskin

D.

,

Provenzale

A.

,

Spiegel

E.A.

,

2015

.

Generation of large-scale winds in horizontally anisotropic convection

,

Phys. Rev. Lett.

115

(

13

),

doi:10.1103/PhysRevLett.115.134501

.

Google Scholar

OpenURL Placeholder Text

WorldCat

Yadav

R.

,

Gastine

T.

,

Christensen

U.

,

Duarte

L.

,

Reiners

A.

,

2016

.

Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells

,

Geophys. J. Int.

204

(

2

),

1120

–

1133

.

Yoshida

S.

Sumita

I.

Kumazawa

M.

,

1996

.

Growth model of the inner core coupled with the outer core dynamics and the resulting elastic anisotropy

,

J. geophys. Res.

101

28 085

–

28 104

.

Zhang

K.

,

1991

.

Convection in a rapidly rotating spherical shell at infinite Prandtl number: steadily drifting rolls

,

Phys. Earth planet. Inter.

68

(

1–2

),

156

–

169

.

Zhang

K.

,

1992

.

Spiralling columnar convection in rapidly rotating spherical fluid shells

,

J. Fluid Mech.

236

535

–

556

.

Download all slides

Month:	Total Views:
July 2017	4
August 2017	33
September 2017	16
October 2017	13
November 2017	18
December 2017	16
January 2018	6
February 2018	2
March 2018	4
April 2018	13
May 2018	4
June 2018	3
July 2018	4
August 2018	9
September 2018	5
October 2018	9
November 2018	12
December 2018	2
January 2019	7
February 2019	3
March 2019	2
April 2019	6
May 2019	4
June 2019	7
July 2019	2
August 2019	2
September 2019	2
October 2019	4
November 2019	2
February 2020	2
March 2020	10
April 2020	3
May 2020	2
June 2020	7
July 2020	4
August 2020	2
October 2020	2
November 2020	15
December 2020	13
January 2021	14
February 2021	4
March 2021	19
April 2021	7
May 2021	7
June 2021	5
July 2021	15
August 2021	10
September 2021	14
October 2021	6
November 2021	8
December 2021	11
January 2022	4
February 2022	4
March 2022	7
April 2022	10
May 2022	2
June 2022	5
July 2022	9
August 2022	11
September 2022	34
October 2022	26
November 2022	6
December 2022	15
January 2023	3
February 2023	2
March 2023	7
April 2023	7
May 2023	2
June 2023	8
July 2023	5
August 2023	4
September 2023	4
October 2023	20
November 2023	8
December 2023	7
January 2024	10
February 2024	10
March 2024	10
April 2024	10
May 2024	9
June 2024	5
July 2024	10
August 2024	1
September 2024	4
October 2024	4
November 2024	7
December 2024	14
January 2025	15
February 2025	10
March 2025	19
April 2025	1
May 2025	10

Article Contents

Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores

Abstract

1 INTRODUCTION

2 MATHEMATICAL FORMULATION

2.1 Governing equation for the non-axisymmetric flow

2.2 Governing equation for the zonal flow

2.3 Governing equation for the temperature

2.4 Numerical method

3 STRUCTURE OF THE CONVECTIVE AND ZONAL FLOWS AND SCALING OF THE VELOCITY

3.1 Radial dependence of the convection

3.2 Zonal jet width and convective lengthscale

3.3 Scaling of the zonal flow velocity

4 DRIFT AND STABILITY OF THE ZONAL JETS

5 DYNAMICAL DIFFERENCE BETWEEN RETROGRADE AND PROGRADE JETS

6 EFFECT OF THE ZONAL FLOWS ON THE HEAT TRANSPORT

7 DISCUSSION

Acknowledgements

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores

Abstract

1 INTRODUCTION

2 MATHEMATICAL FORMULATION

2.1 Governing equation for the non-axisymmetric flow

2.2 Governing equation for the zonal flow

2.3 Governing equation for the temperature

2.4 Numerical method

3 STRUCTURE OF THE CONVECTIVE AND ZONAL FLOWS AND SCALING OF THE VELOCITY

3.1 Radial dependence of the convection

3.2 Zonal jet width and convective lengthscale

3.3 Scaling of the zonal flow velocity

4 DRIFT AND STABILITY OF THE ZONAL JETS

5 DYNAMICAL DIFFERENCE BETWEEN RETROGRADE AND PROGRADE JETS

6 EFFECT OF THE ZONAL FLOWS ON THE HEAT TRANSPORT

7 DISCUSSION

Acknowledgements

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only