Thomas Gastine

Keywords:

We present a new pseudo-spectral open-source code nicknamed pizza. It is dedicated to the study of rapidly-rotating Boussinesq convection under the 2-D spherical quasi-geostrophic approximation, a physical hypothesis that is appropriate to model the turbulent convection that develops in planetary interiors. The code uses a Fourier decomposition in the azimuthal direction and supports both a Chebyshev collocation method and a sparse Chebyshev integration formulation in the cylindrically-radial direction. It supports several temporal discretisation schemes encompassing multi-step time steppers as well as diagonally-implicit Runge-Kutta schemes. The code has been tested and validated by comparing weakly-nonlinear convection with the eigenmodes from a linear solver. The comparison of the two radial discretisation schemes has revealed the superiority of the Chebyshev integration method over the classical collocation approach both in terms of memory requirements and operation counts. The good parallelisation efficiency enables the computation of large problem sizes with O(10 4 × 10 4 ) grid points using several thousands of ranks. This allows the computation of numerical models in the turbulent regime of quasi-geostrophic convection characterised by large Reynolds Re and yet small Rossby numbers Ro. A preliminary result obtained for a strongly supercritical numerical model with a small Ekman number of 10 -9 and a Prandtl number of unity yields Re ≃ 10 5 and Ro ≃ 10 -4 . pizza is hence an efficient tool to study spherical quasi-geostrophic convection in a parameter regime inaccessible to current global 3-D spherical shell models.

pizza: an open-source pseudo-spectral code for spherical quasi-geostrophic convection

INTRODUCTION

Convection under rapid rotation is ubiquitous in astrophysical bodies. The liquid iron cores of terrestrial planets or the atmospheres of the gas giants are selected examples where turbulent convection is strongly influenced by rotational effects (e.g. [START_REF] Aurnou | Rotating convective turbulence in Earth and planetary cores[END_REF]. Such turbulent flows are characterised by very large Reynolds numbers Re > 10 8 and yet small Rossby numbers Ro < 10 -5 , Ro being defined as the ratio between the rotation period and the convective overturn time. This specific combination of Re ≫ 1 and Ro ≪ 1 corresponds to the so-called turbulent quasigeostrophic regime of rotating convection (e.g. [START_REF] Julien | Heat Transport in Low-Rossby-Number Rayleigh-Bénard Convection[END_REF][START_REF] Stellmach | Approaching the Asymptotic Regime of Rapidly Rotating Convection: Boundary Layers versus Interior Dynamics[END_REF]. This implies that, in absence of a magnetic field, the pressure gradients balance the Coriolis force at leading order. As a consequence, the convective flow shows a pronounced invariance along the axis of rotation. At onset of rotating convection for instance, the flow pattern takes the form of quasi-geostrophic elongated columnar structures that have a typical size of E 1/3 , where E = ν/Ωd 2 is the Ekman number with ν the kinematic viscosity, Ω the rotation frequency and d the thickness of the convective layer (e.g. [START_REF] Busse | Thermal instabilities in rapidly rotating systems[END_REF][START_REF] Dormy | The onset of thermal convection in rotating spherical shells[END_REF]. Convection in natural objects corresponds to extremely small Ekman numbers with for instance E ≃ 10 -15 in the Earth core or E ≃ 10 -18 in the gas giants. The quasi-geostrophy of the convective flow is expected to hold as long as the dynamics is dominated by rotation, or in other words as long as the buoyancy force remains relatively small compared to the Coriolis force [START_REF] Gilman | Nonlinear Dynamics of Boussinesq Convection in a Deep Rotating Spherical Shell[END_REF][START_REF] Julien | Heat Transport in Low-Rossby-Number Rayleigh-Bénard Convection[END_REF][START_REF] King | Scaling behaviour in Rayleigh-Bénard convection with and without rotation[END_REF][START_REF] Cheng | Laboratory-numerical models of rapidly rotating convection in planetary cores[END_REF][START_REF] Horn | Toroidal and poloidal energy in rotating Rayleigh-Bénard convection[END_REF][START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF].

Many laboratory experiments of rotating convection in spherical geometry have been carried out, either under micro-gravity conditions (e.g. [START_REF] Hart | Space-laboratory and numerical simulations of thermal convection in a rotating hemispherical shell with radial gravity[END_REF][START_REF] Egbers | The geoflow-experiment on ISS (part I): Experimental preparation and design of laboratory testing hardware[END_REF]; or on the ground using the centrifugal force as a surrogate of the radial distribution of buoyancy (e.g. [START_REF] Busse | Convection induced by centrifugal buoyancy[END_REF][START_REF] Sumita | Experiments on highly supercritical thermal convection in a rapidly rotating hemispherical shell[END_REF][START_REF] Shew | Liquid sodium model of geophysical core convection[END_REF]. Because of their limited size, those experiments could only reach E ≃ 5 × 10 -6 , far from the geophysical/astrophysical regime. In complement to the laboratory experiments, rotating convection in spherical geometry can also be studied by means of three-dimensional global numerical simulations. Because of computational limitations, those numerical models are currently limited to E 10 -7 , Re 10 4 and Ro 10 -3 , hardly scratching into the turbulent quasi-geostrophic (hereafter QG) regime [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF][START_REF] Schaeffer | Turbulent geodynamo simulations: a leap towards Earth's core[END_REF]). Reaching lower Ekman numbers is hence mandatory to further explore this regime with Re ≫ 1 and Ro ≪ 1.

A way to alleviate the computational constraints inherent in global 3-D computations is to consider a spherical QG approximation of the convective flow (e.g. [START_REF] Busse | Convection in a rotating cylindrical annulus -Thermal Rossby waves[END_REF][START_REF] Cardin | Chaotic thermal convection in a rapidly rotating spherical shell: consequences for flow in the outer core[END_REF][START_REF] Plaut | Low-Prandtl-number convection in a rotating cylindrical annulus[END_REF][START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF][START_REF] Morin | Time dependent β-convection in rapidly rotating spherical shells[END_REF][START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF][START_REF] Calkins | The influence of fluid properties on the morphology of core turbulence and the geomagnetic field[END_REF][START_REF] Teed | On the necessary conditions for bursts of convection within the rapidly rotating cylindrical annulus[END_REF][START_REF] Guervilly | Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores[END_REF][START_REF] More | Convectively driven decadal zonal accelerations in Earth's fluid core[END_REF] . The underlying assumption of the spherical QG approximation is that the leading-order cylindrically-radial and azimuthal velocity components are invariant along the axis of rotation z. Under this approximation, the variations of the axial vorticity along the rotation axis are also neglected and an averaging of the continuity equation along the rotation axis implies a linear dependence of the axial velocity on z (Schaeffer & Cardin 2005a;[START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF]. The spherical QG approximation hence restricts the computation of the evolution of the convective velocity to two dimensions only. This is a limitation compared to the 3-D QG convective models developed by [START_REF] Calkins | Three-dimensional quasigeostrophic convection in the rotating cylindrical annulus with steeply sloping endwalls[END_REF] which allow spatial modulations of the convective features along the rotation axis. Because of the radial distribution of the buoyancy forcing in spherical geometry, the temperature is not necessarily well-described by the quasi-geostrophic approximation. Spherical QG models with either a three-dimensional or a two-dimensional treatment of the temperature however yield very similar results [START_REF] Guervilly | Subcritical convection of liquid metals in a rotating sphere using a quasi-geostrophic model[END_REF]. Despite those approximations, the different implementations of the 2-D spherical QG models (e.g. [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF][START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF][START_REF] Calkins | The influence of fluid properties on the morphology of core turbulence and the geomagnetic field[END_REF][START_REF] Teed | On the necessary conditions for bursts of convection within the rapidly rotating cylindrical annulus[END_REF][START_REF] Guervilly | Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores[END_REF] have been found to compare favourably to 3-D direct numerical simulations in spherical geometry (e.g. [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF]Schaeffer & Cardin 2005a;[START_REF] Plaut | Reynolds stresses and mean fields generated by pure waves: applications to shear flows and convection in a rotating shell[END_REF]. This indicates that such 2-D spherical QG models could be efficiently used to explore the turbulent QG regime of convection with E < 10 -8 and Re 10 5 , a parameter regime currently inaccessible to 3-D computations. Quasi-geostrophy is expected to hold as long as the dynamics is dominated by rotation, or in other words as long as the buoyancy force remains relatively small compared to the Coriolis force [START_REF] Gilman | Nonlinear Dynamics of Boussinesq Convection in a Deep Rotating Spherical Shell[END_REF][START_REF] Julien | Heat Transport in Low-Rossby-Number Rayleigh-Bénard Convection[END_REF][START_REF] King | Scaling behaviour in Rayleigh-Bénard convection with and without rotation[END_REF][START_REF] Cheng | Laboratory-numerical models of rapidly rotating convection in planetary cores[END_REF][START_REF] Horn | Toroidal and poloidal energy in rotating Rayleigh-Bénard convection[END_REF][START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF].

The spatial discretisation strategy adopted in spherical QG models usually relies on a hybrid scheme with a truncated Fourier expansion in the azimuthal direction φ and second-order finite differences in the cylindrically-radial direction s (e.g. [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF][START_REF] Calkins | The influence of fluid properties on the morphology of core turbulence and the geomagnetic field[END_REF]). Note that [START_REF] Brummell | High Rayleigh number β-convection[END_REF] and [START_REF] Teed | On the necessary conditions for bursts of convection within the rapidly rotating cylindrical annulus[END_REF] rather employed a spectral Chebyshev collocation technique in s but in the case of a cartesian QG model. The vast majority of those numerical codes adopt a pseudo-spectral approach where the nonlinear terms are treated in the physical space and time-advanced with an explicit AdamsBashforth time scheme, while the linear terms are time-advanced in the Fourier space using a CrankNicolson scheme. In contrast to 3-D models where several codes with active on-going developments are freely accessible to the community (see [START_REF] Matsui | Performance benchmarks for a next generation numerical dynamo model[END_REF], there is a no open-source code for spherical QG convection available to the community.

The purpose of this study is precisely to introduce a new opensource pseudo-spectral spherical QG code, nicknamed pizza. pizza is available at https://github.com/magic-sph/pizza as a free software that can be used, modified, and redistributed under the terms of the GNU GPL v3 license. The package also comes with a suite of python classes to allow a full analysis of the outputs and diagnostics produced by the code during its execution. The code, written in Fortran, uses a Fourier decomposition in φ and either a Chebyshev collocation or a sparse Chebyshev integration method in s (e.g. [START_REF] Stellmach | An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers[END_REF][START_REF] Muite | A numerical comparison of Chebyshev methods for solving fourth order semilinear initial boundary value problems[END_REF][START_REF] Marti | A computationally efficient spectral method for modeling core dynamics[END_REF]. It supports a broad variety of implicit-explicit time schemes encompassing multi-step methods (e.g. [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF] and implicit Runge-Kutta schemes (e.g. [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF]. The parallelisation strategy relies on the Message Passing Interface (MPI) library.

The paper is organised as follows. Section 2 presents the equations for spherical QG convection. Section 3 and 4 are dedicated to the spatial and temporal discretisation schemes implemented in pizza. The parallelisation strategy is described in section 5. The code validation and several examples are discussed in section 6 before concluding in section 7.

A QUASI-GEOSTROPHIC MODEL OF CONVECTION

Because of the strong axial invariance of the flow under rapid rotation, the QG models approximate 3-D convection in spherical geometry by a 2-D fluid domain which corresponds to the equatorial plane of a spherical shell. Using the cylindrical coordinates (s, φ, z), the QG fluid domain hence corresponds to an annulus of inner radius si and outer radius so rotating against the z-axis with an angular frequency Ω. In the following, we adopt a dimensionless formulation of the spherical QG equations using the annulus gap d = sosi as a reference length scale and the viscous diffusion time d 2 /ν as the reference time scale. The temperature contrast ∆T between both boundaries defines the temperature scale. Gravity is assumed to grow linearly with the cylindrical radius s and is non-dimensionalised using its value at the external radius go.

The formulation of the QG model implemented in pizza is based on the spherical QG approximation introduced by [START_REF] Busse | Convection in a rotating cylindrical annulus -Thermal Rossby waves[END_REF] and further expanded by [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF] and [START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF] to include the effects of Ekman pumping. Following Schaeffer & Cardin (2005a) and [START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF] the axial velocity uz is assumed to vary linearly with z. Under this assumption, the Boussinesq continuity equation under the spherical QG approximation yields

1 s ∂(sus) ∂s + 1 s ∂u φ ∂φ + βus = 0 , (1)
where

β = 1 h dh ds = - s h 2 , (2)
and h = (s 2 os 2 ) 1/2 is half the height of the geostrophic cylinder at the cylindrical radius s. We adopt a vorticitystreamfunction formulation to fulfill the QG continuity equation (1). The cylindrically-radial and azimuthal velocity components are hence expanded as follows

us = 1 s ∂ψ ∂φ , u φ = u φ - ∂ψ ∂s -βψ, (3)
where the streamfunction ψ accounts for the non-axisymmetric motions, while u φ corresponds to the axisymmetric zonal flow component, the overbar denoting an azimuthal average. The axial vorticity ω is then expressed by

ω = 1 s ∂(su φ ) ∂s -L β ψ, (4)
where the operator L β is defined by

L β ψ = ∆ψ + 1 s ∂(βsψ) ∂s .
In the above equation, ∆ is the Laplacian operator in cylindrical coordinates. Under the QG approximation, the time evolution of the axial vorticity becomes

∂ω ∂t +∇•(u ω) = 2 E βus - Ra P r 1 so ∂ϑ ∂φ +F(E, u, ω)+∆ω , (5)
where ϑ denotes the temperature perturbation. The reader is referred to [START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF] for a comprehensive derivation of this equation. In the above equation, F(E, u, ω) corresponds to the Ekman-pumping contribution (Schaeffer & Cardin 2005a) to nonaxisymmetric motions expressed by

F(E, u, ω) = -Υ ω - β 2 u φ + β ∂ ∂φ - 5so 2h us . ( 6
)
where

Υ = so E 1/2 1 (s 2 o -s 2 ) 3/4 .
To ensure a correct force balance in the azimuthal direction, the axial vorticity equation ( 5) is supplemented by an equation dedicated to the axisymmetric motions [START_REF] Plaut | Low-Prandtl-number convection in a rotating cylindrical annulus[END_REF]. Taking a φ-average of the azimuthal component of the Navier-Stokes equations yields

∂u φ ∂t + usω = -Υ u φ + ∆u φ - u φ s 2 , (7)
where the first term in the right-hand-side corresponds to the Ekman-pumping contribution for the axisymmetric motions [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF]. The governing equations for the temperature perturbation under the QG approximation is given by

∂ϑ ∂t + ∇ • (u ϑ) + βusϑ + us dTc ds = 1 P r ∆ϑ , ( 8
)
where Tc is the conducting background state [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF][START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF]). In the case of a fixed-temperature contrast between si and so, Tc is given by

Tc = α ln η ln[(1 -η)s], dTc ds = α s ln η ,
where α is a constant coefficient that can be used to rescale the temperature contrast to get a better agreement with the z-average of the conducting temperature of a 3-D spherical shell [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF][START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF]. In the case of fixed temperature boundary conditions,

α = η 1 -η 1 (1 -η 2 ) 1/2 arcsinh (1 -η 2 ) 1/2 η -1 .
The dimensionless equations (4-8) are governed by the Ekman number E, the Rayleigh number Ra and the Prandtl number P r defined by

E = ν Ωd 2 , Ra = αT go∆T d 3 νκ , P r = ν κ , ( 9
)
where αT is the thermal expansion coefficient and κ is the thermal diffusivity.

We assume in the following no-slip and fixed temperature at both boundaries. This yields us = u φ = ϑ = 0 at s = si, so .

(10)

With the definition of the streamfunction (Eq. 3), this corresponds to

ψ = ∂ψ ∂s = ϑ = u φ = 0 at s = si, so . ( 11
)

SPATIAL DISCRETISATION

The unknowns us, u φ , ω and ϑ are expanded in truncated Fourier series in the azimuthal direction up to a maximum order Nm. For each field f = [us, u φ , ω, ϑ], one has

f (s, φ k , t) ≈ Nm m=-Nm fm(s, t) e imφ k ,
where φ k = 2π(k-1)/N φ with k = 1, ..., N φ defines N φ equallyspaced discrete azimuthal grid points. Since all the physical quantities are real, f * -m = fm, where the star denotes a complex conjugate. Complex to real Fast Fourier Transforms (FFTs) can hence be employed to transform each quantity from a spectral representation to a grid representation

f (s, φ k , t) = 2 Nm ′ m=0 ℜ fm(s, t) e imφ k , (12)
where the prime on the summation indicates that the m = 0 coefficient needs to be multiplied by one half. The inverse transforms are handled by real to complex FFTs defined by

fm(s, t) = 1 N φ N φ k=1 f (s, φ k , t) e -imφ k . (13)
Using N φ ≥ 3Nm prevents aliasing errors when treating the nonlinear terms [START_REF] Orszag | On the Elimination of Aliasing in Finite-Difference Schemes by Filtering High-Wavenumber Components[END_REF][START_REF] Boyd | Chebyshev and Fourier Spectral Methods, Second Revised Edition[END_REF]. This implies to discard the Fourier modes with Nm < m ≤ N φ when doing the direct FFT (12) and to pad with zeroes when computing the inverse transforms (13).

In the radial direction, the Fourier coefficients fm are further expanded in truncated Chebyshev series up to degree Nc -1

fm(s k , t) = C Nc-1 ′′ n=0 fmn(t) Tn(x k ) , (14)
where the hat symbols are employed in the following to denote the Chebyshev coefficients. The discrete Chebyshev transform from a spectral representation to a grid representation is given by

fmn(t) = C Nr ′′ k=1 fm(s k , t) Tn(x k ) . (15)
In the above equations C = [2/(Nr -1)] 1/2 is a normalisation factor and the double primes on the summations now indicate that both the first and the last indices are multiplied by one half. Tn(x k ) is the nth-order first-kind Chebyshev polynomial defined by

Tn(x k ) = T kn = cos[n arccos(x k )] = cos πn(k -1) Nr -1 ,
where

x k = cos π(k -1) Nr -1 , k = 1, ..., Nr,
is the kth-point of a Gauss-Lobatto grid with Nr collocation grid points. For an annulus of inner radius si and outer radius so, the Gauss-Lobatto interval that ranges from -1 to 1 is remapped to the interval [si, so] by the following affine mapping

s k = so -si 2 x k + so + si 2 , k = 1, ..., Nr .
The choice of using Gauss-Lobatto grid points also ensures that fast Discrete Cosine Transforms of first kind (DCTs) can be employed to compute the transforms between Chebyshev representation and radial grid space (14-15). pizza relies on the FFTW ⋆ library [START_REF] Frigo | The design and implementation of FFTW3[END_REF] for all the FFTs and DCTs. This ensure that each single spectral transform is computed in O(N ln N ) operations, where N = [Nr, Nm].

Spectral equations using Chebyshev collocation

Several approaches can be employed to approximate the solution of a differential equation using Chebyshev polynomials. The most straightforward choice when dealing with a set of non-constant partial differential equations such as Eqs. (4-8) is to resort to a Chebyshev collocation method (e.g. [START_REF] Canuto | Spectral methods. Fundamentals in Single Domains[END_REF]. In this kind of approach, the unknowns can be either the Chebyshev coefficients fn or the values of the approximate solution at the collocation points f (x k ). Both collocation techniques yield dense matrices with similar condition numbers [START_REF] Peyret | Spectral Methods for Incompressible Viscous Flow[END_REF]). The first one has been widely adopted by the astrophysical and geophysical communities after the seminal work by [START_REF] Glatzmaier | Numerical simulations of stellar convective dynamos. I -The model and method[END_REF].

Semi-discrete formulation

Expanding ω, ψ and ϑ in Fourier and Chebyshev modes yield the following set of coupled semi-discrete equations for the time evolution of ωm and ψm for the non-axisymmetric modes with m > 0

C Nc-1 ′′ n=0 d dt T kn -A C mkn ωmn(t) + B C mkn ψmn(t) = - Ra P r im so ϑm(s k , t) -Nω m (s k , t) C Nc-1 ′′ n=0 T kn ωmn(t) + C C mkn ψmn(t) = 0 , (16)
where the collocation matrices are expressed by

A C mkn =T ′′ kn + 1 s k T ′ kn - m 2 s 2 k + Υ k T kn , B C mkn = Υ k β k 2 T ′ kn + β k β k Υ k 2 + im s k imΥ k - 5soΥ k 2h k - 2 E T kn , C C mkn =T ′′ kn + β k + 1 s k T ′ kn - dβ k ds + β k s + m 2 s 2 k T kn ,
In the above equations, the superscripts C have been introduced to differentiate the collocation matrices from the forthcoming sparse formulation. For clarity, a given function f discretised at the collocation point x k is expressed as f k = f (x k ). T ′ kn and T ′′ kn are the first and second derivative of the nth-order Chebyshev polynomial at the collocation point x k . Nω m (s k , t) corresponds to the Fourier transform (13) of the advection terms that enters Eq. ( 5)

Nω m (s k , t) = 1 N φ N φ j=1 [∇ • (u ω)] e -imφ j .
⋆ http://fftw.org/

where N φ = 3Nm to ensure that the nonlinear terms are alias-free in φ [START_REF] Orszag | On the Elimination of Aliasing in Finite-Difference Schemes by Filtering High-Wavenumber Components[END_REF]. Instead of introducing the intermediate variable ω, we could rather have substituted its definition (4) into Eq. (5) to derive a single time-evolution equation that would depend on ψ only. This would imply to solve an equation of the form

∂ ∂t ∂ 2 ψ ∂s 2 + • • • = ∂ 4 ψ ∂s 4 + • • •
Though appealing this strategy is however not viable since this kind of time-dependent problem has been shown to be unconditionally unstable when using Chebyshev collocation discretisation [START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF][START_REF] Hollerbach | A spectral solution of the magneto-convection equations in spherical geometry[END_REF].

We proceed the same way to discretise the equations for the mean azimuthal flow u φ (7)

C Nc-1 ′′ n=0 d dt T kn -T ′′ kn - 1 s k T ′ kn + Υ k + 1 s 2 k T kn u φ 0n (t) = -Nu φ (s k , t), (17)
where the nonlinear term is expressed by

Nu φ (s k , t) = E 2 Υ k u φ 0 ω0 + 2 Nm 1 ℜ {us m ω * m } .
The first term in the right hand side corresponds to the selfinteraction of the zonal wind [START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF]. Finally, the spatial discretisation of the temperature equation ( 8) yields

C Nc-1 ′′ n=0 d dt T kn - 1 P r T ′′ kn + 1 s k T ′ kn - m 2 s 2 k T kn ϑmn(t) = im s k dTc ds ψm(s k , t) -N ϑm (s k , t) , (18)
where N ϑm (s k , t) corresponds to the FFT of the nonlinear terms that enter Eq. ( 8):

N ϑm (s k , t) = 1 N φ N φ j=1 [∇ • (uϑ) + β k usϑ] e -imφ j .

Boundary conditions

In the collocation method, equations ( 16), ( 17) and ( 18) are prescribed for the Nr -2 internal collocation grid points. The remaining boundary points s = si and s = so are used to impose the boundary conditions (11). This implies that the singularity of β and its derivatives at the outer boundary so is not necessarily an issue when using the collocation method since boundary conditions provide additional constraints there. When a given physical field f = [ψ, ω, ϑ, u φ ] is subject to Dirichlet boundary conditions at both boundaries, the following conditions on the Chebyshev coefficients fn should be fulfilled (e.g. Canuto et al. 2006, Eq. 3.3.19) Nc -1 ′′ n=0 fnm = 0, s = so;

Nc-1 ′′ n=0 (-1) n fnm = 0, s = si , ( 19
)
while for Neumann boundary conditions (e.g. Canuto et al. 2006, Eq. 3.3.23) Nc-1 ′′ n=0 n 2 fnm = 0, s = so;

Nc-1 ′′ n=0 (-1) n+1 n 2 fnm = 0, s = si .

(20) Independently of the subsequent details of the chosen implicitexplicit time scheme employed to time advance the QG equations, Eq. ( 16) forms a complex-type dense matrix operator of size (2Nr × 2Nr) for each Fourier mode m. Figure 1a shows the structure of the matrix that enters the left-hand-side of Eq. ( 16). The top Nr rows corresponds to the time-dependent vorticity equation ( 5), while the bottom Nr rows corresponds to the streamfunction equation ( 4). The four mechanical boundary conditions ( 11) are imposed on the first and last rows of the top-right and bottom-right quadrants of this matrix.

From a numerical implementation standpoint, Chebyshev polynomials at the collocation points T kn and their first and second derivatives T ′ kn and T ′′ kn form dense real matrices of dimensions (Nr ×Nr) that are precalculated and stored in the initialisation procedure of the code. In pizza, the discretised equations (16-18) supplemented by the boundary conditions ( 19) or (20) are solved using LAPACK † . The LU decomposition is handled by the routine dgetrf or its complex-arithmetic counterpart zgetrf and require O(N 3 r ) operations per Fourier mode m. This needs to be done at the initialisation stage of the code or at each iteration where a change in the time-step size occurs (see § 4). During each time step, the routines dgetrs (or zgetrs) are employed for the matrix solve and correspond to O(N 2 r ) operations per Fourier mode m. The amount of memory required to store the dense complex-type matrix that enters the left-hand-side of Eq. ( 16) grows as 64 N 2 r for one single azimuthal wavenumber m for a double-precision calculation. This corresponds to 1 Gigabyte of memory per Fourier mode for Nr = 4096 and hence makes the collocation approach extremely costly when Nr 10 3 .

Spectral equations using a Chebyshev integration method

To circumvent the limitations inherent in the collocation approach, several efficient Chebyshev spectral methods have been developed (e.g. [START_REF] Coutsias | An efficient spectral method for ordinary differential equations with rational function coefficients[END_REF][START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF][START_REF] Olver | A fast and well-conditioned spectral method[END_REF]. They all involve the solve of sparse matrices that are almost banded and can be inverted in O(p Nr) operations, p being the number of bands of the matrices. One approach, first introduced by [START_REF] Clenshaw | The numerical solution of linear differential equations in Chebyshev series[END_REF], consists of integrating q times a set of qth-order ordinary differential equations (ODEs) in Chebyshev space (see also [START_REF] Fox | Chebyshev polynomials in numerical analysis[END_REF][START_REF] Phillips | On the Coefficients of Integrated Expan-sions of Ultraspherical Polynomials[END_REF][START_REF] Greengard | Spectral integration and two-point boundary value problems[END_REF]. First limited to ODEs with constant coefficients, this method has been further extended by [START_REF] Coutsias | An efficient spectral method for ordinary differential equations with rational function coefficients[END_REF] to ODEs with rational function coefficients. The comparison of several Chebyshev methods for fourth-order ODEs carried out by [START_REF] Muite | A numerical comparison of Chebyshev methods for solving fourth order semilinear initial boundary value problems[END_REF] showed the advantages of such a Chebyshev integration method both in terms of matrix condition number and computational cost in the limit of large Nr. This technique has been successfully applied to the problem of rotating convection in both Cartesian [START_REF] Stellmach | An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers[END_REF] and spherical geometry [START_REF] Marti | A computationally efficient spectral method for modeling core dynamics[END_REF]. † http://www.netlib.org/lapack/ , Ra = 3 × 10 4 and P r = 1 and a CNAB2 time scheme with a fixed δt = 10 -4 . (a) corresponds to the collocation method (Eq. 16). T corresponds to the matrix with the coefficients T kn = Tn(x k ). (b) corresponds to the Chebyshev integration method with boundary conditions imposed as the first four tau lines (Eq. 25). (c) corresponds to Chebyshev integration method with boundary conditions enforced via a Galerkin formulation (Eq. 32). For the three panels, the matrix coefficients have been normalised by their maxima such that they share the same color axis. Zero entries are displayed in white.

Semi-discrete formulation

The Chebyshev integration methodology relies on the following indefinite integral identity (e.g. Canuto et al. 2006, Eq. 2.4.23)

Tn(x)dx = 1 2 Tn+1(x) n + 1 - Tn-1(x) n -1 for n > 1, (21)
which in its discrete form corresponds to the following sparse operator

I kn = - 1 2k δ k+1,n + 1 2k δ k-1,n for k > 1,
where δ corresponds to the Kronecker symbol. Identities for multiple integration can then be easily derived by recursive applications of Eq. ( 21).

Because of the singularity of β, we first need to regularise the set of equation (4-8) to make it suitable for a Chebyshev integration method. We hence adopt the following different definition for the streamfunction

Ψ us = 1 s ∂[ζ(s)Ψ ] ∂φ ; u φ = u φ - ∂[ζ(s)Ψ ] ∂s -βζ(s)Ψ . Using ζ(s) = h 2 = s 2 o -s 2 then yields us = h 2 s ∂Ψ ∂φ ; u φ = u φ -h 2 ∂Ψ ∂s + 3s Ψ . ( 22
)
From these definitions, one derives the following expression for the axial vorticity

ω ω = 1 s ∂(su φ ) ∂s -LI Ψ , (23)
where the operator LI is given by

LI Ψ = ∆ h 2 Ψ - 1 s ∂ ∂s s 2 Ψ .
The expansion of Ψ and ϑ in Fourier modes yields the following equation for the time evolution of Ψ for the non-axisymmetric Fourier modes

∂ ∂t -∆ LI - 2 E im Ψm = Ra P r im so ϑm+Nωm-Fǫ(E, Ψm) .
In the above equation, the classical Ekman pumping term (Eq. 6) has been replaced by the approximated form Fǫ defined by

Fǫ = Υǫ LI + s 2 ∂ ∂s - 3s 2 2h 2 ǫ + m 2 + 5imso 2hǫ Ψm (24)
where hǫ = [(so + ǫ) 2s 2 ] 1/2 corresponds to half the height of a geostrophic cylinder that would intersect a sphere with a slightly larger radius so + ǫ, with ǫ ≪ 1. Υǫ is defined accordingly by Υǫ = s

1/2 o /E 1/2 /h 3/2 ǫ
. This implies that Fǫ corresponds to the exact Ekman pumping contribution that would occur in a spherical QG set-up with an outer radius so + ǫ. In other words, the approximated Ekman pumping Fǫ tends to approach the exact contribution F in the limit of vanishing ǫ. This approximation is required when using a Chebyshev integration method to avoid the outer boundary singularity of the exact Ekman pumping term and to get a good spectral representation of this quantity once transformed to Chebyshev space. The error introduced by this approximation will be further assessed in § 6.

In addition, the Ekman pumping term requires special care since it comprises non-rational function coefficients. In contrast to the collocation method where it can be treated implicitly without any additional cost, this term shall hence be treated as yet another non-linear term since its implicit treatment would yield a dense operator with the Chebyshev integration method [START_REF] Hiegemann | Chebyshev matrix operator method for the solution of integrated forms of linear ordinary differential equations[END_REF]).

The equation for the time evolution of Ψ is regularised by a multiplication by s 4 and then integrated four times to yield

s 4 ∂ ∂t -∆ LI - 2 E im Ψm = αr 3 + βr 2 + γr + δ + s 4 Ra P r im so ϑm + Nωm -Fǫ(E, Ψm) , (25)
where α, β, γ and δ are constant of integration that will not be required once this equation has been supplemented by boundary conditions. At this stage, any single term that enters the above equation can be written as the product x q ∂ p f /∂x p , where p and q are positive integers. Following [START_REF] Marti | A computationally efficient spectral method for modeling core dynamics[END_REF], this equation is then integrated by parts until no differential operator remains, such that each term has the following form

4 p=0 • • • p× q x q f (x) dx p .
After expanding f (x) in Chebyshev polynomials using Eq. ( 14), the semi-discrete representation of Eq. ( 25) can be derived by multiple application of the recurrence relation ( 21). This yields

Nc-1 ′′ n=0 d dt A I mkn -B I mkn Ψmn(t) = Nc-1 ′′ n=0 C I kn Ra P r im so ϑmn(t) + Nωmn -Fǫ n(E, Ψm) , (26)
for k > 4. A I mkn , B I mkn , and C I kn are the discrete representations of the following operators

A I m = s 4 LI ; B I m = s 4 ∆LI + 2 E im ; C I = s 4
The internal matrix elements are determined using the freely available python package developed by Marti et al. (2016) ‡ that allows the symbolic computation of those operators § . Excluding boundary conditions, A I m , B I m and C I correspond to band matrices with pu super-diagonals and p ℓ sub-diagonals that have a bandwidth defined by

q = p ℓ + pu + 1 .
The bandwidth of A I m , B I m , and C I is 17, 13 and 17, respectively. We proceed the same way to establish the equations for the axisymmetric zonal flow component and for the temperature perturbation. Eq. ( 7) and Eq. ( 8) are multiplied by s 2 and integrated twice to yield

Nc-1 ′′ n=0 d dt D I kn -E I kn u φ 0n (t) = - Nc-1 ′′ n=0 D I kn Nu φ mn + Υǫu φ 0n , (27)
for the axisymmetric zonal flow component and

Nc -1 ′′ n=0 d dt D I kn - 1 P r F I kmn ϑmn(t) = - Nc-1 ′′ n=0 D kn im h 2 s dTc ds Ψm n + N ϑmn , (28)
for the temperature. Both equations are only valid for k > 2. D I kn , E I kn and F I mkn are the discrete representation of the following operators

D I = s 2 ; E I = s 2 -3 s; F I m = s 2 ∆ .
The bandwidth of D I , E I and F I m is 9, 5 and 5, respectively. In contrast to the semi-discrete equations obtained with the collocation approach, the right-hand-sides of Eq. (26-28) now involve nonlinear terms that are in Chebyshev space. To avoid aliasing errors, the Chebyshev coefficients of nonlinear terms that have n > 2Nr/3 are hence set to zero [START_REF] Orszag | On the Elimination of Aliasing in Finite-Difference Schemes by Filtering High-Wavenumber Components[END_REF].

Boundary conditions

At this stage, the system of equation (26-28) needs to be supplemented by boundary conditions. Given the definition of Ψ , the rigid mechanical boundary conditions that require the cancellation of us and u φ at both boundaries are already ensured by the three following identities: An extra boundary condition on Ψ is thus required. Following [START_REF] Bardsley | Could hydrodynamic Rossby waves explain the westward drift?[END_REF], we make the ansatz

Ψ (s = si) = ∂Ψ ∂s (s = si) = 0, Ψ (s = so) = 0 . ( 29
Ψ ∼ (s 2 o -s 2 ) n when s → so .
This yields the following expression for the viscous term

∆LI Ψ = 1 s 4 (s 2 o -s 2 ) n-3 8n s 8 o (-2n 3 + 3n 2 + 5n -6) ,
when s → so. A finite solution requires either n > 3 or the cancellation of the poynomial on n, which has four roots (-3/2, 0, 1, 2). n = -3/2 is not allowed and n = 0 is redundant with the cancellation of Ψ at s = so. Hence the first possible solution is n = 1 which yields

Ψ ∼ s 2 -s 2 o when s → so .
This corresponds to the following additional boundary condition

∂ 3 Ψ ∂s 3 = 0 for s = so . ( 30
)
When using the Chebyshev integration method, the boundary conditions can be either enforced via the tau-Lanczos method or by setting up an adapted Galerkin basis function [START_REF] Canuto | Spectral methods. Fundamentals in Single Domains[END_REF][START_REF] Boyd | Chebyshev and Fourier Spectral Methods, Second Revised Edition[END_REF]). In the tau-Lanczos formulation, the top rows of the matrices are used to enforce the boundary conditions, which are actually identical to the ones used in the collocation method . The fourth condition on Ψ given in Eq. ( 30) corresponds to the following last tau line (see [START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF])

Nc-1 ′′ n=0 n 2 (n 2 -1)(n 2 -4) Ψn = 0 . ( 31
)
Figure 1b shows the structure of the matrix that enters the left-handside of Eq. ( 26) when the boundary conditions are enforced using a tau-Lanczos formulation. The two first rows of the matrix correspond to the Dirichlet boundary conditions (Eqs. 19 and 29), the third one to the above equation and the fourth one to the Neumann boundary condition (Eqs. 20 and 29). Below those four full lines the matrix has a banded structure with 8 sub-and super-diagonals. This corresponds to a so-called bordered matrix wich can be inverted in O(17 Nr) operations as long as the number of full rows is small compared to the problem size (e.g. [START_REF] Boyd | Chebyshev and Fourier Spectral Methods, Second Revised Edition[END_REF]). Appendix A gives the details of the matrix inversion procedure as implemented in pizza.

We proceed the same way for the boundary conditions on the axisymmetric zonal flow and on the temperature. In those cases the Dirichlet boundary conditions ( 19) are imposed as the two first tau lines of the matrix, while the banded structure below is given by ( 27) and ( 28), respectively.

Alternatively, the boundary conditions can be imposed by introducing a suitable Galerkin basis. The underlying idea is to define basis functions that satisfy the boundary conditions such that the solutions expressed on this set of functions will also directly fulfill the boundary conditions. The Galerkin basis of functions φm is usually defined as a linear combination of a small number nc of Chebyshev polynomials

φn(x) = nc-1 i=0 γ n i Tn+i(x) .
We first construct the Galerkin basis for the four boundary conditions on Ψ (Eqs. 29 and 30). Following [START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF], the tau conditions (19, 20, 31) are used to establish a related Galerkin set. Appendix B gives the details of the calculation of the γ n i coefficients for 0 ≤ i ≤ 4. Ψ is then decomposed on the Galerkin basis as follows

Ψ (s) = Nr -5 n=0 Ψnφn(x) ,
where the tilda notation denotes the Galerkin coefficients. The Galerkin coefficients Ψ relate to the Chebyshev coefficients Ψ via

Ψ = SΨ Ψ ,
where SΨ is the stencil matrix that contains the coefficients γi. For the Galerkin basis employed for the equation on Ψ , SΨ is a band matrix with four sub-diagonals. The Galerkin formulation of Eq. ( 26) can be hence written in its matrix form as

P4 dA I m dt -B I m SΨ Ψm = P4C I Ra P r im so ϑm + Nωm -Fǫ , (32)
where P4 is an operator that removes the top four rows of the matrices, which correspond to the number of boundary conditions [START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF]. Figure 1c shows the structure of the matrix that enters the left-hand-side of Eq. ( 32). Compared to the bordered matrix obtained when using the tau method, the matrix has now a pure banded structure with an increased bandwidth with 8 sub-and 12 super-diagonals. Those matrices could be solved using standard band matrix solvers. In pizza, the LU decomposition is handled by the LAPACK routine dgbtrf or its complex arithmetic counterpart zgbtrf in O(q 2 Nr) operations per Fourier mode m. dgbtrs (or zgbtrs) routines are then employed for the matrix solve in O(q Nr) operations per Fourier mode m.

We proceed the same way for the zonal velocity and the temperature equations by defining a Galerkin basis that ensures Dirichlet boundary conditions at both boundaries. Several different Galerkin basis sets that satisfy this type of boundary conditions have been frequently used in the context of modelling rotating convection (e.g. [START_REF] Pino | Thermal and inertial modes of convection in a rapidly rotating annulus[END_REF][START_REF] Stellmach | An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers[END_REF]. Following [START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF], we decide here to adopt the following set φn(x) = Tn+2(x) -Tn(x), for n < Nr -3 .

(33)

In matrix form, the Galerkin formulations of equations ( 27) and (28) yield

P2 dD I dt -E I SD u φ 0 = -P2D I Nu φ m + Υǫu φ 0 , (34)
for the axisymmetric zonal flow component and

P2 dD I dt - F I m P r SD ϑm = -P2D I im h 2 s dTc ds Ψm + N ϑm , (35)
for the temperature, where SD is the stencil matrix (33) and P2 is an operator that removes the top two rows.

We note that different type of boundary conditions, such as stress-free and/or fixed flux thermal boundary conditions, would necessitate the derivation of dedicated Galerkin bases following a procedure similar to the one discussed in the appendix B.

Previous analysis by [START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF] showed that the Galerkin approach usually yield matrices with a better condition number than the bordered matrices obtained when using the tau-Lanczos method. This is particularly critical when 2-D or 3-D Chebyshev domains are considered but remains acceptable for 1-D problem as considered here (see Table 1 in Julien &Watson 2009). The Galerkin approach should hence be privileged as long as homogeneous boundary conditions are enforced, while inhomogeneous boundary conditions for which a Galerkin description becomes cumbersome are easier to handle with a tau-Lanczos formulation.

TEMPORAL DISCRETISATION

The equations discretised in space can be written as a general ordinary differential equation in time where the right-hand-side is split in two contributions

dy dt = E (y, t) + I(y, t), y(t0) = y0, (36)
where I(y, t) corresponds to the linear terms, while E (y, t) corresponds to the nonlinear advective terms. Temporal stability constraints coming from the linear terms that enter Eqs. (5-8) is usually more stringent that the one coming from the nonlinear terms.

Except for weakly nonlinear calculations, this precludes the usage of purely explicit time schemes such as the popular fourth order Runge-Kutta (e.g. [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF]. Although they offer an enhanced stability, purely implicit schemes are extremely costly since they involve the coupling of all Fourier modes due to the implicit treatment of the nonlinear terms. The potential gain in time step size is hence cancelled by the numerical cost associated with the solve of large matrices. In the following, we hence only consider implicit-explicit schemes (hereafter IMEX) to solve Eq. ( 36) and to produce the numerical approximation yn ≃ y(tn). We first consider the general k-step IMEX linear multistep scheme

yn+1 = k j=1 ajyn+1-j + δt k j=1 b E j En+1-j + k j=0 b I j In+1-j , (37)
where En+1-j = E (yn+1-j, tn+1-j ) and In+1-j = I(yn+1-j , tn+1-j ). The vectors a, b E and b I correspond to the weighting factors of the IMEX multistep scheme. For instance, the commonly-used second-order scheme assembled from the combination of a Crank-Nicolson for the implicit terms and a secondorder Adams-Bashforth for the explicit terms (hereafter CNAB2) corresponds to the following vectors a = (1, 0), b I = (1/2, 1/2) and b E = (3/2, -1/2) for a constant δt. In practice, Eq. ( 37) is rearranged as follows

(I -b I 0 δt I) yn+1 = k j=1 ajyn+1-j + δt k j=1 b E j En+1-j + b I j In+1-j , ( 38
)
where I is the identity matrix. In addition to CNAB2, pizza supports several semi-implicit backward differentiation schemes of second, third and fourth order that are known to have good stability properties (heareafter SBDF2, SBDF3 and SBDF4, see [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF][START_REF] Garcia | A comparison of high-order time integrators for thermal convection in rotating spherical shells[END_REF]. The interested reader is referred to the work by [START_REF] Wang | Variable step-size implicit-explicit linear multistep methods for time-dependent partial differential equations[END_REF] for the derivation of the vectors a, b I and b E when the time step size is variable. Table 1 summarises the main properties of the multistep schemes implemented in pizza.

Multistep schemes suffer from several possible limitations: (i) when the order is larger than two, they are not self-starting and hence require to be initiated with another lower-order starting scheme; (ii) limitations of the time step size to maintain stability is more severe for higher-order schemes (e.g. [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF][START_REF] Carpenter | Fourth-order Runge-Kutta schemes for fluid mechanics applications[END_REF]. In contrast, the multi-stage Runge-Kutta schemes are self-starting and frequently show a stability region that grows with the order of the scheme. To examine their efficiency in the context of spherical QG convection, we have also implemented in pizza several Additive Runge Kutta schemes. For this type of IMEX, we restrict ourself to the so-called Diagonally Implicit Runge Kutta schemes (hereafter DIRK) for which each sub-stage can be solved sequentially. For such schemes, the equation ( 36) is time-advanced from tn to tn+1 by solving ν sub-stages

I -a I ii δt I yi = yn + δt i-1 j=1 a E i,j Ej + a I i,j Ij , 1 ≤ i ≤ ν, (
b I 2 • • • b I ν ,
for the implicit terms, and

c E A E b E = 0 0 c E 2 a E 21 0 . . . . . . . . . . . . c E ν a E ν1 a E ν2 • • • 0 b E 1 b E 2 • • • b E ν ,
for the explicit terms, where zero values above the diagonal have been omitted. In the following, we only consider the stiffly accurate DIRK schemes for which the outcome of the last stage gives the end-result, without needing any assembly stage [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF]. This corresponds to b I j = a I νj and b E j = a E νj for 1 < j < ν. In addition, to minimise the memory storage which is particularly critical in the Chebyshev collocation approach, only the DIRK schemes that involve one single matrix storage in the implicit solve are retained, i.e. a I

ii is independent of i. The latter restriction corresponds to the so-called SDIRK (Singly Diagonally Implicit RungeKutta) schemes. In the following we discuss the convergence and the stability properties of two second order -ARS222 from [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF] and LZ232 from [START_REF] Liu | Some new additive Runge-Kutta methods and their applications[END_REF]-; and two third order SDIRK schemes -ARS443 from [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF] and BPR353 from [START_REF] Boscarino | Implicit-explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit[END_REF]-.

The nonlinear advection terms that enter Eqs. (4-7) are treated explicitly, while the dissipation terms and the vortex streching term in Eq. ( 5) are treated implicitly. As long as the fluid domain is entirely convecting, the buoyancy term that enters the vorticity equation ( 5) can either be treated explicitly or implicitly without a notable change of the stability properties of the IMEX Table 1. Time schemes implemented in pizza. The fifth (I) and the sixth columns (E) correspond to the number of implicit and explicit terms computed for one time step, respectively. The seventh column (Storage) is the number of state vectors that need to be stored to time-advance one physical quantity. The eighth column (Cost) corresponds to the elapsed wall time for one iteration normalised by the cost for one iteration of CNAB2. The last column contains the maximum CFL α obtained for a case with E = 10 -7 , Ra = 2 × 10 11 , P r = 1 and ǫ = 10 -3 which has been computed using the Chebyshev integration and Galerkin methods with (Nr, Nc, Nm) = (1025, 682, 1280). The asterisks corresponds to the models which have been run with an explicit treatment of the buoyancy term. [START_REF] Stellmach | An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers[END_REF]. We can expect more significant differences when some regions of the fluid are stably stratified.

Name

An implicit treatment of the buoyancy term only implies that the temperature equation ( 8) shall be first time-advanced to produce ϑ(tn+1) before time-advancing the vorticity and streamfunction (e.g. [START_REF] Glatzmaier | Numerical simulations of stellar convective dynamos. I -The model and method[END_REF]). The treatment of the Ekman pumping terms depends on the spatial discretisation strategy: while this can be treated implicitly without additional cost in the collocation method, this term has to be treated explicitly when using the Chebyshev integration method.

For an illustrative purpose, we give here the time-stepping equation for Ψm when the Chebyshev integration method (Eq. 26) is used in conjunction with an SDIRK time scheme (Eq. 39)

A I m -a I ii δt B I m Ψm(ti) = A I m Ψm(tn) + δt i-1 j=1 a I i,j B I m Ψm(tj) +δt i-1 j=1 a E i,j C I Ra P r im so ϑm(tj) + Nωm(tj) -Fǫ(tj) ,
where the buoyancy term has been treated explicitly and 1 ≤ i ≤ ν. This equation needs to be solved ν times per time step and the outcome of the final stage produces the time-advanced quantity Ψm(tn+1) for the azimuthal wavenumber m. A summary of the main properties of the SDIRK schemes implemented in pizza is also given in Table 1. Both families of time integrators ( 38) and ( 39) have a very similar structure and can hence be implemented using a shared framework, provided the programming language supports objectoriented implementation [START_REF] Vos | A generic framework for time-stepping partial differential equations (PDEs): general linear methods, object-oriented implementation and application to fluid problems[END_REF]. In pizza we rely on the object-oriented features provided by the Fortran 2003 norm to implement an abstract framework that allows easy switching between different schemes while minimising the number of code lines.

The different time steppers have been validated by running convergence tests. To do so, we consider a physical test problem with E = 3 × 10 -6 , Ra = 10 7 , P r = 0.025 and initiate the numerical experiment with a random temperature perturbation. We then run the numerical model using an SBDF4 time stepper until a statistically steady-state has been reached. This final state serves as the starting conditions of a suite of numerical simulations that use different fixed time step size δt between 10 -9 and 3 × 10 -6 over a fixed physical timespan t = 1.2 × 10 -3 . Following [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF], the error associated with the time stepper is defined as the sum of the relative errors on ϑ, us and u φ , where the relative error for one physical quantity f is expressed by

erel(f ) = (f -fref) 2 f 2 ref 1/2 .
In the above equation, the angular brackets correspond to an integration over the annulus

f = 2π 0 so s i f (s, φ) s ds dφ .
The fourth-order SBDF4 time stepper with the smallest time step size δt = 10 -9 has been used to define the reference solution fref.

Figure 2 shows the error as a function of δt for the time schemes given in Table 1 for both the collocation method (left panel) and the Chebyshev integration method with a Galerkin approach to enforce the boundary conditions (right panel). All schemes converge with their expected theoretical order until a plateau is reached around 3 × 10 -9 for the Chebyshev collocation and 10 -8 for the Chebyshev integration method. This can be attributed to the propagation of rounding errors that occur in the spectral transforms and in the calculation of the radial derivatives [START_REF] Sánchez | Newton-Krylov continuation of periodic orbits for Navier-Stokes flows[END_REF]). In other words, at this level of δt the error becomes dominated by the spatial discretisation errors. For a given order, SDIRK schemes are found to be more accurate than their multistep counterparts for the majority of the cases. This time scheme validation has been carried out with fixed time step sizes on a physical test case that is close to the onset of convection. To examine the efficiency of the different time schemes to model quasi-geostrophic turbulent convection, we also perform a stability analysis on a more turbulent setup. Indeed a precision of a fraction of a percent is usually sufficient when considering parameter studies of turbulent rotating convection (e.g. [START_REF] Gastine | Scaling regimes in spherical shell rotating convection[END_REF]. Hence, the determination of the largest time step size δt is of practical interest to assess the efficiency of a given time scheme. To do so, we consider a problem with E = 10 -7 , Ra = 2 × 10 11 and P r = 1, which is approximately 60 times supercritical. We first time-advance the solution until the nonlinear saturation has been reached using a CNAB2 time scheme. We then use the final state of this computation as the starting conditions of several numerical simulations that use different time schemes. Those simulations are computed over 3 × 10 -4 viscous time, which roughly corresponds to two turnover times. Since the advection terms are treated explicitly, the maximum eligible time step size must satisfy the following

Integration + Galerkin

Figure 2. Relative error as a function of the time step size δt for several multistep and SDIRK time schemes when using Chebyshev collocation (left panel) and Chebyshev integration method with boundary conditions enforced by a Galerkin approach (right panel). For comparison, the expected convergence orders have been denoted by black lines. Those convergence analyses have been carried out at the saturated stage of a numerical setup with E = 3 × 10 -6 , Ra = 10 7 and P r = 0.025. Spatial resolution corresponds to (Nr,Nc,Nm = 193,193,128) for the collocation method and (Nr, Nc, Nm = 193, 128, 128) for the integration method. ǫ = 3 × 10 -3 has been assumed for the Chebyshev integration method.

Courant criterion

δt ≤ α min max s,φ |us| δs -1 , max s,φ |u φ | s δφ -1 , (40)
where δs correspond to the local spacing of the Gauss-Lobatto grid and δφ = 2π/3Nm to the constant spacing in the azimuthal direction. In the above equation, α corresponds to the Courant-Friedrichs-Lewy number (hereafter CFL). To determine the CFL number of each time scheme, we compute series of simulations with different values of α and let the code runs with the maximum allowed δt that fulfills Eq. ( 40). This implies that δt will change at each iteration and hence that the matrices will be rebuilt at each time step. Since LU factorisation is very demanding when using Chebyshev collocation (O(N 3 r ) operations), we restrict the stability analysis to the sparse Chebyshev integration method with a Galerkin approach to enforce the boundary conditions. We use the time evolution of the total enstrophy ω 2 as a diagnostic to estimate the maximum CFL number α. Because of the clustering of the Gauss-Lobatto grid points, the time step size limitation usually occurs in the vicinity of the boundaries. Since ω 2 reaches its maximum value in the viscous boundary layers, any violation of Eq. ( 40) yields spurious spikes in the time evolution of the total enstrophy, well before the code actually crashes. For comparison, we define a reference solution that has been run with an SBDF4 time scheme with the smallest value of α = 0.05.

Figure 3a shows the time-averaged and the standard deviation of ω 2 as a function of α for the time schemes given in Table 1. The curves are comprised of two parts: one horizontal part where the time-averaged total enstrophy remains in close agreement with the reference case and the other featuring a rapid increase of both the time-averaged and the standard deviation of ω 2 . We hence define the largest acceptable α for a given time scheme as the value above which the time-averaged total enstrophy becomes more than 0.3% larger than the reference value. The rightmost column of Table 1 documents the obtained values. All multi-step schemes exhibit comparable CFL numbers with only a weak dependence on the theoretical order of the scheme. This is in agreement with the study by [START_REF] Carpenter | Fourth-order Runge-Kutta schemes for fluid mechanics applications[END_REF] who report comparable time step limitations for several SBDF schemes when the problem becomes numerically stiff. In contrast, the SDIRK schemes allow significantly larger CFL numbers with third-order schemes being more stable than the second-order ones. We quantify the efficiency of a time scheme by the ratio

σ = α cost , ( 41
)
where the cost corresponds to the average wall time of one iteration without LU factorisation (see the before last column in Table 1). Figure 3b shows a comparison of the relative efficiency of the time schemes compared to CNAB2. Although the CFL numbers are larger for the SDIRK schemes, they actually have a similar efficiency to multistep schemes due to their higher numerical cost. CNAB2 and ARS222 are found to be the most efficient second-order schemes, while BPR353 and SBDF3 are the best third-order schemes. The CFL numbers derived here are however only indicative since the stability of the schemes is expected to depend on the stiffness of the physical problem (e.g. [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF][START_REF] Carpenter | Fourth-order Runge-Kutta schemes for fluid mechanics applications[END_REF]. It is yet unclear whether the SDIRK schemes considered here will be able to compete with the multistep methods in the limit of turbulent quasi-geostrophic convection. Addressing this question would necessitate a systematic survey of the limits of stability of the time schemes over a broad range of Reynolds and Rossby numbers.

PARALLELISATION STRATEGY

The implementation of the algorithm presented before in pizza has been designed to run efficiently on massively-parallel architectures.

We rely on a message-passing communication framework based on the MPI (Message Passing Interface) standard. Several approaches have been considered to efficiently parallelise spectral transforms between physical and spectral space (e.g. [START_REF] Foster | Parallel algorithms for the spectral transform method[END_REF]).

Here we decide to resort to a transpose-based approach, such that all the spectral transforms are applied to data that are local to each processor. Whenever needed global transpositions of the data arrays are used to ensure that the dimension that needs to be transformed becomes local.

In pizza the data is distributed in two different configurations. In the first one, the radial level are distributed among MPI ranks while all azimuthal wavenumbers are local to each processor. This allows the computation of the 1D Fourier transforms (Eq. 12), the nonlinear terms in the physical space and the backward inverse transforms (Eq. 13). At this stage the data are rearranged in a second MPI configuration such that the wavenumbers m are distributed, while all radial levels are now in processor. Since each processor can possibly have a different amount of data to be sent to other processors, this parallel transposition is handled by the MPI variant routine mpi alltoallv that offers dedicated arguments to specify the amount of data to be sent and received from each partner. This configuration is used to time-advance the solution either via Chebyshev collocation or via Chebyshev integration method . This implies the solve of linear problems and possibly DCTs (Eq. 14) to transform the data from Chebyshev to radial space. Figure 4 summarises the data distribution used in pizza.

In the following, we examine the scalability performance of pizza using the occigen cluster ¶ . This cluster consists of more than 2000 computational nodes, each node being configured with two Intel 12 cores E5-2690V3 series processor with a clock frequency of 2.6 GHz. To build the executable, we make use of the Intel compiler version 17.0, Intel MPI version 5.1.3, Intel MKL version 17.0 for the linear solve and the matrix vector products and FFTW version 3.3.5 for Fourier and Chebyshev transforms. We first analyse the strong scaling performance of the code by running sequences of numerical simulations with several fixed problem size and an increasing number of MPI ranks. The left panels in Figure 5 show the wall time per iteration as a function of the number of cores for several problem sizes for both Chebyshev collocation and Chebyshev integration methods. The resolution (Nr, Nm) range from (97, 96) to (12289, 12288). Because of the dense complextype matrices of size (2Nr × 2Nr) involved in the time advance of the coupled vorticity-streamfunction equation ( 16), we cannot use the collocation method for the largest problem sizes since it already requires more than 1 GB per rank when Nr = 1537 and Nm = 1536 with 128 MPI ranks. For the spatial resolutions that are sufficiently small to be computed on one single node, we observe an improved performance when the code is running on one single processor (i.e. up to 12 cores) with the Chebyshev collocation. This is not observed in the sparse cases and hence might be attributed to an internal speed-up of the dense matrix solver of the Intel MKL library. Apart from this performance shift, both methods show a scalability performance that improves with the problem size. While the efficiency of the strong scalings are quickly degraded for Nranks > Nm/8 for small problem sizes, pizza shows a very good scalability up to Nranks = Nm/2 for the largest problem sizes. The scalability performance of the collocation method is usually better than the Chebyshev integration method for a given ¶ https://www.cines.fr/calcul/materiels/occigen to the models that use the Chebyshev collocation method, while panels (c) and (d) correspond to the models where the Chebyshev integration is were used in conjunction with a Galerkin approach to enforce the boundary condition. In panels (a) and (c) the dotted black lines correspond to the ideal scalings. The linear fits displayed in panels (b) and (d) have been computed from the cases with (Nr, Nm) = (1537, 1536). All the simulations have been computed using the CNAB2 time scheme.

problem size. This has to do with the larger amount of computational work spent in solving the dense matrices, which comparatively reduces the fraction of the wall time that corresponds to the MPI global transposes.

In complement to the strong scaling analyses, we also examine weak scaling performance tests. This consists of increasing the number of MPI ranks and the problem size accordingly, such that the amount of local data per rank stays constant. The spectral transforms implemented in pizza require O(NrNm ln Nm) operations for the FFTs (Eq. 12) and O(NmNr ln Nr) for the DCTs (Eq. 14). The solve of the linear problems involved in the time advance of the equations (4-8) grows like O(NmN 2 r ) for the collocation method and only O(NmNr) for the Chebyshev integration method. With the 1-D MPI domain decomposition discussed above, this implies that an increase of the spatial resolution while keeping a fixed amount of local data corresponds to an increase of the wall time that should scale with O(Nr) for the collocation method and with O(ln Nr) for the Chebyshev integration method. The right panels of

where the run time is expressed in seconds. For both methods, the normalised wall time per iteration is nearly proportional to the data volume per rank, indicating a good agreement with the expected theoretical scalings. We can make use of those scalings to estimate the minimum theoretical execution time as a function of the problem size. Based on the results of the strong scaling analyses, we assume that pizza shows a good parallel efficiency up to Nranks = Nm/2 when the collocation method is used and up to Nranks = Nm/4 when a sparse Chebyshev formulation is employed. This yields

min(t coll. run ) = 4.4 × 10 -8 N 1.98 r , min(t int. run ) = 1.2 × 10 -6 N 0.97 r ln Nr . (43)
Figure 6 shows a comparison between the actual minimum wall times for different spatial resolutions (see Fig. 5) and the above scalings. A good agreement is found for the sparse Chebyshev for- The lines correspond to the linear fits derived from the weak scaling tests (see Fig. 5b andd) for both radial discretisation strategies assuming N ranks = Nm/2 for the collocation method and N ranks = Nm/4 for the Chebyshev integration method combined with a Galerkin enforcement of boundary conditions. The symbols correspond to the minimum wall times obtained in the strong scaling analyses (Fig. 5a andc). mulation and for the collocation method with NrNm > 10 5 . Since the computational time of FFTs and DCTs still represents a significant fraction of one time step for small problem sizes, this is not surprising that the scaling given in Eq. ( 43) is only approached for sufficiently large problem sizes when the collocation method is employed.

Adopting a Chebyshev integration formulation for the radial scheme provides a significant speed up over the collocation approach, with for instance a factor 10 gain when NrNm ≃ 10 7 . Furthermore, while the collocation method becomes intractable for problem sizes with NrNm > 10 7 because of its intrinsic large memory prerequisite, the sparse formulation can be employed for spatial resolution larger than 10 4 ×10 4 . Global synchronisation and file lock contention can become an issue when reaching this range of problem sizes. In pizza this is remedied by collective calls to MPI-IO write operations to handle the outputting of checkpoints and snapshots.

CODE VALIDATION AND EXAMPLES

Weakly-nonlinear convection

In absence of a documented benchmark of spherical QG convection, we test the numerical implementation by first looking at the onset of convection. The underlying idea being to compare the results coming from a linear eigensolver with the results from pizza. The comparison of the different radial discretisation strategies is of particular interest to quantify the error introduced by the approximation of the Ekman pumping term involved in the sparse formulation (Eq. 24). To determine the onset of spherical QG convection, we linearise the system of equation (4-8) and seek for normal modes with

f (s, φ, t) = ℜ ∞ m=0 fm(s)e imφ+λt ,
where fm = (ψm, ϑm) T and λ = τ + iω d , τ being the growth rate and ω d the angular frequency. Since there is no coupling between the Fourier modes, we can seek for the solution fm of one individual azimuthal wavenumber. This forms the following generalised Table 2. Onset of convection for E = 10 -6 , P r = 0.025 and r i /ro = 4/11 from [START_REF] Gillet | Experimental and numerical studies of convection in a rapidly rotating spherical shell[END_REF] and obtained with the Linear Solver Builder package. Note that the critical Rayleigh number Rac from [START_REF] Gillet | Experimental and numerical studies of convection in a rapidly rotating spherical shell[END_REF]

λL β ψm = Ra P r im so ϑm - 2 E im s βψm -F(E, ψm) + ∆(L β ψm) , λϑm = ∆ϑm - im s dTc ds ψm , (44)
that is supplemented by the boundary conditions (11). We solve this generalised eigenvalue problem using the Linear Solver Builder package (hereafter LSB) developed by [START_REF] Valdettaro | Convergence and round-off errors in a two-dimensional eigenvalue problem using spectral methods and Arnoldi-Chebyshev algorithm[END_REF]. The linear operators that enter Eq. ( 44) are discretised on the Gauss-Lobatto grid using a Chebyshev collocation method in real space (e.g. [START_REF] Canuto | Spectral methods. Fundamentals in Single Domains[END_REF]). The entire spectrum of complex eigenvalues λ is first computed using the QZ algorithm [START_REF] Moler | An Algorithm for Generalized Matrix Eigenvalue Problems[END_REF]. One selected eigenvalue can then be used as a guess to accurately determine the closest eigenpair using the iterative Arnoldi-Chebyshev algorithm (e.g. [START_REF] Saad | Numerical methods for large eigenvalue problems[END_REF]). As indicated in Table 2, the linear solver has been tested and validated against published values of critical Rayleigh numbers for spherical QG convection with or without Ekman pumping [START_REF] Gillet | Experimental and numerical studies of convection in a rapidly rotating spherical shell[END_REF].

In the following we focus on weakly nonlinear QG convection with E = 3 × 10 -6 and P r = 0.025 and a radius ratio ri/ro = 0.35, a physical set up that is quite similar to the one considered by [START_REF] Gillet | Experimental and numerical studies of convection in a rapidly rotating spherical shell[END_REF] for liquid Gallium. Figure 7 shows the critical eigenmode (with τ ≃ 0) computed with LSB for these parameters. The onset of convection takes the form of a thermal Rossby wave that drifts in the retrograde direction with a critical azimuthal wavenumber m = 12, a drifting frequency ω d = -9.42690 × 10 3 and a critical Rayleigh number Rac = 9.55263×10 6 . The numerical convergence of this calculation has been assessed by computing the Chebyshev spectra of the different eigenfunctions as illustrated on Fig. 7c.

To validate the numerical implementation, the growth rate and the drift frequency obtained with pizza are compared to the eigenvalues derived with LSB. This requires a finite growth rate τ , hence we adopt in the following a marginally supercritical Rayleigh number Ra = 10 7 and compute the most critical eigenmodes for this Ra both in absence and in presence of Ekman pumping. The corresponding eigenmodes (ψ, ϑ) T computed with LSB are then used as starting conditions in pizza. A meaningful comparison necessitates that the nonlinear calculation remains in the weakly nonlinear regime. We hence restrict the computation to a short time interval of 10 -2 viscous time, which roughly corresponds to 15 periods of the most unstable drifting thermal Rossby wave. To ensure that the numerical error is dominated by the spatial discretisation rather than by the temporal one, we employ the BPR353 time scheme with a small time step size δt = 10 -7 (see Fig. 2). Figure 8 shows a comparison of the time evolution of the temperature fluc- tuation ℜ(ϑm=12) at mid depth using the linear eigenmode calculated with LSB and using the different radial discretisation schemes implemented in pizza. In absence of Ekman pumping (left panels), the different radial schemes yield almost indiscernible time evolution curves. The zoomed-in inset reveals a 6 significant digits agreement between the eigenmode and the weakly nonlinear calculations. When the Ekman pumping contribution is included (right panels), similar accuracy is recovered between the simulation computed with the collocation method and the eigenmode. The two nonlinear calculations that use the Chebyshev integration approach show a more pronounced deviation due to the approximated Ekman pumping term with ǫ = 3 × 10 -3 .

To determine the growth rate and the drift frequency in the nonlinear calculations, we fit the time evolution of ℜ(ϑm=12) at mid depth with the function a0 cos(ω d t + φ0)e τ t using least squares, the initial amplitude a0 and phase shift φ0 being determined by the starting conditions. Table 3 shows the obtained eigenpairs for the different radial schemes tested with several time integrators and values of ǫ. Overall the best agreement with the eigenvalues are obtained when the third-order BPR353 time scheme is employed. The superiority of the SDIRK scheme likely has to do with the lack of self-starting capabilities of multistep schemes, which hence require a lower-order starting time stepper to complete the first iterations. This procedure introduces errors larger than the theoretical order of the scheme that could account for the slightly larger inaccuracy of those schemes. The approximation of the Ekman pumping contribution when the Chebyshev integration method is used introduces an error that is more pronounced in the growth rate than in the drift frequency. This is expected since dissipation processes usually have a direct impact on the growth rate of an instability. A decrease of ǫ goes along with a proportional drop of the relative error on τ . This is however accompanied by an increase of the number of radial grid points in order to maintain the spectral convergence of the Ekman pumping term (24). This comparison validates the implementation of all the linear terms that enter Eqs. (4-8) for the different radial discretisation schemes. The approximation of the Ekman pumping contribution yields relative error that grow with ǫ. The collocation method should hence be privileged for small problem size. Because of its fastest execution time, the sparse Chebyshev formulation is the recommended approach when dealing with larger problem sizes. A large number of radial grid points indeed permits to accommodate small values of ǫ < 10 -3 , for which the error associated with the approximate Ekman pumping term becomes negligible.

Nonlinear convection

To pursue the code validation procedure, we now examine another physical setup which is not in the weakly nonlinear regime anymore with E = 10 -7 , P r = 1 and Ra = 2 × 10 11 , roughly 60 times the critical Rayleigh number. This corresponds to the setup that has been previously used to determine the Courant number of the different time schemes in § 4. To compare the different radial discretisation schemes, we first compute a simulation until a statistically steady-state has been reached. We then use this physical solution as a starting condition of several numerical simulations that use different radial discretisation schemes and two values of ǫ with the BPR353 time scheme. Since this is now a turbulent convection model, the time step size will change over time to satisfy the Courant condition (Eq. 40). To avoid the costly reconstruction of the matrices at each iteration, we adopt a time step size that is three quarter of the maximum eligible time step. The simulations are then computed over a timespan of roughly 0.03 viscous time, which corresponds to more than 150 turnover times.

Figure 9a shows the time evolution of the total and the zonal kinetic energy defined by

EK = 1 2 u 2 s + u 2 φ = EZ +2π Nm m=1 so s i |u m s | 2 + |u m φ | 2 s ds ,
where the zonal contribution is expressed by

EZ = 1 2 u φ 2 = π s 0 s i u φ 2 s ds .
The three numerical simulations feature a very similar time evolution with roughly 50% of the energy content in the axisymmetric azimuthal motions. They show a quasi-periodic behaviour with quick energy increases followed by slower relaxations. This can be attributed to the time evolution of the zonal jets that slowly drift towards the inner boundary where they become unstable [START_REF] Rotvig | Multiple zonal jets and drifting: Thermal convection in a rapidly rotating spherical shell compared to a quasigeostrophic model[END_REF]. Panels b and c of Fig. 9 show the time-average radial profiles and m spectra of the kinetic energy, respectively. A good agreement is found between the three radial discretisation schemes. Typical of 2-D QG turbulence, an inverse energy cascade with a m -5/3 slope takes place up to a typical lengthscale where the convective features are sheared apart by the zonal jets (here m ≃ 20, see Rhines , Ra = 10 7 and P r = 0.025. Zoomed-in insets highlight the differences between the eigenmode and the three spatial discretisation strategies implemented in pizza. The DNS have been time-advanced using the BPR353 time scheme with a fixed time step size δt = 10 -7 to ensure that the error of the time scheme is negligible (see Fig. 2). The simulations have been initiated with the most unstable m = 12 eigenmode calculated with LSB. Both sparse Chebyshev formulations assume ǫ = 3 × 10 -3 for the cases with Ekman pumping.

Table 3. Growth rate and drift frequency for the m = 12 mode for E = 3 × 10 -6 , Ra = 10 7 and P r = 0.025 with and without Ekman pumping. The first line has been computed with the linear solver LSB, while the others correspond to nonlinear calculations performed with pizza. All the simulations have been computed with a fixed time step size δt = 10 -7 . The setups highlighted with an asterisk have been time-advanced with an explicit treatment of the buoyancy term. The correct digits compared to the eigenmode are underlined for each solution. Table 4. Time-averaged diagnostics of three numerical simulations with E = 10 -7 , Ra = 2 × 10 11 and P r = 1. The simulations have been computed with the BPR353 time scheme. The fourth column corresponds to the average time step size. The fifth and sixth column contain the time-average and the standard deviation of the total and the zonal kinetic energy, respectively. The last column corresponds to the total number of core hours spent to compute the time interval displayed in Fig. 9. For a better quantification of the difference between the three radial schemes, Tab. 4 contains the time-average and the standard deviation of EK and EZ over the entire run time. Since dealiasing is also required in the radial direction when using a sparse Chebyshev formulation, the two cases that have been computed with the Chebyshev integration method have a larger number of radial grid points to ensure a number of Chebyshev modes comparable to the one used with the collocation method. Because of the change of the grid spacing (Eq. 40), this implies a decrease in the average timestep size. The time averages and standard deviation obtained for the three schemes and the two values of ǫ are found to agree within less than 1%. Given the unsteady nature of the solution, the differences in time step size and the limited time span considered for time averaging, it is not clear whether this difference can solely be attributed to the parametrisation of the Ekman pumping contribution. Notwithstanding this possible source of error, this comparison demonstrates that turbulent convection can be accurately modelled by an efficient sparse Chebyshev formulation with an acceptable error introduced by the Ekman pumping term approximation.

r scheme (Nr, Nc, Nm) ǫ δt E K ± σ(E K ) E Z ± σ(E Z )

Turbulent QG convection

To check the ability of the spectral radial discretisation schemes to model turbulent QG convection, we consider a third numerical configuration with E = 10 -9 , Ra = 1.5 × 10 14 and P r = 1. This corresponds to strongly supercritical convection (Ra > 100 Rac) at a very low Ekman number, a prerequisite to ensure that both large Reynolds and small Rossby numbers are reached at the same time. With the dimensionless units adopted in this study,

Re = 2EK π(s 2 o -s 2 i ) 1/2
, Ro = Re E .

For these control parameters, convection develops in the so-called turbulent QG regime (e.g. [START_REF] Julien | Heat Transport in Low-Rossby-Number Rayleigh-Bénard Convection[END_REF] with Re ≃ 10 5 and Ro ≃ 10 -4 . Numerical models that operate at these extreme parameters demand a large number of grid points -here (Nr, Nm) = (6145, 6144)which becomes intractable for the Chebyshev col-location method. We hence only compute this model using the Chebyshev integration method combined with a Galerkin approach to enforce the boundary conditions. For this physical configuration, a time integration of roughly ten convective overturns requires about 10 5 core hours.

Figure 10 shows a snapshot of the vorticity with two zoomedin insets that emphasise the regions close the boundaries. The mixing of the potential vorticity (ω + 2/E)/h by turbulent convective motions generates multiple zonal jets with alternated directions (e.g. [START_REF] Dritschel | Multiple jets as PV staircases: the Phillips effect and the resilience of eddy-transport barriers[END_REF]. This gives rise to a spatial separation of the vortical structures with alternated concentric rings of cyclonic (ω > 0) and anticyclonic (ω < 0) vorticity. The typical size of these zonal jets is usually well-predicted by the Rhines scale defined by (Ro/|β|) 1/2 (e.g. [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF][START_REF] Gastine | Zonal flow scaling in rapidlyrotating compressible convection[END_REF][START_REF] Verhoeven | The compressional beta effect: A source of zonal winds in planets?[END_REF][START_REF] Heimpel | Simulation of deep-seated zonal jets and shallow vortices in gas giant atmospheres[END_REF][START_REF] Guervilly | Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores[END_REF]. This lengthscale marks the separation between Rossby waves at larger scales and turbulent motions at smaller scales. Because of the increase of |β| with the cylindrical radius s in spherical geometry, the zonal jets are getting thinner outward. Close to the outer boundary, the dynamics becomes dominated by tilted vortices elongated in the azimuthal direction, a typical pattern of the propagation of thermal Rossby waves. Because of the steepening of β at large radii, the vortex stretching term becomes the dominant source of vorticity there, such that the propagation of thermal Rossby waves takes over the nonlinear advective processes. This outer region is hence expected to shrink with an increase of the convective forcing (e.g. [START_REF] Guervilly | Multiple zonal jets and convective heat transport barriers in a quasi-geostrophic model of planetary cores[END_REF]. At the interface between jets, the vortical structures are sheared apart into elongated filaments, indicating a direct cascade of enstrophy towards smaller scales.

CONCLUSION

In this study, we have presented a new open-source code, nicknamed pizza, dedicated to the study of rapidly-rotating convection under the 2-D spherical quasi-geostrophic approximation (e.g. [START_REF] Busse | Convection in a rotating cylindrical annulus -Thermal Rossby waves[END_REF][START_REF] Aubert | Quasigeostrophic models of convection in rotating spherical shells[END_REF][START_REF] Gillet | The quasi-geostrophic model for rapidly rotating spherical convection outside the tangent cylinder[END_REF]. The code is available at https://github.com/magic-sph/pizza as a free software that can be used, modified, and redistributed under the terms of the GNU GPL v3 license. The radial discretisa-Figure 10. Snapshot of the axial vorticity for a numerical simulation with E = 10 -9 , Ra = 1.5 × 10 14 and P r = 1. The Chebyshev integration formulation with a Galerkin approach to enforce the boundary condition has been employed to compute this numerical model. We use a spatial resolution (Nr, Nm) = (6145, 6144) and ǫ = 10 -4 for the approximated Ekman pumping term. For a better visualisation, only one quarter of the solution is displayed. tion relies on a decomposition in Fourier series in the azimuthal direction and in Chebyshev polynomials in the radial direction. For the latter, both a classical Chebyshev collocation method (e.g. [START_REF] Glatzmaier | Numerical simulations of stellar convective dynamos. I -The model and method[END_REF][START_REF] Boyd | Chebyshev and Fourier Spectral Methods, Second Revised Edition[END_REF]) and a sparse integration method (e.g. [START_REF] Stellmach | An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers[END_REF][START_REF] Muite | A numerical comparison of Chebyshev methods for solving fourth order semilinear initial boundary value problems[END_REF][START_REF] Marti | A computationally efficient spectral method for modeling core dynamics[END_REF]) are supported. We adopt a pseudo-spectral approach where the nonlinear advective terms are treated in the physical space and transformed to the spectral space using fast discrete Fourier and Chebyshev transforms. pizza supports several implicit-explicit time schemes encompassing multi-step schemes as well as diagonally-implicit Runge-Kutta schemes (e.g. [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF]) that have been validated by convergence tests. The parallelisation strategy relies on a message-passing communication framework based on the MPI standard. The code has been tested and validated against onset of quasi-geostrophic convection.

The comparison of the two radial discretisation schemes has revealed the superiority of the Chebyshev integration method. In contrast to the collocation technique that requires the storage and the inversion of dense matrices, the integration method indeed only involves sparse operators. As a consequence, the memory requirements only grows with O(Nr) and the operation count with O(Nr ln Nr) as compared to O(N 2 r ) when using a collocation approach. Multi-step and diagonally-implicit Runge-Kutta schemes have shown comparable efficiency, defined in this study by the ratio of the maximum CFL number over the numerical cost of one iteration. Additional parameter studies with various Reynolds and Rossby numbers are however required to assess the differences between both families of time integrators. We have found a good parallel scaling up to roughly four radial grid points per MPI task. This implies that large spatial resolution up to O(10 4 × 10 4 ) grid points can be reached with a reasonable wall time if one uses several thousands of MPI tasks. Such large grid resolutions allows the study of turbulent quasi-geostrophic convection at low Ekman numbers. Preliminary results for a numerical model with E = 10 -9 , Ra = 1.5 × 10 14 and P r = 1 shows the formation of multiple zonal jets, when both the Reynolds number is large O(10 5 ) and the Rossby number is small O(10 -4 ). This specific combination of Re ≫ 1 and Ro ≪ 1 is a prerequisite to study the turbulent quasi-geostrophic convection regime [START_REF] Julien | Heat Transport in Low-Rossby-Number Rayleigh-Bénard Convection[END_REF], an important milestone to better understand the internal dynamics of planetary interiors.

Future developments of the code include the implementation of the time-evolution of chemical composition to study doublediffusive convection under the spherical QG framework. On the longer term, the QG flow and temperature computed in the equatorial plane of the spherical shell will be coupled to an induction equation computed in the entire shell using classical 3-D pseudospectral discretisation (e.g. [START_REF] Schaeffer | Quasi-geostrophic kinematic dynamos at low magnetic Prandtl number[END_REF]).

Figure 1 .

 1 Figure1. Representation of the coefficients of the left-hand-side matrices obtained for m = 4 for a setup with E = 10 -3 , Ra = 3 × 10 4 and P r = 1 and a CNAB2 time scheme with a fixed δt = 10 -4 . (a) corresponds to the collocation method (Eq. 16). T corresponds to the matrix with the coefficients T kn = Tn(x k ). (b) corresponds to the Chebyshev integration method with boundary conditions imposed as the first four tau lines (Eq. 25). (c) corresponds to Chebyshev integration method with boundary conditions enforced via a Galerkin formulation (Eq. 32). For the three panels, the matrix coefficients have been normalised by their maxima such that they share the same color axis. Zero entries are displayed in white.

  ) ‡ It can be downloaded as part of the supplementary materials of the study by[START_REF] Marti | A computationally efficient spectral method for modeling core dynamics[END_REF] here. § https://www.sympy.org/

  39)where yi is the intermediate solution at the stage i. Finally the evaluation ofyn+1 = yn + δt ν j=1 b E j Ej + b I j Ij .allows the determination of yn+1. A DIRK scheme with ν stages can be represented in terms of the following so-called Butcher tables c I A I

Figure 3 .Figure 4 .

 34 Figure3. (a) Time-averaged total enstrophy ω 2 as a function of the CFL coefficient α for the time schemes given in Table1. The error bars correspond to one standard deviation and the horizontal dashed line corresponds to the time-averaged enstrophy of a numerical model computed with the SBDF4 time scheme with α = 0.05. (b) Efficiency σ of the time schemes normalised by the efficiency of CNAB2. The vertical dotted line corresponds to a value of one. This efficiency analysis has been carried out at the saturated stage of a numerical model with E = 10 -7 , P r = 1, Ra = 2 × 10 11 . This model was computed using the Chebyshev integration and a Galerkin method with a spatial resolution of (Nr, Nc, Nm = 1025, 682, 1280) and ǫ = 10 -3 .

ttFigure 5 .

 5 Figure 5. Left panels: wall time per iteration as a function of the number of MPI ranks (strong scaling performance) for several spatial resolutions. Right panels: wall time per iteration as a function of the local data volume per MPI task NrNm/N ranks (weal scaling performance). Panels (a) and (b) correspondsto the models that use the Chebyshev collocation method, while panels (c) and (d) correspond to the models where the Chebyshev integration is were used in conjunction with a Galerkin approach to enforce the boundary condition. In panels (a) and (c) the dotted black lines correspond to the ideal scalings. The linear fits displayed in panels (b) and (d) have been computed from the cases with (Nr, Nm) = (1537, 1536). All the simulations have been computed using the CNAB2 time scheme.

  Fig 5 show the wall time per iteration normalised by those theoretical predictions as a function of the data volume per rank expressed by NrNm/Nranks for both Chebyshev methods. Using the simulations with a spatial resolution of (Nr, Nm) = (1537, 1536) we compute the following best fits between the normalised execu-tion time and the local data volume for each radial discretisation scheme

Figure 6 .

 6 Figure 6. Minimum wall time per iteration as a function of the problem size NrNm.The lines correspond to the linear fits derived from the weak scaling tests (see Fig.5b and d) for both radial discretisation strategies assuming N ranks = Nm/2 for the collocation method and N ranks = Nm/4 for the Chebyshev integration method combined with a Galerkin enforcement of boundary conditions. The symbols correspond to the minimum wall times obtained in the strong scaling analyses (Fig.5a and c).

Figure 7 .

 7 Figure 7. Eigenfunction of the first unstable mode for E = 3 × 10 -6 and P r = 0.025. This mode has a critical Rayleigh number Rac = 9.55263 × 10 6 , a drift frequency ω d = -9.42690 × 10 3 and an azimuthal wavenumber m = 12. (a) Temperature fluctuation ϑ in the equatorial plane. (b) Streamfunction ψ in the equatorial plane. (c) Normalised Chebyshev spectra of the eigenfunction.

Figure 9 .

 9 Figure 9. (a) Total E K and zonal E Z kinetic energy as a function of time for three numerical simulations with different radial discretisation schemes. (b) Time and azimuthally averaged kinetic energy as a function of radius. (c) Time-averaged kinetic energy spectra as a function of the wavenumber m. The shaded region correspond to one standard deviation of temporal fluctuations relative to the time averages. The simulations assume the following control parameters: E = 10 -7 , Ra = 2 × 10 11 and P r = 1. The sparse cases have been computed with ǫ = 10 -3 .

  have been normalised to match our definition.

		Rac	m	ω d
		Without Ekman pumping
	LSB	1.3851 × 10 7	13	-1.3028 × 10 4
	Gillet et al. (2007)	1.39 × 10 7	13	-1.300 × 10 4
		With Ekman pumping
	LSB	1.5231 × 10 7	14	-1.2705 × 10 4
	Gillet et al. (2007)	1.53 × 10 7	14	-1.268 × 10 4
	eigenvalue problem

  Real part of ϑ m=12 at mid-depth sm = 0.5(s i + so) as a function of time without Ekman pumping (left panel) and with Ekman pumping (right panel) for a case with E = 3 × 10 -6

	(a)	6	×10 -9	Without Ekman pumping DNS (collocation) DNS (integ.+tau)		(b) 1.0	×10 -10 1.00 ×10 -10 With Ekman pumping
			4		DNS (integ.+Gal.)					0.98
	m )}		2		Eigenmode				0.5	0.00880 0.00885 0.00890 0.96
	ℜ{ϑ m=12 (s = s	-2 0		-5.96310	×10 -9				0.0 -0.5		-1.1987	×10 -10
		-4
					-5.96315					-1.0		-1.1988
		-6			0.009747 0.009748						0.009848 0.009850
			0.000	0.002	0.004	0.006	0.008	0.010	0.000	0.002	0.004	0.006	0.008	0.010
						Time							Time
	Figure 8.

ACKNOWLEDGMENTS

I want to thank Alexandre Fournier for his comments that helped to improve the manuscript. Stephan Stellmach and Benjamin Miquel are acknowledged for their fruitful advices about Galerkin bases and Philippe Marti for his help with the symbolic python package used to assemble the sparse Chebyshev matrices. I also wish to thank Michel Rieutord for sharing the Linear Solver Builder eigensolver. Numerical computations have been carried out on the S-CAPAD platform at IPGP and on the occigen cluster at GENCI-CINES (Grant A0020410095). All the figures have been generated using matplotlib (Hunter 2007). All the post-processing tools that have been used to construct the different figures are part of the source code of pizza and are hence freely accessible. This is IPGP contribution 4015.

APPENDIX A: DIRECT SOLVE OF A BORDERED MATRIX

Suppose one wants to solve the following linear problem which involves a so-called bordered matrix A Aψ = f, where A comprises p full top rows and a banded structure underneath. The matrix problem is sub-divided as follows

where A is a full square matrix of size (p × p), B is a full matrix of size (p × np), C is a sparse matrix of size (np × p) and D is a band matrix of size (np) with a bandwidth q, q being the total number of bands. One first solves the two following banded linear problems

The LU factorisation of the band matrix D requires O(q 2 n) operations, while the solve requires O(q n) operations (e.g. Boyd 2001, Appendix B2). We then assemble the Schur complement of the banded block

before solving the small dense problem of size (p, p)

This requires O(p 3 ) operations for the LU factorisation and O(p 2 ) for the solve. This cost remains negligible as long as p ≪ n, which is the case for the linear problems considered in the Chebyshev integration method. We finally evaluate ψ2 = xy ψ1 , to assemble the final solution given by ψ = (ψ1, ψ2) T .

APPENDIX B: GALERKIN BASIS FOR STREAMFUNCTION BOUNDARY CONDITIONS

In this section, we derive a Galerkin basis function for the following combination of boundary conditions that is used in the Chebyshev integration method for the streamfunction equation

We start by defining the following ansatz for the Galerkin set

Following [START_REF] Mcfadden | Elimination of spurious eigenvalues in the Chebyshev Tau spectral method[END_REF] and [START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF] we then make use of the tau boundary conditions (Eqs. 19,20 and 31) to form the following system of equations

Since there are only four equations for five unknowns, there is a degree of freedom in the determination of the coefficients. We thus choose in following γ n 0 = 1, which yields the following identities for the other coefficients:

, and γ n 4 =

(n + 1) 2n 4 + 12n 3 + 30n 2 + 36n + 15 (n + 3) (2n 4 + 20n 3 + 78n 2 + 140n + 95) .

This paper has been produced using the Blackwell Scientific Publications GJI L A T E X2e class file.