# Large-scale vortices and zonal flows in spherical rotating convection

Yufeng Lin1<EMAIL_ADDRESS>Andrew Jackson2  1 Department of Earth and
Space Sciences, Southern University of Science and Technology, Shenzhen
518055, China 2 Institute of Geophysics, ETH Zurich, Zurich 8092, Switzerland

###### Abstract

Motivated by understanding the dynamics of stellar and planetary interiors, we
have performed a set of direct numerical simulations of Boussinesq convection
in a rotating full sphere. The domain is internally heated with fixed
temperature and stress-free boundary conditions, but fixed heat flux and no-
slip boundary conditions are also briefly considered. We particularly focus on
the large-scale coherent structures and the mean zonal flows that can develop
in the system. At Prandtl number of unity, as the thermal forcing (measured by
the Rayleigh number) is increased above the value for the onset of convection,
we find a relaxation oscillation regime, followed by a geostrophic turbulence
regime. Beyond this we see for the first time the existence of large-scale
coherent vortices that form on the rotation axis. All regime boundaries are
well described by critical values of the convective Rossby number $Ro_{c}$,
with transitions from oscillatory to geostrophic turbulence, and then to the
large-scale vortex regime at values $Ro_{c}\approx 0.2$ and $Ro_{c}\approx
1.5$, respectively. The zonal flow is controlled by the convective Rossby
number and changes its direction when the flow transitions from the
geostrophic turbulence regime to the large-scale vortex regime. While the non-
zonal flow speed and heat transfer can be described by the so-called inertial
scaling in the geostrophic turbulence regime, the formation of large-scale
vortices appears to reduce both the non-zonal flow speed and the efficiency of
convective heat transfer.

##  1 Introduction

Convection is commonplace in the interiors of stars and planets, in which
rotation plays a dominant role (e.g. Zhang & Liao, 2017). Despite a concerted
study of this system, both numerically and theoretically, over many decades,
our understanding of the process remains incomplete. It is of great interest
to properly characterise both the morphology of fluid flow and the mechanism
of heat transfer as a function of the dynamical regime as controlled by the
relevant non-dimensional numbers. The highly supercritical regime, which is
both the most challenging and the most relevant in planetary settings, remains
poorly understood.

Three classes of problem have been tackled in order to make progress. By far
the most studied is the problem of plane layer Rayleigh-Bénard convection, in
both non-rotating and rotating settings (Plumley & Julien, 2019). The second
class is rotating spherical shell convection, naturally motivated by the
geometry of planets that possess solid cores; such a geometry lends itself
naturally to efficient numerical methods (Glatzmaier, 2014; Marti et al.,
2016). This problem has been studied using laboratory experiments (Aubert et
al., 2001) and numerical simulations in both Boussinesq (Christensen, 2002;
Aurnou et al., 2007) and anelastic settings (Gastine et al., 2013), with a
variety of gravity profiles (Gastine et al., 2016) and boundary conditions
(Mound & Davies, 2017; Long et al., 2020). The results of the first two
problems have been recently reviewed by Jones (2015); Aurnou et al. (2015) and
Plumley & Julien (2019), and we refer the reader to these review articles.

The third class of study, which is the subject of the present research, is the
canonical problem of rotating convection in a whole sphere in which an
internal heat source is present; the strength of the heating is measured by
the Rayleigh number $Ra$. Although early theoretical work concentrated on the
whole-sphere geometry (Chandrasekhar, 1961; Roberts, 1968; Busse, 1970; Kida,
1994; Jones et al., 2000; Zhang & Liao, 2004) and the determination of the
critical Rayleigh number for convection ($Ra_{c}$), there has been much less
numerical work focusing on this area. Notwithstanding the studies that reduce
the size of the inner core to negligible size without actually removing it,
notable studies using a true full-sphere geometry are those of Jones et al.
(2000); Sánchez et al. (2016) and Kaplan et al. (2017), which mainly concerned
themselves with the onset of convection or subcritical convection. More
recently, Guervilly et al. (2019) examined the convective length scales in a
rotating sphere using both quasi-geostrophic and fully three-dimensional
methods in the regime of very low Ekman number and small Prandtl number,
albeit with a weak thermal forcing (around the onset of convection).
Convection was studied experimentally in the full sphere by Chamberlain &
Carrigan (1986) using centrifugal gravity, though practical considerations
precluded the exploration of very supercritical Rayleigh numbers: in that
study, $Ra\leq 2.5Ra_{c}$. The highly supercritical regime of convection in a
whole sphere is of great relevance to the dynamics of the Earth’s core prior
to inner core nucleation and to some fully convective stars, yet the problem
remains largely unexplored.

In the present work we concentrate on the problem in which the diffusivities
of heat and momentum are equal, leading to a so-called Prandtl number ($Pr$)
of unity. At fixed viscosity this leaves the rotation period (measured by the
Ekman number $E$) and the strength of internal heating (characterised by $Ra$)
as the only control parameters. We perform a set of direct numerical
simulations with moderate Ekman numbers ($E\sim 10^{-5}$), but over a wide
range of Rayleigh numbers (from near the onset up to $10^{4}Ra_{c}$).

A primary focus of our work has been to delineate the regime boundaries that
occur on increasing the Rayleigh number. We make use of the significant step
forward made by Gilman (1977) who identified the relevance of the convective
Rossby number $Ro_{c}$ in supercritical convection. In terms of the
aforementioned parameters, $Ro_{c}=\sqrt{RaE/Pr}$, and we find that $Ro_{c}$
is indeed the control parameter that defines the regime boundaries. A
relaxation oscillation regime occurs just beyond the onset of convection and
the boundary with the geostrophic turbulence (GT) regime that follows
subsequently is well described by a critical value of $Ro_{c}$. Likewise, we
find a second bifurcation to a regime with a large-scale vortex (LSV), which
is perfectly described by a constant value of $Ro_{c}$. Thus our results
support the significance of this parameter. Conversely, when looking
quantitatively within these regimes, we find that this parameter alone is
insufficient to fully describe the scaling of flow speed and heat transfer. In
particular, we find that the formation of LSVs appears to reduce the
efficiency of the heat transfer and the convective flow speed.

The identification of axisymmetric LSVs in spherical rotating convection is
entirely new. We argue that the omission of an inner core removes the obstacle
to such geostrophic structures that otherwise would be present. Large-scale
coherent vortices have been previously observed in planar rotating convection
in both the anelastic case (Chan, 2007; Chan & Mayr, 2013; Käpylä et al.,
2011) and Boussinesq case (Julien et al., 2012b; Favier et al., 2014;
Guervilly et al., 2014; Stellmach et al., 2014; Rubio et al., 2014). The
formation of LSV is attributed to an inverse cascade, i.e. upscale energy
transfer, in rotating turbulence. In spherical containers, large-scale
coherent structures may be associated with mean zonal flows (Christensen,
2002). We find a reversal of the mean zonal flow when the LSVs are formed.
Such a reversal has been observed in thin spherical shells for both Boussinesq
convection (Aurnou et al., 2007) and anelastic convection (Gastine et al.,
2013) when the Rayleigh number is sufficiently large, but never before in a
full sphere. Effects of boundary conditions on the formation of LSVs are also
briefly discussed.

##  2 Numerical Models

###  2.1 Governing equations

We consider Boussinesq convection in a full sphere with radius of $r_{o}$,
which rotates at $\bm{\Omega}=\Omega\bm{\hat{z}}$. The sphere is filled with
an incompressible fluid of density $\rho$, viscosity $\nu$ and thermal
diffusivity $\kappa$. Convective motions are driven by a homogeneous internal
heat source $S$, under a gravitational field $\bm{g}=-g_{o}\bm{r}/r_{o}$. In
the absence of convection, the basic state temperature is determined by the
conduction alone and is given as

${\widetilde{T_{b}}}=\frac{\beta}{2}(r_{o}^{2}-r^{2}),$ (1)

where $\beta=S/(3k)$ in which $k$ is thermal conductivity.

Using the radius $r_{o}$ as the length scale, the viscous diffusion time
$r_{o}^{2}/\nu$ as the time scale and $\beta r_{o}^{2}$ as the temperature
scale, the system is governed by the following dimensionless equations in the
rotating frame (Marti et al., 2014):

$E\left(\frac{\partial\bm{u}}{\partial t}+\bm{u\cdot\nabla
u}\right)+\bm{\hat{z}}\times\bm{u}=-\bm{\nabla}p+\frac{Ra}{Pr}T\bm{r}+E\nabla^{2}\bm{u},$
(2) $Pr\left(\frac{\partial T}{\partial
t}+\bm{u\cdot\nabla}T\right)=\nabla^{2}(T-T_{b}),$ (3) $\bm{\nabla\cdot u}=0,$
(4)

where $\bm{u}$ is the fluid velocity, $p$ is the reduced pressure, $T$ is the
total temperature and $T_{b}$ is the dimensionless basic temperature given as

$T_{b}=\frac{1}{2}(1-r^{2}).$ (5)

The dimensionless parameters are the Ekman number $E$, the rotationally
modified Rayleigh number $Ra$ and the Prandtl number $Pr$:

$E=\frac{\nu}{2\Omega r_{o}^{2}},\quad Ra=\frac{\alpha\beta
g_{o}r_{o}^{3}}{2\Omega\kappa},\quad\Pran=\frac{\nu}{\kappa},$ (6)

where $\alpha$ is the thermal expansion coefficient. The Prandtl number is
fixed to be unity ($\Pran=1$) in this study. Note that $Ra$ in a full sphere
is naturally flux-based as $\beta$ is related to the total heat flux
$\tilde{Q}$ through the boundary by

$\tilde{Q}=\frac{4\pi r_{0}^{3}}{3}S=4\pi r_{0}^{3}k\beta,$ (7)

regardless of thermal boundary conditions, and the non-dimensional heat flux
$q$ through the boundary is $q=-\nabla T$. Here $Ra$ is related to the
convectional Rayleigh number $Ra_{T}=\alpha g_{0}\Delta Tr_{o}^{3}/\nu\kappa$
based on the temperature difference by $Ra=2NuERa_{T}$, where $Nu$ is the
Nusselt number (see the definition given in equation (16) below).

As pointed out by Gilman (1977), the convective Rossby number, which is a
combination of $E$, $Ra$ and $Pr$, is a key control parameter in supercritical
rotating convection. Based on our definition, the convective Rossby number is
given as

$Ro_{c}=\sqrt{\frac{ERa}{Pr}}=\frac{\sqrt{\alpha\beta g_{0}r_{o}}}{2\Omega},$
(8)

which does not depend on the viscosity $\nu$ and the thermal diffusivity
$\kappa$. The convective Rossby number is the ratio between the rotation time
scale and the free-fall time scale, and characterises the global force balance
between the buoyancy force with respect to the Coriolis force (Gilman, 1977).
Christensen (2002) introduced a similar diffusivity-free parameter
$Ra^{*}=ERa/Pr$, which is the square of the convective Rossby number, i.e.
$Ra^{*}=Ro_{c}^{2}$.

We adopt the stress-free boundary condition for the velocity $\bm{u}$ and a
fixed temperature ($T=0$ for simplicity) at the boundary of the sphere in most
of the numerical simulations. The no-slip boundary condition and a fixed heat
flux, i.e. $\partial T/\partial r=-1.0$, boundary condition are also briefly
considered.

###  2.2 Numerical method

The governing equations (2-4) subject to the boundary conditions are
numerically solved using a spectral method in a whole sphere (Marti & Jackson,
2016). As we consider an incompressible fluid, the velocity $\bm{u}$ is
decomposed into toroidal and poloidal components:

$\bm{u}=\bm{\nabla\times}(\mathcal{T}\bm{r})+\bm{\nabla\times\nabla\times}(\mathcal{P}\bm{r}).$
(9)

The toroidal $\mathcal{T}$ and poloidal $\mathcal{P}$ scalar fields and the
temperature field $T$ are expanded in terms of spherical harmonic expansion on
spherical surfaces and the so-called Jones-Worland polynomials in the radial
direction. The spectral expansion is truncated up to spherical harmonics of
degree $L=255$ and order $M=255$, and up to $N=127$ for the Jones-Worland
polynomials. A critical element of the technique is its graceful handling of
the origin at the centre of the sphere, such that the Courant–Friedrichs–Lewy
condition is not unduly stringent in this region. Time-stepping has been
performed using a second-order Runge-Kutta scheme, with adaptive time steps.
The typical time step is about $10^{-7}$ (viscous time scale) for the most
demanding calculations. We refer the reader to Marti & Jackson (2016) for more
details regarding the numerical scheme. The numerical code has been
benchmarked for several hydrodynamic and magnetohydrodynamic problems in a
full sphere, including the rotating thermal convection problem (Marti et al.,
2014).

###  2.3 Diagnostics

In order to qualitatively describe numerical results, the following quantities
derived from the velocity $\mathbf{u}$ and the temperature $T$ will be used.

The total kinetic energy is given by

$\mathcal{E}=\frac{1}{2}\int_{v}\bm{u}^{2}\mathrm{d}V.$ (10)

The development of strong zonal flow is a distinct feature in the nonlinear
regime, hence we can decompose the total energy into zonal and non-zonal
parts:

$\mathcal{E}_{\mathrm{zonal}}=\frac{1}{2}\int_{v}\left(u_{\phi}^{0}\right)^{2}\mathrm{d}V,\quad\mathcal{E}_{\mathrm{non}}=\mathcal{E}-\mathcal{E}_{\mathrm{zonal}},$
(11)

where $u_{\phi}^{0}$ is the axisymmetric ($m=0$) component of the azimuthal
velocity $u_{\phi}$, i.e. zonal flow. We define the Rossby number as

$Ro=\frac{U_{rms}}{2\Omega r_{0}},$ (12)

where $U_{rms}$ is the dimensional root mean square velocity. Base on the non-
dimensionalisation we used, the Rossby number is related to the non-
dimensional energy as

$Ro=E\sqrt{\frac{2\mathcal{E}}{V}},$ (13)

where $V=4\pi/3$ is the non-dimensional volume of the sphere. Accordingly, we
have the zonal Rossby number

$Ro_{\mathrm{zon}}=E\sqrt{\frac{2\mathcal{E_{\mathrm{zonal}}}}{V}},$ (14)

and the non-zonal Rossby number

$Ro_{\mathrm{non}}=E\sqrt{\frac{2\mathcal{E_{\mathrm{non}}}}{V}}.$ (15)

Heat transport owing to convection can be measured by the Nusselt number. For
the convection of internal heating in a full sphere, the Nusselt number can be
estimated by the drop of the temperature with respect to the basic state at
the centre (Guervilly & Cardin, 2016; Kaplan et al., 2017):

$Nu=\frac{T_{b}(r=0)}{T(r=0)}.$ (16)

Apart from these global quantities, we often use the axial vorticity and the
zonal velocity to visualise flows. The axial vorticity is calculated by

$\omega_{z}=2E\mathbf{\hat{z}}\mathbf{\cdot}(\mathbf{\nabla\times u}).$ (17)

Note that the prefactor of $2E$ leads to the vorticity in the units of the
rotation frequency $\Omega$ given the non-dimensionalisation we have used. The
zonal velocity is also rescaled, i.e. $U_{\phi}^{0}=2Eu_{\phi}^{0}$, such that
$U_{\phi}^{0}$ is in the units of $\Omega r_{0}$. The mean zonal flow
$\overline{U_{\phi}^{0}}$ is time-averaged $U_{\phi}^{0}$.

##  3 Results

Figure 1: Regime diagram in the plane of $(E,Ra)$. Different symbols represent
different flow regimes. Open circles correspond to the steadily drifting
regime near the onset; filled circles correspond to the relaxation oscillation
regime; blue diamonds correspond to the GT regime; filled triangles correspond
to the LSV regime; the open triangle in the top right corresponds to the non-
rotating regime.

Numerical simulations in this study are listed in detail in table 1 in
Appendix A. We fix $\Pran=1$ in this study and vary the control parameters $E$
and $Ra$, leading to a variety of dynamical regimes. Stress-free and fixed
temperature boundary conditions are assumed unless stated otherwise. Figure 1
summarises simulations as a function of $(E,Ra)$ with symbols indicating
different flow regimes, which are discussed in section 3.1.

###  3.1 Flow regimes

(a) (b) (c)

Figure 2: Flows near the onset at $E=10^{-5}$, $Pr=1$, $Ra=220\approx
1.05Ra_{c}$. (a) Axial vorticity in the equatorial plane and (b) in the
meridional plane; (c) zonal flow in the meridional plane.

Convection takes place only when $Ra$ is larger than a critical value
$Ra_{c}$. At $Pr=1$ with the stress-free and the fixed temperature boundary
conditions, the critical value is given as $Ra_{c}\approx 4.1E^{-1/3}+17.8$
(black solid line in figure 1) according to the asymptotic theory of the
linear onset Jones et al. (2000). Our fully nonlinear numerical simulations
are consistent with the linear theory regarding the onset of convection.
Figure 2 shows a snapshot of the axial vorticity and the zonal flow when $Ra$
is marginally above the critical value ($Ra=220\approx 1.05Ra_{c}$) at
$E=10^{-5}$. The flow near the onset is in the form of the so-called thermal
Rossby waves, which are quasi-geostrophic, i.e. nearly invariant along the
rotation axis (figure 2 (b)). The kinetic energy becomes quasi-steady after
the initial growing phase. The azimuthal wavenumber $m=14$ is also in
agreement with the linear theory. We also see that a weak zonal flow develops
(figure 2 (c)) with differential rotation (prograde near the equator and
retrograde in the inner region) due to the nonlinear interactions of the
linear onset mode (Zhang, 1992).

(a)

(b) (c)

Figure 3: An example of the oscillatory flow regime at $E=10^{-5}$,
$Ra=1000\approx 4.8Ra_{c}$. (a) Time evolution of the zonal and non-zonal
kinetic energy. Time is in units of the viscous time scale. (b) and (c) Axial
vorticity in the equatorial plane and in the meridional plane at the instants
indicated by the blue dot and red dot in (a) respectively. (A movie of this
time series is provided; see Movie 1).

When $Ra$ is a few times larger than the critical value, convective motions
become quasi-periodic. Figure 3 shows an example of the relaxation oscillation
regime at $E=10^{-5}$ and $Ra=1000\approx 4.8Ra_{c}$. We can see that the time
evolution of the energy exhibits a quasi-periodic manner (3a). The growth of
the onset mode (3b) leads to subsequent instabilities, resulting in more
vigorous and complex convection (3c). A movie of the axial vorticity for this
case is provided (see Movie 1). As the non-zonal energy approaches the maxima,
there is a burst of the zonal energy, followed by a relaminarisation process
(decreasing of both the zonal and non-zonal energy). During these cycles, the
non-zonal energy can vary over several orders of magnitude between peaks and
troughs. Similar relaxation oscillations have been observed in numerical
simulations of convection in spherical shells, and are also referred to as
convective bursts (Grote & Busse, 2001; Busse, 2002; Christensen, 2002;
Heimpel & Aurnou, 2012). The relaxation oscillations can be attributed to the
competition between the nonlinear effect and the viscous effect (Christensen,
2002). The oscillatory regime exists only at intermediate Rayleigh numbers
(filled circles in figure 1).

(a) (b) (c)

(d) (e)

Figure 4: An example of the GT regime at $E=10^{-5}$ and $Ra=10^{5}$. (a)
Axial vorticity in the equatorial plane; (b) Axial vorticity in the meridional
plane. (c) Zonal flow in the meridional plane. (d-e) Energy spectra as a
function the azimuthal wavenumber $m$ and spherical harmonic degree $l$. Black
dashed lines represent the $-5/3$ power law purely for reference.

As we increase $Ra$, time evolution of the kinetic energy becomes less
periodic and eventually becomes irregular fluctuations, but the flow
structures are still elongated along the rotation axis due to the constraint
of rotation. We refer to this regime as the GT regime. The transition from the
relaxation oscillation to the GT regime seems to be controlled by the
convective Rossby number, i.e. $Ro_{c}\approx 0.2$ (blue dashed line in figure
1). Figure 4 shows a typical case of the GT regime at $E=10^{-5}$ and
$Ra=10^{5}$. In this regime, convective motions are charaterised by space-
filling and sustained turbulence (4a-b). The axial vorticity in the meridional
plane (4b) shows elongated structures along the rotation axis due to the
constraint of rotation in this regime. This is a typical flow regime in
rotating convection when $Ra$ is supercritical, which has been observed in a
wide range of numerical simulations and laboratory experiments (Aurnou et al.,
2015).

Apart from small-scale turbulence, geostrophic zonal flows (figure 4(c)) are
developed in this case, with a strong prograde jet near the equator and mainly
retrograde flows in the bulk. The significance of the zonal flow is clearly
evident in the energy spectra (figure 4 (d-e)) as a function of the azimuthal
wavenumber $m$ and spherical harmonic degree $l$. The kinetic energy is
dominated by the $m=0$ toroidal component, which is basically the zonal flow.
In the $l$ spectra, $l=0$ and $l=2$ toroidal components contribute the most
kinetic energy in the flow.

(a) (b) (c)

(d) (e)

Figure 5: An example of the LSV regime at $E=10^{-5}$ and $Ra=3\times 10^{5}$.
(a) Axial vorticity in the equatorial plane; (b) Axial vorticity in the
meridional plane. (c) Zonal flow in the meridional plane. (d-e) Energy spectra
as a function the azimuthal wavenumber $m$ and spherical harmonic degree $l$.
Black dashed lines represent the $-5/3$ power law purely for reference.

Further increasing $Ra$ leads to a completely different flow regime (red
triangles in figure 1), which is characterised by the formation of LSV. The
transition from the GT regime to the LSV regime is also determined by the
convective Rossby number, i.e. $Ro_{c}\approx 1.5$ (the red dashed line in
figure 1). Figure 5 shows an example of the LSV at $E=10^{-5}$ and $Ra=3\times
10^{5}$. We can see a notable cyclonic vortex located at the centre in the
equatorial plane (figure 5(a)) and extended along the rotation axis (figure
5(b)). The formation process of the LSV is further discussed in section 3.2.
The vortex is accompanied by a strong geostrophic zonal flow (figure 5(c)).
However, the direction of the zonal flow is reversed compared to the previous
regime. The zonal flow is retrograde near the equator and prograde in the
inner region. The kinetic energy is dominated by the $m=0$ toroidal component,
which corresponds to the columnar vortex and the zonal flow. The sawtooth-
shaped curves in the $l$spectra at large $l$ reflect the equatorially
symmetric nature of the LSV and the zonal flow.

(a) (b) (c)

(d) (e)

Figure 6:  An example of LSV breakdown at $E=6.0\times 10^{-5}$ and
$Ra=2.0\times 10^{6}$. (a) Axial vorticity in the equatorial plane; (b) Axial
vorticity in the meridional plane. (c) Zonal flow in the meridional plane.
(d-e) Energy spectra as a function the azimuthal wavenumber $m$ and spherical
harmonic degree $l$. Black dashed lines represent the $-5/3$ power law.

At sufficiently large $Ra$, large-scale coherent structures break down and
convection tends to behaviour similar to that in a non-rotating regime. Figure
6 shows an example in this regime at $E=6.0\times 10^{-5}$ and $Ra=2.0\times
10^{6}\approx 1.6\times 10^{4}Ra_{c}$ (corresponding to the open triangle in
figure 1). Unlike previous cases, axial vorticity in the equatorial and
meridional planes does not exhibit obvious anisotropy figure 6(a-b). The zonal
flow becomes much less geostrophic compared with previous cases (figure 6(c)).
Energy spectra in both $m$ and $l$ are in agreement with the Kolmogorov
scaling of $k^{-5/3}$ over a wide range of wavenumbers (figure 6(d-e)). This
non-rotating regime exists at very large $Ra$ and it is computationally
demanding to simulate this regime at lower $E$. As our focus in this study is
the LSV, we do not explore this regime further.

###  3.2 Formation of the LSV

Figure 7: Time series of energy at $E=10^{-5}$ and $Ra=3.0\times 10^{5}$. (a)
Zonal and non-zonal energy. (b) Energy contained in several different
azimuthal wavenumbers $m$ during the formation stage of the LSV. Time is in
units of the viscous time scale.

We now turn to analyse the formation process of the LSV. Figure 7 shows the
time evolution of the zonal energy and the non-zonal energy (a), as well as
the kinetic energy carried by several different $m$ during the growth phase of
the LSV (b) at $E=10^{-5}$ and $Ra=3.0\times 10^{5}$ (corresponding to the
case in figure 5). The simulation started from a saturated state at lower $Ra$
in the GT regime. We can see that the zonal energy overtakes the non-zonal
energy after the transient stage and then the non-zonal energy drops. The
simulation saturates to a state of which the zonal energy is dominant over the
non-zonal energy.

The energy transfer from a small-scale flows to a large-scale flow is probably
best illustrated in figure 7 (b). While the energy in the $m=0$ component
keeps growing and becomes dominant, the energy in other components initially
grows but subsequently drops. We note that the larger the azimuthal wavenumber
$m$, the earlier the energy drop, suggesting a successive energy transfer from
larger $m$ to smaller $m$. Similar upscale energy transfer has been reported
in rotating Rayleigh-Bénard convection, leading to box-sized vortices
(Guervilly et al., 2014; Rubio et al., 2014). In the full sphere we consider
here, the upscale energy transfer leads to a large cyclonic vortex located at
the centre and strong zonal flow (figure 5).

(a) (b)

(c) (d)

(e) (f)

Figure 8: Snapshots (corresponding to the black dots in figure 7(b)) of the
axial vorticity in the equatorial plane and zonal flow in the meridional plane
at $E=10^{-5}$ and $Ra=3.0\times 10^{5}$. (A movie of this time series is
provided; see Movie 2).

Figure 8 shows a sequence of snapshots (corresponding to the black dots in
figure 7(b)) of the axial vorticity and the zonal flow during the formation
stage of the LSV at $E=10^{-5}$ and $Ra=3.0\times 10^{5}$. We clearly see the
process of vortex merging from small scales (figure 8(a)) to a few prominent
vortices (figure 8 (b-c)), which further merge into one dominant cyclonic
vortex (figure 8 (d)). The eye of the cyclone spirals inwards and eventually
settles at the centre (figure 8 (e-f)). The vortex is elongated along the
rotation axis as we have shown in figure 5 (b). As the vortices merge and
migrate inwards, the prograde zonal flow also shifts inwards, suggesting an
inwards transport of angular momentum. We further discuss the mean zonal flow
in section 3.3.

Figure 9 shows the same snapshots as in figure 8 but for the temperature
field. We also see the merging process from small-scale to larger-scale
structures of the temperature field, culminating with a relatively hot
columnar aggregation around the rotation axis as the LSV is formed. The
instantaneous Nusselt number $Nu$ is indicated at the bottom of each snapshot.
We note that the formation of the LSV results in a significant drop of $Nu$,
e.g. from $Nu=56.9$ (in figure 9e) to $Nu=25.11$ (in figure 9f). This is
perhaps due to the fact that a hot column of fluid is trapped around the
rotation axis by the LSV, which reduces the amount of heat transported by the
convection. We further discuss effects of the LSV on the heat transport in
section 3.4.

(a) (b)

(c) (d)

(e) (f)

Figure 9: As for figure 8 but for the temperature in the equatorial plane and
in the meridional plane. Instantaneous Nusselt number $Nu$ is also indicated.

###  3.3 Mean zonal flows

As we have shown, zonal flows are always developed once the convective motions
set in. In this section, we compare the zonal flow in different regimes and
investigate how the zonal flow varies depending on the control parameters.

(a) (b) (c)

Figure 10: Mean zonal flows in the meridional plane at different $Ra$ for
fixed $E=10^{-5}$. (a) $Ra=1.0\times 10^{3},~{}Ro_{c}=0.1$, in the relaxation
oscillation regime; (b) $Ra=1.0\times 10^{5},~{}Ro_{c}=1.0$, in the GT regime;
(c) $Ra=3.0\times 10^{5},~{}Ro_{c}=\sqrt{3}$, in the LSV regime.

Figure 10 shows the mean zonal flow, i.e. time-averaged zonal flow
$\overline{U_{\phi}^{0}}$, in the meridional plane at various $Ra$ but fixed
$E=10^{-5}$. In all cases, the mean zonal flows are almost invariant along the
rotation axis. However, the direction of the mean zonal flow reverses when the
convective Rossby number is sufficiently large such that the LSV forms, i.e.
$Ro_{c}\gtrsim 1.5$. Similar reversal of the zonal flow has been observed in
thin shells for both Boussinesq convection (Aurnou et al., 2007) and anelastic
convection (Gastine et al., 2013). They also found that the reversal is
controlled by the convective Rossby number $Ro_{c}$, or equivalently by
$Ra^{*}$.

(a) (b)

Figure 11: (a) Mean zonal flow on the surface as a function of the latitude at
various $E$ and $Ra$, but all cases in the GT regime. (b) As for (a) but the
mean zonal flow is rescaled by $Ro_{c}^{3/5}$. Solid (dashed) lines correspond
to cases at $E=10^{-5}$ ($E=3\times 10^{-5}$).

As the mean zonal flow is essentially geostrophic, we use the zonal velocity
profile in the equatorial plane to represent the zonal flow. Figure 11 (a)
shows the zonal profiles at various $E$ and $Ra$ in the GT regime. These
profiles collapse into an invariant profile when the zonal velocity is divided
by $Ro_{c}^{3/5}$. Note that the exponent $3/5$ is merely an empirical
estimate and may be subject to small revisions. Nevertheless, it is clear that
the zonal flow is controlled by the the convective Rossby number in the GT
regime.

(a) (b)

Figure 12: (a) Mean zonal flow on the surface as a function of the latitude at
various $E$ and $Ra$, but all cases in the LSV regime. As for (a) but showing
the specific angular momentum in the inertial frame.

Figure 12(a) shows the zonal profiles in the LSV regime for various $E$ and
$Ra$, but the same $Ro_{c}=\sqrt{3}$. The zonal flow profile in the LSV regime
(figure 12 (a)) is almost, but not exactly, the reverse of the profile in the
GT regime (figure 11 (b)). Aurnou et al. (2007) proposed that the vigorous
convection at large $Ra$ homogenises the angular momentum in the inertial
frame, leading to a retrograde equatorial jet and a prograde zonal flow inside
in the rotating frame. Figure 12 (b) shows the specific angular momentum in
the inertial frame, i.e. $M_{I}=s^{2}+s\overline{U_{\phi}^{0}}$ (non-
dimensional). We see an inward transfer of the angular momentum compared to
solid body rotation, but the angular momentum profiles are far from a
homogeneous distribution, even in the outer region. One may expect a more
homogeneous distribution of the angular momentum on further increasing $Ra$,
but the zonal flow appears to saturate in the LSV regime as we can see from
figure 13 (a).

We have been unable to discover a principle that predicts this profile. Of
interest is the fact that the angular momentum in the LSV regime is clearly
prograde for $s\lessapprox 0.72$ and retrograde for $s\gtrapprox 0.72$, with
the sum being zero, of course, in the rotating frame. But what determines this
changeover point $s^{\ast}$? We merely remark that a fluid in solid body
rotation contains half its angular momentum in $s\leq 0.77$ and half exterior
to this radius. Any connection between $s^{\ast}$ and this value may be
fortuitous.

(a) (b)

Figure 13: Zonal Rossby number $Ro_{zon}$ (a) and rescaled non-zonal Rossby
number $Ro_{non}/E^{2/5}$ (b) as a function of the convective Rossby number
$Ro_{c}$. Different symbols represent different flow regimes as in figure 1.
$Ro_{zon}$ and $Ro_{non}$ are time-averaged over several convective turnover
time scales.

The zonal Rossby number $Ro_{zon}$ increases with increasing $Ro_{c}$ with
varying slopes depending on the flow regimes, and becomes saturated
($Ro_{zon}\sim 0.15$) in the LSV regime (figure 13(a)). The plot is truncated
to $Ro_{zon}\geq 10^{-3}$ to highlight the supercritical cases as $Ro_{zon}$
is very small near the onset. In the oscillatory regime $Ro_{c}<0.2$,
$Ro_{zon}$ also depends on the Ekman number, suggesting that the viscosity
plays a part in the zonal flow. Based on a balance between the Reynolds stress
and the internal viscous stress, Christensen (2002) proposed
$Ro_{zon}\propto\left(Ro_{non}\right)^{2}/E$. The non-zonal Rossby number as a
function of $Ro_{c}$ is shown in figure 13(b), but the above relation does not
hold for our numerical results as both $R_{zon}$ and $Ro_{non}$ follow the
solid black line and roughly scale as $Ro_{c}^{2}$ in the oscillatory regime
$Ro_{c}<0.2$.

In the GT regime $0.2\leq Ro_{c}\leq 1.5$, $Ro_{zon}$ becomes independent of
the viscosity and the empirical scaling $Ro_{zon}\sim Ro_{c}^{3/5}$ (black
dashed line in figure 13(a))is evident, in agreement with figure 11(b). In
this regime, the non-zonal Rossby number is expected to follow the so-called
inertial scaling, which is based on a force balance between Coriolis, inertia
and Archimedean force, if the viscous effect is negligible (Aubert et al.,
2001). The inertial scaling can be expressed as the following based on our
definition of the control parameters (see Appendix B):

$Ro_{non}\sim\left(\frac{E}{Pr}\right)^{2/5}\left(Ro_{c}\right)^{4/5}.$ (18)

In figure 13(b), $Ro_{non}$ is divided by $E^{2/5}$ ($Pr=1$) to compare with
the inertial scaling (black dashed line). While the data points at different
$E$ are well collapsed, we note that numerical results exhibit a steeper trend
(larger exponent) than the inertial scaling. Previous numerical simulations in
spherical shells also found steeper exponents (e.g. Gastine et al., 2016; Long
et al., 2020) and the deviation is usually attributed to the non-negligible
effects of the viscosity in numerical simulations.

###  3.4 Heat transport

(a) (b)

Figure 14: Nusselt number $Nu$ (a) and $NuE^{2/5}$ (b) as a function of the
convective Rossby number $Ro_{c}$. Different symbols represent different flow
regimes as in figure 1. The blue line depicting supercriticality is plotted in
the lower left-hand corner of (a), for $Ro_{c}\leq 0.15$ and $E=10^{-5}$. Here
$Nu$ is time-averaged over several convective turnover time scales.

Figure 14 shows the Nusselt number as a function of the convective Rossby
number $Ro_{c}$. Near the onset of the convection, the convective heat
transfer is proportional to the supercriticality according to a weakly
nonlinear theory (Gillet & Jones, 2006), i.e.

$Nu-1\propto\frac{Ra}{Ra_{c}}-1.$ (19)

We can see from the lower left-hand corner of figure 14(a) that our numerical
results can be well fitted by the aforementioned scaling when $Ra<6Ra_{c}$ for
$E=10^{-5}$ (the blue solid line).

In rotating turbulent convection, the inertial scaling has been promoted to
describe heat transport (Stevenson, 1979; Aubert et al., 2001; Julien et al.,
2012a; Barker et al., 2014; Jones, 2015). Based on our definition of the non-
dimensional parameters, this inertial scaling for the Nusselt number can be
written as (Appendix B):

$Nu\sim\left(\frac{E}{Pr}\right)^{-2/5}\left(Ro_{c}\right)^{6/5}.$ (20)

We can see that our numerical results in the GT regime (blue diamonds in
figure 14) are in very good agreement with the inertial scaling. Note that we
use the stress-free boundary condition in this study. Previous studies using
the no-slip boundary found that the scaling exponent in $Ra$ depends on the
Ekman number (Cheng et al., 2015; Gastine et al., 2016; Long et al., 2020),
i.e. $Nu\sim Ra^{\lambda(E)}$, though it may asymptotically approach the
inertial scaling at low Ekman numbers (Gastine et al., 2016).

There is a sudden drop of $Nu$ on increasing $Ro_{c}$ when the flow transits
from the GT regime to the LSV regime $Ro_{c}>1.5$. As we have mentioned, the
reduction of the heat transport by the LSV is mainly due to the fact that a
relatively hot column is formed around the rotation axis. Nusselt number $Nu$
continues to increase on further increase of $Ra$ in the LSV regime, but the
slope seems to be less steep compared to the GT regime. It is difficult to
deduce any definitive scaling tendency, as higher forcings would be required,
and this is computationally prohibitive.

(a) (b)

Figure 15: Heat flux $q$ on the surface in the LSV regime at $E=10^{-5}$ and
$Ra=3.0\times 10^{5}$. (a) Instantaneous heat flux on the surface (Mollweide
projection). (b) Time and azimuthal average of the heat flux $q$ on the
surface as a function of latitude.

Figure 15 (a) shows the heat flux $q$ on the surface in the LSV regime at
$E=10^{-5}$ and $Ra=3.0\times 10^{5}$. The heat flux exhibits both small-scale
perturbations and large-scale coherent variation in latitude, in line with the
formation of LSV. The time and azimuthal average in figure 15 (b) shows that
the heat flux at high latitudes is about three times larger than that around
the equator. The latitudal variations of the heat flux again suggest that the
formation of LSV provides a barrier for the horizontal heat transport in the
system.

###  3.5 Effects of boundary conditions

Thus far, we have made use of stress-free and fixed temperature boundary
conditions. In this section, we briefly discuss effects of boundary conditions
on the formation of LSV.

(a) (b) (c)

Figure 16: Axial vorticity and zonal flow with the no-slip boundary condition
at $E=1.0\times 10^{-5}$ and $Ra=5.0\times 10^{5}$. (a-b) Axial vorticity in
the equatorial plane and in the meridional plane; (c) zonal flow in the
meridional plane.

(a) (b)

(c) (d)

Figure 17: Comparison of the fixed temperature boundary condition (a-b) and
the fixed flux boundary condition (c-d) for $E=3.0\times 10^{-5}$ and
$Ra=10^{5}$. (a) and (c) Axial vorticity in the equatorial plane and in the
meridional plane. (b) and (d) Temperature in the equatorial plane and in the
meridional plane.

Previous studies (Guervilly et al., 2014; Stellmach et al., 2014), studying
planar rotating convection, found that the LSV exists only when a stress-free
boundary condition is imposed. In the present study, we briefly explore the
possibility of formation of LSV with a no-slip boundary, but do not observe
any LSV being formed, at least in the parameter range we could access. Figure
16 shows an example with a no-slip boundary condition at $E=10^{-5}$ and
$Ra=5\times 10^{5}$, at which the LSV would form with the stress-free boundary
condition. With the no-slip boundary condition, we see that the flow is
characterised by small-scale structures elongated along the rotation axis
(figure 16(a-b)). We note that a zonal flow is also developed with a prograde
jet near the equator and a retrograde flow in the inner region (figure 16(c)),
but the amplitude is much smaller compared to the zonal flow with a stress-
free boundary condition. In fact, the zonal Rossby number is smaller than the
non-zonal Rossby number when no-slip boundary conditions are employed (table
1), in contrast to the case with the stress-free boundary condition.

Figure 17 compares different thermal boundary conditions, i.e. fixed
temperature (17a-b) and fixed flux (17c-d) conditions at $E=1.0\times 10^{-5}$
and $Ra=5.0\times 10^{5}$. Both cases use the stress-free boundary condition.
It can be seen that in this regime the thermal boundary conditions play a
minor role. A major difference is the temperature field in the polar regions,
where the hot columnar structure associated with the LSV needs to adapt
different thermal boundary conditions. The hot columnar region can extend to
the polar surface when the fixed-flux condition is imposed.

##  4 Conclusions

We have observed the previously unseen presence of LSVs in full-sphere
rotating convection, thus complementing the previous discoveries in Boussinesq
rotating boxes (Favier et al., 2014; Guervilly et al., 2014; Stellmach et al.,
2014; Rubio et al., 2014). Relatively simple large-scale zonal flow profiles
are developed in this regime, with a prograde axial vortex and a retrograde
equatorial zonal flow. While the geostrophic character of these flows is
pervasive, the ubiquitous radial structure of the zonal flow in the LSV regime
remains a theoretical challenge to comprehend. We have been able to fully
document the transition between regimes as the Rayleigh number is increased
from that required to first initiate convection, finding that the convective
Rossby number is a completely prognostic control parameter: $Ro_{c}\approx
0.2$ and $Ro_{c}\approx 1.5$ delineate the boundaries between the oscillatory
regime, the geostrophic turbulence regime and the LSV regime.

In the geostrophic turbulence regime, the zonal Rossby number $Ro_{zon}$
solely depends on the convective Rossby number and follows an empirical
scaling law $Ro_{zon}\sim Ro_{c}^{3/5}$, while the non-zonal Rossby number and
the Nusselt number can be described by the well-known inertial scaling. The
direction of the zonal flow is reversed when the flow transitions from the
geostrophic turbulence regime to the LSV regime, while the non-zonal Rossby
number and the Nusselt number are reduced by the formation of LSV. In the LSV
regime we witness a saturation of the zonal flow speed, while the non-zonal
Rossby number and the Nusselt number show an abrupt drop followed by a gradual
increase, whose slopes are reduced compared to their systematics in the
geostrophic turbulence regime.

The effects of both the mechanical and thermal boundary conditions on the
formation of LSVs are briefly considered. While the thermal boundary condition
plays a minor role, the no-slip boundary condition prohibits the formation of
LSVs, at least for the parameter regimes we could access, in line with
previous studies in planar rotating convection (Guervilly et al., 2014;
Stellmach et al., 2014). It is possible that the dynamics with the stress-free
and the no-slip conditions would converge at very low Ekman numbers, but it is
computationally demanding to study the highly supercrtical convection at very
low Ekman numbers via direct numerical simulations. We note the prospect of
better experimental delineation of zonal flows in the future using the ZoRo
experiment (Su et al., 2020). This rapidly rotating spheroid has the
capability of zonal flow determination using acoustic normal mode velocimetry,
and will be subject to differential heating in the future.

## Declaration of Interests

The authors report no conflict of interest.

## Acknowledgements

We thank Philippe Marti for the provision of the computer code EPMDynamo, and
much implementation advice. Simulations were performed at the Swiss National
Supercomputer Centre (CSCS) under account s872, on the Euler cluster of ETH
Zurich and on the Taiyi cluster supported by the Center for Computational
Science and Engineering of Southern University of Science and Technology. The
research was partially supported by a European Research Council Advanced grant
833848 (UEMHP) to AJ under the European Union’s Horizon 2020 research and
innovation programme. YL was supported by the B-type Strategic Priority
Program of the Chinese Academy of Sciences (XDB41000000), the National Natural
Science Foundation of China (grant No.41904066), the pre-research project on
Civil Aerospace Technologies of China National Space Administration (D020308)
and the Macau Foundation. YL was partly supported by a Swiss NSF Advanced
PostDoc Mobility Fellowship when this study was initiated.

##  Appendix A List of numerical simulations

All of the numerical simulations are listed in table 1 including the control
parameters, time-averaged diagnostic parameters, numerical resolutions
($N,L,M$) and flow regimes. Here SD represents the steadily drifting regime
near the onset; RO represents the relaxation oscillation regime; GT represents
the geostrophic turbulence regime; LSV represents the large scale vortices; NR
represents the non-rotating regime. The case marked with # used the stress-
free and fixed flux boundary condition, while the one marked with * used the
no-slip and fixed temperature boundary conditions. All other cases used the
stress-free and fixed temperature boundary conditions.

$E$ | $Ra$ | $Ro_{zon}$ | $Ro_{non}$ | $Nu$ | $(N,L,M)$ | Regime
---|---|---|---|---|---|---
$6.0\times 10^{-5}$ | $Ra_{c}=126.04$ |  |  |  |  |
$6.0\times 10^{-5}$ | $1.3\times 10^{2}$ | $2.85\times 10^{-5}$ | $1.55\times 10^{-4}$ | 1.002 | (63,63,63) | SD
$6.0\times 10^{-5}$ | $2.0\times 10^{2}$ | $5.81\times 10^{-4}$ | $6.94\times 10^{-4}$ | 1.041 | (63,63,63) | SD
$6.0\times 10^{-5}$ | $3.0\times 10^{2}$ | $1.56\times 10^{-3}$ | $1.20\times 10^{-3}$ | 1.085 | (63,63,63) | RO
$6.0\times 10^{-5}$ | $5.0\times 10^{2}$ | $3.77\times 10^{-3}$ | $2.02\times 10^{-3}$ | 1.774 | (63,63,63) | RO
$6.0\times 10^{-5}$ | $1.0\times 10^{3}$ | $9.08\times 10^{-3}$ | $3.93\times 10^{-3}$ | 1.774 | (63,63,63) | GT
$6.0\times 10^{-5}$ | $1.5\times 10^{3}$ | $1.26\times 10^{-2}$ | $4.90\times 10^{-3}$ | 2.158 | (63,63,63) | GT
$6.0\times 10^{-5}$ | $2.0\times 10^{3}$ | $1.49\times 10^{-2}$ | $5.80\times 10^{-3}$ | 2.683 | (63,63,63) | GT
$6.0\times 10^{-5}$ | $5.0\times 10^{3}$ | $2.23\times 10^{-2}$ | $1.07\times 10^{-2}$ | 4.731 | (63,127,127) | GT
$6.0\times 10^{-5}$ | $1.0\times 10^{4}$ | $2.87\times 10^{-2}$ | $2.27\times 10^{-2}$ | 7.162 | (63,127,127) | GT
$6.0\times 10^{-5}$ | $2.0\times 10^{4}$ | $3.65\times 10^{-2}$ | $3.15\times 10^{-2}$ | 10.57 | (63,127,127) | GT
$6.0\times 10^{-5}$ | $3.3\times 10^{4}$ | $4.11\times 10^{-2}$ | $3.99\times 10^{-2}$ | 13.55 | (63,127,127) | GT
$6.0\times 10^{-5}$ | $5.0\times 10^{4}$ | $1.05\times 10^{-1}$ | $3.80\times 10^{-2}$ | 11.88 | (63,127,127) | LSV
$6.0\times 10^{-5}$ | $7.0\times 10^{4}$ | $1.20\times 10^{-1}$ | $4.37\times 10^{-2}$ | 13.01 | (63,127,127) | LSV
$6.0\times 10^{-5}$ | $1.0\times 10^{5}$ | $1.33\times 10^{-1}$ | $6.06\times 10^{-2}$ | 14.37 | (63,127,127) | LSV
$6.0\times 10^{-5}$ | $5.0\times 10^{5}$ | $1.54\times 10^{-1}$ | $1.26\times 10^{-1}$ | 29.43 | (127,255,255) | LSV
$6.0\times 10^{-5}$ | $1.0\times 10^{6}$ | $1.50\times 10^{-1}$ | $1.32\times 10^{-1}$ | 38.49 | (127,255,255) | LSV
$6.0\times 10^{-5}$ | $2.0\times 10^{6}$ | $1.46\times 10^{-1}$ | $1.81\times 10^{-1}$ | 59.41 | (127,255,255) | NR
$3.0\times 10^{-5}$ | $Ra_{c}=152.95$ |  |  |  |  |
$3.0\times 10^{-5}$ | $1.6\times 10^{2}$ | $2.43\times 10^{-5}$ | $9.62\times 10^{-5}$ | 1.002 | (63,63,63) | SD
$3.0\times 10^{-5}$ | $2.0\times 10^{2}$ | $2.06\times 10^{-4}$ | $2.93\times 10^{-4}$ | 1.019 | (63,63,63) | SD
$3.0\times 10^{-5}$ | $3.0\times 10^{2}$ | $7.54\times 10^{-4}$ | $5.36\times 10^{-4}$ | 1.054 | (63,63,63) | RO
$3.0\times 10^{-5}$ | $5.0\times 10^{2}$ | $2.09\times 10^{-3}$ | $9.18\times 10^{-4}$ | 1.159 | (63,127,127) | RO
$3.0\times 10^{-5}$ | $1.0\times 10^{3}$ | $5.83\times 10^{-3}$ | $2.22\times 10^{-3}$ | 1.694 | (63,127,127) | RO
$3.0\times 10^{-5}$ | $2.0\times 10^{3}$ | $1.10\times 10^{-2}$ | $3.38\times 10^{-3}$ | 2.361 | (63,127,127) | GT
$3.0\times 10^{-5}$ | $3.0\times 10^{3}$ | $1.38\times 10^{-2}$ | $4.23\times 10^{-3}$ | 2.982 | (63,127,127) | GT
$3.0\times 10^{-5}$ | $5.0\times 10^{3}$ | $1.71\times 10^{-2}$ | $5.84\times 10^{-3}$ | 4.152 | (63,127,127) | GT
$3.0\times 10^{-5}$ | $8.0\times 10^{3}$ | $2.04\times 10^{-2}$ | $7.75\times 10^{-3}$ | 5.577 | (63,127,127) | GT
$3.0\times 10^{-5}$ | $1.0\times 10^{4}$ | $2.21\times 10^{-2}$ | $8.79\times 10^{-3}$ | 6.615 | (63,127,127) | GT
$3.0\times 10^{-5}$ | $2.0\times 10^{4}$ | $2.86\times 10^{-2}$ | $1.33\times 10^{-2}$ | 9.684 | (63,127,127) | GT
$3.0\times 10^{-5}$ | $5.0\times 10^{4}$ | $3.78\times 10^{-2}$ | $2.93\times 10^{-2}$ | 15.25 | (63,127,127) | GT
$3.0\times 10^{-5}$ | $7.0\times 10^{4}$ | $3.92\times 10^{-2}$ | $3.56\times 10^{-2}$ | 17.89 | (63,127,127) | GT
$3.0\times 10^{-5}$ | $1.0\times 10^{5}$ | $1.04\times 10^{-1}$ | $2.46\times 10^{-2}$ | 15.08 | (127,255,255) | LSV
$3.0\times 10^{-5}$ | $1.0\times 10^{5}$ | $1.05\times 10^{-1}$ | $2.58\times 10^{-2}$ | – | (127,255,255) | LSV #
$3.0\times 10^{-5}$ | $2.0\times 10^{5}$ | $1.23\times 10^{-1}$ | $3.76\times 10^{-2}$ | 18.25 | (127,255,255) | LSV
$3.0\times 10^{-5}$ | $5.0\times 10^{5}$ | $1.54\times 10^{-1}$ | $8.13\times 10^{-2}$ | 24.27 | (127,255,255) | LSV
$3.0\times 10^{-5}$ | $1.0\times 10^{6}$ | $1.67\times 10^{-1}$ | $9.65\times 10^{-2}$ | 31.92 | (127,255,255) | LSV
$1.0\times 10^{-5}$ | $Ra_{c}=210.52$ |  |  |  |  |
$1.0\times 10^{-5}$ | $2.2\times 10^{2}$ | $1.01\times 10^{-5}$ | $3.42\times 10^{-5}$ | 1.001 | (63,63,63) | SD
$1.0\times 10^{-5}$ | $3.0\times 10^{2}$ | $1.75\times 10^{-4}$ | $1.52\times 10^{-4}$ | 1.020 | (63,127,127) | SD
$1.0\times 10^{-5}$ | $5.0\times 10^{2}$ | $7.25\times 10^{-4}$ | $2.17\times 10^{-4}$ | 1.116 | (63,127,127) | RO
$1.0\times 10^{-5}$ | $1.0\times 10^{3}$ | $2.56\times 10^{-3}$ | $4.52\times 10^{-4}$ | 1.366 | (63,127,127) | RO
$1.0\times 10^{-5}$ | $2.0\times 10^{3}$ | $6.19\times 10^{-3}$ | $1.17\times 10^{-3}$ | 2.046 | (63,127,127) | RO
$1.0\times 10^{-5}$ | $5.0\times 10^{3}$ | $1.11\times 10^{-2}$ | $2.25\times 10^{-3}$ | 3.326 | (63,127,127) | RO
$1.0\times 10^{-5}$ | $1.0\times 10^{4}$ | $1.44\times 10^{-2}$ | $3.46\times 10^{-3}$ | 5.120 | (127,255,255) | GT
$1.0\times 10^{-5}$ | $3.0\times 10^{4}$ | $2.12\times 10^{-2}$ | $6.60\times 10^{-3}$ | 10.37 | (127,255,255) | GT
$1.0\times 10^{-5}$ | $3.0\times 10^{4}$ | $2.12\times 10^{-2}$ | $6.60\times 10^{-3}$ | 10.37 | (127,255,255) | GT
$1.0\times 10^{-5}$ | $1.0\times 10^{5}$ | $3.23\times 10^{-2}$ | $1.18\times 10^{-2}$ | 19.76 | (127,255,255) | GT
$1.0\times 10^{-5}$ | $2.0\times 10^{5}$ | $3.84\times 10^{-2}$ | $1.74\times 10^{-2}$ | 27.03 | (127,255,255) | GT
$1.0\times 10^{-5}$ | $3.0\times 10^{5}$ | $1.08\times 10^{-1}$ | $1.56\times 10^{-2}$ | 20.03 | (127,255,255) | LSV
$1.0\times 10^{-5}$ | $5.0\times 10^{5}$ | $1.19\times 10^{-1}$ | $1.92\times 10^{-2}$ | 29.57 | (127,255,255) | LSV
$1.0\times 10^{-5}$ | $5.0\times 10^{5}$ | $2.14\times 10^{-2}$ | $2.60\times 10^{-2}$ | 43.17 | (127,255,255) | GT *
Table 1: Summary of numerical simulations. Symbols # and * denote the cases
using different boundary boundary conditions (see text for details).

##  Appendix B Inertial scaling expressed in our parameters

In rapidly rotating convection, the so-called inertial scaling is proposed to
describe the convective flow speed, typical length scale and heat transport
(Stevenson, 1979; Aubert et al., 2001). The inertial scaling is based on a
balance of Coriolis-Inertial-Archimedean (CIA) terms, incorporated with the
mixing length theory, leading to a scaling independent of the viscosity. Here
we follow the review article of Jones (2015) and express the scaling in terms
of dimensionless parameters defined in this study. Note that the Rayleigh
number $Ra$ we defined is rotationally modified and heat-flux based, which is
loosely related to the $Ra_{Q}$ in Jones (2015) by $Ra\approx ERa_{Q}$.
According to equation (119) in Jones (2015), the inertial scaling predicts the
non-zonal Rossby number (representing the convective flow velocity) as

$Ro_{non}\sim\left(\frac{E}{Pr}\right)(EPr)^{1/5}Ra_{Q}^{2/5}\sim\left(\frac{E}{Pr}\right)^{4/5}Ra^{2/5}.$
(21)

We can write the above scaling in the convective Rossby number
$Ro_{c}=\sqrt{RaE/Pr}$ as

$Ro_{non}\sim\left(\frac{E}{Pr}\right)^{2/5}Ro_{c}^{4/5}.$ (22)

Note that the inertial scaling does not provide direct prediction of the zonal
Rossby number.

The inertial scaling predicts the Nusselt number as (equation (119) in Jones,
2015)

$Nu\sim
E^{4/5}Pr^{-1/5}Ra_{Q}^{3/5}\sim\left(\frac{E}{Pr}\right)^{1/5}Ra^{3/5}.$ (23)

Note that we have assumed $Nu-1\sim Nu$ in the strongly supercritical regime.
Using $Ra^{*}=RaE/Pr$, we obtain

$Nu\sim\left(\frac{E}{Pr}\right)^{-2/5}\left(Ra^{*}\right)^{3/5}.$ (24)

The above scaling is equivalent to

$Nu\sim E^{2}Pr^{-1/2}Ra_{T}^{3/2},$ (25)

where $Ra_{T}$ is the convectional Rayleigh number based on the temperature
difference, i.e. $Ra_{T}=\alpha g_{0}\Delta Tr_{o}^{3}/{\nu\kappa}$, and is
related to $Ra$ by

$Ra=2NuERa_{T}.$ (26)

Equation (24) can be also represented in terms of the convective Rossby number
$Ro_{c}$ using $Ra^{*}=Ro_{c}^{2}$:

$Nu\sim\left(\frac{E}{Pr}\right)^{-2/5}\left(Ro_{c}\right)^{6/5}.$ (27)

## References

  * Aubert et al. (2001) Aubert, J, Brito, D, Nataf, H.-C., Cardin, P & Masson, J.-P. 2001 Asystematic experimental study of rapidly rotating spherical convection in water and liquid gallium. Phys. Earth Planet. Inter. 128, 51–74.
  * Aurnou et al. (2015) Aurnou, J.M., Calkins, M.A., Cheng, J.S., Julien, K., King, E.M., Nieves, D., Soderlund, K.M. & Stellmach, S. 2015 Rotating convective turbulence in Earth and planetary cores. Phys. Earth Planet. Int. 246, 52–71.
  * Aurnou et al. (2007) Aurnou, Jonathan, Heimpel, Moritz & Wicht, Johannes 2007 The effects of vigorous mixing in a convective model of zonal flow on the ice giants. Icarus 190 (1), 110–126.
  * Barker et al. (2014) Barker, Adrian J., Dempsey, Adam M. & Lithwick, Yoram 2014 Theory and simulations of rotating convection. Astrophys. J. 791 (1).
  * Busse (1970) Busse, F. H. 1970 Thermal instabilities in rapidly rotating systems. J. Fluid Mech. 44, 441–460.
  * Busse (2002) Busse, F. H. 2002 Convective flows in rapidly rotating spheres and their dynamo action. Phys. Fluids 14 (4), 1301–1314.
  * Chamberlain & Carrigan (1986) Chamberlain, J. A. & Carrigan, C. R. 1986 An experimental investigation of convection in a rotating sphere subject to time varying thermal boundary conditions. Geophys. Astrophys. Fluid Dyn. 35 (1-4), 303–327.
  * Chan (2007) Chan, K. L. 2007 Rotating convection in f-boxes: Faster rotation. Astron. Nachr. 328 (10), 1059–1061.
  * Chan & Mayr (2013) Chan, Kwing L. & Mayr, Hans G. 2013 Numerical simulation of convectively generated vortices: Application to the Jovian planets. Earth Planet. Sci. Lett. 371-372, 212–219.
  * Chandrasekhar (1961) Chandrasekhar, S. 1961 Hydrodynamic and hydromagnetic stability. Clarendon Press.
  * Cheng et al. (2015) Cheng, J. S., Stellmach, S., Ribeiro, A., Grannan, A., King, E. M. & Aurnou, J. M. 2015 Laboratory-numerical models of rapidly rotating convection in planetary cores. Geophys. J. Int. 201 (1), 1–17.
  * Christensen (2002) Christensen, U.~R. 2002 Zonal flow driven by strongly supercritical convection in rotating spherical shells. J. Fluid Mech. 470, 115–133.
  * Favier et al. (2014) Favier, B., Silvers, L. J. & Proctor, M. R. E. 2014 Inverse cascade and symmetry breaking in rapidly rotating Boussinesq convection. Phys. Fluids 26 (9), 096605.
  * Gastine et al. (2016) Gastine, Thomas, Wicht, Johannes & Aubert, Julien 2016 Scaling regimes in spherical shell rotating convection. J. Fluid Mech. 808, 690–732.
  * Gastine et al. (2013) Gastine, T., Wicht, J. & Aurnou, J. M. 2013 Zonal flow regimes in rotating anelastic spherical shells: An application to giant planets. Icarus 225 (1), 156–172.
  * Gillet & Jones (2006) Gillet, N. & Jones, C. A. 2006 The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder. J. Fluid Mech. 554, 343.
  * Gilman (1977) Gilman, Peter A. 1977 Nonlinear dynamics of Boussinesq convection in a deep rotating spherical shell-i. Geophys. Astrophys. Fluid Dyn. 8 (1), 93–135.
  * Glatzmaier (2014) Glatzmaier, Gary A. 2014 Introduction to modeling convection in planets and stars: magnetic field, density stratification, rotation. Princeton series in astrophysics . Princeton University Press.
  * Grote & Busse (2001) Grote, E & Busse, F H 2001 Dynamics of convection and dynamos in rotating spherical fluid shells. Fluid Dyn. Res. 28 (5), 349–368.
  * Guervilly & Cardin (2016) Guervilly, Céline & Cardin, Philippe 2016 Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model. J. Fluid Mech. 808, 61–89.
  * Guervilly et al. (2019) Guervilly, Céline, Cardin, Philippe & Schaeffer, Nathanaël 2019 Turbulent convective length scale in planetary cores. Nature 570 (7761), 368–371.
  * Guervilly et al. (2014) Guervilly, Céline, Hughes, David W. & Jones, Chris A. 2014 Large-scale vortices in rapidly rotating Rayleigh–Bénard convection. J. Fluid Mech. 758, 407–435.
  * Heimpel & Aurnou (2012) Heimpel, Moritz & Aurnou, Jonathan M. 2012 Convective bursts and the coupling of Saturn’S equatorial storms and interior rotation. Astrophys. J. 746 (1).
  * Jones (2015) Jones, C.A. 2015 Thermal and Compositional Convection in the Outer Core. In Treatise on Geophysics 2nd Edition (ed. G. Schubert), , vol. 8, pp. 115–159. Elsevier.
  * Jones et al. (2000) Jones, Chris A., Soward, Andrew M. & Mussa, Ali I. 2000 The onset of thermal convection in a rapidly rotating sphere. J. Fluid Mech. 405 (2000), 157–179.
  * Julien et al. (2012a) Julien, Keith, Knobloch, Edgar, Rubio, Antonio M. & Vasil, Geoffrey M. 2012a Heat transport in low-Rossby-number Rayleigh-Bénard convection. Phys. Rev. Lett. 109 (25), 1–5.
  * Julien et al. (2012b) Julien, K., Rubio, A. M., Grooms, I. & Knobloch, E. 2012b Statistical and physical balances in low Rossby number Rayleigh-Bénard convection. Geophys. Astrophys. Fluid Dyn. 106 (4-5), 392–428.
  * Kaplan et al. (2017) Kaplan, E. J., Schaeffer, N., Vidal, J. & Cardin, P. 2017 Subcritical Thermal Convection of Liquid Metals in a Rapidly Rotating Sphere. Phys. Rev. Lett. 119 (9), 094501.
  * Käpylä et al. (2011) Käpylä, Petri J., Mantere, Maarit J. & Hackman, Thomas 2011 Starspots due to large-scale vortices in rotating turbulent convection. Astrophys. J. 742 (1), 34.
  * Kida (1994) Kida, Shigeo 1994 Stability of thermal convection in a rapidly rotating sphere. J. Phys. Soc. Japan 63 (8), 2964–2973.
  * Long et al. (2020) Long, R. S., Mound, J. E., Davies, C. J. & Tobias, S. M. 2020 Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions. J. Fluid Mech. 889, A7.
  * Marti et al. (2016) Marti, P., Calkins, M. A. & Julien, K. 2016 A computationally efficient spectral method for modeling core dynamics. Geochem. Geophys. Geosyst. 17 (8), 3031–3053.
  * Marti & Jackson (2016) Marti, P. & Jackson, A. 2016 A fully spectral methodology for magnetohydrodynamic calculations in a whole sphere. J. Comp. Phys. 305, 403–422.
  * Marti et al. (2014) Marti, P., Schaeffer, N., Hollerbach, R., Cebron, D., Nore, C., Luddens, F., Guermond, J.-L., Aubert, J., Takehiro, S., Sasaki, Y., Hayashi, Y.-Y., Simitev, R., Busse, F., Vantieghem, S. & Jackson, A. 2014 Full sphere hydrodynamic and dynamo benchmarks. Geophys. J. Int. 197 (1), 119–134.
  * Mound & Davies (2017) Mound, JE & Davies, CJ 2017 Heat transfer in rapidly rotating convection with heterogeneous thermal boundary conditions. J. Fluid Mech. .
  * Plumley & Julien (2019) Plumley, Meredith & Julien, Keith 2019 Scaling laws in rayleigh-bénard convection. Earth Space Sci. 6 (9), 1580–1592.
  * Roberts (1968) Roberts, P. H. 1968 On the thermal instability of a rapidly rotating fluid sphere containing heat sources. Phil. Trans. R. Soc. Lond. A 263, 93–117.
  * Rubio et al. (2014) Rubio, Antonio M., Julien, Keith, Knobloch, Edgar & Weiss, Jeffrey B. 2014 Upscale Energy Transfer in Three-Dimensional Rapidly Rotating Turbulent Convection. Phys. Rev. Lett. 112 (14), 144501.
  * Stellmach et al. (2014) Stellmach, S., Lischper, M., Julien, K., Vasil, G., Cheng, J., Ribeiro, A., King, E. & Aurnou, J. 2014 Approaching the Asymptotic Regime of Rapidly Rotating Convection: Boundary Layers versus Interior Dynamics. Phys. Rev. Lett. 113 (25), 254501.
  * Stevenson (1979) Stevenson, David J. 1979 Turbulent Thermal Convection in the Presence of Rotation and a Magnetic Field: A Heuristic Theory. Geophys. Astrophys. Fluid Dyn. 12 (1), 139–169.
  * Su et al. (2020) Su, Sylvie, Cébron, David, Nataf, Henri-Claude, Cardin, Philippe, Vidal, Jérémie, Solazzo, Max & Do, Yann 2020 Acoustic spectra of a gas-filled rotating spheroid. Eur. J. Mech. B Fluids 84, 302 – 310.
  * Sánchez et al. (2016) Sánchez, Juan, Garcia, Ferran & Net, Marta 2016 Critical torsional modes of convection in rotating fluid spheres at high Taylor numbers. J. Fluid Mech. 791, R1.
  * Zhang (1992) Zhang, K. 1992 Spiralling columnar convection in rapidly rotating spherical fluid shells. J. Fluid Mech. 236, 535–556.
  * Zhang & Liao (2004) Zhang, Keke & Liao, Xinhao 2004 A new asymptotic method for the analysis of convection in a rapidly rotating sphere. J. Fluid Mech. 518, 319–346.
  * Zhang & Liao (2017) Zhang, Keke & Liao, Xinhao 2017 Theory and Modeling of Rotating Fluids: Convection, Inertial Waves and Precession, 1st edn. Cambridge Monographs on Mechanics . Cambridge University Press.