Comparison of Quasi-Geostrophic, Hybrid and 3D models of planetary core convection

We present investigations of rapidly-rotating convection in a thick spherical shell geometry relevant to planetary cores, comparing results from Quasi-Geostrophic, 3D 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}$. 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 specific geometry studied here the hybrid method is best suited for studying convection at modest forcing, $Ra \leq 10 \, Ra_c$ when $Pr=1$, and departs from the 3D model results at higher $Ra$, displaying systematically lower heat transport. We find that the lack of equatorially anti-symmetric 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 3D models, for the specific range of parameters explored here. We cannot point to major disagreements between these two datasets at $Pr \geq 0.1$, although the QG model is effectively more strongly driven than the hybrid case. When $Pr$ is decreased, the range of agreement between the Hybrid and 3D models expands, indicating the hybrid method may be better suited to study convection in the regime $Pr \ll 1$. Previously proposed scaling laws for rapidly-rotating convection are retrieved: 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.

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 ≪ 1 and Ro ≪ 1 (Julien et al. 2012), and when the typical time-scale of the convection is much longer that the rotation period, where Ek = ν/Ωd 2 measures the viscous effects compared to the Coriolis force and Ro = ReEk measures the nonlinear inertial effects compared to the Coriolis force.Here Re = U d/ν is the Reynolds number with ν the kinematic viscosity of the fluid, Ω 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 ≫ 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 nonlinearities 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 ∼ 10 −15 , Re ∼ 10 9 and Ra/Rac ≫ 10 3 where Ra is the Rayleigh number and Rac is the critical Rayleigh number for the onset of convection.The most ambitious 3D simulations and experimental studies are only able to reach Ek ∼ 10 −7 ; Ra/Rac ∼ 10 3 and Re ∼ 5 × 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 2D approximation of the real 3D situation.The 2D 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 3D 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 3D 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 3D temperature field in a spherical shell geometry.Here we explore advantages and limitations of QG and hybrid QG-3D models compared with full 3D 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 heatflux 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 ≥ 2 × 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 ≥ 10 −6 with our hybrid approach and Ek ≥ 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 P r < 1.Here our main goal is to complement these studies, using QG, 3D and hybrid QG-3D simulations in a thick spherical shell geometry, focusing on relatively strongly-driven cases (high Ra/Rac) 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 article is organised as follows: Section 2 presents the equations and the methodology of our QG, hybrid QG-3D, and 3D models.Section 3 presents results obtained using our models focusing on comparisons between QG, Hybrid and 3D calculations, and including cases with inhomogoneous 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.

Quasi-Geostrophic model formulation
In this section we first describe the basic QG model employed before moving on to the new 3D 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, ϕ, z) in a spherical shell between the inner radius si and the outer radius so, rotating about the z-axis with a constant angular velocity Ω.The horizontal components of the velocity field u ⊥ , perpendicular to the rotation axis, are assumed to be invariant along the rotation axis, i.e. u ⊥ = (us, u ϕ , 0), where us and u ϕ 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 = so − si as the reference length-scale, the viscous diffusion time d 2 /ν as the reference time-scale, and the temperature contrast between the boundaries ∆T = Ti − To = T (si) − T (so) as the reference for temperature.Throughout this study, we adopt η = si/so = 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) ∝ s and it is non-dimensionalized based on its value at the outer boundary go = g(so).
Following Schaeffer & Cardin (2005) and Gastine (2019) it is assumed that the axial velocity uz varies linearly with z in the direction of the rotation axis, including contributions from the radial velocity us and the Ekman pumping P, i.e.
with P(Ek, u ⊥ , ωz) the Ekman pumping term deduced from Greenspan's formula in a rigid sphere (see Eq. 8), and h ≡ √ s 2 o − s 2 , the half-height of a cylinder aligned with the rotation axis at a radius s.The spherical-QG continuity equation then reads We represent the non-axisymmetric QG-velocity by a streamfunction ψ such that which ensures that (3) is satisfied.u ϕ is the remaining axisymmetric zonal flow component.The overbar x denotes the azimuthal average of any quantity x, such that x ≡ 1 2π 2π 0 x dϕ . (5) The axial vorticity ωz = ez • ∇ × u is then 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 2D form of the contributing fields, may be written where the subscript ⊥ corresponds to the horizontal part of the operators -e.g., ∇ ⊥ = es • ∂s + r −1 e ϕ • ∂ ϕ , u ⊥ = (us, u ϕ , 0)and P(Ek, u ⊥ , ωz) corresponds to the Ekman-pumping contribution (Schaeffer & Cardin 2005) for the non-axisymmetric motions with The non-dimensional control parameters, the Ekman number, the Rayleigh number and the Prandtl number are respectively defined by where αT is the thermal expansion coefficient, ν is the kinematic viscosity, and κ is the thermal diffusivity.The z-averaged axial vorticity equation ( 7) has to be supplemented by an equation to account for the axisymmetric motions.This is obtained by taking the ϕ-average of the azimuthal component of the Navier-Stokes equation and reads where the last term of 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 §2.4.The other coupled prognostic equation used to complete the system is the QG-temperature perturbation equation, where the temperature is written as a perturbation about a mean 2D conducting state, i.e.T2D = T cond 2D + ϑ2D, where T cond 2D is the conducting background state, a solution of ∇ 2 T cond 2D = 0 subject to the chosen boundary conditions.For fixed-temperature boundary conditions at si and so this yields Further details on the boundary conditions for the temperature field, including the possibility of heat flux boundary conditions are given in §2.4.

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 equation ( 12) is replaced by the full 3D temperature equation where r is the spherical radius and u3D = (ur, u θ , u ϕ3D ) is the 3D-velocity in spherical coordinates.Similarly to the QG case, we write the temperature as a perturbation temperature about the conducting background state, i.e.T3D = T cond 3D + ϑ3D.We compute the conducting temperature profile, or dimensionless radial temperature profile, T cond 3D , as the solution of ∇ 2 T cond 3D = 0, which for a fixed temperature contrast between ri and ro without internal heating, yields This is the full 3D version of Eq. ( 13) in spherical geometry.We reconstruct the 3D-velocity field, u3D using the conversion between cylindrical and spherical coordinate systems as needed in this equation, from the QG velocity field, such that where uz is proportional to z and incorporates the effects of the Ekman pumping (see Eq. 1) consistent with our initial quasigeostrophic assumption (Eq.3), and where the cylindrical quantities are cast onto the 3D-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 occurs in that region.
Considering a 3D temperature field, Eq. ( 7) becomes where the angular brackets ⟨x⟩ refer to the axial or z-average of any quantity x defined by The equation for the zonal motions ( 11) is not modified as it does not involve the temperature.
Since we treat the temperature in 3D, 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 ϕ3D , and which satisfies Integrating between the position z and the height of the column above the equator, h, we obtain where g(r) = r/ro is the 3D 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 3D 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 3D case, it is natural to expect the hybrid QG-3D model to behave better than the classical QG model when compared with a full-3D model, a hypothesis that will be further assessed in the Results Section.

3D model formulation
In order to compare the results of our QG and hybrid QG-3D models, we also consider a purely 3D model.Similarly to the two previous setups, we consider convection of a Boussinesq fluid enclosed in a spherical shell of inner radius ri and outer radius ro rotating about the z axis.The same scales and dimensionless parameters are used and thus the 3D Navier-Stokes equations read where P is the pressure, and er,z are respectively the unit vectors in the radial and the axial directions.

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 so and si in the QG or Hybrid models and at ri and ro in the 3D calculations.
Imposing fixed temperature at the boundaries yields, respectively for the QG-and the 3D-temperature field ϑ3D = 0, at r = {ri , ro} .
In this study, the majority of our simulations are conducted under these boundary conditions, but we also consider an other set of thermal boundary conditions with a fixed temperature at the inner radius and an imposed flux at the outer boundary.The latter involves ∂ϑ ∂r = 0, at r = ro, ϑ = 0, at r = ri . ( where ϑ can either be ϑ2D or ϑ3D.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 §3.3).With heat-flux boundary conditions at the outer boundary the radial conductive profiles become where Ti and Qo 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, i.e.
Lateral variations in the amplitude in the heat flux are then defined by

Numerics
The calculations presented here were carried out using an extension of the open-source pseudo-spectral spherical QG code pizza (Gastine 2019) -freely 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 ϕ 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 multistep 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 3D simulations were computed with the opensource 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, 3D fields in pizza and MagIC are expanded in Spherical Harmonics up to the degree and order ℓmax in the angular (θ, ϕ) directions and in Chebyshev polynomials with Nr collocation grid points in the radial direction.The 2D quantities in pizza are expanded in Fourier series up to the degree Nm in the azimuthal direction and in Chebyshev polynomials up to degree Ns in the radial direction.The open-source SHTns ⋆ library is employed in both codes to handle the Spherical Harmonic Transforms (Schaeffer 2013).Parallelisation 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 employed finite differences in radius while we resort to a Chebyshev collocation method.In the hybrid setup, the z-extrapolation of the variables from the 2D-grid to the 3D-grid is computed using Eq. ( 16).Reduction of the quantities from the 3D grid back onto the 2D grid, and the computation of the thermal wind, relies on two zintegration functions described in the Appendix D1.For clarity, all ⋆ https://bitbucket.org/nschaeff/shtns3D-quantities are labelled with a subscript 3D (such as u ϕ3D ), QG quantities have no subscripts.

Posterior diagnostics
We next introduce the various diagnostics and notations that are used for the analysis of the simulations.The Nusselt number N u, which characterises 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, i.e.
where, the temperature perturbation ϑ can be either ϑ2D or ϑ3D and T cond is either T cond 2D (13) or T cond 3D (15) and in QG calculations r and ro are replaced by s, so.Note that in the heat flux boundary case, this definition of N u leads to N u = 1 because ∂rϑ(ro) = 0, so following Goluskin (2016) we instead use In the above expressions, x corresponds to the time average of any quantity x, such that with τ the averaging time.
The dimensionless kinetic energy, Ekin per unit volume, is defined by where V corresponds to the full spherical shell volume in the 3D 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ϕ.
From this expression, we define a diagnostic for the fluid velocity which characterises the average flow speed, based on the rootmean-square (r.m.s.) of the velocity, and which is denoted by the Reynolds number 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), where Ezon is the dimensionless axisymmetric kinetic energy per unit volume, similar to (34), that is defined by Ezon = 1/2V V u ϕ 2 dV ′ and is associated with the time-averaged zonal Reynolds number, Finally, for the typical flow length-scale we use the typical cylindrical radial velocity length-scale L −1 us , determined from the time-averaged us energy spectrum

RESULTS
We present here results of rapidly-rotating convection, focusing on a regime well above the onset of convection, i.e.Ra ⩾ 5 Rac.
We explore the Ekman numbers from Ek = 10 −4 down to Ek = 10 −10 , and consider Prandtl numbers from P r = 10 down to P r = 10 ).In total 144 runs have been performed and a list of their key diagnostics is given in Table A1 in Appendix A.
Figure 1 summarises all the runs we have carried out with fixed temperature contrast for this study in terms of their heat transfer N u 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 3D runs are marked with hatched and empty symbols respectively.

Heat transfer
Previous work by Guervilly et al. (2019) have 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 P r 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 ≤ 10 −10 ) with a focus around Ek = 10 −6 in order to facilitate comparisons with full-3D simulations.All our runs have Ro ≪ 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 §3.5).
Figure 2 displays all the Nusselt numbers of our data set versus the supercriticality, i.e.Ra/Rac.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 3D configurations, although we have used the asymptotic expansion by Dormy et al. (2004) for P r = 1 when Ek < 10 −7 in this latter case.Concerning the hybrid QG-3D model, we have determined the onset for 3 configurations -at P r = 1 and Ek = {10 −4 , 10 −5 , 10 −6 } -by time-integrating the nonlinear equations (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 ≤ 10 −7 or P r ̸ = 1 we have assumed the same critical Rayleigh as for 3D configurations.In all cases, a simple extrapolation employing the asymptotic scaling for rotating convection Rac ∼ Ek −4/3 has been used whenever the aforementioned techniques could not be applied for practical reasons.Concerning P r = 1, the Rac values obtained with LSB and SINGE methods agree within ∼ 6% at all Ek and follow the expected converging trend (Dormy et al. 2004).The Rac value of the hybrid model is ∼ 2% lower than that of the 3D at Ek = 10 −4 and ∼ 13% lower at Ek = 10 −6 .The mc 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 Rac.
We can observe that for the lowest supercriticalities (Ra ≤ 10 Rac) all the points in the weakly non-linear regime follow a power law of the form N u − 1 ∼ Ra/Rac − 1.For stronger forcing with Ra > 10 Rac, the numerical models seem to approach an asymptotic behaviour of the form N u ∝ (Ra/Rac) α (black and blue dotted-line for the QG and 3D 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 ri/ro = 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 N u ∼ Ra 3/2 Ek 2 P r −1 put forward by Julien et al. (2012) and retrieved in the 3D spherical shell computations by Gastine et al. (2016) when ri/ro = 0.6.The QG runs are slightly offset compared to 3D cases towards larger Nusselt numbers for the same supercriticality (Ra/Rac).Strikingly however, this asymptotic scaling is followed only by the QG and 3D simulations, while the hybrid runs follow a much shallower slope N u ∼ (Ra/Rac) 1/5 .Several outliers also appear in the QG and 3D configurations: the series of QG points at P r = 10 seem to follow a different slope with N u values considerably higher than the values obtained at P r ̸ = 10 for the same Ra/Rac ratio.All the QG and 3D runs at P r = 0.1 (filled and empty diamonds) lie below the mean trend, suggesting that the heat transport is less efficient for the same supercriticality when P r < 1.Thus, the purely 3D and QG cases are generally in agreement and we find two different behaviours (splitting around Ra/Rac ∼ 10) with 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 3D configuration for Ra > 10 Rac (at P r = 1) and follows a much shallower scaling behaviour.In constrast, when P r = 10 −2 the QG and 3D models disagree above Ra/Rac = 10, while the range of accordance between the hybrid QG-3D and the 3D configurations seems to extend to higher Ra, up to Ra > 15 − 20 Rac at P r = 0.1 − 0.01, suggesting the range of agreement between the two latter models may be larger at low P r and low Ra/Rac (as suggested in previous studies, e.g. in Guervilly & Cardin 2016, 2017;Guervilly et al. 2019).

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 3D and QG simulations at modest driving, which we define here to be the parameter regime where Ra ≤ 20 Rac, for a case at Ek = 10 −6 , P r = 1, and Ra = 2 × 10 9 (= 10.6 Rac for the 3D, = 12.3 Rac for the Hybrid, and = 11.3Rac for the QG setups).
Fig. 4 additionally displays the ϕ-averaged temperature field T3D (top panel) in the purely 3D (a), Hybrid (b) and QG (c) cases as well as equatorial sections of the z-averaged vorticity ωz (middle panel) in 3D (d), hybrid (e) cases, and in the purely QG case (bottom -f).Note that in the 3D configuration, this involves a zaveraging 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 3D 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 3D and Hybrid cases, these two jets display a similar columnar structure that span 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 Rezon ≪ Rec in all three cases (the exact values are given at the end of this 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 3D 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 3D 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 3D case, although there are discrepancies towards the inner and outer boundaries.
Turning to the azimuthally-averaged temperature fields (Fig. 4(a-b-c)), we find that the profiles in the 3D 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 3D 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 uz.For   For the three temperature and the three vorticity plots respectively, the same colorscales are used.Note that the colorscale for the vorticity is saturated to highlight the fine structure of the flows.
these parameters, the 3D 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.
Figure 4(c-d-e) shows a comparison of the axial vorticity ωz between the purely 3D, 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 ∼ so − 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 3D 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 |β| (2) increases with the cylindrical radius.This treatment impedes radial motions and favors the propagation of thermal Rossby waves over the advective processes as the stretching term due to βus 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 3D and hybrid cases are rather similar with a convective Reynolds Rec, a zonal Reynolds Rezon and a Rossby number Ro of respectively 448.1, 97.1 and 4.58 × 10 −4 in the 3D-case, and 404.4,86.9 and 4.14 × 10 −4 in the hybrid case.The Nusselt number N u differs more strongly with values of 2.38 in the 3D 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×10 −4 , respectively for Rec, Rezon and Ro while N u = 5.05.This example indicates how the hybrid approach is capable of accounting the 3D convective dynamics happening in the fluid bulk at modest driving, i.e. here with Ra = 10.6 Rac.

Limitations of the hybrid approach
We now compare results from our hybrid QG-3D model with 3D and QG simulations for a more strongly driven case at Ek = 10 −6 , P r = 1, and Ra = 10 10 (= 53.2 Rac for the 3D, = 61.3Rac for the Hybrid, and = 63.4Rac for the QG setup).Below we use the term 'strong driving' to refer to the parameter regime where Ra > 20 Rac.
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 × 10 9 , while in the purely 3D 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 3D 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 3D cases now show significant differences.The vorticity in the 3D 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 3D configuration and is already vigorous at these parameters.Interestingly, the QG case seems to be closer to the 3D case than the Hybrid, with filaments of vorticity strongly sheared in the ϕ-direction and a vigorous convective pattern degenerating into thin Rossby waves towards the outer boundary.Overall, the purely 3D 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 Rec, Rezon, Ro and N u that are respectively equal to 3339.8, 2405.2, 4.12 × 10 −3 , and 19.5 in the 3D case; 1455.9, 829.5, 1.68 × 10 −3 , and 2.20 in the hybrid case; and 2584.4,1845.5, 3.18 × 10 −3 , and 36.3 in the QG case.Between the N u and Rec numbers of the Hybrid and 3D models, we have observed a relation of the form N u Hyb /N u 3D ∝ (Re Hyb c /Re 3D c ) 2/5 .Examining the kinetic energy spectra for the horizontal velocity, shown in Fig. 6, we find that the QG, 3D and hybrid models show similar decreasing slopes up to m ∼ 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 T3D(r) (a) and example snapshots of z-averaged cylindrical temperature profiles ⟨T3D⟩ (s) (b) obtained with the purely 3D and the hybrid setups at Ek = 10 −6 and P r = 1 for a series of increasing supercriticalities, ranging from Ra = 10 9 = 5.3 Rac up to Ra = 2.66 × 10 10 = 141.5 Rac.We retrieve the fact that the profiles are fairly similar when Ra ≤ 2 × 10 9 ∼ 10 Rac (green and cyan curves), while the Hybrid and purely 3D 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 3D case while the hybrid temperature profiles do not vary much (red and yellow curves) and stay close to the radial conducting state (blue dashedcurve).The same conclusions can be drawn when looking at the zaveraged cylindrical radial temperature profiles (Fig. 7(b)) although here we can observe that the inner boundary temperature decreases with increasing Ra in both setups with a slightly larger decrease of ⟨T3D⟩ (si) in the 3D case compared to the hybrid case.Note the increased activity inside the tangent cylinder in the 3D case when Ra ≥ 10 10 which is not accounted for in our hybrid QG-3D model.The evolution of ⟨T3D⟩ (si) in both the 3D and hybrid cases is important: the temperature at the ICB drops significantly when Ra  is increased, suggesting that approximating the fixed 3D temperature at ri by a fixed temperature at si in the QG-approach (about T2D(si) ∼ 0.445 at all Ra) is rather crude (again see Fig. 5(c)).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 3D setup at very large Rayleigh numbers before the 3D convection reaches the non-rotating regime.
The fundamental difference between the 3D 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 anti-symmetric parts and the z-component of the buoyancy force are missing in the QG and hybrid models compared to the full 3D 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 3D 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 Ra P r urgT3D, which we can decompose into its equatorially-symmetric (ES) and equatorially anti-symmetric (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 3D and the hybrid models for the two cases of §3.2.1 and §3.2.2, i.e. at Ek = 10 −6 , P r = 1, Ra = 2 × 10 9 and at Ek = 10 −6 , P r = 1, Ra = 1 × 10 10 , respectively.At Ra = 2 × 10 9 , we retrieve a similar buoyancy power profile for the hybrid and 3D models, although the hybrid model has less energy especially towards the CMB.EAS modes become noticeable between Ra = 2 × 10 9 (where they are almost zero) and Ra = 5.5 × 10 9 (where they account for 9% of the total power) and grow increasingly strong, reaching 23% of the convective power by Ra = 10 10 .In addition, we find that 10% of the power is driven by the Ra P r uzgzT3D at Ra = 2 × 10 9 which grows to 30% at Ra = 10 10 .
At Ra = 2 × 10 9 , we find that the peak-to-peak ratio between the convective power in the 3D and hybrid models is around 1.5 from which (based on the IAC scaling, see Sec. 3.5 and Fig. 13)) Figure 7. Top: spherical time-averaged radial temperature profiles T3D (r) of 4 cases conducted at the same parameters at Ek = 10 −6 and P r = 1 using our Hybrid method (dotted lines) or a 3D model (full lines) varying Ra from Ra = 10 9 (green curves) up to Ra = 2.66 × 10 10 (red curves).Bottom: Same as above but for the cylindrical radial z-averaged temperature profiles ⟨T 3D ⟩ (s).Note that the latter profiles are not time-averaged and were derived from snapshots.
we expect a ratio between the velocities of 1.5 2/5 ∼ 1.17, close to the actual ratio of Re 3D c /Re Hyb c = 448.1/404.4∼ 1.11.At Ra = 1 × 10 10 , we have Re 3D c /Re Hyb c = 3339.9/1455.9∼ 2.29 which requires a ratio of 2.29 5/2 ∼ 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% of the total volume of the shell) and the enforced linearity of uz (see Eq. 1) may also contribute to missing power at high Ra/Rac, 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 3D, QG and Hybrid come from differences in the underlying convective power.The Hybrid and QG models lack the equatorially anti-symmetric 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 uz 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 3D cases.In the hybrid con- figurations, 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 3D computations which is limited to Ra ≤ 10 Rac, at P r = 1, and has a decreasing predictive capacity with increasing Ra.On the contrary, we cannot point to important disagreements between the QG and 3D datasets, suggesting that the QG model retains its predictive power, even at high Ra/Rac, at least in this particular configuration.The impact of that conclusion on the whole dataset will be further discussed in §3.5.

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 A2 in Appendix B. We first present in detail a comparison of results obtained from QG, 3D and hybrid QG-3D calculations carried out at the same parameters Ek = 10 −6 , P r = 1, and RaQ = 8 × 10 9 with a ℓ = m = 2, Q * = 3 lateral variation about the imposed heatflux condition.The resolutions in the 3D, Hybrid and QG cases are respectively (Nr, ℓmax) = (321, 341), (Ns, Nm)/(Nr, ℓmax) = (385, 416)/(385, 416) and (Ns, Nm) = (513, 512).Fig. 9 displays example snapshots of the 3D temperature field at the CMB T3D(ro) using heat-flux boundary conditions for the 3D model in (a) and the hybrid model in (b).Both cases show the expected ℓ = 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 3D cases spanning −1.49≤ T3D(ro) ≤ 0.38 and 0.23 ≤ T3D(ro) ≤ 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 3D configuration, as discussed in section §3.2.2.This is also reflected in global diagnostics, with lower N u∆ values of 2.03 in the hybrid case compared with 7.15 in the 3D 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 (Fig. 9(c-d-e)) we observe significant differences between the three cases.The vorticity structures up to r ∼ 2/3ro are qualitatively similar in the 3D 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 3D 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 center 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 3D case, while the purely QG is slightly over-driven compared with the 3D case.Note that standard global diagnostics can be misleading here as they involve averages over the entire shell.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 ∼ so − 1/3 and displays only a wide spiraling arm covering this region, a result similar to previous 3D 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 RaQ = 8 × 10 9 .The boundary perturbation does not penetrate very deep in the shell at these parameters (s ∼ so − 1/3 in the m = 1 case, and s ∼ so − 1/4 in the m = 2).The final case (d) has been conducted at a larger forcing (RaQ = 3.6 × 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 s ∼ so − 1/6.It may be that in the limit of very large RaQ convection, the region affected by the inhomogeneous bound-ary 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 3D computations (Mound & Davies 2017).

Influence of the Prandtl number at low Rossby number
Focusing on the Hybrid and 3D series at Ek = 10 −6 , P r = 0.1, we find the same limitations as in Sec.3.2.2, with an apparent lack of energy in the hybrid configuration when Ra is increased (see Fig. 2).However as P r is decreased, larger velocities are attained at smaller Nusselt numbers, and the range of agreement across the Hybrid and 3D configurations expands (up to Ra ∼ 15 × Rac at P r = 0.1).We observed similar convective patterns between the Hybrid and 3D models at all P r (1, 10 −1 , 10 −2 ) when Ra ≤ 15×Rac.Unfortunately, lowest Prandtl runs are extremely computationally costly to run, especially with the 3D approach, because powerful zonal flows and velocities are triggered even at low N u inducing difficulties to reach a converged power balance.For example, for a run at Ek = 10 −6 , P r = 0.01, Ra = 5 × 10 8 = 19.3Rac, the convective Reynolds numbers reach 11968.4 for a N u of 2.07 in the 3D case; Rec values that were not even reached by Ra = 140 Rac at P r = 1.The two 3D runs at P r = 10 −2 however appear closer to the Hybrid trend of Fig. 2 with comparable N u up to Ra/Rac = 20, whereas the QG model departs from the 3D trend around Ra ∼ 10 Rac with much higher N u in the QG cases.Our results do thus indicate better agreement between the Hybrid and the 3D configurations at low P r (as observed in 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 3D when Ra is sufficiently increased, as observed for our P r = 0.1 series.The exact Ra/Rac range of agreement when P r is decreased further below 0.1 requires more 3D and hybrid runs in this challenging regime in order to be determined.
Thus, we now take advantage of the computationally more efficient purely QG setup to study the convective flows at more extreme parameters, that is 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 ωz (top panel) and of azimuthal velocity u ϕ (bottom panel) in the equatorial plane for two cases: Ek = 10 −8 , P r = 1, Ra = 8.99 × 10 12 = 142 Rac (left column) and Ek = 10 −8 , P r = 10 −2 , Ra = 3.37 × 10 10 = 6.1 Rac (right column).Despite the much higher supercriticality attained in the P r = 1 case, these two runs have comparable Rossby numbers: 4.10×10 −4 for the P r = 1 case and 5.97 × 10 −4 for the P r = 10 −2 case.Both cases are purely QG calculations, and 3D temperature effects have not been included.
In both cases, the azimuthal velocity (c-d) displays multiple zonal jets of alternating direction (blue is retrograde and red is prograde flows), which directly translate in the axial vorticity (a-b) into alternating rings of cyclonic (ωz > 0 in red) and anticyclonic (ω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, (ω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/β) 1/2 (Rhines 1975).Closer to the boundary, we see the influence of the slope and the β-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 motions typical of thermal Rossby waves, as already observed in §3.2.1.Figure 12 additionally shows the time-averaged radial profiles of potential vorticity along with the time-averaged zonal flows.Retrograde zonal jets where potential vorticity gradients are slightly stronger (marked by white stripes in Fig. 12), seem narrower than the regions where the gradients are weaker (corresponding to prograde jets), a result already observed by Guervilly & Cardin (2017).Since large supercriticalities are required to obtain well-formed potential vorticity staircases, it also appears on this figure that the case with P r = 10 −2 does not show a comparable degree of homogeni-sation, due to the significantly lower Ra/Rac reached in that case, despite similar values of Ro.
Besides the two cases having comparable Rossby numbers and time-averaged kinetic energy spectra, the P r = 1 case (Fig. 11(ac)) displays smaller eddies and thinner jets than the P r = 10 −2 case with a larger number of coherent jets (8 in the P r = 1 case compared to 3 in the P r = 10 −2 case).Near the outer boundary the transition of the dynamics into thermal Rossby waves happens deeper in the shell in the P r = 10 −2 compared to the P r = 1 case.This transition is also visible in the ϕ-velocity where the jets lose their coherence around s ∼ so − 1/3; a direct consequence of the lower supercriticality attained in the P r = 10 −2 case (Ra = 6.1 Rac in that latter case, compared to Ra = 142 Rac for the P r = 1 case).
In the low Rossby regime explored here (Ro < 6 × 10 −4 ), changing the Prandtl number by a factor 100 drastically modifies the form of the convective pattern.We find that decreasing P r results in wider and fewer 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 P r (see Table A1 or the next section).

Scaling laws for rapidly-rotating convection
We now finally explore the scaling behaviour of rapidly-rotating turbulent convection in our three different model setups.Theoretical scaling laws of rotating convection can be derived by considering the following dimensional 3D vorticity equation In the limit of rapid rotation, the Proudman-Taylor theorem promotes z invariant flows with l ⊥ ≪ l // , where l ⊥ and l // correspond to the convective flow length-scale perpendicular and parallel to the rotation axis.Assuming that l // ∼ d, this implies that the gradients orthogonal to the axis of rotation ∇ ⊥ can be approximated by 1/l ⊥ , while the axial gradients ∂/∂z simply scale as 1/d.It also follows that ω ∼ Uc/l ⊥ , where Uc 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) or in terms of scaling quantities where Θ is a typical temperature perturbation.The balance between Coriolis and Inertia yields This diffusion-free scaling is commonly known as the Rhines scaling (Rhines 1975) and it is expected to hold in the limit of Ek ≪ 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 αT gUcΘ ∼ ν 3 d 4 Ra(N u − 1)P r −2 (see Jones 2015) yields in its dimensionless form 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, Rec ∼ [ Ra P r 2 (N u − 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 quasigeostrophic 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.Figure 13 shows all our numerical simulations compared with the CIA scaling laws for convective velocity and length-scale.On the top panel, the typical nondimensional length-scale of the convection Lu s is plotted against the Rossby number Ro = RecEk, corresponding to (42), while the bottom figure shows Rec as a function of RaQP r −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 dot-dashed line), and departures to the theory are confined to the highest Ekman numbers (Ek ≥ 10 −5 ).Introducing a local Rossby number RoL = Ro d/Lu s we found that all our runs have RoL < 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 RoL (associated with the highest Ek) of our dataset.There is no additional P r dependence since all the simulations with P r ̸ = 1 are close to the average scaling behaviour.The best-fit for the whole data-set yields a power law with an exponent of 0.31 for the Ro dependency but considering only for the runs with Ro ≤ 10 −3 , we find a steeper slope with an exponent of α ∼ 0.41, a value comparable to that obtained in 3D 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 ⊥ /d ∼ 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 ≤ 10 −6 .Cases at higher Ekman number (Ek ≥ 10 −5 ; red and orange symbols) depart from the theory, following instead a power  43).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 3D runs with empty symbols.The theoretically predicted scaling corresponding to the CIA balance (Eq.( 43) is shown by the dot-dashed line.
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 P r = 1 series but shifted upward for P r < 1 and downward for P r > 1.The clear separation between the series at different Prandtl numbers suggests that there is a dependency on P r in the Rec 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 § 3.2.2 with an effective lack of buoyancy power when Ra is increased, yielding lower velocities -i.e.Rec -and lower heat transport effectiveness -i.e.N u, hence weaker Ra(N u − 1) -for the same control parameters {Ek , P r , Ra}.Note that the purely 3D 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 3D 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 3D setup (see Fig. 7(b)).Parameter studies with either QG or 3D 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 Rec < 10 4 .Discrepancies arise at high Ek and low Ra/Rac where the VAC balance is probably more suitable and at high Ra/Rac 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 P r = 1 yields Rec = 0.53 [Ra(N u − 1)Ek 1/2 ] 0.43 in reasonable agreement with the theory, but the P r 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/Rac 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 P r 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).

SUMMARY AND DISCUSSION
We have used QG, 3D and Hybrid models, the latter involving a QG velocity field and a 3D temperature field, to explore in a thick spherical shell the regime of strongly-driven, rapidly rotating convection, focusing on low Ekman numbers 10 −10 ≤ Ek ≤ 10 −4 , reaching supercriticalities up to Ra ∼ 160 Rac and considering a range of Prandtl numbers close to and below unity 10 −2 ≤ P r ≤ 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 laterallyvarying heat flux boundary conditions, a 3D-temperature field, and a thermal wind.
Using the hybrid QG-3D approach at parameters Ek ≤ 10 −6 and P r ≤ 1 we are able to reproduce important aspects of convective dynamics seen in 3D models for weak to moderate supercriticalities (Ra/Rac ≤ 10 − 15).In that regime, the meridional temperature profile and ϕ-averaged azimuthal velocity are well retrieved although, as for purely QG models, the dynamical behaviour in the hybrid model deviates from the 3D models close to the outer boundary.When Ra is further increased, we find our hybrid model develops much weaker convective flows compared with the 3D configuration in all cases.In contrast, the range of agreement between the hybrid QG-3D and the 3D configurations increases when P r = 10 −2 (as suggested by Guervilly et al. 2019) while the QG model departs from the 3D around Ra ∼ 10 Rac.We expect that the hybrid model will eventually diverge from the full 3D model when Ra/Rac is sufficiently high, but the exact value of the diverging point remains to be determined at P r ≤ 10 −2 and has not been yet numerically reached despite reaching large Re values.This is to some extent expected since 3D effects become important when the thermal forcing increases (Calkins et al. 2013): the appearance of non-QG, equatorially anti-symmetric, axial flows that do not vary linearly along the rotation axis -breaking the underlying classical QG assumption -and the associated missing correlations between uz and T , are part of this discrepancy.We found these effects account for up to 46% of the missing convective power for a case at Ek = 10 −6 , P r = 1, Ra = 53 Rac.
The thick shell geometry studied and the omission of the dynamics inside the tangent cylinder, as well as the enforced linearity of uz, 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 3−5 times faster than a 3D 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 P r = 1 mainly because the z-interpolation scheme requires the resolution to be high.At P r < 1, the hybrid model becomes more advantageous; since the 3D grid involves only the temperature field it can be much coarser than the 2D 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/Rac ≤ 15) at all P r and we envision that the hybrid method could become even more attractive in terms of computational resources, even at P r = 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 characterised 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 3D (e.g., Mound & Davies 2017) and are evident in our 3D comparison calculations carried out in a similar regime.
When comparing the 3D, QG and hybrid models using inhomogneous thermal boundary conditions, we find that the QG and 3D configurations are qualitatively similar whereas the hybrid model again seems to lack buoyancy power.Despite the absence of smallscale 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 3D cases, suggesting that the dynamics is relatively well captured at mid-to-low latitudes.For Q * in the range from 2 − 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 ≤ 10 −7 , Re ≫ 1 and Ro ≪ 1, as reported in previous QG studies (e.g., Guervilly & Cardin 2017).When decreasing the ratio of diffusivities such that P r ≤ 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 P r < 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 ⊥ /d gradually approaches the asymptotic Rhines scaling l ⊥ ∼ 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 P r = 1 follow reasonably well the theoretical CIA scaling Rec ∼ [Ra(N u − 1)P r −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 3D 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 si, T2D(si) ∼ 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/Rac ≥ 10 − 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 3D magnetic field on the QG-convection.

APPENDIX A: NUMERICAL SIMULATIONS WITH FIXED TEMPERATURE CONTRAST
Table A1: Summary of the numerical simulations with fixed temperature contrast computed in this study.All models have been computed with η = ri/ro = 0.35.Ra is the Rayleigh number (the supercriticality Ra = • × Rac is also provided), P r is the Prandtl number, N u is the Nusselt number, Rec is the convective Reynolds number, Rezon is the zonal Reynolds number, Lu s is the typical length-scale for the cylindrical radial velocity field and (Ns, Nm)/(Nr, ℓmax) are the grid-size for the run.We compute the onset of rotating convection using the open source software SINGE † (Vidal & Schaeffer 2015) for the 3D 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 3D configuration, except for 3 cases at Ek = {10 −4 , 10 −5 , 10 −6 } and P r = 1, for which, we have determined the onset by time-integrating the nonlinear equations (11-14-17) using pizza with an initial sectorial temperature perturbation and by bracketing the Rayleigh number until the critical value is attained.Figure A1 displays the Rac values obtained for the 3 methods at P r = 1 and 10 −4 ≤ Ek ≤ 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, θ, ϕ) = F(r, θ) exp(λt + imϕ) , in the 3D configuration and of the form g(s, ϕ) = G(s) exp(λt + imϕ) , in the QG setup.Starting at a given Ra, the critical Rayleigh number Rac for a given azimuthal wavenumber mc is attained when ℜ(λ) = 0. Note that for the 3D configuration, it becomes numerically demanding to determine the onset of convection using a linear solver for Ek < 10 −7 for P r ≥ 1 and for Ek < 10 −6 ; P r < 1.For the cases with P r = 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 Rac ∼ Ek −4/3 is employed.Table A3 summarises the critical Rayleigh numbers Rac and azimuthal wavenumbers mc for the different setups.Ek −1 at P r = 1 with fixed temperature difference across the shell, for our QG (using LSB, in red), Hybrid (time-stepped with pizza for Ek ≤ 10 −6 , in green) and 3D (using SINGE for Ek ≤ 10 −7 , in cyan) models, and compared with an analytical solution in the Ek → 0 limit (Dormy et al. 2004, , in black).

Method
Table A3: Summary of the critical Rayleigh numbers Rac and critical azimuthal wavenumbers mc computed for our different setups.We have computed the onset of the hybrid method for 3 configurations but have otherwise assumed that hybrid QG-3D and purely 3D runs have the same Rac.

Figure 1 .
Figure 1.Nusselt number, N u 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 symbols and colors.The runs performed using the hybrid approach are represented with hatched symbols and the full 3D runs with empty symbols.

Figure 2 .
Figure 2. Nusselt number, N u as a function of the supercriticality, Ra/Rac.Various Ekman and Prandtl numbers displayed with respectively different symbols and colors.The runs performed using the hybrid approach are represented with hatched symbols and the full 3D runs with empty symbols.

Figure 4 .
Figure 4. Top panel: Comparison of the meridional section of the ϕ-averaged of the temperature field T 3D for the 3D (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 3D simulation (d), and equatorial section of the axial vorticity ω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 , P r = 1, and Ra = 2 × 10 9 .The spatial resolution in the 3D case is (Nr, ℓmax) = (129, 341); in the hybrid case is (Ns, Nm)/(Nr, ℓmax) = (257, 256)/(257, 256); and in the QG case is (Ns, Nm) = (385, 384) (bottom).For the three temperature and the three vorticity plots respectively, the same colorscales are used.Note that the colorscale for the vorticity is saturated to highlight the fine structure of the flows.

Figure 5 .
Figure 5. Top panel: comparison of the meridional section of the ϕ-averaged of the temperature field T 3D for the 3D (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 3D simulation (d), and equatorial section of the axial vorticity ω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 , P r = 1, and Ra = 1×10 10 (= 53.2 Rac for the 3D, = 61.3Rac for the Hybrid, and = 63.4Rac for the QG setup).The resolution in the 3D, hybrid and QG cases is respectively (Nr, ℓmax) = (321, 682), (Ns, Nm)/(Nr, ℓmax) = (513, 512)/(513, 341) and (Ns, Nm) = (577, 576).For the three temperature and the three vorticity plots respectively, the same colorscales are used.Note that the colorscale for the vorticity is saturated to highlight the fine structure of the flows.
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 ωz in four cases with top heat flux/bottomtemperature imposed boundary conditions: (a) Ek = 10 −6 , P r =

Figure 9 .
Figure 9. Top panel: Comparison of the temperature field at the CMB T 3D (ro) for (a) the 3D case and (b) the Hybrid case.Bottom panel: z-averaged vorticity for the 3D simulation (c), and equatorial section of the axial vorticity ω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 , P r = 1, and Ra Q = 8 × 10 9 imposing a ℓ = m = 2, Q * = 3 lateral variation around the imposed heat-flux condition.The resolutions in the 3D and hybrid cases are respectively (Nr, ℓmax) = (321, 682), and (Ns, Nm)/(Nr, ℓmax) = (385, 416)/(385, 416) and is (Ns, Nm) = (513, 512) in the QG case.For the three vorticity plots, the same colorscale is used and is saturated to highlight the fine structure of the flows.

Figure 10 .
Figure 10.Equatorial snapshots of axial vorticity, ωz, purely QG with heat-flux imposed at the top and a fixed temperature imposed at the bottom: Ek = 10 −6 , Ra Q = 4 × 10 9 , P r = 1 with no lateral variations (a); same parameters (Ek = 10 −6 , Ra Q = 4 × 10 9 , P r = 1) but with a m = 2, Q * = 3 lateral variation (b); same parameters (Ek = 10 −6 , Ra Q = 4 × 10 9 , P r = 1) but with a m = 1, Q * = 3 lateral variation and (c); and Ek = 10 −6 , Ra Q = 3.6 × 10 10 , P r = 1 but with a m = 2, Q * = 3 lateral variation (d).The cases (a-b-c) have the same colorscale displayed on the top right and the case (d) has its own colorscale displayed on the bottom right.Note that the colorscales are saturated to highlight the fine structure of the flows.The four cases are purely QG and their spatial resolution is (Ns, Nm) = (513, 512) in all cases.

Figure 13 .
Figure 13.Top: Typical length-scale of the radial velocity, Lu s , as a function of the Rossby number Rec E. A best-fit power law using only the runs with Ro ≤ 10 −3 is shown by the dot-dashed line.Bottom: Convective Reynolds number, Rec, as a function of Ra(N u − 1)P r −2 Ek 1/2 (Eq.(43).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 3D runs with empty symbols.The theoretically predicted scaling corresponding to the CIA balance (Eq.(43) is shown by the dot-dashed line.

Figure A1 .
Figure A1.Evolution of the critical Rayleigh number Rac as a function of the inverse Ekman number Ek −1 at P r = 1 with fixed temperature difference across the shell, for our QG (using LSB, in red), Hybrid (time-stepped with pizza for Ek ≤ 10 −6 , in green) and 3D (using SINGE for Ek ≤ 10 −7 , in cyan) models, and compared with an analytical solution in the Ek → 0 limit(Dormy et al. 2004, , in black).

Table A2 :
Continued on next page . . .Summary of the numerical simulations with inhomogeneous heat flux boundary conditions computed in this study.All models have been computed with η = ri/ro = 0.35.RaQ is the flux-based Rayleigh number, Y (m)/(ℓ,m) is the mode m (QG) or (ℓ, m) (Hybrid or 3D) of the imposed lateral flux variations (Y0 or Y0,0 indicates no lateral variations), Q * is the relative amplitude of the lateral flux variations, P r is the Prandtl number, N u∆ is the Nusselt number based on the temperature contrast in the shell, Rec is the convective Reynolds number, Rezon is the zonal Reynolds number, Lu s is the typical length-scale for the cylindrical radial velocity field and (Ns, Nm)/(Nr, ℓmax) are the grid-size for the run.Continued on next page . . .