# Comparison of quasi-geostrophic, hybrid and 3D models of planetary core convection

O. Barrois, T. Gastine, C. C. Finlay

► To cite this version:

O. Barrois, T. Gastine, C. C. Finlay. Comparison of quasi-geostrophic, hybrid and 3D models of planetary core convection. Geophysical Journal International, 2022, 10.1093/gji/ggac141. insu-03748541

HAL Id: insu-03748541

https://insu.hal.science/insu-03748541

Submitted on 10 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Comparison of quasi-geostrophic, hybrid and 3-D models of planetary core convection

O. Barrois,<fnref n="1" /> T. Gastine Ⓔ<fnref n="2" /> and C.C. Finlay Ⓒ<fnref n="1" />

Accepted 2022 April 20. Received 2022 March 29; in original form 2021 June 23

## SUMMARY

We present investigations of rapidly rotating convection in a thick spherical shell geometry relevant to planetary cores, comparing results from quasi-geostrophic (QG), 3-D and hybrid QG-3D models. The 170 reported calculations span Ekman numbers, *Ek*, between $10^{-4}$ and $10^{-10}$, Rayleigh numbers, *Ra*, between 2 and 150 times supercritical and Prandtl numbers, *Pr*, between 10 and $10^{-2}$. The default boundary conditions are no-slip at both the ICB and the CMB for the velocity field, with fixed temperatures at the ICB and the CMB. Cases driven by both homogeneous and inhomogeneous CMB heat flux patterns are also explored, the latter including lateral variations, as measured by *Q*\*, the peak-to-peak amplitude of the pattern divided by its mean, taking values up to 5. The QG model is based on the open-source *pizza* code. We extend this in a hybrid approach to include the temperature field on a 3-D grid. In general, we find convection is dominated by zonal jets at mid-depths in the shell, with thermal Rossby waves prominent close to the outer boundary when the driving is weaker. For the thick spherical shell geometry studied here the hybrid method is best suited for studying convection at modest forcing, $Ra \le 10 Ra_c$ when $Pr = 1$, and departs from the 3-D model results at higher *Ra*, displaying systematically lower heat transport characterized by lower Nusselt and Reynolds numbers. We find that the lack of equatorially-antisymmetric motions and z-correlations between temperature and velocity in the buoyancy force contributes to the weaker flows in the hybrid formulation. On the other hand, the QG models yield broadly similar results to the 3-D models, for the specific aspect ratio and range of Rayleigh numbers explored here. We cannot point to major disagreements between these two data sets at $Pr \ge 0.1$, with the QG model effectively more strongly driven than the hybrid case due to its cylindrically averaged thermal boundary conditions. When *Pr* is decreased, the range of agreement between the hybrid and 3-D models expands, for example up to $Ra \le 15 Ra_c$ at $Pr = 0.1$, indicating the hybrid method may be better suited to study convection in the low *Pr* regime. We thus observe a transition between two regimes: (i) at $Pr \ge 0.1$ the QG and 3-D models agree in the studied range of *Ra/Ra*<sub>c</sub> while the hybrid model fails when $Ra > 15 Ra_c$ and (ii) at $Pr = 0.01$ the QG and 3-D models disagree for $Ra > 10 Ra_c$ while the hybrid and 3-D models agree fairly well up to $Ra \sim 20 Ra_c$. Models that include laterally varying heat flux at the outer boundary reproduce regional convection patterns that compare well with those found in similarly forced 3-D models. Previously proposed scaling laws for rapidly rotating convection are tested; our simulations are overall well described by a triple balance between Coriolis, inertia and Archimedean forces with the length-scale of the convection following the diffusion-free Rhines-scaling. The magnitude of *Pr* affects the number and the size of the jets with larger structures obtained at lower *Pr*. Higher velocities and lower heat transport are seen on decreasing *Pr* with the scaling behaviour of the convective velocity displaying a strong dependence on *Pr*. This study is an intermediate step towards a hybrid model of core convection also including 3-D magnetic effects.

**Key words:** Earth Core; thermal convection; theories and simulations; Numerical modelling.

<footnote n="1">¹Division of Geomagnetism, DTU Space, Technical University of Denmark, Lyngby DK-2800, Denmark. E-mail: ovrbr@space.dtu.dk</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, Université de Paris, CNRS, F-75005 Paris, France</footnote>
# 1 INTRODUCTION

Many celestial bodies such as rocky and gas planets of the Solar system are rapidly rotating. The effects of rotation on fluid systems have been widely studied—it impedes the onset of the convection (Chandrasekhar 1961), constrains heat transport (Rossby 1969) and shapes the convection into the form of thin columns nearly invariant along the rotation axis (Busse 1970). Such convective flows are subject to a zeroth order Geostrophic balance between the Coriolis force and the pressure gradient that arises when $Ek \ll 1$ and $Ro \ll 1$ (Julien et al. 2012), and when the typical timescale of the convection is much longer than the rotation period, where $Ek = v/\Omega d^2$ measures the viscous effects compared to the Coriolis force and $Ro = ReEk$ measures the non-linear inertial effects compared to the Coriolis force. Here $Re = Ud/v$ is the Reynolds number with $v$ the kinematic viscosity of the fluid, $\Omega$ the angular velocity of the planet, $d$ the typical size of the fluid container and $U$ the typical fluid velocity.

In this study we focus on the convective dynamics relevant for the Earth's outer core. It is expected to be in a strongly driven state, with $Re \gg 1$ and turbulent convection (see Roberts & King 2013, for a review). Such a regime is extremely challenging to explore both experimentally and numerically, because of the important non-linearities and the necessity of resolving fast rotational dynamics while wishing to track the evolution of long-lived jets and vortices (e.g., Stellmach et al. 2014; Aurnou et al. 2015; Gastine et al. 2016). The relevant dimensionless parameters for the Earth's core are thought to be $Ek \sim 10^{-15}$, $Re \sim 10^9$ and $Ra/Ra_c \gg 10^3$ where $Ra$ is the Rayleigh number and $Ra_c$ is the critical Rayleigh number for the onset of convection. The most ambitious 3-D simulations and experimental studies are only able to reach $Ek \sim 10^{-7}$; $Ra/Ra_c \sim 10^3$ and $Re \sim 5 \times 10^3$ (Aubert 2015; Schaeffer et al. 2017; Sheyko et al. 2018).

One alternative avenue for studying this challenging regime is to use reduced quasi-geostrophic (QG) convection models. In their classical form QG models consider perturbations about a leading order balance between Coriolis and Pressure gradient forces, whose axial vorticity is invariant along the rotation axis. The dynamics is then essentially confined to the equatorial plane. Busse (1970) initially developed QG models in an annulus geometry assuming a small boundary slope. The QG framework was later modified and extended to better account for spherical geometry and to include phenomenon such as Ekman pumping by Cardin & Olson (1994); Aubert et al. (2003); Schaeffer & Cardin (2005); Gillet & Jones (2006); Calkins et al. (2012). With such QG models it has been possible to investigate rotating convection for Ekman numbers as low as $Ek = 10^{-11}$ close to the onset of convection (Guervilly et al. 2019).

Such QG models are essentially a 2-D approximation of the real 3-D situation. The 2-D treatment of temperature often used in QG models is not rigorously justified (see, e.g. Gillet & Jones 2006) and fails to capture thermal wind contributions. Furthermore such classical QG models focus on an axially invariant approximation of the axial vorticity and on related horizontal flows in the equatorial plane. In spherical geometry they perform worst close to the outer boundary where the boundary slope becomes large, or when the forcing becomes large enough that the vertical velocity becomes significant, and 3-D motions set in (Calkins et al. 2013). More advanced extensions of the QG framework have recently been proposed where the full velocity field is better accounted for by projecting onto a QG basis (Labbé et al. 2015; Maffei et al. 2017; Gerick et al. 2020) or by z-averaging before taking the curl (Jackson & Maffei 2020). In this work, we follow Gastine (2019) and use the QG formulation proposed by Schaeffer & Cardin (2006) that was expanded in a hybrid approach by Guervilly & Cardin (2016, 2017) to also include a 3-D temperature field. Our numerical implementation of this hybrid QG-3D method (or simply hybrid) is an extension of the pizza code by Gastine (2019) to include a 3-D temperature field in a spherical shell geometry. Here we explore advantages and limitations of QG and hybrid QG-3D models compared with full 3-D core convection models.

Thermal boundary conditions may play an important role in convection in planetary cores. Strictly speaking these boundary conditions are not fixed but time-dependent, and coupling to compositional effects should be considered (Glatzmaier & Roberts 1996). In practice however, when considering Earth's core, fixed heat-flux conditions at the core-mantle boundary and fixed temperature conditions at the inner core boundary are often argued to be relevant (e.g., Gubbins et al. 2003). Early studies with heat-flux boundary conditions suggested these might promote slightly longer wavelengths and larger convective flows (Gibbons et al. 2007; Sakuraba & Roberts 2009) although such discrepancies have more recently been attributed to different levels of forcing (Yadav et al. 2016; Schwaiger et al. 2021). This is consistent with an asymptotic equivalence between heat flux and temperature boundary conditions (Calkins et al. 2015) with standard universal scaling laws retrieved far from the onset in both cases (Clarté et al. 2021). More dramatic effects are possible when the heat flux boundary conditions vary laterally, which in some locations will enhance heat transport and can result in a preservation of large-scale downwelling systems (Mound & Davies 2017; Long et al. 2020; Sahoo & Sreenivasan 2020). The above statements are primarily based on studies carried out at moderate Ekman numbers ($Ek \ge 2 \times 10^{-6}$) that were often weakly or moderately driven. Here, although the majority of our simulations use fixed temperature boundary conditions, we report results of a number of calculations with imposed heat flux outer boundary condition, including inhomogeneous cases where this varies laterally. We examine whether inhomogeneous boundary conditions continue to impact the convective pattern in more strongly driven cases at Ekman numbers slightly smaller than those considered in previous studies ($Ek \ge 10^{-6}$ with our hybrid approach and $Ek \ge 10^{-7}$ with our QG approach) and investigate whether QG and hybrid QG-3D models can capture relevant aspects of convection in such cases.

Scaling laws describe how global quantities characterizing the convection, such as the convective length scale, flow speed and heat transport, vary with the control parameters based on the underlying dynamics (e.g., Gillet & Jones 2006; King & Buffett 2013; Gastine et al. 2016). Two theoretical scaling laws have attracted much attention for describing the properties of rapidly rotating convective flows: one based on a triple balance between the Coriolis, Inertia and Archimedean forces—called the CIA-scaling (Ingersoll & Pollard 1982; Cardin & Olson 1994)—and the other based on a triple balance between Viscous, Archimedean and Coriolis forces—called the VAC-scaling (King & Buffett 2013). Early studies at modest rotation rates had difficulty in distinguishing between the two scalings (Aubert et al. 2001; Gillet & Jones 2006;
King & Buffett 2013), but more recent investigations have shown a preference for the CIA balance in the bulk of the fluid, away from viscous boundary layers (Gastine et al. 2016; Long et al. 2020; Schwaiger et al. 2021). An impressive convergence towards the viscous-free scaling in the limit of low viscosity and close to the onset of convection has also recently been described by Guervilly et al. (2019) in the context of fluids with Prandtl numbers $Pr < 1$. Here our main goal is to complement these studies, using QG, 3-D and hybrid QG-3D simulations in a thick spherical shell geometry, focusing on relatively strongly driven cases (high Ra/Ra$_c$) and exploring the role of the Prandtl number, which may influence the typical size of the convective pattern (Calkins et al. 2012; King & Aurnou 2013; Guervilly & Cardin 2016).

The paper is organized as follows: Section 2 presents the equations and the methodology of our QG, hybrid QG-3D and 3-D models. Section 3 presents results obtained using our models focusing on comparisons between QG, hybrid and 3-D calculations, and including cases with inhomogeneous heat flux boundary conditions. We also describe the impact of Prandtl number on the form of convection at low Rossby number and examine how well our results satisfy convective scaling laws. We conclude with a discussion and a summary of our findings in Section 4.

# 2 METHODOLOGY

## 2.1 QG model formulation

In this section we first describe the basic QG model used before moving on to the new 3-D modifications we have implemented. We use the same QG model formulation and notation as Gastine (2019), who followed closely the approach set out by Schaeffer & Cardin (2005) and Gillet & Jones (2006). We work in cylindrical coordinates $(s, \phi, z)$ in a spherical shell between the inner radius $s_i$ and the outer radius $s_o$, rotating about the z-axis with a constant angular velocity $\Omega$. The horizontal components of the velocity field $\mathbf{u}_{\perp}$, perpendicular to the rotation axis, are assumed to be invariant along the rotation axis, that is $\mathbf{u}_{\perp} = (u_s, u_{\phi}, 0)$, where $u_s$ and $u_{\phi}$ are, respectively, the radial and azimuthal velocities. It is further assumed that the dynamics is encapsulated by the evolution of the axial vorticity averaged in the z direction, such that the dynamics is restricted to that in the equatorial plane of the spherical shell (Maffei et al. 2017). Below we refer to this as the classical QG model in order to distinguish it from recently developed variants (Labbé et al. 2015; Gerick et al. 2020; Jackson & Maffei 2020).

Non-dimensionalization is carried out using the shell thickness $d = s_o - s_i$ as the reference length-scale, the viscous diffusion time $d^2/v$ as the reference timescale, and the temperature contrast between the boundaries $\Delta T = T_i - T_o = T(s_i) - T(s_o)$ as the reference for temperature. Throughout this study, we adopt $\eta = s_i/s_o = 0.35$ suitable for a thick shell such as the Earth's outer core. The gravity $g$ is assumed to be linear with respect to the cylindrical radius such that $g(s) \propto s$ and it is non-dimensionalized based on its value at the outer boundary $g_o = g(s_o)$.

Following Schaeffer & Cardin (2005) and Gastine (2019) it is assumed that the axial velocity $u_z$ varies linearly with $z$ in the direction of the rotation axis, including contributions from the radial velocity $u_r$ and the Ekman pumping $\mathcal{P}$, that is

$$u_z(s, \phi, z) = z[\beta u_s + \mathcal{P}(Ek, \mathbf{u}_{\perp}, \omega_z)], \quad (1)$$

with $\mathcal{P}(Ek, \mathbf{u}_{\perp}, \omega_z)$ the Ekman pumping term deduced from Greenspan's formula in a rigid sphere (see eq. 8),

$$\beta = \frac{1}{h} \frac{dh}{ds} = -\frac{s}{h^2}, \quad (2)$$

and $h \equiv \sqrt{s_o^2 - s^2}$, the half-height of a cylinder aligned with the rotation axis at a radius $s$. The spherical-QG continuity equation then reads

$$\frac{\partial(su_s)}{\partial s} + \frac{\partial u_{\phi}}{\partial \phi} + \beta s u_s = 0. \quad (3)$$

We represent the non-axisymmetric QG-velocity by a streamfunction $\psi$ such that

$$u_s = \frac{1}{s} \frac{\partial \psi}{\partial \phi}, \quad u_{\phi} = \overline{u_{\phi}} - \frac{\partial \psi}{\partial s} - \beta \psi, \quad (4)$$

which ensures that (3) is satisfied. $\overline{u_{\phi}}$ is the remaining axisymmetric zonal flow component. The overbar $\bar{x}$ denotes the azimuthal average of any quantity $x$, such that

$$\bar{x} \equiv \frac{1}{2\pi} \int_{0}^{2\pi} x \, d\phi. \quad (5)$$

The axial vorticity $\omega_z = \mathbf{e}_z \cdot \nabla \times \mathbf{u}$ is then

$$\omega_z = \frac{1}{s} \frac{\partial(s\overline{u_{\phi}})}{\partial s} - \nabla^2 \psi - \frac{1}{s} \frac{\partial(\beta s \psi)}{\partial s}. \quad (6)$$

The dynamics of the axial vorticity can then be described by the axial component of the curl of the momentum equation in cylindrical coordinates, averaged in the z-axis direction, which due to the assumed 2-D form of the contributing fields, may be written

$$\frac{\partial \omega_z}{\partial t} + \nabla_{\perp} \cdot (\mathbf{u}_{\perp} \omega_z) - \frac{2}{Ek} \beta u_s = \nabla_{\perp}^2 \omega_z - \frac{Ra}{Pr} \frac{1}{s_o} \frac{\partial \vartheta_{2-D}}{\partial \phi} + \mathcal{P}(Ek, \mathbf{u}_{\perp}, \omega_z), \quad (7)$$
where the subscript $\perp$ corresponds to the horizontal part of the operators—for example $\nabla_{\perp} = (e_s \cdot \partial_s + r^{-1}e_{\phi} \cdot \partial_{\phi})$, $u_{\perp} = (u_s, u_{\phi}, 0)$—and $P(Ek, u_{\perp}, \omega_z)$ corresponds to the Ekman-pumping contribution (Schaeffer & Cardin 2005) for the non-axisymmetric motions

$$P(Ek, u_{\perp}, \omega_z) = -\Pi \left[ \omega_z - \frac{\beta}{2} u_{\phi} + \beta \frac{\partial}{\partial \phi} u_s - \frac{5s_o}{2h} u_s \right], \quad (8)$$

$$\Pi = \left( \frac{s_o}{Ek} \right)^{1/2} \frac{1}{h^{3/2}}. \quad (9)$$

The non-dimensional control parameters, the Ekman number, the Rayleigh number and the Prandtl number are, respectively, defined by

$$Ek = \frac{\nu}{\Omega d^2}, \quad Ra = \frac{\alpha_T g_o \Delta T d^3}{\kappa \nu}, \quad Pr = \frac{\nu}{\kappa}, \quad (10)$$

where $\alpha_T$ is the thermal expansion coefficient, $\nu$ is the kinematic viscosity and $\kappa$ is the thermal diffusivity.

The z-averaged axial vorticity eq. (7) has to be supplemented by an equation to account for the axisymmetric motions. This is obtained by taking the $\phi$-average of the azimuthal component of the Navier–Stokes equation and reads

$$\frac{\partial \overline{u_{\phi}}}{\partial t} + u_s \frac{\partial \overline{u_{\phi}}}{\partial s} + \frac{\overline{u_s u_{\phi}}}{s} = \nabla_{\perp}^2 \overline{u_{\phi}} - \frac{1}{s^2} \overline{u_{\phi}} - \Pi \overline{u_{\phi}}, \quad (11)$$

where the last term on the right-hand side is the Ekman-pumping contribution for the axisymmetric motions. The boundary conditions for the velocity field are described in detail in Section 2.4.

The other coupled prognostic equation used to complete the system is the QG-temperature perturbation equation,

$$\frac{\partial \vartheta_{2-D}}{\partial t} + \nabla_{\perp} \cdot (\mathbf{u}_{\perp} \vartheta_{2-D}) + \beta u_s \vartheta_{2-D} + u_s \frac{dT_{2-D}^{\text{cond}}}{ds} = \frac{1}{Pr} \nabla_{\perp}^2 \vartheta_{2-D}, \quad (12)$$

where the temperature is written as a perturbation about a mean 2-D conducting state, that is $T_{2-D} = T_{2-D}^{\text{cond}} + \vartheta_{2-D}$, where $T_{2-D}^{\text{cond}}$ is the conducting background state, a solution of $\nabla^2 T_{2-D}^{\text{cond}} = 0$ subject to the chosen boundary conditions. For fixed-temperature boundary conditions at $s_i$ and $s_o$ this yields

$$T_{2-D}^{\text{cond}} = \frac{1}{\ln \eta} \ln[(1 - \eta)s], \quad \frac{dT_{2-D}^{\text{cond}}}{ds} = \frac{1}{s \ln \eta}. \quad (13)$$

Further details on the boundary conditions for the temperature field, including the possibility of heat flux boundary conditions are given in Section 2.4.

## 2.2 Extended hybrid QG-3D model

In this study we follow Guervilly (2010) and Guervilly & Cardin (2016) and go beyond the classical QG model presented in the previous section to develop a hybrid QG-3D model in which the QG perturbation temperature eq. (12) is replaced by the full 3-D temperature equation

$$\frac{\partial \vartheta_{3-D}}{\partial t} + \mathbf{u}_{3-D} \cdot \nabla \vartheta_{3-D} + u_r \frac{dT_{3-D}^{\text{cond}}}{dr} = \frac{1}{Pr} \nabla^2 \vartheta_{3-D}, \quad (14)$$

where $r$ is the spherical radius and $\mathbf{u}_{3-D} = (u_r, u_{\theta}, u_{\phi_{3-D}})$ is the 3-D velocity in spherical coordinates. Similarly to the QG case, we write the temperature as a perturbation temperature about the conducting background state, that is $T_{3-D} = T_{3-D}^{\text{cond}} + \vartheta_{3-D}$. We compute the conducting temperature profile, or dimensionless radial temperature profile, $T_{3-D}^{\text{cond}}$, as the solution of $\nabla^2 T_{3-D}^{\text{cond}} = 0$, which for a fixed temperature contrast between $r_i$ and $r_o$ without internal heating, yields

$$T_{3-D}^{\text{cond}}(r) = \frac{r_i r_i}{r} - r_i, \quad \frac{dT_{3-D}^{\text{cond}}}{dr} = - \frac{r_i r_o}{r^2}. \quad (15)$$

This is the full 3-D version of eq. (13) in spherical geometry. We reconstruct the 3-D velocity field, $\mathbf{u}_{3-D}$ using the conversion between cylindrical and spherical coordinate systems as needed in this equation, from the QG velocity field, such that

$$\left\{ \begin{array}{l} u_r(r, \theta, \phi_{3-D}) = \sin \theta u_s(s, \phi) + \cos \theta u_z(s, \phi, z), \\ u_{\theta}(r, \theta, \phi_{3-D}) = \cos \theta u_s(s, \phi) - \sin \theta u_z(s, \phi, z), \\ u_{\phi_{3-D}}(r, \theta, \phi_{3-D}) = u_{\phi}(s, \phi), \end{array} \right. \quad (16)$$

where $u_z$ is proportional to $z$ and incorporates the effects of the Ekman pumping (see eq. 1) consistent with our initial QG assumption (eq. 3), and where the cylindrical quantities are cast onto the 3-D grid using a bi-linear extrapolation (see Appendix D1 for more details). Inside the Tangent cylinder, the velocities are set to zero and thus only temperature diffusion and thermal wind effects occurs in that region.

Considering a 3-D temperature field, eq. (7) becomes

$$\frac{\partial \omega_z}{\partial t} + \nabla_{\perp} \cdot (\mathbf{u}_{\perp} \omega_z) - \frac{2}{Ek} \beta u_s = \nabla_{\perp}^2 \omega_z - \frac{Ra}{Pr} \left\langle \frac{1}{r_o} \frac{\partial T_{3-D}}{\partial \phi_{3-D}} \right\rangle + P(Ek, \mathbf{u}_{\perp}, \omega_z), \quad (17)$$
where the angular brackets $\langle x \rangle$ refer to the axial or z-average of any quantity $x$ defined by

$$ \langle x \rangle = \frac{1}{2h} \int_{-h}^{h} x \, dz. \qquad (18) $$

The equation for the zonal motions (11) is not modified as it does not involve the temperature.

Since we treat the temperature in 3-D, we can now also take into account the thermal wind contribution to the velocity field which results in an extra term being added to $u_{\phi 3-D}$, and which satisfies

$$ \frac{\partial u_{\phi 3-D}}{\partial z} = \frac{RaEk}{2Pr} \frac{g(r)}{r} \frac{\partial T_{3-D}}{\partial \theta}. \qquad (19) $$

Integrating between the position $z$ and the height of the column above the equator, $h$, we obtain

$$ u_{\phi 3-D}(r, \theta, \phi_{3-D}) = u_{\phi 3-D}(h, \theta, \phi_{3-D}) - \frac{RaEk}{2Pr} \int_{z}^{h} \frac{g(r)}{r} \frac{\partial T_{3-D}}{\partial \theta} dz', \qquad (20) $$

where $g(r) = r/r_o$ is the 3-D gravity field. Here, for efficiency, the thermal wind contribution is assumed to be symmetric about the equatorial plane, although this condition can be relaxed if needed depending on the chosen boundary conditions.

Because of the full 3-D treatment of the heat equation, the consideration of the thermal wind effects and the fact that the thermal boundary conditions are the same as in the 3-D case, it is natural to expect the hybrid QG-3D model to behave better than the classical QG model when compared with a full 3-D model, a hypothesis that will be further assessed in Section 3.

## 2.3 3-D model formulation

In order to compare the results of our QG and hybrid QG-3D models, we also consider a purely 3-D model. Similarly to the two previous setups, we consider convection of a Boussinesq fluid enclosed in a spherical shell of inner radius $r_i$ and outer radius $r_o$ rotating about the z-axis. The same scales and dimensionless parameters are used and thus the 3-D Navier-Stokes equations read

$$ \nabla \cdot \mathbf{u}_{3-D} = 0, \qquad (21) $$

$$ \begin{aligned} \frac{\partial \mathbf{u}_{3-D}}{\partial t} + \mathbf{u}_{3-D} \cdot \nabla \mathbf{u}_{3-D} + \frac{2}{Ek} \mathbf{e}_z \times \mathbf{u}_{3-D} = -\nabla P \\ + \frac{Ra}{Pr} \frac{r}{r_o} \vartheta_{3-D} \mathbf{e}_r + \nabla^2 \mathbf{u}_{3-D}, \end{aligned} \qquad (22) $$

$$ \frac{\partial \vartheta_{3-D}}{\partial t} + \mathbf{u}_{3-D} \cdot \nabla \vartheta_{3-D} + u_r \frac{\mathrm{d} T_{3-D}^{\mathrm{cond}}}{\mathrm{d} r} = \frac{1}{Pr} \nabla^2 \vartheta_{3-D}, \qquad (23) $$

where $P$ is the pressure, and $e_{r,z}$ are, respectively, the unit vectors in the radial and the axial directions.

The velocity field, $\mathbf{u}_{3-D}$, is decomposed into poloidal $W$ and toroidal $Z$ potentials following:

$$ \mathbf{u}_{3-D} = \nabla \times (\nabla \times W \mathbf{e}_r) + \nabla \times \mathcal{Z} \mathbf{e}_r. \qquad (24) $$

## 2.4 Boundary conditions

Since our focus is on modelling the dynamics of the Earth's outer core, we treat the fluid shell as a container with rigid, impenetrable, and co-rotating walls. This implies that in the rotating frame of reference all velocity components should vanish at $s_o$ and $s_i$ in the QG or hybrid models and at $r_i$ and $r_o$ in the 3-D calculations.

Imposing fixed temperature at the boundaries yields, respectively, for the QG- and the 3-D temperature field

$$ \vartheta_{2-D} = 0, \text{ at } s = \{s_i, s_o\}, \qquad (25) $$

$$ \vartheta_{3-D} = 0, \text{ at } r = \{r_i, r_o\}. $$

In this study, the majority of our simulations are conducted under these boundary conditions, but we also consider another set of thermal boundary conditions with a fixed temperature at the inner radius and an imposed flux at the outer boundary. The latter involves

$$ \frac{\partial \vartheta}{\partial r} = 0, \text{ at } r = r_o, \quad \vartheta = 0, \text{ at } r = r_i. \qquad (26) $$

where $\vartheta$ can either be $\vartheta_{2-D}$ or $\vartheta_{3-D}$. The heat flux (or the temperature) may be spatially variable, and any combination of fixed temperature and fixed heat flux at either the inner or the outer boundary can be applied in our model.

Below we present a number of examples performed using a fixed heat flux at the outer boundary and a fixed temperature at the inner boundary (see Section 3.3). With heat-flux boundary conditions at the outer boundary the radial conductive profiles become

$$ T_{2-D}^{\mathrm{cond}-Q} = T_i + Q_o s_o \ln \left( \frac{s}{s_i} \right), \qquad (27) $$
$$T_{3-D}^{\text{cond}-Q} = T_i + Q_o r_o^2 \ln\left(\frac{1}{r_i} - \frac{1}{r}\right), \qquad (28)$$

where $T_i$ and $Q_o$ are, respectively, the temperature at the inner boundary and the heat flux at the outer boundary. The Rayleigh number should then be understood as a flux-based Rayleigh number, that is

$$Ra_Q = \frac{\alpha_T g_o Q_o d^4}{\nu \kappa}. \qquad (29)$$

Lateral variations in the amplitude in the heat flux are then defined by

$$Q^* = \frac{|Q_{\text{max}} - Q_{\text{min}}|}{|Q_o|}. \qquad (30)$$

## 2.5 Numerics

The calculations presented here were carried out using an extension of the open-source pseudo-spectral spherical QG code pizza (Gastine 2019)—available at https://github.org/magic-sph/pizza under the GNU GPL v3 license. The pizza code is written in Fortran, uses a Fourier decomposition in $\phi$ and either Chebyshev collocation (Glatzmaier 1984), or a sparse Chebyshev integration method in s (e.g., Stellmach & Hansen 2008; Muite 2010; Marti et al. 2016). It also supports a number of implicit-explicit (IMEX) time-stepping schemes including multi-step methods (e.g., Ascher et al. 1995) and semi-implicit Runge-Kutta schemes (e.g., Ascher et al. 1997). The reader is invited to consult (Gastine 2019) for further details about the original implementation of pizza and its parallelization.

The purely 3-D simulations were computed with the open-source magnetohydrodynamics code MagIC (Wicht 2002; Gastine et al. 2016) – available at http://www.github.com/magic-sph/magic under the GNU GPL v3 license. Similarly to pizza, MagIC supports various multistep and Runge-Kutta IMEX time schemes.

In this study, 3-D fields in pizza and MagIC are expanded in Spherical Harmonics up to the degree and order $l_{\text{max}}$ in the angular $(\theta, \phi)$ directions and in Chebyshev polynomials with $N_r$ collocation grid points in the radial direction. The 2-D quantities in pizza are expanded in Fourier series up to the degree $N_m$ in the azimuthal direction and in Chebyshev polynomials up to degree $N_s$ in the radial direction. The open-source SHTns<fnref n="1" /> library is used in both codes to handle the Spherical Harmonic Transforms (Schaeffer 2013). Parallelization of the hybrid QG-3D code relies on the Message Passing Interface (MPI) library.

Our numerical implementation follows closely the approach of Schaeffer & Cardin (2005); Guervilly (2010); Guervilly & Cardin (2016), with an important difference that they used finite differences in radius while we resort to a Chebyshev collocation method. In the hybrid setup, the z-extrapolation of the variables from the 2-D grid to the 3-D grid is computed using eq. (16). Reduction of the quantities from the 3-D grid back onto the 2-D grid, and the computation of the thermal wind, relies on two z-integration functions described in the Appendix D1. For clarity, all 3-D quantities are labelled with a subscript 3-D (such as $u_{\phi 3-D}$), QG quantities have no subscripts.

## 2.6 Posterior diagnostics

We next introduce the various diagnostics and notations that are used for the analysis of the simulations. The Nusselt number $Nu$, which characterizes the heat transport of the system, is defined here in the fixed temperature configuration as the ratio between the total heat flux and the heat carried by conduction, that is

$$Nu = 1 + \frac{\frac{d\hat{\vartheta}}{dr}|_{r=r_o}}{\frac{dT^{\text{cond}}}{dr}|_{r=r_o}}, \qquad (31)$$

where the temperature perturbation $\theta$ can be either $\theta_{2-D}$ or $\theta_{3-D}$ and $T^{\text{cond}}$ is either $T_{2-D}^{\text{cond}}$ (13) or $T_{3-D}^{\text{cond}}$ (15) and in QG calculations $r$ and $r_o$ are replaced by $s$, $s_o$. Note that in the heat flux boundary case, this definition of $Nu$ leads to $Nu = 1$ because $\partial_r\theta(r_o) = 0$, so following Goluskin (2016) we instead use

$$Nu_{\Delta} = \frac{T^{\text{cond}}(r_i) - T^{\text{cond}}(r_o)}{\hat{\vartheta}(r_i) - \hat{\vartheta}(r_o)}. \qquad (32)$$

In the above expressions, $\hat{x}$ corresponds to the time average of any quantity $x$, such that

$$\hat{x} = \frac{1}{\tau} \int_{t_0}^{t_0+\tau} x dt, \qquad (33)$$

with $\tau$ the averaging time.

<footnote n="1">¹https://bitbucket.org/nschaeff/shtns</footnote>
The dimensionless kinetic energy, $E_{kin}$ per unit volume, is defined by

$$E_{kin} = \frac{1}{2} \frac{1}{V} \int_V \mathbf{u}_{\perp}^2 dV', \qquad (34)$$

where $V$ corresponds to the full spherical shell volume in the 3-D configurations and the volume outside the tangent cylinder in the QG setups. In the QG case we thus have $dV' = h(s)s ds d\phi$. From this expression, we define a diagnostic for the fluid velocity which characterizes the average flow speed, based on the root-mean-square (r.m.s.) of the velocity, and which is denoted by the Reynolds number

$$Re = \sqrt{2E_{kin}}. \qquad (35)$$

We also define the time-averaged convective Reynolds number, where the axisymmetric zonal flow contribution has been removed, since it can represent a significant fraction of the total kinetic energy without directly contributing to the heat transfer (Gastine et al. 2016),

$$Re_c = \sqrt{2(E_{kin} - E_{zon})}, \qquad (36)$$

where $E_{zon}$ is the dimensionless axisymmetric kinetic energy per unit volume, similar to (34), that is defined by $E_{zon} = 1/2V \int_V \overline{u_{\phi}^2} dV'$ and is associated with the time-averaged zonal Reynolds number,

$$Re_{zon} = \sqrt{2E_{zon}}. \qquad (37)$$

Finally, for the typical flow length-scale we use the typical cylindrical radial velocity length-scale $L_{us}^{-1}$, determined from the time-averaged $u_s$ energy spectrum

$$L_{us}^{-1} = \left( \frac{d \sum_{m=0}^{m_{\max}} m u_{sm}}{\pi \sum_{m=0}^{m_{\max}} u_{sm}} \right). \qquad (38)$$

# 3 RESULTS

We present here results of rapidly rotating convection, focusing on a regime well above the onset of convection, that is $Ra \ge 5 Ra_c$. We explore Ekman numbers from $Ek = 10^{-4}$ down to $Ek = 10^{-10}$, and consider Prandtl numbers from $Pr = 10$ down to $Pr = 10^{-2}$, and moderate-to-high supercriticalities ranging from $Ra = 1.7 Ra_c$ up to $157.3 Ra_c$, reaching $Ra$ as large as $4.83 \times 10^{14}$. As well as hybrid QG-3D simulations we present a large number of purely QG cases; these are more computationally efficient to run and allow a more comprehensive exploration of the parameter space. We also present a collection of fully 3-D runs computed with *Magic* (Wicht 2002) and which have been either specifically computed for this study or taken from Schwaiger et al. (2021). The temporal convergence of the runs has been ensured by running each simulation long enough to obtain a statistical equilibrium of the diagnostics. The numerical truncation ranges from $(N_s, N_m)/(N_r, \ell_{max}) = (97, 96)/(97, 96)$ for the highest/lowest $Ek/Ra$ of $(10^{-4}/2 \times 10^6)$ up to $(N_s, N_m)/(N_r, \ell_{max}) = (9217, 9216)/( -, -)$ for the lowest/highest $Ek/Ra$ numbers $(10^{-10}/4.83 \times 10^{14})$. In total 144 runs have been performed and a list of their key diagnostics is given in Table A1.

Fig. 1 summarizes all the runs we have carried out with fixed temperature contrast for this study in terms of their heat transfer $Nu$ as a function of the applied forcing $Ra$. The colours indicate the Ekman numbers while the different Prandtl numbers explored are indicated with different symbol shapes and transparency. Hybrid QG-3D and purely 3-D runs are marked with hatched and empty symbols, respectively.

## 3.1 Heat transfer

Previous work by Guervilly et al. (2019) has explored the parameters at low Ekman numbers—reaching down to $Ek = 10^{-11}$—and have validated the hybrid approach for a weakly supercritical $Ra$ and low $Pr$ setup (e.g., Guervilly & Cardin 2017; Guervilly et al. 2019) in a full sphere geometry.

Here, we extend the Rayleigh number range to reach higher supercriticalities and restrict ourselves to higher Ekman numbers ($Ek \le 10^{-10}$) with a focus around $Ek = 10^{-6}$ in order to facilitate comparisons with full-3-D simulations. All our runs have $Ro \ll 1$ as is appropriate for QG convection studies. A small number of cases conducted at the highest Ekman numbers have a local Rossby numbers based on the length-scale of the flow which are up to 0.1 (this is discussed further in Section 3.5).

Fig. 2 displays all the Nusselt numbers of our data set versus the supercriticality, $Ra/Ra_c$. The critical Rayleigh numbers have been computed for each configuration using either the Linear Solver Builder package (LSB, Valdettaro et al. 2007) for QG models, or the open-source linear solver SINGE (https://bitbucket.org/vidalje/singe, see Vidal & Schaeffer 2015) for 3-D configurations, although we have used the asymptotic expression of Dormy et al. (2004) for $Pr = 1$ when $Ek < 10^{-7}$ in this latter case. Concerning the hybrid QG-3D model, we have determined the onset for three configurations – at $Pr = 1$ and $Ek = \{10^{-4}, 10^{-5}, 10^{-6}\}$ – by time-integrating the nonlinear eqs (11–14–17) using the pizza code with an initial sectorial temperature perturbation and by bracketing the Rayleigh number until the critical value is attained. When $Ek \le 10^{-7}$ or $Pr \ne 1$ we have assumed the same critical Rayleigh as for 3-D configurations. In all cases, a
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Nusselt number, *Nu* as a function of the Rayleigh number, *Ra*. Summary of all the runs with fixed temperature contrast considered in this study with the various Ekman and Prandtl numbers displayed with, respectively, different colours and symbols. The runs performed using the hybrid approach are represented with hatched symbols and the full 3-D runs with empty symbols.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Nusselt number, *Nu* as a function of the supercriticality, *Ra/Ra*<sub>c</sub>. Various Ekman and Prandtl numbers displayed with, respectively, different colours and symbols. The runs performed using the hybrid approach are represented with hatched symbols and the full 3-D runs with empty symbols.</figcaption></figure>

simple extrapolation using the asymptotic scaling for rotating convection Ra_c ~ Ek<sup>-4/3</sup> has been used whenever the aforementioned techniques could not be applied for practical reasons. Concerning Pr = 1, the Ra_c values obtained with LSB and SINGE methods agree within ~6 per cent at all Ek and follow the expected converging trend (Dormy et al. 2004). The Ra_c value of the hybrid model is ~2 per cent lower than that of the 3-D at Ek = 10<sup>-4</sup> and ~13 per cent lower at Ek = 10<sup>-6</sup>. The m_c values obtained with all methods agree with each other within a range of m ± 2 in all configurations. See Appendix C for more details about our estimates of Ra_c.

We can observe that for the lowest supercriticalities (Ra ≤ 10 Ra_c) all the points in the weakly non-linear regime follow a power law of the form Nu − 1 ~ Ra/Ra_c − 1. For stronger forcing with Ra > 10 Ra_c, the numerical models seem to approach an asymptotic behaviour of the form Nu ∝ (Ra/Ra_c)α (black and blue dotted-line for the QG and 3-D runs, respectively) with no additional dependence on the Ekman number. A simple polynomial fit suggests a power law with a slope of about α ≃ 1.1, an exponent in line with previous findings of rotating convection in spherical shells with r_i/r_o = 0.35 and fixed temperature boundary conditions (e.g., Yadav et al. 2016). This is somewhat lower than the theoretical asymptotic scaling for rapidly rotating convection Nu ~ Ra<sup>3/2</sup>Ek<sup>2</sup>Pr<sup>-1</sup> put forward by Julien et al. (2012) and retrieved in the 3-D spherical shell computations by Gastine et al. (2016) when r_i/r_o = 0.6. The QG runs are slightly offset compared to 3-D cases towards larger Nusselt numbers for the same supercriticality (Ra/Ra_c). Strikingly however, this asymptotic scaling is followed only by the QG and 3-D simulations, while the hybrid runs follow a much shallower slope Nu ~ (Ra/Ra_c)<sup>1/5</sup>. Several outliers also appear in the QG
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Meridional section of the φ-averaged azimuthal velocity $\overline{u}_{\phi 3-D}$ contribution including the effect of the 3-D thermal wind for the 3-D (a) and the hybrid (b) cases. A 3-D extension of the purely QG φ-averaged azimuthal velocity [using eq. (16)] is also presented (c). The three computations have been carried out at the same parameters $Ek = 10^{-6}$, $Pr = 1$, and $Ra = 2 \times 10^9$ (= 10.6 $Ra_c$ for the 3-D, = 12.3 $Ra_c$ for the hybrid, and = 11.3 $Ra_c$ for the QG setups]. The resolution is ($N_r, \ell_{max}$) = (129, 341) in the 3-D case, ($N_s, N_m$)/($N_r, \ell_{max}$) = (257, 256)/(257, 256) in the hybrid case and ($N_s, N_m$) = (385, 384) in the QG case. The same colourbar is used in all cases and is saturated to highlight finer structures.</figcaption></figure>

and 3-D configurations: the series of QG points at $Pr = 10$ seem to follow a different slope with $Nu$ values considerably higher than the values obtained at $Pr \neq 10$ for the same $Ra/Ra_c$ ratio. All the QG and 3-D runs at $Pr = 0.1$ (filled and empty diamonds) lie below the mean trend, suggesting that the heat transport is less efficient for the same supercriticality when $Pr < 1$. Overall, the purely 3-D and QG cases are generally in agreement and we find two different behaviours (diverging around $Ra/Ra_c \sim 10$), a weakly non-linear regime and a regime with steep scaling at higher forcing levels, while the hybrid QG-3D setup starts to significantly depart from the 3-D configuration for $Ra > 10 Ra_c$ (at $Pr = 1$) and follows a much shallower scaling behaviour.

In constrast, when $Pr = 10^{-2}$ the QG and 3-D models disagree above $Ra/Ra_c = 10$, while the range of accordance between the hybrid QG-3D and the 3-D configurations seems to extend to higher $Ra$, up to $Ra > 15 - 20 Ra_c$ at $Pr = 0.1-0.01$, suggesting the range of agreement between the two latter models may be larger at low $Pr$ and low $Ra/Ra_c$ (as suggested in previous studies, e.g., in Guervilly & Cardin 2016, 2017; Guervilly et al. 2019).

## 3.2 Comparison of convective planforms

### 3.2.1 Comparison at modest driving

To further investigate the features observed in Fig. 2, we begin by comparing results from our hybrid QG-3D model with 3-D and QG simulations at modest driving, which we define here to be the parameter regime where $Ra \le 20 Ra_c$, for a case at $Ek = 10^{-6}$, $Pr = 1$ and $Ra = 2 \times 10^9$ (= 10.6 $Ra_c$ for the 3-D, = 12.3 $Ra_c$ for the hybrid, and = 11.3 $Ra_c$ for the QG setups).

In Fig. 3 we present a comparison of meridional sections of the φ-averaged azimuthal velocity $\overline{u}_{\phi 3-D}$ obtained in a 3-D simulation (a), from our hybrid QG-3D model (b) and from the QG model (c). The resolution of the QG case is ($N_s, N_m$) = (385, 384), the hybrid case uses a spatial resolution of ($N_s, N_m$)/($N_r, \ell_{max}$) = (257, 256)/(257, 256), and the resolution in the 3-D case is ($N_r, \ell_{max}$) = (129, 341).

Fig. 4 additionally displays the φ-averaged temperature field $\overline{T}_{3-D}$ (top panel) in the purely 3-D (a), hybrid (b) and QG (c) cases as well as equatorial sections of the z-averaged vorticity $\omega_z$ (lower panel) in 3-D (d), hybrid (e) cases, and in the purely QG case (bottom – f). Note that in the 3-D configuration, this involves a z-averaging procedure—we average over the North Hemisphere only inside the TC—while this is straightforward in QG and hybrid cases. As previously stated the control parameters are strictly the same.

Considering the meridional profiles of the zonal velocities (Fig. 3) we observe that the innermost retrograde jet near the tangent cylinder is slightly offset outwards in the 3-D case while it is very close to the inner boundary in the hybrid case and QG case, creating an artificially strong shear at the tangent cylinder. In both the 3-D and hybrid cases, these two jets display a similar columnar structure that spans the entire height of the shell with the strongest velocity amplitude (compared to the other jets) and that do not vary much with z. In the bulk, beside this fairly geostrophic jet, we find several thinner and weaker jets which are ageostrophic and demonstrate that the thermal wind has an important effect here; these features are reproduced in the hybrid case but not in the QG case.

In the QG case, we retrieve the strongest jet near the tangent cylinder, followed by perfectly geostrophic jets of alternating sign, with prograde jets dominating near the CMB. It is worth noting that $Re_{zon} \ll Re_c$ in all three cases (the exact values are given at the end of this
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top panel: comparison of the meridional section of the $\phi$-averaged of the temperature field $\overline{T}_{3-D}$ for the 3-D (a), the hybrid (b) and the QG case (c). The QG temperature field has been extended in $z$ using the conversion between cylindrical and spherical coordinate systems. Bottom panel: $z$-averaged vorticity for the 3-D simulation (d), and equatorial section of the axial vorticity $\omega_z$ for the hybrid QG-3D simulations (e), and the QG-simulation (f). The three computations have been carried out at the same parameters $Ek = 10^{-6}$, $Pr = 1$ and $Ra = 2 \times 10^9$. The spatial resolution in the 3-D case is $(N_r, \ell_{max}) = (129, 341)$; in the hybrid case is $(N_s, N_m)/(N_r, \ell_{max}) = (257, 256)/(257, 256)$; and in the QG case is $(N_s, N_m) = (385, 384)$ (bottom). For the three temperature and the three vorticity plots, respectively, the same colourscales are used. Note that the colourscale for the vorticity is saturated to highlight the fine structure of the flows.</figcaption></figure>
section) which is consistent with the relatively weak zonal jets found in the bulk. Near the equator, the amplitude of the azimuthal velocity is slightly larger in the 3-D case compared with the hybrid case. Since the velocity inside the tangent cylinder is set to zero before applying the thermal wind approximation in our hybrid approach, significant differences with the 3-D models are visible in that region. Overall however, the hybrid case qualitatively reproduces much of the zonal flow dynamics that happens in the bulk of the 3-D case, although there are discrepancies towards the inner and outer boundaries.

Turning to the azimuthally averaged temperature fields (Figs 4a–c), we find that the profiles in the 3-D and hybrid QG-3D cases are very similar with isothermal lines that are bent across the tangent cylinder and that extend in the equatorial plane. These isothermal lines are slightly more squeezed towards the equatorial plane in the 3-D case compared to the hybrid case and there is a difference in the spacing of the isotherms in the z-direction, likely due to the simple relationship we have used to reconstruct $u_z$. For these parameters, the 3-D temperature profile is rather well retrieved in the hybrid case, in contrast to the QG case which does not have the correct temperature profile and displays a largely homogeneous temperature in the bulk and a sharp drop toward the outer boundary.

Figs 4(c)–(e) shows a comparison of the axial vorticity $\omega_z$ between the purely 3-D, the hybrid, and the purely QG cases. Considering the shell from the inner core to mid-depths, we find it hard to distinguish the three planforms of convection which all display filaments of vorticity of similar amplitude and length-scales, wider near the inner boundary, and sheared in the azimuthal direction with a gradual reduction of the convective cells size with increasing $s$. Closer to the outer boundary, the convective pattern changes in all cases with the filaments becoming more radially elongated. This transition occurs at about the same radius $s \sim s_o - 1/3$ in each case. Obvious differences are seen approaching the outer boundary of the container. In both the hybrid and the purely QG cases the velocity field transitions into elongated azimuthal structures typical of thermal Rossby waves. The vorticity in the 3-D case, on the other hand, becomes almost perpendicular to the outer boundary with very thin and radially elongated filaments. The discrepancy may reflect a fundamental difference in the boundary geometry between the different configurations: in both the QG and the hybrid model the slope of the container $|\beta|$ (2) increases with the cylindrical radius. This treatment impedes radial motions and favors the propagation of thermal Rossby waves over the advective processes such as the stretching term due to $\beta u_s$ which becomes the dominant source of axial vorticity because of the steepening of the slope at large radii; a phenomenon expected to weaken with an increasing forcing (Guervilly & Cardin 2017). Other QG implementations that also incorporate the horizontal components of vorticity have recently been developed and may perform better in this low latitude region (see, e.g. Labbé et al. 2015; Maffei et al. 2017; Gerick et al. 2020).

Global diagnostics in the 3-D and hybrid cases are rather similar with a convective Reynolds $Re_c$, a zonal Reynolds $Re_{zon}$ and a Rossby number $Ro$ of, respectively, 448.1, 97.1 and $4.58 \times 10^{-4}$ in the 3-D case, and 404.4, 86.9 and $4.14 \times 10^{-4}$ in the hybrid case. The Nusselt number $Nu$ differs more strongly with values of 2.38 in the 3-D case and of 1.67 in the hybrid case. The same diagnostics obtained for the purely QG simulation are $550.4, 126.6$ and $5.65 \times 10^{-4}$, respectively for $Re_c, Re_{zon}$ and $Ro$ while $Nu = 5.05$. This example indicates how the hybrid approach is capable of accounting the 3-D convective dynamics happening in the fluid bulk at modest driving (here with $Ra = 10.6 Ra_c$).

### 3.2.2 Limitations of the hybrid approach

We now compare results from our hybrid QG-3D model with 3-D and QG simulations for a more strongly driven case at $Ek = 10^{-6}$, $Pr = 1$ and $Ra = 10^{10}$ (= 53.2 $Ra_c$ for the 3-D, = 61.3 $Ra_c$ for the hybrid and = 63.4 $Ra_c$ for the QG setup). Below we use the term 'strong driving' to refer to the parameter regime where $Ra > 20 Ra_c$.

In Fig. 5 we present the $\phi$-averaged temperature field $\overline{T}_{3-D}$ (top panel) in the purely 3-D (a), hybrid (b) and QG (c) cases as well as the z-averaged vorticity $\omega_z$ in the 3-D case (d), and equatorial sections of the axial vorticity in the hybrid (e), and the QG cases (f). The resolution in the 3-D, hybrid and QG cases is, respectively $(N_r, \ell_{max}) = (321, 682)$, $(N_s, N_m)/(N_r, \ell_{max}) = (513, 512)/(513, 341)$ and $(N_s, N_m) = (577, 576)$.

Compared to the previously moderately forced case, the meridional sections of the temperature field now significantly differ. The hybrid temperature stays similar to the previous case at $Ra = 2 \times 10^9$, while in the purely 3-D case we find the temperature is better mixed with isotherms further away from each others and less contrast in the fluid bulk. The QG temperature profile still does not present the correct temperature variation across the bulk but is more homogeneous and displays a sharper contrast toward the CMB when compared with the lower forced case. Rapid variations of the isotherms close to the boundaries in the 3-D and QG cases indicate the formation of thermal boundary layers, whereas the hybrid model has not developed such layers. Similarly, looking at sections of the axial vorticity, the hybrid and 3-D cases now show significant differences. The vorticity in the 3-D case is much stronger than in the hybrid case, with filaments that are more sheared in the azimuthal direction and with significantly perturbed thin Rossby waves near the outer boundary. The hybrid case has in contrast not departed far from the previously moderately driven case, the main difference being that the convective motions now span the entire shell with larger convective cells. We also observe that convection has started inside the tangent cylinder in the 3-D configuration and is already vigorous at these parameters. Interestingly, the QG case seems to be closer to the 3-D case than the hybrid, with filaments of vorticity strongly sheared in the $\phi$-direction and a vigorous convective pattern degenerating into thin Rossby waves towards the outer boundary. Overall, the purely 3-D and QG cases have reached a regime of vigorous convection with a well-mixed temperature background while the hybrid case displays much weaker convection, comparable to the modest driving regime. These differences are also observed in the global diagnostics with $Re_c, Re_{zon}, Ro$ and $Nu$ that are respectively equal to 3339.8, 2405.2, $4.12 \times 10^{-3}$ and 19.5 in the 3-D case; 1455.9, 829.5, $1.68 \times 10^{-3}$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top panel: comparison of the meridional section of the φ-averaged of the temperature field $\overline{T}_{3-D}$ for the 3-D (a), the hybrid (b) and the QG case (c). The QG temperature field has been extended in z using the conversion between cylindrical and spherical coordinate systems. Bottom panel: z-averaged vorticity for the 3-D simulation (d), and equatorial section of the axial vorticity $\omega_z$ for the hybrid QG-3D simulations (e), and the QG-simulation (f). The three computations have been carried out at the same parameters $Ek = 10^{-6}$, $Pr = 1$ and $Ra = 1 \times 10^{10}$ (= 53.2 $Ra_c$ for the 3-D, = 61.3 $Ra_c$ for the hybrid and = 63.4 $Ra_c$ for the QG setup). The resolution in the 3-D, hybrid and QG cases is, respectively ($N_r, \ell_{max}$) = (321, 682), ($N_s, N_m$)/($N_r, \ell_{max}$) = (513, 512)/(513, 341) and ($N_s, N_m$) = (577, 576). For the three temperature and the three vorticity plots, respectively, the same colourscales are used. Note that the colourscale for the vorticity is saturated to highlight the fine structure of the flows.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Kinetic energy spectra of the vertically averaged total velocity, $\frac{1}{2}u_{\perp}^2 = \frac{1}{2}(u_s^2 + v_{\phi}^2)$, as a function of the order $m + 1$, for the 3-D (cyan plain line), hybrid (green dotted-dashed line) and QG (orange dotted line) models, for a configuration conducted with the same parameters in the three cases $Ek = 10^{-6}$, $Pr = 1$ and $Ra = 1 \times 10^{10}$. The resolution in the 3-D, hybrid and QG cases is, respectively $(N_r, \ell_{\max}) = (321, 682)$, $(N_s, N_m)/(N_r, \ell_{\max}) = (513, 512)/(513, 341)$ and $(N_s, N_m) = (577, 576)$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top panel: spherical time-averaged radial temperature profiles $\hat{T}_{3-D}(r)$ for four cases conducted at the same parameters at $Ek = 10^{-6}$ and $Pr = 1$ using our hybrid method (dotted lines) or a 3-D model (full lines) with varying $Ra$ from $Ra = 10^9$ (green curves) up to $Ra = 2.66 \times 10^{10}$ (red curves). Bottom panel: same as above but for the cylindrical radial z-averaged temperature profiles $\langle T_{3-D} \rangle (s)$. Note that the latter profiles are not time-averaged and were derived from snapshots.</figcaption></figure>

and 2.20 in the hybrid case; and 2584.4, 1845.5, 3.18 × 10⁻³ and 36.3 in the QG case. Between the $Nu$ and $Re_c$ numbers of the hybrid and 3-D models, we have observed a relation of the form $Nu^{Hyb}/Nu^{3-D} \propto (Re_c^{Hyb}/Re_c^{3-D})^{2/5}$.

Examining the kinetic energy spectra for the horizontal velocity, shown in Fig. 6, we find that the QG, 3-D and hybrid models show similar decreasing slopes up to $m \sim 200$, although there is less energy in the hybrid model. The spectrum for the hybrid configuration also shows more steeply decreasing slope at large $m$, confirming the lack of power at small length-scales already seen in the convective planforms.

Fig. 7 shows the time-averaged radial temperature profiles $\hat{T}_{3-D}(r)$ (a) and example snapshots of z-averaged cylindrical temperature profiles $\langle T_{3-D} \rangle (s)$ (b) obtained with the purely 3-D and the hybrid setups at $Ek = 10^{-6}$ and $Pr = 1$ for a series of increasing supercriticalities, ranging from $Ra = 10^9 = 5.3 Ra_c$ up to $Ra = 2.66 \times 10^{10} = 141.5 Ra_c$. We retrieve the fact that the profiles are fairly similar when $Ra \le 2 \times 10^9 \sim 10 Ra_c$ (green and cyan curves), while the hybrid and purely 3-D temperature profiles diverge significantly as $Ra$ is further increased. We observe the formation of thermal boundary layers at both spherical shell boundaries with a well-mixed interior in the 3-D case while the hybrid temperature profiles do not vary much (red and yellow curves) and stay close to the radial conducting state (blue dashed-curve). The same conclusions can be drawn when looking at the z-averaged cylindrical radial temperature profiles (Fig. 7b) although
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Buoyancy power $\frac{Ra}{Pr} u_r g T_{3-D}$, decomposed into its equatorially symmetric (ES, $u_r^{ES} g T_{3-D}^{ES}$) and equatorially antisymmetric (EAS, $u_r^{EAS} g T_{3-D}^{EAS}$) for the 3-D (cyan curves) and the hybrid (green curves) models for two different simulations carried out at $Ek = 10^{-6}$, $Pr = 1$ and $Ra = 2 \times 10^9$ (a), and at $Ek = 10^{-6}$, $Pr = 1$ and $Ra = 1 \times 10^{10}$ (b).</figcaption></figure>

here we can observe that the inner boundary temperature decreases with increasing Ra in both setups with a slightly larger decrease of $<T_{3-D}>$ ($s_i$) in the 3-D case compared to the hybrid case. Note the increased activity inside the tangent cylinder in the 3-D case when $Ra \ge 10^{10}$ which is not accounted for in our hybrid QG-3D model. The evolution of $<T_{3-D}>$ ($s_i$) in both the 3-D and hybrid cases is important: the temperature at the ICB drops significantly when Ra is increased, suggesting that approximating the fixed 3-D temperature at $r_i$ by a fixed temperature at $s_i$ in the QG-approach [about $T_{2-D}(s_i) \sim 0.445$ at all Ra] is rather crude (again see Fig. 5c). Because of this reduction of temperature with increasing Ra at the inner boundary in cylindrical coordinates, we expect that the driving in the purely QG configuration will eventually disagree with the 3-D setup at very large Rayleigh numbers before the 3-D convection reaches the non-rotating regime.

The fundamental difference between the 3-D and the QG/hybrid models lies in the assumed QG nature of the velocity field and in the use of z-averaging to represent the convective dynamics. This implies that both the equatorially antisymmetric parts and the z-component of the buoyancy force are missing in the QG and hybrid models compared to the full 3-D setup. An important difference between the hybrid and the QG setups is that the boundary conditions at the inner core are applied over the entire tangent cylinder in the QG case, artificially providing more power to the QG setup compared to the hybrid or the 3-D configurations. This enables the QG model to transition more quickly towards a turbulent state as Ra is increased.

Convective power is locally given by the quantity $\frac{Ra}{Pr} u_r g T_{3-D}$, which we can decompose into its equatorially symmetric (ES) and equatorially antisymmetric (EAS) components. Fig. 8 presents results concerning the integrated ES and EAS convective power profiles as a function of the spherical radius obtained with the 3-D and the hybrid models for the two cases of Section 3.2.1 and Section 3.2.2, that is at $Ek = 10^{-6}$, $Pr = 1$, $Ra = 2 \times 10^9$ and at $Ek = 10^{-6}$, $Pr = 1$, $Ra = 1 \times 10^{10}$, respectively. At $Ra = 2 \times 10^9$, we retrieve a similar buoyancy power profile for the hybrid and 3-D models, although the hybrid model has less energy especially towards the CMB. EAS modes become noticeable between $Ra = 2 \times 10^9$ (where they are almost zero) and $Ra = 5.5 \times 10^9$ (where they account for 9 per cent of the total power) and grow increasingly strong, reaching 23 per cent of the convective power by $Ra = 10^{10}$. In addition, we find that 10 per cent of the power is driven by the $\frac{Ra}{Pr} u_z g_z T_{3-D}$ at $Ra = 2 \times 10^9$ which grows to 30 per cent at $Ra = 10^{10}$.

At $Ra = 2 \times 10^9$, we find that the peak-to-peak ratio between the convective power in the 3-D and hybrid models is around 1.5 from which (based on the IAC scaling, see Section 3.5 and Fig. 13) we expect a ratio between the velocities of $1.5^{2/5} \sim 1.17$, close to the actual ratio of $Re_c^{3-D}/Re_c^{Hyb} = 448.1/404.4 \sim 1.11$. At $Ra = 1 \times 10^{10}$, we have $Re_c^{3-D}/Re_c^{Hyb} = 3339.9/1455.9 \sim 2.29$ which requires a ratio of $2.29^{5/2} \sim 7.97$ in terms of the convective power. The missing power due to the EAS modes and the z-component of the buoyancy force are alone not sufficient to explain all of this difference. In the strongly driven regime, it is possible that the lack of convection inside the tangent cylinder (although this represents only 15–20 per cent of the total volume of the shell) and the enforced linearity of $u_z$ (see eq. 1) may also contribute to missing power at high $Ra/Ra_c$, but it is difficult to separate these contributions given the models in our database.

To summarize, our results suggest that the differences in the transition to turbulence in the 3-D, QG and hybrid come from differences in the underlying convective power. The Hybrid and QG models lack the equatorially antisymmetric and the z-component of the convective power. This leads to a delayed transition to turbulent flow in the hybrid model. Differences are less noticeable in the QG case, likely because the inner thermal boundary condition is applied over the entire tangent cylinder. At strong forcing the lack of convection inside the tangent cylinder and assumed linearity of $u_z$ in the hybrid model may also play a role.

The above limitations result in the hybrid setup remaining in the weakly non-linear regime with only a small increase of the heat transport and of the velocity with increasing Ra. This is in line with the previously observed discrepancies in the heat transport (see Fig. 2) with Nusselt numbers that stay on a lower slope in the hybrid case than in the purely QG and the 3-D cases. In the hybrid configurations, the thermal-boundary layers do not fully develop, and the temperature profiles do not significantly depart from the background conducting
state, which translates into weaker convection when the forcing is increased. In the QG configuration the missing buoyancy power is partly offset by additional power provided by the cylindrical boundary conditions. In practice, this means the hybrid QG-3D method has a range of good agreement with full 3-D computations which is limited to $Ra \le 10 \, Ra_c$, at $Pr = 1$, and has a decreasing predictive capacity with increasing $Ra$. On the contrary, we cannot point to important disagreements between the QG and 3-D data sets, suggesting that the QG model retains its predictive power, even at high $Ra/Ra_c$, at least in this particular configuration. The impact of that conclusion on the whole data set will be further discussed in Section 3.5.

### 3.3 Influence of laterally varying heat flux boundary conditions

In this section we present examples of rapidly rotating convection with an imposed laterally varying heat flux condition at the outer boundary and a fixed temperature condition at the inner boundary. Under these conditions, 26 additional runs have been performed and key control parameters are given in Table B1.

We first present in detail a comparison of results obtained from QG, 3-D and hybrid QG-3D calculations carried out at the same parameters $Ek = 10^{-6}$, $Pr = 1$ and $Ra_Q = 8 \times 10^9$ with a $\ell = m = 2$, $Q^* = 3$ lateral variation about the imposed heat-flux condition. The resolutions in the 3-D, hybrid and QG cases are, respectively ($N_r$, $\ell_{max}$) = (321, 341), ($N_s$, $N_m$)/($N_r$, $\ell_{max}$) = (385, 416)/(385, 416) and ($N_s$, $N_m$) = (513, 512).

Fig. 9 displays example snapshots of the 3-D temperature field at the CMB $T_{3-D}(r_o)$ using heat-flux boundary conditions for the 3-D model in (a) and the hybrid model in (b). Both cases show the expected $\ell = m = 2$ variation due to regions of higher and lower heat flux, and we also see the imprint of the underlying convection linked to the regions of enhanced heat-flux at the equator alternating with large quiescent regions of high temperature associated with lower heat-flux. The amplitude of the temperature anomalies in the hybrid case is however larger with temperature variations in the hybrid and 3-D cases spanning $-1.49 \le T_{3-D}(r_o) \le 0.38$ and $0.23 \le T_{3-D}(r_o) \le 1.11$, respectively. The imprint of the underlying convection is more clearly seen in the larger eddies evident in the hybrid case, especially at mid-to-high latitudes. This is consistent with the weaker convective forcing that occurs in the hybrid compared with the 3-D configuration, as discussed in Section 3.2.2. This is also reflected in global diagnostics, with lower $Nu_\Delta$ values of 2.03 in the hybrid case compared with 7.15 in the 3-D case. The sharp transition from convective to diffusive-only dynamics inside the tangent cylinder due to our hybrid implementation is again obvious in Fig. 9(b).

Turning to the convective dynamics of the column-averaged axial vorticity (Figs 9c-e) we observe significant differences between the three cases. The vorticity structures up to $r \sim 2/3r_o$ are qualitatively similar in the 3-D and the QG cases (c and e) while the hybrid case (d) displays larger scale vortices and a less turbulent structure, consistent with it being less strongly driven. Close to the outer boundary, all cases show the expected $m = 2$ lateral variation with alternating regions of weak and enhanced convection, although there are differences in the exact locations and morphologies of these regions in the presented snapshots.

At these parameters the hybrid approach fails to retrieve the small length-scale convective structures at high latitudes found in the 3-D case and has a noticeable temperature offset at the outer boundary. On the other hand it does reproduce similar, albeit larger scale, structures at mid-to-low latitudes. This is evident for example in the plots of the temperature anomaly in the top panels of Fig. 9, towards the centre of the images and close to the equator, where signatures of vigorous convection and related wave structures are seen. The observed differences in the bottom panels of Fig. 9 can largely be attributed to differences of the buoyancy power in the three cases. At these parameters our hybrid approach is unable to drive convection which is as turbulent as that seen in the 3-D case, while the purely QG is slightly over-driven compared with the 3-D case. Note that standard global diagnostics can be misleading here as they involve averages over the entire shell.

Fig. 10 explores further the impact of laterally varying heat flux boundary conditions, showing results from a more extreme convective regime using the QG model which is computationally least expensive. It presents examples of equatorial snapshots of the axial vorticity $\omega_z$ in four cases with top heat flux/bottom-temperature imposed boundary conditions: (a) $Ek = 10^{-6}$, $Pr = 1$, $Ra_Q = 4 \times 10^9$ with no lateral variations; (b) same parameters ($Ek = 10^{-6}$, $Pr = 1$, $Ra_Q = 4 \times 10^9$) but with a $m = 2$, $Q^* = 3$ lateral variation; (c) same parameters ($Ek = 10^{-6}$, $Pr = 1$, $Ra_Q = 4 \times 10^9$) but with a $m = 1$, $Q^* = 3$ lateral variation; and $Ek = 10^{-6}$, $Pr = 1$, $Ra_Q = 3.6 \times 10^{10}$ with a $m = 2$, $Q^* = 3$ lateral variation (d). Cases (a-c) have thus the same parameters but different lateral heat-flux conditions applied while case (d) has the same lateral heat-flux conditions as case (b) but is approximately ten times more supercritical. The case without lateral variations in (a) features very similar convective patterns compared with fixed temperature boundary conditions setups (see e.g., Fig. 4). Namely, filaments of axial vorticity are sheared in the bulk by zonal jets of alternating direction, becoming thinner, and degenerating into thermal Rossby waves towards the outer boundary.

The bottom panel demonstrates that the laterally varying heat flux has been successfully imposed and can drastically modify the convective planform when these lateral variations are sufficiently large. In the case with a $m = 1$ and $Q^* = 3$ pattern (c), we observe that the right hemisphere is not convecting above $s \sim s_o - 1/3$ and displays only a wide spiraling arm covering this region, a result similar to previous 3-D studies (see, e.g. the fig. 4 of Mound & Davies 2017). For the case with a $m = 2$ pattern with $Q^* = 3$ (b), we observe a similar behaviour with regions of weak convection dominated exclusively by azimuthal motions near the outer boundary, as was also seen in Fig. 9(e) with $Ra_Q = 8 \times 10^9$. The boundary perturbation does not penetrate very deep in the shell at these parameters ($s \sim s_o - 1/3$ in the $m = 1$ case, and $s \sim s_o - 1/4$ in the $m = 2$). The final case (d) has been conducted at a larger forcing ($Ra_Q = 3.6 \times 10^{10}$) and displays similar, although more turbulent, features compared to the previous cases. The region of weak convection is however smaller and limited to fluid regions above
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Top panel: comparison of the temperature field at the CMB $T_{3-D}(r_o)$ for (a) the 3-D case and (b) the hybrid case. Bottom panel: z-averaged vorticity for the 3-D simulation (c), and equatorial section of the axial vorticity $\omega_z$ for the hybrid QG-3D simulations (d) and the QG-simulation (e). The three computations have been carried out at the same parameters $Ek = 10^{-6}$, $Pr = 1$ and $Ra_Q = 8 \times 10^9$ imposing a $\ell = m = 2$, $Q^* = 3$ lateral variation about the imposed heat-flux condition. The resolutions in the 3-D and hybrid cases are, respectively $(N_r, \ell_{max}) = (321, 682)$ and $(N_s, N_m)/(N_r, \ell_{max}) = (385, 416)/(385, 416)$ and is $(N_s, N_m) = (513, 512)$ in the QG case. For the three vorticity plots, the same colourscale is used and is saturated to highlight the fine structure of the flows.</figcaption></figure>

$s \sim s_o - 1/6$. It may be that in the limit of very large $Ra_Q$ convection, the region affected by the inhomogeneous boundary conditions could shrink to a very thin layer close to the outer boundary.

Our results demonstrate that imposing a fixed heat flux at the outer boundary does not drastically change the QG-convection compared to a fixed temperature boundary condition. However, imposing a lateral variation of heat flux at the CMB certainly can inhibit the convective motions in a region near the surface whose size depends on $Q^*$ and the supercriticality, consistent with the findings in 3-D computations (Mound & Davies 2017).

## 3.4 Influence of the Prandtl number at low Rossby number

Focusing on the hybrid and 3-D series at $Ek = 10^{-6}$, $Pr = 0.1$, we find the same limitations as in Section 3.2.2, with an apparent lack of energy in the hybrid configuration when $Ra$ is increased (see Fig. 2). However as $Pr$ is decreased, larger velocities are attained at smaller Nusselt numbers, and the range of agreement across the hybrid and 3-D configurations expands (up to $Ra \sim 15 \times Ra_c$ at $Pr = 0.1$). We observed similar convective patterns between the hybrid and 3-D models at all $Pr$ (1, $10^{-1}$, $10^{-2}$) when $Ra \le 15 \times Ra_c$. Unfortunately, the lowest Prandtl runs are extremely computationally costly to run, especially with the 3-D approach, because powerful zonal flows and velocities are triggered even at low $Nu$ inducing difficulties to reach a converged power balance. For example, for a run at $Ek = 10^{-6}$, $Pr = 0.01$, $Ra = 5 \times 10^8 = 19.3 Ra_c$,
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Equatorial snapshots of axial vorticity, $\omega_z$, for purely QG runs with heat-flux imposed at the top and a fixed temperature imposed at the bottom: $Ek = 10^{-6}$, $Ra_Q = 4 \times 10^9$, $Pr = 1$ with no lateral variations (a); same parameters ($Ek = 10^{-6}$, $Ra_Q = 4 \times 10^9$, $Pr = 1$) but with a $m = 2$, $Q^* = 3$ lateral variation (b); same parameters ($Ek = 10^{-6}$, $Ra_Q = 4 \times 10^9$, $Pr = 1$) but with a $m = 1$, $Q^* = 3$ lateral variation and (c) and $Ek = 10^{-6}$, $Ra_Q = 3.6 \times 10^{10}$, $Pr = 1$ with a $m = 2$, $Q^* = 3$ lateral variation (d). The cases (a–c) have the same colourscale displayed on the top right and the case (d) has its own colourscale displayed on the bottom right. Note that the colourscales are saturated to highlight the fine structure of the flows. The four cases are purely QG and their spatial resolution is ($N_s$, $N_m$) = (513, 512) in all cases.</figcaption></figure>

the convective Reynolds number reaches 11 968.4 for a $Nu$ of 2.07 in the 3-D case; $Re_c$ values that were not even reached by $Ra = 140 Ra_c$ at $Pr = 1$. The two 3-D runs at $Pr = 10^{-2}$ however appear closer to the hybrid trend of Fig. 2 with comparable $Nu$ up to $Ra/Ra_c = 20$, whereas the QG model departs from the 3-D trend around $Ra \sim 10 Ra_c$ with much higher $Nu$ in the QG cases. Our results do thus indicate better agreement between the hybrid and the 3-D configurations at low $Pr$ (as also observed by Guervilly et al. 2019); in contrast the QG configuration performs less well in this regime. We expect nonetheless the hybrid model to depart from the 3-D when $Ra$ is sufficiently increased, as observed for our $Pr = 0.1$ series. The exact $Ra/Ra_c$ range of agreement when $Pr$ is decreased further below 0.1 requires more 3-D and hybrid runs in this challenging regime in order to be determined.

We next take advantage of the computationally more efficient purely QG setup to study convective flows at more extreme parameters, that is at lower $Ek$ and higher $Ra$. In particular, we investigate here the impact of the Prandtl number in this regime. Fig. 11 shows example snapshots of axial vorticity $\omega_z$ (top panel) and of azimuthal velocity $u_\phi$ (bottom panel) in the equatorial plane for two cases: $Ek = 10^{-8}$, $Pr = 1$, $Ra = 8.99 \times 10^{12} = 142 Ra_c$ (left-hand column) and $Ek = 10^{-8}$, $Pr = 10^{-2}$, $Ra = 3.37 \times 10^{10} = 6.1 Ra_c$ (right-hand column). Despite the much higher supercriticality attained in the $Pr = 1$ case, these two runs have comparable Rossby numbers: $4.10 \times 10^{-4}$ for the $Pr = 1$ case and $5.97 \times 10^{-4}$ for the $Pr = 10^{-2}$ case. Both cases are purely QG calculations, and 3-D temperature effects have not been included.

In both cases, the azimuthal velocity (c and d) displays multiple zonal jets of alternating direction (blue is retrograde and red is prograde flows), which directly translate in the axial vorticity (a and b) into alternating rings of cyclonic ($\omega_z > 0$ in red) and anticyclonic ($\omega_z < 0$ in blue) vortices. Between two alternating jets, the vortices are streched out and sheared into azimuthally elongated filaments, which involve a direct cascade of energy from the large to the small length-scales (Rhines 1975; Gastine 2019). Potential vorticity, $(\omega_z + 2/Ek)/h$, is mixed due to stirring by the turbulent motions, and creates these characteristic concentric jets with a typical size that is approximately predicted by the Rhines-scale ($Ro/\beta)^{1/2}$ (Rhines 1975). Closer to the boundary, we see the influence of the slope and the $\beta$-effect where the steepening of the curvature of the container impedes the radial advection of the vortices and causes the dynamics to degenerate into azimuthally elongated
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Left: Equatorial snapshot of the axial vorticity $\omega_z$ (a) and of the azimuthal velocity $u_\phi$ (c), for a numerical simulation with parameters $Ek = 10^{-8}$, $Pr = 1$ and $Ra = 8.99 \times 10^{12} = 142 Ra_c$. For this case the spatial resolution is $(N_s, N_m) = (3073, 3072)$. Right: Equatorial snapshot of the axial vorticity $\omega_z$ (b) and of the azimuthal velocity $u_\phi$ (d), for a numerical simulation with parameters $Ek = 10^{-8}$, $Pr = 10^{-2}$ and $Ra = 3.37 \times 10^{10} = 4.1 Ra_c$. The spatial resolution in that case is $(N_s, N_m) = (3457, 3456)$. Respectively, for the $\omega_z$ and $u_\phi$, the same colourbars are used for both cases and are saturated.</figcaption></figure>

motions typical of thermal Rossby waves, as previously observed in Section 3.2.1. Fig. 12 additionally shows the time-averaged cylindrical radial profiles of potential vorticity along with the time-averaged zonal flow. Retrograde zonal jets where potential vorticity gradients are slightly stronger (marked by white stripes in Fig. 12), seem narrower than in the regions where the gradients are weaker (corresponding to prograde jets), a result already found by Guervilly & Cardin (2017). Since large supercriticalities are required to obtain well-formed potential vorticity staircases, the case with $Pr = 10^{-2}$ does not show a comparable degree of homogenization, due to the significantly lower $Ra/Ra_c$ reached in that case, despite similar values of $Ro$.

Besides the two cases having comparable Rossby numbers and time-averaged kinetic energy spectra, the $Pr = 1$ case (Figs 11a-c) displays smaller eddies and thinner jets than the $Pr = 10^{-2}$ case with a larger number of coherent jets (8 in the $Pr = 1$ case compared to 3 in the $Pr = 10^{-2}$ case). Near the outer boundary the transition of the dynamics into thermal Rossby waves happens deeper in the shell in the $Pr = 10^{-2}$ compared to the $Pr = 1$ case. This transition is also visible in the $\phi$-velocity where the jets lose their coherence around $s \sim s_o - 1/3$; a direct consequence of the lower supercriticality attained in the $Pr = 10^{-2}$ case ($Ra = 6.1 Ra_c$ in that case, compared to $Ra = 142 Ra_c$ for the $Pr = 1$ case).

In the low Rossby regime explored here ($Ro < 6 \times 10^{-4}$), changing the Prandtl number by a factor 100 drastically modifies the form of the convective pattern. We find that decreasing $Pr$ results in fewer and wider jets as well as larger convective structures that are maintained at a much lower supercriticality. This effect has previously been reported by Guervilly & Cardin (2017) who suggested a weak dependence of $Ro$ on $Pr$ (see Table A1 or the next section).

## 3.5 Scaling laws for rapidly rotating convection

We now finally explore the scaling behaviour of rapidly rotating turbulent convection in our three model setups. Theoretical scaling laws of rotating convection can be derived by considering the following dimensional 3-D vorticity equation

$$ \frac{\partial \omega}{\partial t} + (\mathbf{u} \cdot \nabla) \omega + 2\Omega \cdot \nabla \mathbf{u} = \nabla \times (\alpha_T \vartheta_{3,D} \mathbf{g}) + \nu \nabla^2 \omega. \quad (39) $$

In the limit of rapid rotation, the Proudman-Taylor theorem promotes $z$ invariant flows with $l_\perp \ll l_\parallel$, where $l_\perp$ and $l_\parallel$ correspond to the convective flow length-scale perpendicular and parallel to the rotation axis. Assuming that $l_\parallel \sim d$, this implies that the gradients orthogonal
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Time averaged radial profiles of the the planetary vorticity $2Ek^{-1}/h$ (black-dotted curve), potential vorticity $(\omega_z + 2/Ek)/h(s)$ (orange curve) and the azimuthal velocity $\dot{u}_\phi(s)$ (blue curve); for the numerical simulation with parameters $Ek = 10^{-8}$ $Pr = 1$, $Ra = 8.99 \times 10^{12} = 142 Ra_c$ and a spatial resolution of $(N_s, N_m) = (3073, 3072)$ (a), and the numerical simulation with parameters $Ek = 10^{-8}$, $Pr = 10^{-2}$, $Ra = 3.37 \times 10^{10} = 6.1 Ra_c$ and a spatial resolution of $(N_s, N_m) = (3457, 3456)$ (b). The shaded area and the white stripes correspond, respectively, to the fluid regions where the azimuthal velocity is positive (prograde jets) and negative (retrograde jets). The planetary vorticity (black-dotted curve) indicates the asymptotic behaviour of the potential vorticity in a perfectly mixed shell.</figcaption></figure>

to the axis of rotation $\nabla_{\perp}$ can be approximated by $1/l_{\perp}$, while the axial gradients $\partial/\partial z$ simply scale as $1/d$. It also follows that $\omega \sim U_c / l_{\perp}$, where $U_c$ is a typical convective velocity.

In the diffusivity-free limit relevant for planetary convective cores, the dominant terms entering eq. (39) involve a triple balance between Coriolis, Inertia and Archimedean forces (Hide 1974; Ingersoll & Pollard 1982; Cardin & Olson 1994; Gillet & Jones 2006)

$$2\Omega \cdot \nabla u \sim u \cdot \nabla \omega \sim \nabla \times (\alpha_T \vartheta_{3-D} g). \qquad (40)$$

or in terms of scaling quantities

$$\Omega \frac{U_c}{d} \sim U_c \frac{\omega_c}{l_{\perp}} \sim \alpha_T g \frac{\Theta}{l_{\perp}}, \qquad (41)$$

where $\Theta$ is a typical temperature perturbation. The balance between Coriolis and Inertia terms yields

$$\frac{l_{\perp}}{d} \sim (Re_c Ek)^{1/2} \sim Ro^{1/2}. \qquad (42)$$

This diffusion-free scaling is commonly known as the Rhines scaling (Rhines 1975) and it is expected to hold in the limit of $Ek \ll 1$ when viscous effects become negligible in the bulk of the fluid (e.g., Gastine et al. 2016; Guervilly et al. 2019). The other equality which enters eq. (41) coupled with the additional assumption that $\alpha_T g U_c \Theta \sim \frac{v^3}{d^4} Ra(Nu - 1) Pr^{-2}$ (see Jones 2015) yields in its dimensionless form

$$Re_c \sim \left[ \frac{Ra}{Pr^2} (Nu - 1) \right]^{2/5} Ek^{1/5}. \qquad (43)$$
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Top panel: typical length-scale of the radial velocity, $L_{u_s}$, as a function of the Rossby number $Re_c E$. A best-fitting power law using only the runs with $Ro \le 10^{-3}$ is shown by the dotted-dashed line. Bottom panel: convective Reynolds number, $Re_c$, as a function of $Ra(Nu - 1)Pr^{-2}Ek^{1/2}$ (eq. 43). The Ekman and Prandtl numbers are indicated with colours and shapes, respectively, and the hybrid runs are indicated with the corresponding hatched symbols and the full 3-D runs with empty symbols. The theoretically predicted scaling corresponding to the CIA balance (eq. 43) is shown by the dotted-dashed line.</figcaption></figure>

This equation is known as the inertial scaling of the convective velocity for rotating convection (or the CIA scaling, e.g Gillet & Jones 2006; King & Buffett 2013).

Note that another equilibrium would hold if viscous effects would replace inertia in the vorticity balance (40). This equilibrium is sometimes referred to as the VAC scaling, $Re_c \sim [\frac{Ra}{Pr^2}(Nu - 1)]^{1/2}Ek^{1/3}$, where Viscous, Archimedean and Coriolis effects are the dominant terms (Aubert et al. 2001; King & Buffett 2013). We will not discuss this scaling since it does not provide a suitable interpretation of the numerical simulations in the turbulent QG regime (Gastine et al. 2016; Guervilly et al. 2019; Schwaiger et al. 2021), as is also found with our simulations.

We now analyse the relevance of the asymptotic scaling laws eqs (42-43) in the context of our ensemble of numerical simulations with fixed temperature contrast. Fig. 13 shows all our numerical simulations compared with the CIA scaling laws for convective velocity and length-scale. On the top panel, the typical non-dimensional length-scale of the convection $L_{u_s}$ is plotted against the Rossby number $Ro = Re_c Ek$, corresponding to (42), while the bottom figure shows $Re_c$ as a function of $Ra(Nu-1)Pr^{-2}Ek^{1/2}$, corresponding to (43). We observe that the Rhines scale captures well the behaviour observed in our simulations. The majority of the points are aligned (black dotted-dashed line), and departures from the theory are confined to the highest Ekman numbers ($Ek \ge 10^{-5}$). Introducing a local Rossby number $Ro_L = Ro d/L_{u_s}$ we find that all our runs have $Ro_L < 0.1$ (not shown) which indicates that our primary assumption holds at the local level, even if the geostrophic constraint can be weaker for the highest values of $Ro_L$ (associated with the highest Ek) of our data set. There is no additional Pr
dependence since all the simulations with $Pr \neq 1$ are close to the average scaling behaviour. The best-fitting for the whole data set yields a power law with an exponent of 0.31 for the $Ro$ dependency but considering only the runs with $Ro \le 10^{-3}$, we find a steeper slope with an exponent of $\alpha \sim 0.41$, a value comparable to that obtained in 3-D parameters studies (Gastine et al. 2016; Long et al. 2020) but shallower than the theoretical 1/2 scaling. Using QG models at Ekman numbers as low as $Ek = 10^{-11}$, Guervilly et al. (2019) showed that $l_{\perp}/d \sim Ro^{1/2}$ is gradually approached in the low viscosity limit appropriate for planetary cores.

The bottom panel of Fig. 13 shows that most of our simulations agree well with the CIA theoretical scaling law—with the power exponent 2/5 (black dotted line)—especially for $Ek \le 10^{-6}$. Cases at higher Ekman number ($Ek \ge 10^{-5}$; red and orange symbols) depart from the theory, following instead a power law with a lower exponent consistent with a deviation observed from the Rhines scaling at higher $Ek$. The simulations conducted at different Prandtl numbers are also well aligned with the CIA power law and are parallel to the $Pr = 1$ series but shifted upward for $Pr < 1$ and downward for $Pr > 1$. The clear separation between the series at different Prandtl numbers suggests that there is a dependency on $Pr$ in the $Re_c$ scaling which affects the prefactor of the scaling law and is not accounted for in the CIA scaling law (43). The hybrid runs (hatched symbols) are offset leftwards and downwards compared to the purely QG runs, reflecting that lower Nusselt numbers and velocities are reached for the same parameters. This shift can be understood in terms of the limitations discussed in Section 3.2.2 with an effective lack of buoyancy power when $Ra$ is increased, yielding lower velocities—that is $Re_c$—and lower heat transport effectiveness—that is $Nu$, hence weaker $Ra(Nu - 1)$—for the same control parameters ($Ek$, $Pr$, $Ra$). Note that the purely 3-D runs (empty symbols) are also shifted towards the left and stand between the hybrid and the QG cases, indicating that the QG runs are, in contrast, overpowered compared to the 3-D cases. This is likely due to the cylindrical boundary conditions in the QG case: the temperature imposed at the inner boundary is fixed for the whole column (and not only for the inner core surface) at all $Ra$, artificially supplying more thermal power to the bulk compared to the purely 3-D setup (see Fig. 7b). Parameter studies with either QG or 3-D models (e.g., Gillet & Jones 2006; Guervilly 2010; King & Buffett 2013) have reported exponents steeper than 2/5. This has been attributed by Gastine et al. (2016) to the sizeable role played by viscous dissipation in the boundary layers for $Re_c < 10^4$. Discrepancies arise at high $Ek$ and low $Ra/Re_c$ where the VAC balance is probably more suitable and at high $Ra/Re_c$ where the QG approximation no longer holds (as has been observed in, e.g. Gastine et al. 2016). The CIA scaling law hence partly captures the actual scaling behaviour of the convective velocity: fitting all the data with $Pr = 1$ yields $Re_c = 0.53 [Ra(Nu - 1)Ek^{1/2}]^{0.43}$ in reasonable agreement with the theory, although the $Pr$ dependence is not well accounted for.

Despite the well-known limitations of the QG approximation, and the limitations of our hybrid method at high $Ra/Re_c$ as documented above, we find that across our entire suite of calculations, results are broadly consistent with the Rhines and CIA scaling laws with a remaining dependence on $Pr$ for the latter. This lends additional support to findings of previous studies that also favoured a CIA balance but focused on a weaker forcing regime in a full sphere geometry (Guervilly et al. 2019), used a thinner spherical shell geometry and a different gravity profile (Gastine et al. 2016) or used heat flux boundary conditions (Long et al. 2020).

# 4 DISCUSSION AND CONCLUSIONS

We have used QG, 3-D and hybrid models, the latter involving a QG velocity field and a 3-D temperature field, to explore in a thick spherical shell the regime of strongly driven, rapidly rotating convection, focusing on low Ekman numbers $10^{-10} \le Ek \le 10^{-4}$, reaching supercriticalities up to $Ra \sim 160 Re_c$ and considering a range of Prandtl numbers close to and below unity $10^{-2} \le Pr \le 10$, also exploring the impact of laterally varying heat flux boundary conditions. This work involved extending the QG convection code pizza (Gastine 2019) to include the possibility to work with laterally varying heat flux boundary conditions, a 3-D temperature field, and a thermal wind.

Using the hybrid QG-3D approach at parameters $Ek \le 10^{-6}$ and $Pr \le 1$ we are able to reproduce important aspects of convective dynamics seen in 3-D models for weak to moderate supercriticalities ($Ra/Re_c \le 10-15$). In that regime, the meridional temperature profile and $\phi$-averaged azimuthal velocity are well retrieved although, as for purely QG models, the dynamical behaviour in the hybrid model deviates from the 3-D models close to the outer boundary. When $Ra$ is further increased, we find our hybrid model develops much weaker convective flows compared with the 3-D configuration in all cases. In contrast, the range of agreement between the hybrid QG-3D and the 3-D configurations increases when $Pr = 10^{-2}$ (as suggested by Guervilly et al. 2019) while the QG model departs from the 3-D around $Ra \sim 10 Re_c$. We expect that the hybrid model will eventually diverge from the full 3-D model when $Ra/Re_c$ is sufficiently high, but the exact value of the diverging point remains to be determined at $Pr \le 10^{-2}$ and has not been yet numerically reached despite reaching large $Re$ values. This is to some extent expected since 3-D effects become important when the thermal forcing increases (Calkins et al. 2013): the appearance of non-QG, equatorially antisymmetric, axial flows that do not vary linearly along the rotation axis—breaking the underlying classical QG assumption—and the associated missing correlations between $u_z$ and $T$, are part of this discrepancy. We found these effects account for up to 46 per cent of the missing convective power for a case at $Ek = 10^{-6}$, $Pr = 1$, $Ra = 53 Re_c$. The thick shell geometry studied and the omission of the dynamics inside the tangent cylinder, as well as the enforced linearity of $u_z$, may also play a role in the less vigorous convection found in the hybrid model. By construction in our hybrid setup the vertical motions remain weak compared to the horizontal motions. Alternative formulations of QG-type models have recently been proposed that aim to better represent all flow components (Gerick et al. 2020; Jackson & Maffei 2020); these could perhaps provide a means to improve on the results presented here.

In theory the hybrid method can be three to five times faster than a 3-D model when using the same resolution—because compute time for Legendre transforms associated with the velocity field are saved. However, we did not find major computational advantages in using the hybrid method at $Pr = 1$ mainly because the z-interpolation scheme requires the resolution to be high. At $Pr < 1$, the hybrid
model becomes more advantageous; since the 3-D grid involves only the temperature field it can be much coarser than the 2-D grid associated with the velocity field. The hybrid QG-3D approach studied here is nevertheless suitable for studying the rapidly rotating regime of convection at moderate forcing (i.e. $Ra/Ra_c \le 15$) at all $Pr$ and we envision that the hybrid method could become more attractive in terms of computational resources, even at $Pr = 1$, if the accuracy of the interpolation methods can be improved without sacrificing too much speed.

At the relatively low Ekman numbers considered here imposing a fixed heat flux condition at the outer boundary had little impact on the convective dynamics, compared to using a fixed temperature condition. However imposing lateral variations in the heat flux at the outer boundary can result in regions close to the boundary where the convection is inhibited or suppressed, and characterized by spiral arms where only azimuthal motions are possible, even at high supercriticalities, provided $Q^*$—the peak-to-peak relative amplitude of the flux perturbation—is sufficiently large. We also observe enhanced convective patterns in our QG models attached to the fluid regions with higher heat-flux. Such alternating regions of inhibited and enhanced convection are well known from previous studies in 3-D (e.g., Mound & Davies 2017) and are evident in our 3-D comparison calculations carried out in a similar regime. When comparing the 3-D, QG and hybrid models using inhomogeneous thermal boundary conditions, we find that the QG and 3-D configurations are qualitatively similar whereas the hybrid model again seems to lack buoyancy power. Despite the absence of small-scale convection found in the hybrid setup, the basic heat anomaly pattern at the outer boundary and the upwelling/downwelling system under the enhanced flux regions are similar to those found in the 3-D cases, suggesting that the dynamics is relatively well captured at mid-to-low latitudes. For $Q^*$ in the range from 2 to 5, and for the relatively high supercriticalities explored here the underlying convection deep in the shell is not greatly affected.

In general, we find that azimuthal shearing of axial vorticity dominates the convective dynamics in the bulk and leads to the formation of multiple zonal jets of alternating sign when $Ek \le 10^{-7}$, $Re \gg 1$ and $Ro \ll 1$, as reported in previous QG studies (e.g., Guervilly & Cardin 2017). When decreasing the ratio of diffusivities such that $Pr \le 10^{-1}$, we find the QG-dynamics is not fundamentally modified: multiple zonal jets still dominate in the bulk but a lower supercriticality is required for the same $Ro$ number, leading to the formation of fewer and wider zonal jets at $Pr < 1$.

Regarding scaling laws, our data set follows reasonably well the inertial scaling of rotating convection, which relies on a triple force balance between buoyancy, Coriolis force and inertia (e.g., Cardin & Olson 1994; Gillet & Jones 2006). The convective flow length-scale $l_\perp/d$ gradually approaches the asymptotic Rhines scaling $l_\perp \sim Ro^{1/2}$ at low $Ek$ (Rhines 1975), albeit with an exponent lower than 1/2 for the $Ek$ considered here. For the velocity scaling, we find that the simulations with $Pr = 1$ follow reasonably well the theoretical CIA scaling $Re_c \sim [Ra(Nu - 1)Pr^{-2}Ek^{1/2}]^{2/5}$ with a retrieved exponent equal to 0.43. We find a clear dependence of the velocity scaling behaviour on the Prandtl number, which is not well described by the classical inertial scaling of rotating convection. The 3-D runs stand between the results of the QG and the hybrid setups, suggesting that the QG configuration produces too much convective power. This is most likely a consequence of applying temperature boundary conditions in cylindrical geometry which involves the crude approximation of a fixed temperature on the whole tangent cylinder at $s_i$, $T_{2-D}(s_i) \sim 0.445$ for all $Ra$ and may partly compensate for the lack of ageostrophic components in the QG models. On the other hand, the hybrid configuration clearly lacks convective power at $Ra/Ra_c \ge 15$. Our results are overall consistent with other recent parameter studies in different geometries and studying different ranges of the control parameters (e.g., Gastine et al. 2016; Guervilly et al. 2019). The Prandtl number dependence seems to mainly affect the prefactors of the scaling laws, suggesting there is no fundamental change in the dynamics, at least for the parameter range explored here.

This study is a first step towards a more general hybrid QG-3D approach to Earth's core dynamics that will include the crucial effects of a 3-D magnetic field on the QG-convection.

## ACKNOWLEDGMENTS

We thank two anonymous reviewers for constructive comments that helped improve the manuscript. Nathanaël Schaeffer and Céline Guervilly are thanked for helpful discussions and for sharing their code for the z-averaging functions and for the z-integration of the thermal wind. We also thank Tobias Schwaiger for providing us with his non-magnetic data base of simulations. We additionally thank Jérémie Vidal for helping with the SINGE code and the convergence of the onset of convection for 3-D runs at $Ek \le 10^{-6}$. 3-D global numerical computations using the MagIC code were performed at S-CAPAD, IPGP and using HPC resources from the IDRIS Jean Zay CPU cluster (Grant 2021-A0070410095). QG and hybrid QG-3D numerical computations using the pizza code were performed at DTU Space, using the Humboldt and Larmor CPU clusters. OB and CCF were supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 772561).

## DATA AVAILABILITY

All our data are available upon reasonable request to the corresponding author and some key parameters from our whole data-set are already included in this paper. The codes used (pizza and MagIC) are both freely available (at http://www.github.com/magic-sph/pizza and http://www.github.com/magic-sph/magic, respectively) under the GNU GPL v3 license.
## REFERENCES

Ascher, U.M., Ruuth, S.J. & Wetton, B.T., 1995. Implicit-explicit methods for time-dependent partial differential equations, *SIAM J. Numer. Anal.*, **32**(3), 797-823.

Ascher, U.M., Ruuth, S.J. & Spiteri, R.J., 1997. Implicit-explicit runge-kutta methods for time-dependent partial differential equations, *Appl. Numer. Math.*, **25**(2-3), 151-167.

Aubert, J., 2015. Geomagnetic forecasts driven by thermal wind dynamics in the Earth's core, *Geophys. J. Int.*, **203**(3), 1738-1751.

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**(1-4), 51-74.

Aubert, J., Gillet, N. & Cardin, P., 2003. Quasigeostrophic models of convection in rotating spherical shells, *Geochem. Geophys. Geosyst.*, **4**(7), 1-19.

Aurnou, J., Calkins, M., Cheng, J., Julien, K., King, E., Nieves, D., Soderlund, K.M. & Stellmach, S., 2015. Rotating convective turbulence in earth and planetary cores, *Phys. Earth planet. Inter.*, **246**, 52-71.

Busse, F.H., 1970. Thermal instabilities in rapidly rotating systems, *J. Fluid Mech.*, **44**(3), 441-460.

Calkins, M.A., Aurnou, J.M., Eldredge, J.D. & 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.

Calkins, M.A., Hale, K., Julien, K., Nieves, D., Driggs, D. & Marti, P., 2015. The asymptotic equivalence of fixed heat flux and fixed temperature thermal boundary conditions for rapidly rotating convection, *J. Fluid Mech.*, **784**, R2, doi:10.1017/jfm.2015.606.

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**(3-4), 235-259.

Chandrasekhar, S., 1961. *Hydrodynamic and Hydrodynamic Stability*, Oxford Univ. Press.

Clarté, T.T., Schaeffer, N., Labrosse, S. & Vidal, J., 2021. The effects of a robin boundary condition on thermal convection in a rotating spherical shell, *J. Fluid Mech.*, **918**, doi:10.1017/jfm.2021.356.

Dormy, E., Soward, A., Jones, C., Jault, D. & Cardin, P., 2004. The onset of thermal convection in rotating spherical shells, *J. Fluid Mech.*, **501**, 43-70.

Gastine, T., 2019. pizza: an open-source pseudo-spectral code for spherical quasi-geostrophic convection, *Geophys. J. Int.*, **217**(3), 1558-1576.

Gastine, T., Wicht, J. & Aubert, J., 2016. Scaling regimes in spherical shell rotating convection, *J. Fluid Mech.*, **808**, 690-732.

Gerick, F., Jault, D., Noir, J. & Vidal, J., 2020. Pressure torque of torsional alfvén modes acting on an ellipsoidal mantle, *Geophys. J. Int.*, **222**(1), 338-351.

Gibbons, S., Gubbins, D. & Zhang, K., 2007. Convection in rotating spherical fluid shells with inhomogeneous heat flux at the outer boundary, *Geophys. Astrophys. Fluid Dyn.*, **101**(5-6), 347-370.

Gillet, N. & Jones, C., 2006. The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder, *J. Fluid Mech.*, **554**, 343-369.

Glatzmaier, G.A., 1984. Numerical simulations of stellar convective dynamos. i. the model and method, *J. Comput. Phys.*, **55**(3), 461-484.

Glatzmaier, G.A. & Roberts, P.H., 1996. An anelastic evolutionary geodynamo simulation driven by compositional and thermal convection, *Phys. D: Nonlin. Phenom.*, **97**(1), 81-94.

Goluskin, D., 2016. *Internally Heated Convection and Rayleigh-Bénard Convection*, Springer.

Gubbins, D., Dehant, V., Creager, K., Karato, S. & Zatman, S., 2003. Thermal core-mantle interactions: theory and observations, in *Earth's Core: Dynamics, Structure, Rotation*, pp. 163-179, eds Dehant, V., Creager, K. C., Karato, S.-i. & Zatman, S., John Wiley & Sons, Inc.

Guervilly, C., 2010. *Dynamos numériques planétaires générées par cisaillement en surface ou chauffage interne*, Université de Grenoble.

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., 2017. Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores, *Geophys. J. Int.*, **211**(1), 455-471.

Guervilly, C., Cardin, P. & Schaeffer, N., 2019. Turbulent convective length scale in planetary cores, *Nature*, **570**(7761), 368-371.

Hide, R., 1974. Jupiter and Saturn, *Proc. R. Soc. Lond.*, A, **336**(1604), 63-84.

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**(1), 62-80.

Jackson, A. & Maffei, S., 2020. Plesio-geostrophy for Earth's core: I. Basic equations, inertial modes and induction, *Proc. R. Soc. A*, **476**(2243), doi:10.1098/rspa.2020.0513.

Jones, C.A., 2015. 8.05 Thermal and compositional convection in the outer core, in *Treatise on Geophysics*, 2nd edn, pp. 115-159, ed. Schubert, G., Elsevier.

Julien, K., Knobloch, E., Rubio, A.M. & Vasil, G.M., 2012. Heat transport in low-Rossby-number Rayleigh-Bénard convection, *Phys. Rev. Lett.*, **109**(25), 254503.

King, E.M. & Aurnou, J.M., 2013. Turbulent convection in liquid metal with and without rotation, *Proc. Natl. Acad. Sci.*, **110**(17), 6688-6693.

King, E.M. & Buffett, B.A., 2013. Flow speeds and length scales in geodynamo models: the role of viscosity, *Earth planet. Sci. Lett.*, **371**, 156-162.

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.

Long, R., Mound, J., Davies, C. & Tobias, S., 2020. Scaling behaviour in spherical shell rotating convection with fixed-flux thermal boundary conditions, *J. Fluid Mech.*, **889**, doi:10.1017/jfm.2020.67.

Maffei, S., Jackson, A. & Livermore, P.W., 2017. Characterization of columnar inertial modes in rapidly rotating spheres and spheroids, *Proc. R. Soc. A*, **473**(2204), doi:10.1098/rspa.2017.0181.

Marti, P., Calkins, M. & Julien, K., 2016. A computationally efficient spectral method for modeling core dynamics, *Geochem. Geophys. Geosyst.*, **17**(8), 3031-3053.

Mound, J.E. & Davies, C.J., 2017. Heat transfer in rapidly rotating convection with heterogeneous thermal boundary conditions, *J. Fluid Mech.*, **828**, 601-629.

Muite, B., 2010. A numerical comparison of Chebyshev methods for solving fourth order semilinear initial boundary value problems, *J. Comput. Appl. Math.*, **234**(2), 317-342.

Rhines, P.B., 1975. Waves and turbulence on a beta-plane, *J. Fluid Mech.*, **69**(3), 417-443.

Roberts, P.H. & King, E.M., 2013. On the genesis of the Earth's magnetism, *Rep. Prog. Phys.*, **76**(9), 096801, doi:10.1088/0034-4885/76/9/096801.

Rossby, H., 1969. A study of Bénard convection with and without rotation, *J. Fluid Mech.*, **36**(2), 309-335.

Sahoo, S. & Sreenivasan, B., 2020. Convection in a rapidly rotating cylindrical annulus with laterally varying boundary heat flux, *J. Fluid Mech.*, **883**, doi:10.1017/jfm.2019.803.

Sakuraba, A. & Roberts, P.H., 2009. Generation of a strong magnetic field using uniform heat flux at the surface of the core, *Nat. Geosci.*, **2**(11), 802-805.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, *Geochem. Geophys. Geosyst.*, **14**(3), 751-758.

Schaeffer, N. & Cardin, P., 2005. Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers, *Phys. Fluids*, **17**(10), doi:10.1063/1.2073547.

Schaeffer, N. & Cardin, P., 2006. Quasi-geostrophic kinematic dynamics at low magnetic Prandtl number, *Earth planet. Sci. Lett.*, **245**(3-4), 595-604.

Schaeffer, N., Jault, D., Nataf, H.-C. & Fournier, A., 2017. Turbulent geodynamo simulations: a leap towards Earth's core, *Geophys. J. Int.*, **211**(1), 1-29.
Schwaiger, T., Gastine, T. & Aubert, J., 2021. Relating force balances and flow length scales in geodynamo simulations, *Geophys. J. Int.*, **224**(3), 1890–1904.

Sheyko, A., Finlay, C., Favre, J. & Jackson, A., 2018. Scale separated low viscosity dynamos and dissipation within the earth's core, *Nat. Sci. Rep.*, **8**(1), 1–7.

Stellmach, S. & Hansen, U., 2008. An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers, *Geochem. Geophys. Geosyst.*, **9**(5), doi:10.1029/2007GC001778.

Stellmach, S., Lischper, M., Julien, K., Vasil, G., Cheng, J.S., Ribeiro, A., King, E.M. & Aurnou, J.M., 2014. Approaching the asymptotic regime of rapidly rotating convection: boundary layers versus interior dynamics, *Phys. Rev. Lett.*, **113**(25), 254501, doi:10.1103/PhysRevLett.113.254501.

Valdettaro, L., Rieutord, M., Braconnier, T. & Frayssé, V., 2007. Convergence and round-off errors in a two-dimensional eigenvalue problem using spectral methods and Arnoldi-Chebyshev algorithm, *J. Comput. Appl. Math.*, **205**(1), 382–393.

Vidal, J. & Schaeffer, N., 2015. Quasi-geostrophic modes in the Earth's fluid core with an outer stably stratified layer, *Geophys. J. Int.*, **202**(3), 2182–2193.

Wicht, J., 2002. Inner-core conductivity in numerical dynamo simulations, *Phys. Earth planet. Inter.*, **132**(4), 281–302.

Yadav, R.K., Gastine, T., Christensen, U.R., Duarte, L.D.V. & 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.

## APPENDIX A: NUMERICAL SIMULATIONS WITH FIXED TEMPERATURE CONTRAST

<table><caption>**Table A1.** Summary of the numerical simulations with fixed temperature contrast computed in this study. All models have been computed with $\eta = r_i/r_o = 0.35$. Ra is the Rayleigh number (the supercriticality $Ra = \cdot \times Ra_c$ is also provided), Pr is the Prandtl number, Nu is the Nusselt number, $Re_c$ is the convective Reynolds number, $Re_{zon}$ is the zonal Reynolds number, $L_{us}$ is the typical length-scale for the cylindrical radial velocity field and $(N_s, N_m)/(N_r, \ell_{max})$ are the grid size for the run.</caption><thead><tr><th>Method</th><th>Ra</th><th>= × Ra<sub>c</sub></th><th>Pr</th><th>Nu</th><th>Re<sub>c</sub></th><th>Re<sub>zon</sub></th><th>L<sub>us</sub></th><th>(N<sub>s</sub>, N<sub>m</sub>)/(N<sub>r</sub>, &ell;<sub>max</sub>)</th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-weight:bold;">Ek = 1 &times; 10<sup>-4</sup></td></tr><tr><td>QG</td><td>4.83 &times; 10<sup>6</sup></td><td>=4.8 &times;</td><td>10</td><td>2.62</td><td>7.0</td><td>0.5</td><td>3.651 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>QG</td><td>1.12 &times; 10<sup>7</sup></td><td>=10.2 &times;</td><td>10</td><td>7.76</td><td>14.7</td><td>1.3</td><td>3.296 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>QG</td><td>3.87 &times; 10<sup>7</sup></td><td>=35.1 &times;</td><td>10</td><td>18.02</td><td>36.0</td><td>6.6</td><td>2.993 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>QG</td><td>2.43 &times; 10<sup>6</sup></td><td>=3.4 &times;</td><td>1</td><td>1.91</td><td>35.9</td><td>14.7</td><td>4.481 &times; 10<sup>-1</sup></td><td>(97, 96)/(-, -)</td></tr><tr><td>Hybrid</td><td>2.00 &times; 10<sup>6</sup></td><td>=6.6 &times;</td><td>1</td><td>1.46</td><td>54.4</td><td>27.2</td><td>4.479 &times; 10<sup>-1</sup></td><td>(97, 96)/(97, 96)</td></tr><tr><td>QG</td><td>4.83 &times; 10<sup>6</sup></td><td>=6.8 &times;</td><td>1</td><td>3.63</td><td>72.2</td><td>40.0</td><td>3.769 &times; 10<sup>-1</sup></td><td>(97, 96)/(-, -)</td></tr><tr><td>QG</td><td>1.93 &times; 10<sup>7</sup></td><td>=27.1 &times;</td><td>1</td><td>11.95</td><td>211.5</td><td>191.7</td><td>5.006 &times; 10<sup>-1</sup></td><td>(97, 96)/(-, -)</td></tr><tr><td>QG</td><td>3.87 &times; 10<sup>7</sup></td><td>=54.3 &times;</td><td>1</td><td>17.69</td><td>335.2</td><td>339.4</td><td>5.884 &times; 10<sup>-1</sup></td><td>(97, 96)/(-, -)</td></tr><tr><td>QG</td><td>4.30 &times; 10<sup>5</sup></td><td>=1.7 &times;</td><td>0.1</td><td>1.05</td><td>44.6</td><td>28.3</td><td>7.412 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>QG</td><td>9.66 &times; 10<sup>5</sup></td><td>=3.3 &times;</td><td>0.1</td><td>1.25</td><td>114.4</td><td>120.5</td><td>6.699 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>QG</td><td>4.83 &times; 10<sup>6</sup></td><td>=16.6 &times;</td><td>0.1</td><td>3.97</td><td>633.0</td><td>912.2</td><td>6.823 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">Ek = 1 &times; 10<sup>-5</sup></td></tr><tr><td>QG</td><td>1.12 &times; 10<sup>8</sup></td><td>=6.2 &times;</td><td>10</td><td>3.97</td><td>19.0</td><td>0.9</td><td>1.559 &times; 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>QG</td><td>2.25 &times; 10<sup>8</sup></td><td>=12.4 &times;</td><td>10</td><td>11.86</td><td>38.8</td><td>2.7</td><td>1.424 &times; 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>3-D</td><td>5.00 &times; 10<sup>8</sup></td><td>=30.8 &times;</td><td>10</td><td>15.03</td><td>97.1</td><td>13.1</td><td>-</td><td>(-, -)/(193, 170)</td></tr><tr><td>QG</td><td>1.12 &times; 10<sup>9</sup></td><td>=61.8 &times;</td><td>10</td><td>50.66</td><td>131.0</td><td>16.1</td><td>1.849 &times; 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>3-D</td><td>2.00 &times; 10<sup>9</sup></td><td>=123.4 &times;</td><td>10</td><td>36.63</td><td>319.8</td><td>53.0</td><td>-</td><td>(-, -)/(256, 213)</td></tr><tr><td>QG</td><td>2.25 &times; 10<sup>9</sup></td><td>=123.5 &times;</td><td>10</td><td>72.15</td><td>210.2</td><td>38.3</td><td>2.215 &times; 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>3-D</td><td>4.00 &times; 10<sup>7</sup></td><td>=3.8 &times;</td><td>1</td><td>1.39</td><td>62.9</td><td>21.9</td><td>-</td><td>(-, -)/(97, 170)</td></tr><tr><td>Hybrid</td><td>4.64 &times; 10<sup>7</sup></td><td>=4.7 &times;</td><td>1</td><td>1.37</td><td>78.8</td><td>27.5</td><td>2.172 &times; 10<sup>-1</sup></td><td>(193, 128)/(193, 128)</td></tr><tr><td>QG</td><td>5.21 &times; 10<sup>7</sup></td><td>=5.0 &times;</td><td>1</td><td>2.48</td><td>111.5</td><td>31.2</td><td>2.107 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>Hybrid</td><td>9.28 &times; 10<sup>7</sup></td><td>=9.4 &times;</td><td>1</td><td>1.60</td><td>154.9</td><td>57.4</td><td>2.385 &times; 10<sup>-1</sup></td><td>(193, 128)/(193, 128)</td></tr><tr><td>3-D</td><td>1.00 &times; 10<sup>8</sup></td><td>=9.5 &times;</td><td>1</td><td>2.45</td><td>194.5</td><td>53.8</td><td>-</td><td>(-, -)/(97, 213)</td></tr><tr><td>QG</td><td>1.04 &times; 10<sup>8</sup></td><td>=9.9 &times;</td><td>1</td><td>4.82</td><td>232.8</td><td>81.5</td><td>2.578 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>Hybrid</td><td>1.86 &times; 10<sup>8</sup></td><td>=18.9 &times;</td><td>1</td><td>1.82</td><td>272.9</td><td>131.2</td><td>2.899 &times; 10<sup>-1</sup></td><td>(193, 128)/(193, 128)</td></tr><tr><td>3-D</td><td>2.00 &times; 10<sup>8</sup></td><td>=19.0 &times;</td><td>1</td><td>5.42</td><td>471.3</td><td>161.0</td><td>-</td><td>(-, -)/(97, 256)</td></tr><tr><td>QG</td><td>2.09 &times; 10<sup>8</sup></td><td>=19.9 &times;</td><td>1</td><td>10.88</td><td>436.8</td><td>214.1</td><td>3.225 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>3-D</td><td>3.00 &times; 10<sup>8</sup></td><td>=28.6 &times;</td><td>1</td><td>8.99</td><td>741.1</td><td>338.3</td><td>-</td><td>(-, -)/(121, 288)</td></tr><tr><td>3-D</td><td>4.00 &times; 10<sup>8</sup></td><td>=38.1 &times;</td><td>1</td><td>11.99</td><td>959.3</td><td>499.5</td><td>-</td><td>(-, -)/(121, 288)</td></tr><tr><td>QG</td><td>4.18 &times; 10<sup>8</sup></td><td>=39.8 &times;</td><td>1</td><td>21.99</td><td>775.1</td><td>528.3</td><td>4.107 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>QG</td><td>5.21 &times; 10<sup>8</sup></td><td>=49.7 &times;</td><td>1</td><td>26.52</td><td>885.7</td><td>687.5</td><td>4.765 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>3-D</td><td>7.00 &times; 10<sup>8</sup></td><td>=66.7 &times;</td><td>1</td><td>18.09</td><td>1485.4</td><td>1031.8</td><td>-</td><td>(-, -)/(161, 426)</td></tr><tr><td>QG</td><td>8.36 &times; 10<sup>8</sup></td><td>=79.6 &times;</td><td>1</td><td>30.82</td><td>1216.9</td><td>1206.1</td><td>4.729 &times; 10<sup>-1</sup></td><td>(385, 384)/(-, -)</td></tr><tr><td>QG</td><td>1.08 &times; 10<sup>7</sup></td><td>=3.0 &times;</td><td>0.1</td><td>1.31</td><td>243.8</td><td>230.9</td><td>3.568 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>Hybrid</td><td>4.00 &times; 10<sup>7</sup></td><td>=11.3 &times;</td><td>0.1</td><td>1.41</td><td>789.1</td><td>1067.7</td><td>5.062 &times; 10<sup>-1</sup></td><td>(193, 240)/(129, 144)</td></tr><tr><td>QG</td><td>2.16 &times; 10<sup>7</sup></td><td>=6.0 &times;</td><td>0.1</td><td>1.99</td><td>554.4</td><td>665.8</td><td>4.432 &times; 10<sup>-1</sup></td><td>(193, 192)/(-, -)</td></tr><tr><td>QG</td><td>4.31 &times; 10<sup>7</sup></td><td>=11.9 &times;</td><td>0.1</td><td>3.47</td><td>1059.5</td><td>1493.5</td><td>5.969 &times; 10<sup>-1</sup></td><td>(385, 384)/(-, -)</td></tr><tr><td>QG</td><td>5.39 &times; 10<sup>7</sup></td><td>=14.9 &times;</td><td>0.1</td><td>4.27</td><td>1280.3</td><td>1876.1</td><td>6.310 &times; 10<sup>-1</sup></td><td>(385, 384)/(-, -)</td></tr><tr><td>QG</td><td>2.25 &times; 10<sup>8</sup></td><td>=62.0 &times;</td><td>0.1</td><td>15.52</td><td>3841.3</td><td>6703.9</td><td>6.159 &times; 10<sup>-1</sup></td><td>(769, 768)/(-, -)</td></tr><tr><td>3-D</td><td>2.50 &times; 10<sup>8</sup></td><td>=70.8 &times;</td><td>0.1</td><td>8.05</td><td>4426.5</td><td>5711.5</td><td>-</td><td>(-, -)/(145, 426)</td></tr></tbody></table>
<table><caption>Table A1. Continued</caption><thead><tr><th>Method</th><th>Ra</th><th>= × Ra<sub>c</sub></th><th>Pr</th><th>Nu</th><th>Re<sub>c</sub></th><th>Re<sub>zon</sub></th><th>L<sub>us</sub></th><th>(N<sub>s</sub>, N<sub>m</sub>)/(N<sub>r</sub>, l<sub>max</sub>)</th></tr></thead><tbody><tr><td>Hybrid</td><td>2.50 × 10<sup>8</sup></td><td>=70.8 ×</td><td>0.1</td><td>1.82</td><td>2584.0</td><td>4611.2</td><td>7.549 × 10<sup>-1</sup></td><td>(449, 480)/(321, 352)</td></tr><tr><td>QG</td><td>1.12 × 10<sup>7</sup></td><td>=5.7 ×</td><td>0.01</td><td>1.87</td><td>3074.0</td><td>5093.7</td><td>8.048 × 10<sup>-1</sup></td><td>(769, 768)/(−, −)</td></tr><tr><td>Hybrid</td><td>1.00 × 10<sup>7</sup></td><td>=5.1 ×</td><td>0.01</td><td>1.22</td><td>2614.7</td><td>3995.6</td><td>7.580 × 10<sup>-1</sup></td><td>(513, 512)/(321, 352)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Ek = 1 × 10<sup>-6</sup></td><td></td><td></td><td></td></tr><tr><td>QG</td><td>2.25 × 10<sup>10</sup></td><td>=67.6 ×</td><td>10</td><td>50.47</td><td>493.7</td><td>89.0</td><td>1.075 × 10<sup>-1</sup></td><td>(4097, 4096)/(−, −)</td></tr><tr><td>3-D</td><td>8.00 × 10<sup>8</sup></td><td>=4.3 ×</td><td>1</td><td>1.42</td><td>146.7</td><td>33.1</td><td>−</td><td>(−, −)/(129, 256)</td></tr><tr><td>3-D</td><td>1.00 × 10<sup>9</sup></td><td>=5.3 ×</td><td>1</td><td>1.55</td><td>195.9</td><td>41.7</td><td>−</td><td>(−, −)/(129, 256)</td></tr><tr><td>Hybrid</td><td>1.00 × 10<sup>9</sup></td><td>=6.1 ×</td><td>1</td><td>1.44</td><td>206.8</td><td>35.4</td><td>1.164 × 10<sup>-1</sup></td><td>(257, 256)/(257, 256)</td></tr><tr><td>QG</td><td>1.12 × 10<sup>9</sup></td><td>=6.3 ×</td><td>1</td><td>2.81</td><td>292.2</td><td>48.5</td><td>1.208 × 10<sup>-1</sup></td><td>(385, 384)/(−, −)</td></tr><tr><td>3-D</td><td>2.00 × 10<sup>9</sup></td><td>=10.6 ×</td><td>1</td><td>2.38</td><td>448.1</td><td>97.1</td><td>−</td><td>(−, −)/(129, 341)</td></tr><tr><td>Hybrid</td><td>2.00 × 10<sup>9</sup></td><td>=12.3 ×</td><td>1</td><td>1.67</td><td>404.4</td><td>86.9</td><td>1.456 × 10<sup>-1</sup></td><td>(257, 256)/(257, 256)</td></tr><tr><td>QG</td><td>2.00 × 10<sup>9</sup></td><td>=11.3 ×</td><td>1</td><td>5.05</td><td>550.4</td><td>126.6</td><td>1.547 × 10<sup>-1</sup></td><td>(385, 384)/(−, −)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>9</sup></td><td>=12.7 ×</td><td>1</td><td>5.75</td><td>620.7</td><td>153.6</td><td>1.598 × 10<sup>-1</sup></td><td>(385, 384)/(−, −)</td></tr><tr><td>3-D</td><td>2.80 × 10<sup>9</sup></td><td>=14.9 ×</td><td>1</td><td>3.44</td><td>700.6</td><td>176.9</td><td>−</td><td>(−, −)/(161, 426)</td></tr><tr><td>3-D</td><td>3.50 × 10<sup>9</sup></td><td>=18.6 ×</td><td>1</td><td>4.77</td><td>997.7</td><td>266.6</td><td>−</td><td>(−, −)/(181, 426)</td></tr><tr><td>Hybrid</td><td>4.00 × 10<sup>9</sup></td><td>=24.5 ×</td><td>1</td><td>2.01</td><td>740.7</td><td>256.0</td><td>1.879 × 10<sup>-1</sup></td><td>(257, 256)/(257, 256)</td></tr><tr><td>QG</td><td>4.49 × 10<sup>9</sup></td><td>=25.3 ×</td><td>1</td><td>12.79</td><td>1190.0</td><td>539.7</td><td>2.121 × 10<sup>-1</sup></td><td>(385, 384)/(−, −)</td></tr><tr><td>3-D</td><td>5.50 × 10<sup>9</sup></td><td>=29.3 ×</td><td>1</td><td>9.43</td><td>1803.9</td><td>1014.4</td><td>−</td><td>(−, −)/(201, 426)</td></tr><tr><td>Hybrid</td><td>7.80 × 10<sup>9</sup></td><td>=47.9 ×</td><td>1</td><td>2.14</td><td>1239.3</td><td>571.8</td><td>2.401 × 10<sup>-1</sup></td><td>(385, 384)/(385, 384)</td></tr><tr><td>QG</td><td>8.99 × 10<sup>9</sup></td><td>=50.7 ×</td><td>1</td><td>28.31</td><td>2145.5</td><td>1419.2</td><td>2.902 × 10<sup>-1</sup></td><td>(577, 576)/(−, −)</td></tr><tr><td>3-D</td><td>1.00 × 10<sup>10</sup></td><td>=53.2 ×</td><td>1</td><td>19.46</td><td>3339.9</td><td>2405.2</td><td>−</td><td>(−, −)/(321, 682)</td></tr><tr><td>Hybrid</td><td>1.00 × 10<sup>10</sup></td><td>=61.3 ×</td><td>1</td><td>2.20</td><td>1455.9</td><td>829.5</td><td>2.622 × 10<sup>-1</sup></td><td>(513, 512)/(513, 341)</td></tr><tr><td>QG</td><td>1.12 × 10<sup>10</sup></td><td>=63.4 ×</td><td>1</td><td>36.27</td><td>2584.4</td><td>1845.5</td><td>3.192 × 10<sup>-1</sup></td><td>(577, 576)/(−, −)</td></tr><tr><td>3-D</td><td>1.60 × 10<sup>10</sup></td><td>=85.1 ×</td><td>1</td><td>31.36</td><td>5272.2</td><td>4102.7</td><td>−</td><td>(−, −)/(433, 682)</td></tr><tr><td>QG</td><td>1.80 × 10<sup>10</sup></td><td>=101.4 ×</td><td>1</td><td>58.65</td><td>3793.5</td><td>3109.7</td><td>3.882 × 10<sup>-1</sup></td><td>(577, 576)/(−, −)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>10</sup></td><td>=126.7 ×</td><td>1</td><td>72.09</td><td>4520.3</td><td>3957.2</td><td>4.265 × 10<sup>-1</sup></td><td>(769, 768)/(−, −)</td></tr><tr><td>Hybrid</td><td>2.60 × 10<sup>10</sup></td><td>=159.5 ×</td><td>1</td><td>2.90</td><td>2963.4</td><td>2090.2</td><td>3.829 × 10<sup>-1</sup></td><td>(513, 512)/(513, 341)</td></tr><tr><td>3-D</td><td>2.66 × 10<sup>10</sup></td><td>=141.5 ×</td><td>1</td><td>45.94</td><td>7285.0</td><td>7945.6</td><td>−</td><td>(−, −)/(577, 853)</td></tr><tr><td>Hybrid</td><td>1.03 × 10<sup>8</sup></td><td>=2.0 ×</td><td>0.1</td><td>1.09</td><td>317.5</td><td>227.5</td><td>1.671 × 10<sup>-1</sup></td><td>(385, 384)/(257, 171)</td></tr><tr><td>QG</td><td>1.16 × 10<sup>8</sup></td><td>=2.2 ×</td><td>0.1</td><td>1.22</td><td>396.8</td><td>285.3</td><td>1.790 × 10<sup>-1</sup></td><td>(385, 384)/(−, −)</td></tr><tr><td>Hybrid</td><td>2.06 × 10<sup>8</sup></td><td>=4.1 ×</td><td>0.1</td><td>1.23</td><td>821.6</td><td>684.4</td><td>2.438 × 10<sup>-1</sup></td><td>(385, 384)/(257, 171)</td></tr><tr><td>QG</td><td>2.31 × 10<sup>8</sup></td><td>=4.4 ×</td><td>0.1</td><td>1.73</td><td>1049.2</td><td>855.3</td><td>2.564 × 10<sup>-1</sup></td><td>(385, 384)/(−, −)</td></tr><tr><td>QG</td><td>4.63 × 10<sup>8</sup></td><td>=8.8 ×</td><td>0.1</td><td>2.84</td><td>2161.3</td><td>2044.0</td><td>3.378 × 10<sup>-1</sup></td><td>(385, 384)/(−, −)</td></tr><tr><td>QG</td><td>9.26 × 10<sup>8</sup></td><td>=17.7 ×</td><td>0.1</td><td>5.43</td><td>3990.9</td><td>4472.9</td><td>4.296 × 10<sup>-1</sup></td><td>(577, 576)/(−, −)</td></tr><tr><td>3-D</td><td>1.00 × 10<sup>9</sup></td><td>=19.7 ×</td><td>0.1</td><td>2.77</td><td>3790.7</td><td>3458.6</td><td>−</td><td>(−, −)/(257, 309)</td></tr><tr><td>Hybrid</td><td>1.00 × 10<sup>9</sup></td><td>=19.7 ×</td><td>0.1</td><td>1.68</td><td>3196.6</td><td>3663.4</td><td>3.678 × 10<sup>-1</sup></td><td>(385, 384)/(257, 171)</td></tr><tr><td>QG</td><td>1.15 × 10<sup>9</sup></td><td>=22.1 ×</td><td>0.1</td><td>6.87</td><td>4768.8</td><td>5706.2</td><td>4.569 × 10<sup>-1</sup></td><td>(577, 576)/(−, −)</td></tr><tr><td>QG</td><td>1.85 × 10<sup>9</sup></td><td>=35.3 ×</td><td>0.1</td><td>11.72</td><td>6791.4</td><td>9383.2</td><td>4.970 × 10<sup>-1</sup></td><td>(769, 768)/(−, −)</td></tr><tr><td>3-D</td><td>2.00 × 10<sup>9</sup></td><td>=39.4 ×</td><td>0.1</td><td>5.95</td><td>7575.3</td><td>8782.8</td><td>−</td><td>(−, −)/(321, 426)</td></tr><tr><td>Hybrid</td><td>2.00 × 10<sup>9</sup></td><td>=39.4 ×</td><td>0.1</td><td>1.89</td><td>4975.7</td><td>6699.7</td><td>4.566 × 10<sup>-1</sup></td><td>(385, 384)/(257, 256)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>9</sup></td><td>=42.9 ×</td><td>0.1</td><td>13.95</td><td>7845.2</td><td>11884.2</td><td>4.495 × 10<sup>-1</sup></td><td>(769, 768)/(−, −)</td></tr><tr><td>3-D</td><td>5.00 × 10<sup>9</sup></td><td>=98.5 ×</td><td>0.1</td><td>12.77</td><td>12743.8</td><td>22752.8</td><td>−</td><td>(−, −)/(385, 1024)</td></tr><tr><td>Hybrid</td><td>5.00 × 10<sup>9</sup></td><td>=98.5 ×</td><td>0.1</td><td>2.14</td><td>8225.4</td><td>13536.3</td><td>5.201 × 10<sup>-1</sup></td><td>(769, 768)/(385, 384)</td></tr><tr><td>QG</td><td>4.85 × 10<sup>7</sup></td><td>=1.8 ×</td><td>0.01</td><td>1.20</td><td>2274.7</td><td>2943.8</td><td>4.144 × 10<sup>-1</sup></td><td>(577, 576)/(−, −)</td></tr><tr><td>QG</td><td>9.71 × 10<sup>7</sup></td><td>=3.6 ×</td><td>0.01</td><td>1.65</td><td>4930.7</td><td>7871.4</td><td>4.566 × 10<sup>-1</sup></td><td>(577, 576)/(−, −)</td></tr><tr><td>Hybrid</td><td>1.00 × 10<sup>8</sup></td><td>=3.9 ×</td><td>0.01</td><td>1.23</td><td>4749.7</td><td>7399.9</td><td>4.566 × 10<sup>-1</sup></td><td>(353, 384)/(225, 256)</td></tr><tr><td>Hybrid</td><td>1.30 × 10<sup>8</sup></td><td>=5.0 ×</td><td>0.01</td><td>1.29</td><td>6020.8</td><td>9797.3</td><td>4.908 × 10<sup>-1</sup></td><td>(513, 512)/(225, 256)</td></tr><tr><td>QG</td><td>1.94 × 10<sup>8</sup></td><td>=7.1 ×</td><td>0.01</td><td>2.54</td><td>8888.8</td><td>17077.8</td><td>4.993 × 10<sup>-1</sup></td><td>(1025, 1024)/(−, −)</td></tr><tr><td>3-D</td><td>2.00 × 10<sup>8</sup></td><td>=7.7 ×</td><td>0.01</td><td>1.45</td><td>6538.2</td><td>11456.1</td><td>−</td><td>(−, −)/(321, 682)</td></tr><tr><td>Hybrid</td><td>2.00 × 10<sup>8</sup></td><td>=7.7 ×</td><td>0.01</td><td>1.38</td><td>8233.1</td><td>14835.2</td><td>4.583 × 10<sup>-1</sup></td><td>(961, 1008)/(241, 288)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>8</sup></td><td>=8.9 ×</td><td>0.01</td><td>2.97</td><td>9761.0</td><td>19614.7</td><td>4.339 × 10<sup>-1</sup></td><td>(1025, 1024)/(−, −)</td></tr><tr><td>QG</td><td>2.43 × 10<sup>8</sup></td><td>=8.2 ×</td><td>0.01</td><td>3.00</td><td>10636.1</td><td>21312.7</td><td>5.072 × 10<sup>-1</sup></td><td>(1537, 1536)/(−, −)</td></tr><tr><td>Hybrid</td><td>4.00 × 10<sup>8</sup></td><td>=15.5 ×</td><td>0.01</td><td>1.52</td><td>13556.1</td><td>24142.9</td><td>6.494 × 10<sup>-1</sup></td><td>(1249, 1344)/(289, 384)</td></tr><tr><td>3-D</td><td>5.00 × 10<sup>8</sup></td><td>=19.3 ×</td><td>0.01</td><td>2.07</td><td>11968.4</td><td>28074.3</td><td>−</td><td>(−, −)/(321, 1024)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Ek = 1 × 10<sup>-7</sup></td><td></td><td></td><td></td></tr><tr><td>QG</td><td>2.43 × 10<sup>10</sup></td><td>=7.4 ×</td><td>1</td><td>2.97</td><td>714.3</td><td>115.3</td><td>6.781 × 10<sup>-2</sup></td><td>(577, 576)/(−, −)</td></tr><tr><td>Hybrid</td><td>5.00 × 10<sup>10</sup></td><td>=15.1 ×</td><td>1</td><td>1.73</td><td>1038.9</td><td>580.6</td><td>8.735 × 10<sup>-2</sup></td><td>(577, 672)/(577, 384)</td></tr><tr><td>QG</td><td>4.83 × 10<sup>10</sup></td><td>=14.8 ×</td><td>1</td><td>5.74</td><td>1459.0</td><td>848.3</td><td>1.148 × 10<sup>-1</sup></td><td>(769, 768)/(−, −)</td></tr><tr><td>Hybrid</td><td>6.45 × 10<sup>10</sup></td><td>=19.4 ×</td><td>1</td><td>1.81</td><td>1315.6</td><td>887.7</td><td>9.837 × 10<sup>-2</sup></td><td>(769, 768)/(769, 384)</td></tr><tr><td>3-D</td><td>6.50 × 10<sup>10</sup></td><td>=19.6 ×</td><td>1</td><td>3.85</td><td>1853.8</td><td>910.8</td><td>−</td><td>(−, −)/(433, 682)</td></tr><tr><td>QG</td><td>6.50 × 10<sup>10</sup></td><td>=19.9 ×</td><td>1</td><td>8.55</td><td>1879.9</td><td>1362.0</td><td>1.073 × 10<sup>-1</sup></td><td>(769, 768)/(−, −)</td></tr><tr><td>Hybrid</td><td>8.00 × 10<sup>10</sup></td><td>=24.1 ×</td><td>1</td><td>1.82</td><td>1501.6</td><td>1048.2</td><td>1.1068 × 10<sup>-1</sup></td><td>(769, 768)/(769, 512)</td></tr><tr><td>QG</td><td>9.66 × 10<sup>10</sup></td><td>=29.6 ×</td><td>1</td><td>13.91</td><td>2685.4</td><td>2343.9</td><td>1.290 × 10<sup>-1</sup></td><td>(769, 768)/(−, −)</td></tr></tbody></table>
<table><caption>Table A1. Continued</caption><thead><tr><th>Method</th><th>Ra</th><th>= × Ra<sub>c</sub></th><th>Pr</th><th>Nu</th><th>Re<sub>c</sub></th><th>Re<sub>zon</sub></th><th>L<sub>us</sub></th><th>(N<sub>s</sub>, N<sub>m</sub>)/(N<sub>r</sub>, l<sub>max</sub>)</th></tr></thead><tbody><tr><td>3-D</td><td>1.00 × 10<sup>11</sup></td><td>=30.1 ×</td><td>1</td><td>7.79</td><td>3325.4</td><td>2798.3</td><td>—</td><td>(—, —)/(513, 896)</td></tr><tr><td>Hybrid</td><td>1.00 × 10<sup>11</sup></td><td>=30.1 ×</td><td>1</td><td>1.94</td><td>1821.0</td><td>1382.9</td><td>1.112 × 10<sup>-1</sup></td><td>(769, 512)/(769, 512)</td></tr><tr><td>QG</td><td>1.93 × 10<sup>11</sup></td><td>=59.3 ×</td><td>1</td><td>31.71</td><td>4833.0</td><td>5217.6</td><td>1.711 × 10<sup>-1</sup></td><td>(769, 768)/(—, —)</td></tr><tr><td>QG</td><td>2.43 × 10<sup>11</sup></td><td>=74.5 ×</td><td>1</td><td>40.75</td><td>5865.3</td><td>6659.6</td><td>1.857 × 10<sup>-1</sup></td><td>(1153, 1152)/(—, —)</td></tr><tr><td>QG</td><td>3.87 × 10<sup>11</sup></td><td>=118.6 ×</td><td>1</td><td>69.78</td><td>8834.0</td><td>10483.3</td><td>2.202 × 10<sup>-1</sup></td><td>(1537, 1536)/(—, —)</td></tr><tr><td>QG</td><td>4.83 × 10<sup>11</sup></td><td>=148.2 ×</td><td>1</td><td>92.81</td><td>10796.3</td><td>12937.1</td><td>2.390 × 10<sup>-1</sup></td><td>(2049, 2048)/(—, —)</td></tr><tr><td>Hybrid</td><td>2.20 × 10<sup>9</sup></td><td>=2.0 ×</td><td>0.1</td><td>1.16</td><td>1182.5</td><td>703.5</td><td>1.208 × 10<sup>-1</sup></td><td>(769, 768)/(513, 256)</td></tr><tr><td>Hybrid</td><td>2.70 × 10<sup>9</sup></td><td>=2.5 ×</td><td>0.1</td><td>1.19</td><td>1517.0</td><td>925.8</td><td>1.259 × 10<sup>-1</sup></td><td>(769, 768)/(513, 256)</td></tr><tr><td>QG</td><td>4.99 × 10<sup>9</sup></td><td>=6.0 ×</td><td>0.1</td><td>1.82</td><td>3005.5</td><td>1994.4</td><td>1.434 × 10<sup>-1</sup></td><td>(769, 768)/(—, —)</td></tr><tr><td>QG</td><td>1.12 × 10<sup>10</sup></td><td>=13.6 ×</td><td>0.1</td><td>3.23</td><td>6091.6</td><td>5945.9</td><td>1.704 × 10<sup>-1</sup></td><td>(769, 768)/(—, —)</td></tr><tr><td>QG</td><td>2.00 × 10<sup>10</sup></td><td>=24.1 ×</td><td>0.1</td><td>5.67</td><td>9833.7</td><td>11422.3</td><td>2.066 × 10<sup>-1</sup></td><td>(769, 768)/(—, —)</td></tr><tr><td>QG</td><td>2.49 × 10<sup>10</sup></td><td>=30.1 ×</td><td>0.1</td><td>7.28</td><td>11883.6</td><td>14389.5</td><td>2.254 × 10<sup>-1</sup></td><td>(1025, 1024)/(—, —)</td></tr><tr><td>QG</td><td>4.00 × 10<sup>10</sup></td><td>=48.3 ×</td><td>0.1</td><td>12.71</td><td>17857.6</td><td>23170.1</td><td>2.555 × 10<sup>-1</sup></td><td>(1537, 1536)/(—, —)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>9</sup></td><td>=5.8 ×</td><td>0.01</td><td>2.54</td><td>20527.5</td><td>30576.4</td><td>2.840 × 10<sup>-1</sup></td><td>(1537, 1536)/(—, —)</td></tr><tr><td>QG</td><td>5.62 × 10<sup>9</sup></td><td>=14.5 ×</td><td>0.01</td><td>4.73</td><td>43765.4</td><td>65138.1</td><td>4.050 × 10<sup>-1</sup></td><td>(4609, 4608)/(—, —)</td></tr><tr><td>Hybrid</td><td>1.07 × 10<sup>11</sup></td><td>=6.8 ×</td><td>1</td><td>1.44</td><td>748.5</td><td>259.5</td><td>4.613 × 10<sup>-2</sup></td><td>(769, 768)/(769, 256)</td></tr><tr><td>QG</td><td>1.20 × 10<sup>11</sup></td><td>=7.9 ×</td><td>1</td><td>2.98</td><td>1088.2</td><td>364.3</td><td>5.043 × 10<sup>-2</sup></td><td>(769, 768)/(—, —)</td></tr><tr><td>QG</td><td>2.40 × 10<sup>11</sup></td><td>=15.7 ×</td><td>1</td><td>6.05</td><td>2144.8</td><td>1667.4</td><td>7.173 × 10<sup>-2</sup></td><td>(1025, 1024)/(—, —)</td></tr><tr><td>QG</td><td>4.81 × 10<sup>11</sup></td><td>=31.5 ×</td><td>1</td><td>13.88</td><td>4052.7</td><td>4268.4</td><td>9.658 × 10<sup>-2</sup></td><td>(1025, 1024)/(—, —)</td></tr><tr><td>QG</td><td>9.62 × 10<sup>11</sup></td><td>=62.9 ×</td><td>1</td><td>33.00</td><td>7528.6</td><td>9158.5</td><td>1.249 × 10<sup>-1</sup></td><td>(1025, 1024)/(—, —)</td></tr><tr><td>QG</td><td>1.20 × 10<sup>12</sup></td><td>=78.6 ×</td><td>1</td><td>42.72</td><td>9179.8</td><td>11560.0</td><td>1.361 × 10<sup>-1</sup></td><td>(1025, 1024)/(—, —)</td></tr><tr><td>QG</td><td>1.92 × 10<sup>12</sup></td><td>=125.8 ×</td><td>1</td><td>71.35</td><td>13560.0</td><td>18168.4</td><td>1.589 × 10<sup>-1</sup></td><td>(1537, 1536)/(—, —)</td></tr><tr><td>QG</td><td>2.40 × 10<sup>12</sup></td><td>=157.3 ×</td><td>1</td><td>90.78</td><td>18102.3</td><td>24309.7</td><td>1.836 × 10<sup>-1</sup></td><td>(2305, 2304)/(—, —)</td></tr><tr><td>QG</td><td>5.21 × 10<sup>11</sup></td><td>=8.2 ×</td><td>1</td><td>2.92</td><td>1523.8</td><td>889.8</td><td>3.807 × 10<sup>-2</sup></td><td>(1025, 1024)/(—, —)</td></tr><tr><td>QG</td><td>1.04 × 10<sup>12</sup></td><td>=16.5 ×</td><td>1</td><td>5.95</td><td>3073.2</td><td>2900.9</td><td>5.469 × 10<sup>-2</sup></td><td>(1537, 1536)/(—, —)</td></tr><tr><td>QG</td><td>2.09 × 10<sup>12</sup></td><td>=32.9 ×</td><td>1</td><td>13.99</td><td>5994.2</td><td>7201.6</td><td>7.423 × 10<sup>-2</sup></td><td>(1537, 1536)/(—, —)</td></tr><tr><td>QG</td><td>4.18 × 10<sup>12</sup></td><td>=66.0 ×</td><td>1</td><td>32.65</td><td>11052.2</td><td>15298.2</td><td>9.662 × 10<sup>-2</sup></td><td>(2049, 2048)/(—, —)</td></tr><tr><td>QG</td><td>4.49 × 10<sup>12</sup></td><td>=71.0 ×</td><td>1</td><td>35.20</td><td>12059.8</td><td>16988.9</td><td>9.931 × 10<sup>-2</sup></td><td>(3073, 3072)/(—, —)</td></tr><tr><td>QG</td><td>8.99 × 10<sup>12</sup></td><td>=142.0 ×</td><td>1</td><td>77.22</td><td>22869.3</td><td>34026.5</td><td>1.308 × 10<sup>-1</sup></td><td>(3073, 3072)/(—, —)</td></tr><tr><td>QG</td><td>5.62 × 10<sup>10</sup></td><td>=4.1 ×</td><td>0.1</td><td>1.36</td><td>3747.8</td><td>1887.6</td><td>6.280 × 10<sup>-2</sup></td><td>(1537, 1536)/(—, —)</td></tr><tr><td>Hybrid</td><td>1.00 × 10<sup>11</sup></td><td>=4.2 ×</td><td>0.1</td><td>1.58</td><td>4358.6</td><td>2942.0</td><td>8.157 × 10<sup>-2</sup></td><td>(1537, 1536)/(513, 171)</td></tr><tr><td>QG</td><td>1.12 × 10<sup>11</sup></td><td>=8.1 ×</td><td>0.1</td><td>1.83</td><td>7519.1</td><td>5751.9</td><td>8.310 × 10<sup>-2</sup></td><td>(1537, 1536)/(—, —)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>11</sup></td><td>=16.2 ×</td><td>0.1</td><td>2.93</td><td>13375.8</td><td>14941.6</td><td>1.030 × 10<sup>-1</sup></td><td>(2049, 2048)/(—, —)</td></tr><tr><td>QG</td><td>1.12 × 10<sup>10</sup></td><td>=2.0 ×</td><td>0.01</td><td>1.33</td><td>15817.8</td><td>12622.1</td><td>1.220 × 10<sup>-1</sup></td><td>(2049, 1536)/(—, —)</td></tr><tr><td>QG</td><td>2.00 × 10<sup>10</sup></td><td>=3.6 ×</td><td>0.01</td><td>1.60</td><td>27367.1</td><td>24838.9</td><td>1.399 × 10<sup>-1</sup></td><td>(2049, 2048)/(—, —)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>10</sup></td><td>=4.1 ×</td><td>0.01</td><td>1.96</td><td>28161.1</td><td>24492.0</td><td>1.448 × 10<sup>-1</sup></td><td>(3073, 2048)/(—, —)</td></tr><tr><td>QG</td><td>3.37 × 10<sup>10</sup></td><td>=6.1 ×</td><td>0.01</td><td>2.05</td><td>40509.8</td><td>43890.5</td><td>1.470 × 10<sup>-1</sup></td><td>(3457, 3456)/(—, —)</td></tr><tr><td>QG</td><td>3.93 × 10<sup>10</sup></td><td>=7.1 ×</td><td>0.01</td><td>2.24</td><td>44887.0</td><td>54466.6</td><td>1.319 × 10<sup>-1</sup></td><td>(3457, 3456)/(—, —)</td></tr><tr><td>QG</td><td>4.49 × 10<sup>10</sup></td><td>=8.1 ×</td><td>0.01</td><td>2.47</td><td>49801.2</td><td>61299.8</td><td>1.647 × 10<sup>-1</sup></td><td>(4097, 4096)/(—, —)</td></tr><tr><td>QG</td><td>1.12 × 10<sup>11</sup></td><td>=20.3 ×</td><td>0.01</td><td>4.43</td><td>120663.5</td><td>206311.7</td><td>1.496 × 10<sup>-1</sup></td><td>(9217, 9216)/(—, —)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>13</sup></td><td>=17.6 ×</td><td>1</td><td>5.99</td><td>6366.9</td><td>7671.7</td><td>3.002 × 10<sup>-2</sup></td><td>(3073, 2048)/(—, —)</td></tr><tr><td>QG</td><td>4.49 × 10<sup>13</sup></td><td>=35.2 ×</td><td>1</td><td>11.75</td><td>16231.2</td><td>24728.6</td><td>5.514 × 10<sup>-2</sup></td><td>(3073, 3072)/(—, —)</td></tr><tr><td>QG</td><td>2.25 × 10<sup>12</sup></td><td>=9.4 ×</td><td>0.1</td><td>2.20</td><td>16253.8</td><td>7816.3</td><td>4.624 × 10<sup>-2</sup></td><td>(4609, 4096)/(—, —)</td></tr><tr><td>QG</td><td>4.83 × 10<sup>14</sup></td><td>=17.6 ×</td><td>1</td><td>5.99</td><td>15563.3</td><td>17876.4</td><td>1.964 × 10<sup>-2</sup></td><td>(6145, 6144)/(—, —)</td></tr></tbody></table>

APPENDIX B: NUMERICAL SIMULATIONS WITH INHOMOGENEOUS HEAT FLUX BOUNDARY CONDITIONS
<table><caption>**Table B1.** Summary of the numerical simulations with inhomogeneous heat flux boundary conditions computed in this study. All models have been computed with $\eta = r_i/r_o = 0.35$. $Ra_Q$ is the flux-based Rayleigh number, $Y_{(m)}(\ell, m)$ is the mode m (QG) or $(\ell, m)$ (hybrid or 3-D) of the imposed lateral flux variations ($Y_0$ or $Y_{0,0}$ indicates no lateral variations), $Q^*$ is the relative amplitude of the lateral flux variations, $Pr$ is the Prandtl number, $Nu_\Delta$ is the Nusselt number based on the temperature contrast in the shell, $Re_c$ is the convective Reynolds number, $Re_{zon}$ is the zonal Reynolds number, $L_{us}$ is the typical length-scale for the cylindrical radial velocity field and $(N_s, N_m)/(N_r, \ell_{max})$ are the grid size for the run.</caption><thead><tr><th>Method</th><th>Ra<sub>Q</sub></th><th>Y<sub>(m)</sub>(ℓ, m)</th><th>Q<sup>*</sup></th><th>Pr</th><th>Nu<sub>Δ</sub></th><th>Re<sub>c</sub></th><th>Re<sub>zon</sub></th><th>L<sub>us</sub></th><th>(N<sub>s</sub>, N<sub>m</sub>)/(N<sub>r</sub>, ℓ<sub>max</sub>)</th></tr></thead><tbody><tr><td rowspan="2">Hybrid</td><td rowspan="2">1.00 × 10<sup>7</sup></td><td rowspan="2">Y<sub>2,2</sub></td><td rowspan="2">3</td><td rowspan="2">1</td><td>Ek = 1 × 10<sup>-4</sup></td><td rowspan="2">135.1</td><td rowspan="2">68.3</td><td rowspan="2">3.142 × 10<sup>0</sup></td><td rowspan="2">(65, 96)/(65, 96)</td></tr><tr><td>Ek = 1 × 10<sup>-5</sup></td></tr><tr><td>Hybrid</td><td>2.00 × 10<sup>8</sup></td><td>Y<sub>2,2</sub></td><td>3</td><td>1</td><td>1.67</td><td>425.4</td><td>87.7</td><td>4.488 × 10<sup>-1</sup></td><td>(129, 128)/(129, 128)</td></tr><tr><td>Hybrid</td><td>4.00 × 10<sup>8</sup></td><td>Y<sub>2,2</sub></td><td>3</td><td>1</td><td>1.71</td><td>543.1</td><td>110.0</td><td>6.283 × 10<sup>-1</sup></td><td>(129, 128)/(129, 128)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Ek = 1 × 10<sup>-6</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>QG</td><td>2.00 × 10<sup>9</sup></td><td>Y<sub>0</sub></td><td>—</td><td>1</td><td>4.20</td><td>460.2</td><td>119.5</td><td>2.417 × 10<sup>-1</sup></td><td>(385, 384)/(-, -)</td></tr><tr><td>Hybrid</td><td>2.00 × 10<sup>9</sup></td><td>Y<sub>2,2</sub></td><td>3</td><td>1</td><td>2.11</td><td>484.6</td><td>195.7</td><td>1.963 × 10<sup>-1</sup></td><td>(289, 320)/(289, 320)</td></tr><tr><td>QG</td><td>2.00 × 10<sup>9</sup></td><td>Y<sub>2</sub></td><td>3</td><td>1</td><td>30.32</td><td>144.4</td><td>192.5</td><td>8.491 × 10<sup>-2</sup></td><td>(385, 384)/(-, -)</td></tr><tr><td>QG</td><td>4.00 × 10<sup>9</sup></td><td>Y<sub>0</sub></td><td>—</td><td>1</td><td>4.25</td><td>927.8</td><td>300.8</td><td>1.745 × 10<sup>-1</sup></td><td>(385, 384)/(-, -)</td></tr><tr><td>QG</td><td>4.00 × 10<sup>9</sup></td><td>Y<sub>1</sub></td><td>3</td><td>1</td><td>13.16</td><td>611.9</td><td>135.4</td><td>1.366 × 10<sup>-1</sup></td><td>(513, 512)/(-, -)</td></tr><tr><td>3-D</td><td>4.00 × 10<sup>9</sup></td><td>Y<sub>2,2</sub></td><td>3</td><td>1</td><td>4.33</td><td>735.6</td><td>256.8</td><td>—</td><td>(-, -)/(257, 341)</td></tr><tr><td>Hybrid</td><td>4.00 × 10<sup>9</sup></td><td>Y<sub>2,2</sub></td><td>3</td><td>1</td><td>1.59</td><td>773.5</td><td>429.4</td><td>1.496 × 10<sup>-1</sup></td><td>(337, 384)/(337, 384)</td></tr><tr><td>QG</td><td>4.00 × 10<sup>9</sup></td><td>Y<sub>2</sub></td><td>3</td><td>1</td><td>8.86</td><td>640.3</td><td>120.2</td><td>1.745 × 10<sup>-1</sup></td><td>(513, 512)/(-, -)</td></tr><tr><td>QG</td><td>8.00 × 10<sup>9</sup></td><td>Y<sub>0</sub></td><td>—</td><td>1</td><td>5.88</td><td>1257.9</td><td>656.8</td><td>2.417 × 10<sup>-1</sup></td><td>(513, 512)/(-, -)</td></tr><tr><td>3-D</td><td>8.00 × 10<sup>9</sup></td><td>Y<sub>2,2</sub></td><td>3</td><td>1</td><td>7.15</td><td>974.7</td><td>352.4</td><td>—</td><td>(-, -)/(321, 341)</td></tr><tr><td>Hybrid</td><td>8.00 × 10<sup>9</sup></td><td>Y<sub>2,2</sub></td><td>3</td><td>1</td><td>2.03</td><td>1414.7</td><td>489.6</td><td>2.618 × 10<sup>-1</sup></td><td>(385, 416)/(385, 416)</td></tr><tr><td>QG</td><td>8.00 × 10<sup>9</sup></td><td>Y<sub>2</sub></td><td>3</td><td>1</td><td>13.92</td><td>893.8</td><td>238.4</td><td>1.496 × 10<sup>-1</sup></td><td>(513, 512)/(-, -)</td></tr><tr><td>Hybrid</td><td>1.26 × 10<sup>10</sup></td><td>Y<sub>2,2</sub></td><td>3</td><td>1</td><td>1.90</td><td>1306.3</td><td>1162.2</td><td>2.417 × 10<sup>-1</sup></td><td>(417, 448)/(417, 448)</td></tr><tr><td>QG</td><td>2.60 × 10<sup>10</sup></td><td>Y<sub>2</sub></td><td>3</td><td>1</td><td>21.59</td><td>1456.0</td><td>835.3</td><td>1.745 × 10<sup>-1</sup></td><td>(513, 512)/(-, -)</td></tr><tr><td>QG</td><td>3.60 × 10<sup>10</sup></td><td>Y<sub>0</sub></td><td>—</td><td>1</td><td>10.26</td><td>231.0</td><td>186.2</td><td>1.653 × 10<sup>-1</sup></td><td>(513, 512)/(-, -)</td></tr><tr><td>QG</td><td>3.60 × 10<sup>10</sup></td><td>Y<sub>2</sub></td><td>3</td><td>1</td><td>22.61</td><td>1914.1</td><td>957.0</td><td>3.142 × 10<sup>-1</sup></td><td>(513, 512)/(-, -)</td></tr><tr><td>QG</td><td>3.60 × 10<sup>10</sup></td><td>Y<sub>2</sub></td><td>5</td><td>1</td><td>59.32</td><td>720.5</td><td>185.5</td><td>1.963 × 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>Ek = 1 × 10<sup>-7</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>QG</td><td>1.62 × 10<sup>11</sup></td><td>Y<sub>0</sub></td><td>—</td><td>1</td><td>5.88</td><td>2736.8</td><td>2336.3</td><td>1.745 × 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>QG</td><td>1.62 × 10<sup>11</sup></td><td>Y<sub>2</sub></td><td>2</td><td>1</td><td>7.18</td><td>2392.5</td><td>2067.1</td><td>1.257 × 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>QG</td><td>1.62 × 10<sup>11</sup></td><td>Y<sub>2</sub></td><td>5</td><td>1</td><td>17.85</td><td>1538.2</td><td>822.0</td><td>8.055 × 10<sup>-2</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>QG</td><td>3.20 × 10<sup>11</sup></td><td>Y<sub>0</sub></td><td>—</td><td>1</td><td>10.86</td><td>2699.1</td><td>2280.1</td><td>1.653 × 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>QG</td><td>3.20 × 10<sup>11</sup></td><td>Y<sub>1</sub></td><td>3</td><td>1</td><td>119.83</td><td>1150.9</td><td>385.8</td><td>9.240 × 10<sup>-2</sup></td><td>(1537, 1536)/(-, -)</td></tr><tr><td>QG</td><td>3.20 × 10<sup>11</sup></td><td>Y<sub>2</sub></td><td>3</td><td>1</td><td>24.42</td><td>1752.4</td><td>1117.0</td><td>1.527 × 10<sup>-1</sup></td><td>(1537, 1536)/(-, -)</td></tr></tbody></table>

## APPENDIX C: ONSET OF CONVECTION

We compute the onset of rotating convection using the open source software SINGE<fnref n="2" /> (Vidal & Schaeffer 2015) for the 3-D configuration and the Linear Solver Builder package (LSB, Valdettaro et al. 2007) for the QG setup. In absence of a dedicated linear solver for the hybrid QG-3D configuration, we make the assumption that the critical Rayleigh number for this setup is the same as in the 3-D configuration, except for 3 cases at $Ek = \{10^{-4}, 10^{-5}, 10^{-6}\}$ and $Pr = 1$, for which, we have determined the onset by time-integrating the non-linear eqs (11–14–17) using pizza with an initial sectorial temperature perturbation and by bracketing the Rayleigh number until the critical value is attained. Fig. C1 displays the $Ra_c$ values obtained for the three methods at $Pr = 1$ and $10^{-4} \le Ek \le 10^{-8}$.

Both SINGE and LSB codes solve for the generalized eigenvalue problems formed by the linearized Navier–Stokes and temperature equations. They seek normal modes of the form

$$f(r, \theta, \phi) = \mathcal{F}(r, \theta) \exp(\lambda t + im\phi),$$

in the 3-D configuration and of the form

$$g(s, \phi) = \mathcal{G}(s) \exp(\lambda t + im\phi),$$

in the QG setup. Starting at a given $Ra$, the critical Rayleigh number $Ra_c$ for a given azimuthal wavenumber $m_c$ is attained when $\Re(\lambda) = 0$. Note that for the 3-D configuration, it becomes numerically demanding to determine the onset of convection using a linear solver for $Ek < 10^{-7}$ for $Pr \ge 1$ and for $Ek < 10^{-6}$; $Pr < 1$. For the cases with $Pr = 1$, we then resort to using the asymptotic expansion derived by Dormy et al. (2004) for spherical shells with differential heating (see their eq. 3.25a). For the remaining configurations the leading-order asymptotic scaling for the onset of rotating convection $Ra_c \sim Ek^{-4/3}$ is used. Table C1 summarizes the critical Rayleigh numbers $Ra_c$ and azimuthal wavenumbers $m_c$ for the different setups.

<footnote n="2">²https://bitbucket.org/vidalje/singe</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure C1.** Evolution of the critical Rayleigh number $Ra_c$ as a function of the inverse Ekman number $Ek^{-1}$ at $Pr = 1$ with fixed temperature difference across the shell, for our QG (using LSB, in red), hybrid (time-stepped with pizza for $Ek \le 10^{-6}$, in green) and 3-D (using SINGE for $Ek \le 10^{-7}$, in cyan) models, and compared with an analytical solution in the $Ek \to 0$ limit (Dormy *et al. 2004*, in black).</figcaption></figure>

<table><caption>**Table C1.** Summary of the critical Rayleigh numbers $Ra_c$ and critical azimuthal wavenumbers $m_c$ computed for our different setups. We have computed the onset of the hybrid method for three configurations but have otherwise assumed that hybrid QG-3D and purely 3-D runs have the same $Ra_c$.</caption><thead><tr><th>Ek</th><th>Pr</th><th>Setup</th><th>Ra<sub>c</sub></th><th>m<sub>c</sub></th><th>Computation method</th></tr></thead><tbody><tr><td>1 × 10<sup>-4</sup></td><td>10</td><td>QG</td><td>1.101 × 10<sup>6</sup></td><td>8</td><td>LSB</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>1</td><td>QG</td><td>7.121 × 10<sup>5</sup></td><td>8</td><td>LSB</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>1</td><td>Hybrid</td><td>6.82 × 10<sup>5</sup></td><td>8</td><td>Time-integrated with pizza</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>1</td><td>3-D</td><td>6.94 × 10<sup>5</sup></td><td>7</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>0.1</td><td>QG</td><td>2.911 × 10<sup>5</sup></td><td>6</td><td>LSB</td></tr><tr><td>1 × 10<sup>-4</sup></td><td>0.01</td><td>QG</td><td>1.663 × 10<sup>5</sup></td><td>4</td><td>LSB</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>10</td><td>QG</td><td>1.819 × 10<sup>7</sup></td><td>18</td><td>LSB</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>10</td><td>3-D and Hybrid</td><td>1.621 × 10<sup>7</sup></td><td>17</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>1</td><td>QG</td><td>1.050 × 10<sup>7</sup></td><td>16</td><td>LSB</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>1</td><td>Hybrid</td><td>9.83 × 10<sup>6</sup></td><td>16</td><td>Time-integrated with pizza</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>1</td><td>3-D</td><td>1.05 × 10<sup>7</sup></td><td>15</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>0.1</td><td>QG</td><td>3.623 × 10<sup>6</sup></td><td>11</td><td>LSB</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>0.1</td><td>3-D and Hybrid</td><td>3.53 × 10<sup>6</sup></td><td>11</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>0.01</td><td>QG</td><td>1.967 × 10<sup>6</sup></td><td>6</td><td>LSB</td></tr><tr><td>1 × 10<sup>-5</sup></td><td>0.01</td><td>3-D and Hybrid</td><td>1.952 × 10<sup>6</sup></td><td>6</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>10</td><td>QG</td><td>3.322 × 10<sup>8</sup></td><td>38</td><td>LSB</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>1</td><td>QG</td><td>1.773 × 10<sup>8</sup></td><td>33</td><td>LSB</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>1</td><td>Hybrid</td><td>1.63 × 10<sup>6</sup></td><td>33</td><td>Time-integrated with pizza</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>1</td><td>3-D</td><td>1.88 × 10<sup>8</sup></td><td>31</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>0.1</td><td>QG</td><td>5.24 × 10<sup>7</sup></td><td>23</td><td>LSB</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>0.1</td><td>3-D and Hybrid</td><td>5.076 × 10<sup>7</sup></td><td>23</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>0.01</td><td>QG</td><td>2.728 × 10<sup>7</sup></td><td>12</td><td>LSB</td></tr><tr><td>1 × 10<sup>-6</sup></td><td>0.01</td><td>3-D and Hybrid</td><td>2.587 × 10<sup>7</sup></td><td>12</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>10</td><td>QG</td><td>6.446 × 10<sup>9</sup></td><td>80</td><td>LSB</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>1</td><td>QG</td><td>3.259 × 10<sup>9</sup></td><td>69</td><td>LSB</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>1</td><td>3-D and Hybrid</td><td>3.321 × 10<sup>9</sup></td><td>67</td><td>SINGE</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>0.1</td><td>QG</td><td>8.287 × 10<sup>8</sup></td><td>47</td><td>LSB</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>0.1</td><td>3-D and Hybrid</td><td>2.219 × 10<sup>8</sup></td><td>—</td><td>Ek<sup>4/3</sup> extrapolation from Ek = 1 × 10<sup>-6</sup>; Pr = 0.1</td></tr><tr><td>1 × 10<sup>-7</sup></td><td>0.01</td><td>QG</td><td>3.868 × 10<sup>8</sup></td><td>26</td><td>LSB</td></tr><tr><td>3 × 10<sup>-8</sup></td><td>1</td><td>QG</td><td>1.529 × 10<sup>10</sup></td><td>102</td><td>LSB</td></tr><tr><td>3 × 10<sup>-8</sup></td><td>1</td><td>3-D and Hybrid</td><td>1.568 × 10<sup>10</sup></td><td>—</td><td>Dormy *et al. (2004)*</td></tr><tr><td>1 × 10<sup>-8</sup></td><td>1</td><td>QG</td><td>6.332 × 10<sup>10</sup></td><td>145</td><td>LSB</td></tr><tr><td>1 × 10<sup>-8</sup></td><td>1</td><td>3-D and Hybrid</td><td>6.492 × 10<sup>10</sup></td><td>—</td><td>Dormy *et al. (2004)*</td></tr><tr><td>1 × 10<sup>-8</sup></td><td>0.1</td><td>QG</td><td>1.383 × 10<sup>10</sup></td><td>98</td><td>LSB</td></tr><tr><td>1 × 10<sup>-8</sup></td><td>0.01</td><td>QG</td><td>5.548 × 10<sup>9</sup></td><td>56</td><td>LSB</td></tr><tr><td>1 × 10<sup>-8</sup></td><td>0.1</td><td>3-D and Hybrid</td><td>4.781 × 10<sup>9</sup></td><td>—</td><td>Ek<sup>4/3</sup> extrapolation from Ek = 1 × 10<sup>-6</sup>; Pr = 0.1</td></tr><tr><td>1 × 10<sup>-9</sup></td><td>1</td><td>QG</td><td>1.276 × 10<sup>12</sup></td><td>310</td><td>LSB</td></tr><tr><td>1 × 10<sup>-9</sup></td><td>0.1</td><td>QG</td><td>2.396 × 10<sup>11</sup></td><td>205</td><td>LSB</td></tr><tr><td>1 × 10<sup>-10</sup></td><td>1</td><td>QG</td><td>2.749 × 10<sup>13</sup></td><td>—</td><td>Ek<sup>4/3</sup> extrapolation from Ek = 1 × 10<sup>-9</sup>; Pr = 1</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>**Figure D1.** Convergence of the relative error for the two z-integration schemes applied to an analytical field as a function of the resolution: the scheme used for the z-average of the buoyancy (blue curve) and the z-integration scheme used for the thermal wind (green curve). The dotted grey-curve displays the order of convergence of each scheme. The analytical function used for this test is $z \sin(\pi z)$ (eq. D4).</figcaption></figure>

# APPENDIX D: CODE VALIDATION, BENCHMARKS AND SCALING

## D1 Benchmark of the z-integral functions

In this section, we describe numerical tests we have carried out in order to validate our hybrid QG-3D extension of pizza. In linking the QG and 3-D parts of the code there are two key aspects that require validation: (i) the z-averaging (18) of the 3-D buoyancy to obtain its contribution on the QG-grid and (ii) the z-integration of the 3-D temperature field to compute the thermal wind contribution (20) on the 3-D grid.

In order to discuss the validation and the accuracy of our numerical schemes, we define a relative error estimate, $e_{rel}$, as

$$e_{rel}(f) = \left[ \frac{\{ (f_{\text{ref}} - f)^2 \}_s}{\{ f_{\text{ref}}^2 \}_s} \right]^{1/2}, \qquad (\text{D1})$$

where the brackets in the above equation correspond to an average over the annulus

$$\{f\}_s = \frac{1}{V} \int_0^{2\pi} \int_{s_l}^{s_o} f h(s) s \, ds \, d\phi. \qquad (\text{D2})$$

### D1.1 Analytical benchmark of z-averaging

Our computation of the z-averaged buoyancy term in eq. (17) relies on a simple bilinear interpolation and a basic averaging: $N_z = 2N_s$ points are regularly distributed between $h$ and $-h$ for each cylindrical radius $s$ on the QG-grid and the 4 nearest points of the 3-D grid are used to interpolate the value of the field at the corresponding $z$ points; then a simple summation divided by the number of points is performed along the $z$-direction at each location in the equatorial plane to obtain the average. The bilinear interpolation is expected to have an accuracy of order 2 and the basic summation an accuracy of order 1. We choose a simple and fast order 1–2 scheme because retaining the speed and efficiency of the QG approach is our priority and it is readily parallelized. It has also been chosen in previous studies (Guervilly 2010; Guervilly & Cardin 2016).

We tested our parallel implementation of this z-averaging scheme by comparison with an analytical solution. We considered the following 3-D field

$$f'(r, \theta, \phi) = z \sin(\pi z), \quad \text{with } z = r \cos \theta, \qquad (\text{D3})$$

whose z-integral is

$$f(s, \phi) = \int_{-h}^{h} f'(r, \theta, \phi) dz = \left[ \frac{\sin(\pi z)}{\pi^2} - \frac{z \cos(\pi z)}{\pi} \right]_{-h}^{h}. \qquad (\text{D4})$$

Then we compute the relative error value, $e_{rel}$, obtained between our scheme and this analytical solution while testing different grid sizes, to test the accuracy and the convergence of our scheme. We set $N_r = N_s = N_z/2$, $N_m = N_{\phi,3-D} = 2 \times N_\theta$, and we vary the grid size from $(N_s, N_m) = (16, 16)$ to $(N_s, N_m) = (512, 512)$. The results are displayed as a function of the grid resolution in Fig. D1 (blue curve). The maximum value of the averaged function $f$ is 0.338396. As expected, we find that the accuracy varies approximately as 1.2–1.4 times the grid resolution.
The accuracy ($\sim \{1 \times 10^{-1}, 4 \times 10^{-3}\}$) and the convergence of our scheme (order 1.2–1.4) validates our implementation of eq. (17).

## D1.2 Analytical benchmark of z-integration for thermal wind

Turning to the implementation of the thermal wind (20), our implementation relies on a 2-neighbours interpolation scheme and performs a z-integration between any point, $(r, \theta)$, on the 3-D grid and the half-height of the column of fluid above this point, $h$, by interpolating and summing the $N_z$ points across the $\theta$ lines directly above the position of interest. Because only 2 points are involved in the interpolation and the integration is again a basic sum, we expect the scheme to have an accuracy of order 1. In order to test the accuracy of this procedure we compare against the integral of the same analytical function considered in the previous test (eq. D4) but here evaluate the analytic integral between $z$ and $h$. We again compute $e_{rel}$ compared to the analytical solution (D4), fix $N_r = N_s = N_z/2$, $N_m = N_{\phi 3-D} = 2 \times N_{\theta}$ and vary the grid size from $(N_s, N_m) = (16, 16)$ to $(N_s, N_m) = (1024, 1024)$.

The results are displayed in Fig. D1 (green curve) and we can observe that the accuracy at the lowest resolution is $\sim 2 \times 10^{-1}$ and $7 \times 10^{-3}$ at the highest resolution. The maximum value of the integrated function $f$ is 0.462040. We find that the scheme for the z-integration of the thermal-wind contribution has the expected accuracy of order 1 (see the slope), validating our implementation of the scheme.