Venkatesh Gopinath

Alexandre Fournier
email: fournier@ipgp.fr

Thomas Gastine

An assessment of implicit-explicit time integrators for the pseudo-spectral approximation of Boussinesq thermal convection in an annulus

Keywords: Boussinesq convection, pseudo-spectral methods, stiff ODE/PDE/DAE time marching, IMEX time integrators, stability, turbulence

We analyze the behaviour of an ensemble of time integrators applied to the semi-discrete problem resulting from the spectral discretization of the equations describing Boussinesq thermal convection in a cylindrical annulus.

The equations are cast in their vorticity-streamfunction formulation that yields a differential algebraic equation (DAE). The ensemble comprises 28 members: 4 implicit-explicit multistep schemes, 22 implicit-explicit Runge-Kutta (IMEX-RK) schemes, and 2 fully explicit schemes used for reference. The schemes whose theoretical order varies from 2 to 5 are assessed for 11 different physical setups that cover laminar and turbulent regimes.

Multistep and order 2 IMEX-RK methods exhibit their expected order of convergence under all circumstances. IMEX-RK methods of higher-order show occasional order reduction that impacts both algebraic and differential field variables. We ascribe the order reduction to the stiffness of the problem at hand and, to a larger extent, the presence of the DAE. Using the popular Crank-Nicolson Adams-Bashforth of order 2 (CNAB2) integrator as reference, performance is defined by the ratio of maximum admissible time step to the cost of performing one iteration; the maximum admissible time step is determined by inspection of the time series of viscous dissipation within the system, which guarantees a physically acceptable solution. Relative performance is bounded between 0.5 and 1.5 across all studied configurations. Considering accuracy jointly with performance, we find that 6 schemes consistently outperform CNAB2, meaning that in addition to allowing for a more efficient calculation, the accuracy that they achieve at their operational, dissipation-based limit of stability yields a lower error. In our most turbulent setup, where the behaviour of the methods is almost entirely dictated by their explicit component, 13 IMEX-RK integrators outperform CNAB2 in terms of accuracy and efficiency.

Introduction

This study is concerned with the evaluation of implicit-explicit (IMEX) time marching methods applied to the numerical simulation of thermal convection for geophysical or astrophysical bodies. Thermal convection is an ubiquitous process in natural systems of large size; it drives the internal and external evolution of planets and stars as they receive and shed heat from and to outer space. In the case of Earth's internal dynamics, two envelopes undergo convection: the rocky mantle and the liquid outer core underneath it, which is essentially composed of an Iron-Nickel alloy. Both systems host vigorous time-dependent convective currents, over vastly different time scales, since the mantle turnover time, O(10 8 ) years, is a million times larger than that of the core, O(10 2 ) years.

Numerical models of mantle convection appeared in the nineteen-seventies (e.g. [START_REF] Mckenzie | Convection in the Earth's mantle: towards a numerical simulation[END_REF][START_REF] Sato | Finite element models for creeping convection[END_REF][START_REF] Kopitzke | Finite element convection models: comparison of shallow and deep mantle convection, and temperatures in the mantle[END_REF][START_REF] Jarvis | Time-dependent convection in the Earth's mantle[END_REF] and have grown steadily in size and complexity since (see e.g. [START_REF] Zhong | Numerical methods for mantle convection[END_REF], for a review). Efforts have been carried out in view of handling both complex geometries and rheologies, in order for instance to simulate plate tectonics in an inertialess framework that necessitates the solve of a modified Stokes problem for the flow. This led to the development of multilevel elliptic solvers and the implementation of adaptive mesh refinement techniques. The design of high-order integration schemes has logically not been the focus of attention, since the priority was to design methods able in particular to cope with viscosity contrasts spanning several decades. Two freely available codes whose development continues today and used to simulate mantle convection in two or three dimensions are fluidity [START_REF] Davies | Fluidity: A fully unstructured anisotropic adaptive mesh computational modeling framework for geodynamics[END_REF] and ASPECT [START_REF] Kronbichler | High accuracy mantle convection simulation through modern numerical methods[END_REF]. With regard to the advection-diffusion of temperature anomalies, fluidity resorts to an implicit θ-method (Canuto et al., 2006, §D.2.3) for the diffusive and advective terms, and has θ = 1/2, which corresponds to the second-order Crank-Nicolson (CN henceforth) method. The ASPECT code opted for a compromise between stability and accuracy in choosing a backward-difference formula of second order (BDF2, e.g., Canuto et al., 2006, §D.2.4). BDF2 was also praised by these authors on the account of its "efficiency of implementation (higher-order schemes often become unwieldy as they require complicated initialization during the first few time steps, and require the storage of many solution vectors from previous time steps)" [START_REF] Kronbichler | High accuracy mantle convection simulation through modern numerical methods[END_REF]. This statement does not consider high-order self-restart time integration methods of the implicit-explicit Runge -Kutta (IMEX-RK henceforth) type, to be discussed below.

In contrast, self-consistent models of core dynamics, which comprise in their full form an electromagnetic component responsible for the generation of the geomagnetic field by dynamo action, involve a Newtonian fluid of constant viscosity whose dynamics is strongly affected by planetary rotation. Their complexity lies in their inherent three-dimensional, global character and in the nonlinear coupling between field variables (velocity, pressure, temperature, magnetic field). These models came to the fore in the nineteen-nineties [START_REF] Glatzmaier | A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle[END_REF][START_REF] Kageyama | Computer simulation of a magnetohydrodynamic dynamo[END_REF], in the wake of the pioneering work of [START_REF] Glatzmaier | Numerical simulations of stellar convective dynamos. I-The model and method[END_REF] on the solar dynamo.

Most codes to date are part of Glatzmaier's legacy; they rely in the horizontal directions on a spectral representation of field variables using spherical harmonics. Nonlinear terms are computed in physical space, which requires forward and inverse spherical harmonic transforms to be performed at each time step. This step represents the most expensive computation of three-dimensional spherical simulations. Efforts to improve code performances focused on increasing the efficiency of spectral transforms by parallelization, using strategies based either on distributedmemory [START_REF] Clune | Computational aspects of a code to study rotating turbulent convection in spherical shells[END_REF] or shared-memory [START_REF] Schaeffer | Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations[END_REF]. In comparison, little work was devoted to reducing the time-to-solution using efficient time integration schemes. The performance benchmark by [START_REF] Matsui | Performance benchmarks for a next generation numerical dynamo model[END_REF] offers an interesting perspective on this state of affairs, as it reports the performance of 13 spectral codes on various laminar test problems (the study also includes two finite element codes). All spectral codes rely on a implicit-explicit scheme that treats linear terms, at the exception of the Coriolis force, implicitly, and nonlinear terms explicitly. The majority of codes resort to a θ-method for the linear term and have the Crank-Nicolson value θ = 1/2, (although some may use the stabler first-order θ = 0.6, see e.g. Hollerbach, 2000, around Eq. ( 21)).

Nonlinear and Coriolis terms are evaluated in 9 instances out of 13 with a second-order Adams-Bashforth (AB2 henceforth) method. Most codes combine CN for the implicit part with AB2 for the explicit part, forming what will we refer to in the following the traditional CNAB2 method. This popularity appears at odds with the flaws of the CNAB2 method, reported for instance by [START_REF] Tilgner | Spectral methods for the simulation of incompressible flows in spherical shells[END_REF] in his analysis of time integrators for fluid flow problems in spherical shell geometry: there may exist circumstances under which "the unreasonably popular CNAB2 " [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF] may not damp oscillatory modes at the expected physical rate (meaning they are not damped enough), which can then lead to a global instability if non-linearities are present (see also the discussions in [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF]Canuto et al., 2006, , §D.2.2.). An alternative consists of a second-order predictor-corrector approach, with the predictor and corrector stages based on AB2 and CN, respectively (as used in the Leeds code, see [START_REF] Willis | Thermal core-mantle interaction: Exploring regimes for "locked" dynamo action[END_REF]. The focus of the study by [START_REF] Matsui | Performance benchmarks for a next generation numerical dynamo model[END_REF] is the scalability of codes towards using petascale computers, based on the assessment of the efficiency of the various spatial parallelization strategies followed by the various contributors to the performance benchmark. It is intriguing to note that the paper does not even mention the benefit (in terms of efficiency and also accuracy, in view of performing turbulent simulations using spectral methods in space) that one could potentially gain from using more accurate and stabler time-schemes, on the condition that the gain compensates the extra cost.

A survey of literature reveals that alternatives to CNAB2 were considered occasionally in the past, for threedimensional Boussinesq thermal convection in axisymmetric (cylindrical or spherical) domains [START_REF] Fournier | A Fourier-spectral element algorithm for thermal convection in rotating axisymmetric containers[END_REF], where a BDF2/AB3 IMEX scheme was used), convection-driven dynamo action in Cartesian domains [START_REF] Stellmach | An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers[END_REF], where the SBDF2, sometimes called extended Gear of order 2, combination was employed), or rotating convection under the anelastic approximation in 2-D and 3-D Cartesian domains (Verhoeven and Stellmach, 2014, where SBDF3 was used). The influence of the chosen time integrator on the accuracy of the solution was also recently put forward by [START_REF] Lecoanet | Tensor calculus in spherical coordinates using Jacobi polynomials. part-II: Implementation and examples[END_REF]. Comparing CNAB2 and SBDF4, they showed that using the latter fourth order scheme allowed to assess a refined convergence of the reference values for benchmarks of convection and dynamo action in full spheres [START_REF] Marti | Full sphere hydrodynamic and dynamo benchmarks[END_REF] up to ten decimal places.

Following up on the investigation of [START_REF] Tilgner | Spectral methods for the simulation of incompressible flows in spherical shells[END_REF], [START_REF] Livermore | An implementation of the exponential time differencing scheme to the magnetohydrodynamic equations in a spherical shell[END_REF] analyzed the performance of second-order exponential time differencing (ETD) applied to dynamo action in spherical geometry, finding that for order 2, ETD methods were not to be preferred over the traditional CNAB2, given that their performance was similar, and the implementation of the latter much easier. [START_REF] Livermore | An implementation of the exponential time differencing scheme to the magnetohydrodynamic equations in a spherical shell[END_REF] also noted that the situation may become different would one resort to higher-order schemes. Such an endeavor was undertaken by [START_REF] Garcia | Exponential versus IMEX high-order time integrators for thermal convection in rotating spherical shells[END_REF] in the context of three-dimensional rotating convection. A comparison was made between exponential integrators and multistep IMEX schemes that had been previously investigated by [START_REF] Garcia | A comparison of high-order time integrators for thermal convection in rotating spherical shells[END_REF]. Based on the investigation of moderately supercritical configurations, the authors concluded on the superior accuracy of exponential integrators, at the expense of a larger cost.

In this study, the emphasis is set on the accuracy and efficiency of multistage IMEX-RK methods, which are compared with some multistep IMEX methods of order 2, 3 and 4. IMEX-RK methods have been almost ignored by the core dynamics community. Interestingly, however, [START_REF] Glatzmaier | An anelastic evolutionary geodynamo simulation driven by compositional and thermal convection[END_REF] implemented early on the IMEX-RK method proposed by [START_REF] Spalart | Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions[END_REF] to three-dimensional dynamo modeling in spherical geometry. This scheme, which will be evaluated in this study, combines a third-order explicit component with a secondorder implicit component. To our knowledge, no subsequent usage of that specific scheme was reported in the planetary core dynamics / dynamo community until a few years ago. It recently resurfaced in the study of magneto-convection in Cartesian geometry by [START_REF] Yan | Heat transfer and flow regimes in quasi-static magnetoconvection with a vertical magnetic field[END_REF]. Aside from the scheme by [START_REF] Spalart | Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions[END_REF], we note that [START_REF] Hollerbach | A spectral solution of the magneto-convection equations in spherical geometry[END_REF] defined a second order IMEX-RK scheme assembled from an explicit RK2 for its explicit part and a Crank-Nicolson for its implicit component. More recently, [START_REF] Marti | A computationally efficient spectral method for modeling core dynamics[END_REF] tested several of the IMEX-RK schemes proposed by [START_REF] Cavaglieri | Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff highdimensional ODE systems[END_REF] of second, third and fourth order applied to standard core dynamics problems in a spherical shell geometry, with a focus on the impact of the implicit or explicit treatment of the linear Coriolis force on the efficiency of their code. IMEX-RK methods were also investigated by one of us in the context of two-dimensional quasi-geostrophic spherical convection [START_REF] Gastine | pizza: an open-source pseudo-spectral code for spherical quasi-geostrophic convection[END_REF]. It was noted that of the three third-order schemes that were tested, the multistep SBDF3 by [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF] and the IMEX-RK BPR353 by [START_REF] Boscarino | Implicit-explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit[END_REF] (both of which will also be evaluated in this paper) had a similar efficiency in a rapidly-rotating, moderately turbulent configuration. The latter IMEX-RK scheme was recently employed by [START_REF] Tassin | Geomagnetic semblance and dipolar-multipolar transition in top-heavy doublediffusive geodynamo models[END_REF] to obtain turbulent double diffusive geodynamo models of higher accuracy. [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF] performed a thorough and inspiring comparison of IMEX-BDF, IMEX-RK and exponential integrators applied to a variety a problems, including the two-dimensional stratified Boussinesq equations and the quasigeostrophic equation , both in a periodic Cartesian domain. (The 6 IMEX-RK schemes that they analyzed are part of the 22 IMEX-RK schemes analyzed in this work.) Their conclusions can be tentatively summarized as follows: in the setups that they considered, exponential integrators are vastly superiors in terms of accuracy to multistep and multistage methods, even if some IMEX-RK schemes, such as the BHR553 scheme of [START_REF] Boscarino | On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF], may display a convergence rate better than the nominal third order. The exact treatment of linear diffusive terms enabled by exponential integrators is key in the moderately nonlinear setups they considered, and [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF] stressed that IMEX scheme could prove superior to exponential integrators when nonlinearities play a more sizeable role in the dynamics. In passing, [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF] also noted that for the 2D stratified Boussinesq equations, IMEX-RK methods "all displayed the disturbing ability to produce stable but inaccurate results at large step sizes" (their Fig. 9). We shall come back to this observation in our own analysis.

Owing in part to the ease of their implementation and interchangeability in modern computing frameworks (e.g. [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] IMEX-RK recently received attention in atmosphere and climate modeling [START_REF] Giraldo | Implicit-explicit formulations of a three-dimensional nonhydrostatic unified model of the atmosphere (NUMA)[END_REF][START_REF] Gardner | Implicit-explicit (IMEX) Runge-Kutta methods for non-hydrostatic atmospheric models[END_REF][START_REF] Vogl | Evaluation of implicit-explicit additive Runge-Kutta integrators for the HOMME-NH dynamical core[END_REF]. These studies looked into the possibility of using such schemes to overcome the severe limitations of pure explicit time marching that arise in nonhydrostatic models of the atmosphere which host acoustic waves propagating in the vertical direction. The smallness of the propagation time of those waves compared with the characteristic time of convective transport makes the problem stiff, and calls for an implicit treatment of those terms responsible for wave propagation. Using a testbed consisting of a gravity wave test and a baroclinic instability test from the 2012 dynamical core model intercomparison project [START_REF] Ullrich | Dynamical core model intercomparison project (DCMIP) test case document[END_REF], [START_REF] Gardner | Implicit-explicit (IMEX) Runge-Kutta methods for non-hydrostatic atmospheric models[END_REF] implemented an anisotropic splitting strategy, termed HEVI, for horizontally explicit -vertically implicit, implemented on 21 IMEX-RK schemes of order 2, 3, 4 and 5 (see their section 3.2.1 for their description). Likewise, a comprehensive comparative study of 27 IMEX-RK schemes led [START_REF] Vogl | Evaluation of implicit-explicit additive Runge-Kutta integrators for the HOMME-NH dynamical core[END_REF] to recommend 5 schemes that consistently perform better than the rest of the pack for the same two test cases (their Table 6).

With the aim of applying this systematic methodology to a problem relevant to the dynamics of planetary interiors, we focus here on 2D Boussinesq convection in a cylindrical annulus, in the absence of background rotation. This setup is admittedly simpler than some problems recently reported in the literature. Yet, its modest size allows us to cover a relatively broad range of behaviors including turbulent solutions and to perform a comprehensive investigation of 22 IMEX-RK schemes, in addition to 4 IMEX multistep schemes and two fully explicit schemes.

The outline of the paper is the following: we present the physical model and its numerical approximation in section 2. Results follow in section 3, where accuracy, order reduction, and computational efficiency are investigated and discussed for 11 different dynamical setups. We next summarize our findings and conclude in section 4 with some tentative recommendations for subsequent use of IMEX-RK schemes in the context of three-dimensional dynamo simulations in spherical geometry.

The model and its numerical approximation

Governing equations

Let us operate in cylindrical coordinates (s, ϕ, z) with local unit vectors ŝ, φ and ẑ. We consider a Newtonian fluid contained in a flat annulus of outer radius s o and inner radius s i . The fluid is subjected to a uniform inward radial gravity field of amplitude g 0 . The outer and inner cylindrical walls are maintained at uniform temperatures T o and T i , respectively. The temperature contrast ∆T ≡ T i -T o is positive and gives rise to convective flow when it exceeds a critical value.

The primitive state variables describing the fluid are its velocity u = u s , u ϕ , pressure p and temperature T .

The material properties of the fluid relevant for the problem of interest here are its density ρ, kinematic viscosity ν, thermal diffusivity κ and thermal expansion coefficient α. The equation of state that relates changes in temperature δT to changes in density δρ is δρ = -αρδT.

Under the Boussinesq approximation, the properties of the fluid are homogeneous, save for density that can vary with the local temperature according to the previous law when (and only when) the gravitational force is computed. The basic state about which convection can take place is the motionless hydrostatic conducting state. The conducting temperature profile is axisymmetric,

T c (s, ϕ) = T c (s) = ∆T log(s/s i ) log(s i /s o ) + T i . (1)
We scale length, time, and temperature by the gap width (D = s os i ), the viscous diffusion time D 2 /ν, and ∆T , respectively. We also choose pressure, p, to be scaled by ρ o ν 2 /D 2 , where, ρ o is the background density.

Conservation of mass, momentum and energy results in the following set of equations

∇ • u = 0, (2) ∂u ∂t = -∇ • (u ⊗ u) -∇p + Ra Pr T ŝ + ∇ 2 u, (3)
∂T ∂t = -∇ • (uT ) + 1 Pr ∇ 2 T, (4)
to be complemented with initial and boundary conditions (see below).

The dimensionless control numbers are the Rayleigh number Ra and the Prandtl number Pr, defined by

Ra = g o α∆T D 3 νκ , (5)
and

Pr = ν κ . (6)
The two-dimensional nature of the problem and the incompressibility constraint prompt us to resort to a vorticitystreamfunction formulation, see e.g. [START_REF] Glatzmaier | Introduction to Modeling Convection in Planets and Stars: Magnetic Field, Density Stratification[END_REF], §2.1 or [START_REF] Peyret | Spectral Methods for Incompressible Viscous Flow[END_REF], §II.6. Let the symbol overbar denote the azimuthal averaging operator,

f = 1 2π 2π 0 f (ϕ)dϕ.
We introduce a streamfunction ψ such that

u s = 1 s ∂ψ ∂ϕ , (7a)
u ϕ = u ϕ - ∂ψ ∂s . ( 7b
)
Given the periodicity of the domain in the azimuthal direction, the decomposition of the azimuthal flow into a mean component u ϕ and a non-zonal component -∂ψ/∂s ensures the periodicity of pressure in the azimuthal direction (e. g. Plaut and Busse, 2002, §2). The evolution of the mean component is governed by the azimuthal average of Eq. ( 3). The axial vorticity ω = ωẑ is given by

ω = 1 s ∂(su ϕ ) ∂s -∇ 2 ψ, (8)
and its time-dependency is controlled by the axial component of the curl of the momentum equation Eq. (3). In summary, the set of dimensionless equations to solve in the vorticity-streamfunction formulation reads

∂u ϕ ∂t = -u s ω + ∆u ϕ , (9a)
∂ω ∂t = -∇ • (uω) + ∇ 2 ω - Ra Pr 1 s ∂T ∂ϕ , (9b)
∂T ∂t = -∇ • (uT ) + 1 Pr ∇ 2 T, (9c)
ω = 1 s ∂(su ϕ ) ∂s -∇ 2 ψ, (9d)
u = u ϕ + ∇ × (ψẑ), (9e)
where the modified Laplacian operator ∆ = ∂ s (∂ s + 1/s).

Boundary conditions

Regarding the mechanical boundary conditions for the annulus, we shall assume throughout a no-slip boundary condition u ϕ = 0 along the curved inner and outer boundary walls, together with with the impermeable condition u s = 0. The latter condition implies that

ψ = 0 at s = s i , s o . ( 10
)
The no-slip boundary condition implies

∂ψ ∂s = u ϕ = 0 at s = s i , s o , (11)
which results in vorticity at the boundaries to be

ω = - ∂ 2 ψ ∂s 2 at s = s i , s o . (12)
Thus, we have four boundary conditions on ψ and two for u ϕ . For the temperature, the dimensionless boundary conditions are

T = 1 at s = s i , (13a)
T = 0 at s = s o .

(13b)

Diagnostics

Before dealing with the numerical approximation of the problem, let us introduce a few diagnostics to obtain useful information and to check solution correctness. In the following, we denote the spatial average over the area A of the annulus by a double overbar and time average by angular brackets • • • . For a field f (s, ϕ),

f = 1 A A f (s, ϕ)sdsdϕ. ( 14
)
where, A = π(s 2 os 2 i ) is the area of the annulus. We compute the kinetic energy E k at specified times of the simulation. At a given instant in time, it is given by

E k (t) = 1 2 u 2 s + u 2 ϕ (t). ( 15
)
The Nusselt number Nu quantifies the ratio of total heat flux to the reference heat flux carried by conduction alone.

We define it at the inner and outer boundaries by considering the time and azimuthal averages of the temperature, e.g.

Nu o =       d T ds       s=s o / dT c ds s=s o =       d T ds       s=s o s o log s o s i , (16)
for the outer boundary; the expression for Nu i being obtained upon substituting the s o factor by s i in this equation.

The balance of inner and outer Nusselt numbers (incoming and outgoing heat fluxes) indicates that thermal relaxation has been reached and it is a good indicator for convergence of the solution. Now, we define the Reynolds number Re which measures the ratio of inertial to viscous forces. With our choice of scales, it is given as the time-averaged root-mean-square magnitude of the velocity

Re = u 2 s + u 2 ϕ 1/2 = 2E k . (17)
The next diagnostic quantity we compute from the solution at specified time intervals is the power balance. From the solution, we check if heat loss by viscous dissipation balances on average the buoyancy input power (e.g. [START_REF] King | Heat transfer by rapidly rotating Rayleigh-Bénard convection[END_REF]. The expression for viscous dissipation is given as

D ν (t) = u • ∇ 2 u(t). ( 18
)
Using the vector identity ∇ × ∇ × u = ∇(∇ • u) -∇ 2 u, the definition of vorticity ω and the incompressibility i constraint, the viscous dissipation term becomes

D ν (t) = -[u • (∇ × ω)](t). (19)
When no-slip boundary conditions are prescribed, it further simplifies to

D ν (t) = -ω 2 (t). ( 20
)
The buoyancy input power P reads

P(t) = Ra Pr u s T (t). ( 21
)
On time average, we expect the solution to satisfy

P = -D ν . (22)

Spatial discretization

We apply a Fourier-collocation approach to discretize Eqs. (9a-9e) in space. The Fourier expansion is performed along the azimuthal direction which is naturally periodic and a Chebyshev collocation method is employed along the radial direction, see e.g. [START_REF] Glatzmaier | Introduction to Modeling Convection in Planets and Stars: Magnetic Field, Density Stratification[END_REF]. The truncated Fourier expansions of field variables with a dependency on the azimuthal angle read

ω(s, ϕ, t) ≈ 2 N m m=0 ω m (s, t)e imϕ , (23a)
T (s, ϕ, t) ≈ 2 N m m=0 T m (s, t)e imϕ , (23b)
ψ(s, ϕ, t) ≈ 2 N m m=1 ψ m (s, t)e imϕ , (23c)
u s (s, ϕ, t) ≈ 2 N m m=1 u sm (s, t)e imϕ , (23d)
u ϕ (s, ϕ, t) ≈ u ϕ (s, t) + 2 N m m=1 u ϕm (s, t)e imϕ . ( 23e
)
where N m is the maximum order of the truncation, and consequently the number of Fourier modes is N m + 1, f is the real part of a complex-valued function f and the single prime on the summation symbol means that the m = 0 term in the series is multiplied by 1/2.

Substituting these expressions into Eq. ( 9) results in the following set of nondimensional equations

∂u ϕ ∂t (s, t) = -u s ω(s, t) + ∆u ϕ (s, t), (24a)
∂ω m ∂t (s, t) = -∇ m • (uω) m (s, t) + ∇ 2 m ω m (s, t) - Ra Pr im s T m (s, t), for m > 0, ( 24b
)
∂T m ∂t (s, t) = -∇ m • (uT ) m (s, t) + 1 Pr ∇ 2 m T m (s, t), for m ≥ 0, ( 24c
)
ω 0 (s, t) = 1 s ∂(su ϕ ) ∂s (s, t), ω m (s, t) = -∇ 2 m ψ m (s, t), for m > 0, ( 24d
)
u sm (s, t) = im s ψ m (s, t), for m > 0, u ϕ0 (s, t) = u ϕ (s, t), u ϕm (s, t) = - ∂ψ m ∂s (s, t), for m > 0, ( 24e
)
where the Fourier mode-dependent divergence and Laplacian operators read

∇ m • a = (1/s)∂ s (sa s ) + (im/s)a ϕ and ∇ 2 m = ∂ 2 s + (1/s)∂ s -(m 2 /s 2 )
, respectively. The notation (• • • ) m refers to the m th Fourier mode of the term inside the brackets. In the following, we shall refer to the previous formulation as the sm formulation.

We proceed with a radial approximation based on Chebyshev polynomials C n up to degree N s -1. Each field variable g m (s, t) appearing in the previous system is expanded according to

g m (s, t) ≈ 2 N s -1 1/2 N s -1 n=0 g mn (t)C n [x(s)], (25)
where the double quote implies that the first and last terms are multiplied by 1/2. The cylindrical radius s is mapped into coordinate x by

x = 2s -s o -s i s o -s i = 2s -s o -s i (26)
in order to use the Chebyshev-Gauss-Lobatto points defined by

x k = cos kπ N s -1 (27)
with k = 0 to N s -1. The discrete Chebyshev expansion evaluates g m (s = s k , t) such that

s k = s o -s i 2 x k + s o + s i 2 . ( 28
)
Conversely,

g mn (t) = 2 N s -1 1/2 N s -1 k=0 g(s k , t)C n (x k ), (29)
in which

C n (x k ) = cos(n arccos x k ) = cos nπk N s -1 .
In practice, our unknowns consist of the g mn . The radial approximation of Eqs. 24 leads to the following semi-discrete problem

d dt M u ϕ = N u ϕ + L u ϕ u ϕ , (30a)
d dt M ω m = N ω,m + L ω,m ω m + B m T m for m > 0, (30b)
d dt M T m = N T,m + L T,m T m for m ≥ 0, (30c)
M ω 0 = L 0 u ϕ , M ω m = L ωψ,m ψ m for m ≥ 0, ( 30d
) M u s,m = L u s ψ,m ψ m , for m > 0, M u ϕ,0 = u ϕ , M u ϕ,m = L u ϕ ψ,m ψ m , for m > 0. ( 30e
)
We have omitted the remaining dependency to time for the sake of conciseness. Bold capital letters refer to matrices acting upon column vectors of the kind g m , that contain the N s coefficients of the expansion of g m (s) on Chebyshev polynomials as given by Eq. ( 25) above,

g m = g m0 , . . . , g m N s -1 T , (31)
where T means transposition without conjugation. The M matrix on the left-hand side of Eqs. 30a-30e is a N s × N s square matrix that converts modal values to gridpoint values, since the equalities 30a-30e are prescribed at the N s collocation points s k , with k ∈ 0, . . . , N s -1. Its k-th row reads

γ C 0 (x k )/2, C 1 (x k ), . . . , C N s -2 (x k ), C N s -1 (x k )/2 , ( 32
)
where the constant γ = [2/(N s -1)] for each collocation grid point s k . For instance,

N T,m = -∇ m • (uT ) m (s 0 ), ∇ m • (uT ) m (s 1 ), . . . , ∇ m • (uT ) m (s N s -2 ), ∇ m • (uT ) m (s N s -1 ) T .
We resort to a pseudo-spectral approach: Nonlinear terms are evaluated on the physical grid, prior to being transformed back to spectral space. For instance, in the previous equation, the product uT is computed on the s -ϕ grid and transformed back to Fourier space using a Fast Fourier transform; a discrete cosine transform (Press et al., 2007, 12.4.2) is next used to transform quantities from the sm space to the nm space. The divergence is computed using the appropriate three-term recurrence relation for Chebyshev polynomials (e.g. Canuto et al., 2006, §2.4.2).

Linear terms in the system 30 are cast in terms of matrix-vector products; they appear in the algebraic equations Eqs. 30d-30e, in addition to the right-hand side of differential Eqs. 30a-30c. The matrices are all N s × N s , and they possibly involve derivatives of the Chebyshev polynomials, C , C , etc. For instance, the k-th row of matrix

L ω,m reads γ 1 2       C 0 (x k ) + 1 s k C 0 (x k ) - m 2 s 2 k C 0 (x k )       , C 1 (x k ) + 1 s k C 1 (x k ) - m 2 s 2 k C 1 (x k ), . . . , C N s -2 (x k ) + 1 s k C N s -2 (x k ) - m 2 s 2 k C N s -2 (x k ), 1 2       C N s -1 (x k ) + 1 s k C N s -1 (x k ) - m 2 s 2 k C N s -1 (x k )       ,
while the k-th row of the B m buoyancy matrix reads

- Ra Pr imγ s k C 0 (x k )/2, C 1 (x k ), . . . , C N s -2 (x k ), C N s -1 (x k )/2 .
Boundary conditions are prescribed through the appropriate modification of some of the matrices entering Eq. ( 30).

We shall get back to this shortly.

Time discretization

The Chebyshev-Fourier collocation method leads to a semi-discrete set of differential algebraic equations (DAE) of the generic form

dx dt = N(x, z) + Lx, (33a)
0 = G(x, z), ( 33b
)
where N is the nonlinear operator acting on the differential state vector x = u ϕ , ω, T and the algebraic state vector z = ψ, u , L is the linear operator acting upon x and G is the linear operator relating x and z. According to [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF] this defines a semi-explicit DAE of order 1. We will solve it in time using methods developed for ordinary differential equations (ODEs), ensuring that constraint (33b) is satisfied through a standard solution technique to be detailed below.

We resort to implicit-explicit methods that treat Lx and N(x, z) implicitly and explicitly, respectively. There are different families of IMEX methods, and here we are interested in both IMEX multistep and IMEX multistage methods [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF][START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II[END_REF][START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF][START_REF] Kennedy | Additive Runge-Kutta schemes for convection-diffusion-reaction equations[END_REF].

Solution technique

Before we get to the details of the multistep and multistage methods, let us describe the backbone of our solution technique, in particular how we deal with boundary conditions and the DAE when advancing in time. The former are enforced through the implicit solves, while the latter is taken care of by means of a block-matrix solve.

To update the field values and advance in time, our approach is the following: the equations for mean flow (30a)

and temperature (30c) are solved first, as their implicit components are not coupled with the vorticity or the streamfunction equations (Eq. (30b) and Eq. ( 30d)). This comes down to inverting a linear system of the form (αM -βL) y = r.h.s.

where the coefficients α and β are both positive real numbers that depend on the time integrator, y stands for the vector containing the N s coefficients of the Chebyshev expansion of mean flow or temperature, and the right-hand side r.h.s contains a mix of linear and nonlinear components. The first and last rows of the N s × N s matrix to invert, αM -βL, are modified in order to enforce the two boundary conditions that apply either on temperature or on the mean flow at s = s i and s = s o , respectively (see e.g. [START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF], Table 4). The corresponding entries of the r.h.s vector are modified accordingly and set to 0 or 1, depending on the boundary condition.

Next, we update the vorticity and the non-zonal streamfunction, Eq. (30b) and Eq. ( 30d), for each Fourier mode m > 0. We do so by inverting the following 2N s × 2N s block-matrix,

          αM -βL ω,m 0 M -L ωψ,m                     y ω,m y ψ,m           =           r.h.s. ω 0          
The first half of the vector [y ω,m , y ψ,m ] T , y ω,m , is the updated vorticity, while its second half is the updated stream- for the streamfunction, by modifying the first, N s -th, N s + 1-th and last row of the block matrix (Glatzmaier, 2013, §10.3.2). The first N s entries of these rows are set to 0, while the remaining N s entries are modified in order enforce the homogeneous Dirichlet and Neumann conditions for the streamfunction, see e.g. Table 4 in [START_REF] Julien | Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods[END_REF]. Accordingly, the first, N s -th, N s + 1-th and last entry of the right-hand side vector [r.h.s. ω , 0] T are set to 0.

See also (Gastine, 2019, his Fig. 1a,) for a graphical illustration of this implementation.

Finally, the updated mean flow and non-zonal streamfunction (which satisfy their respective boundary conditions, given by Eqs. (10-11), allow us to evaluate the velocity components (Eq. ( 30e)), thereby permitting the evaluation of the nonlinear terms N(x) necessary for the next update.

We shall now describe the multistep and multistage time discretization methods.

IMEX multistep methods

Multistep methods rely on a polynomial interpolation in time. Let K be the number of steps of an IMEX multistep method, with K ≥ 1. Let ∆t denote the timestep size and x i denote the approximate solution for the differential state vector at time t i = i∆t, leaving aside the algebraic z, which is updated alongside x following the solution technique detailed above. Then, for fixed ∆t, following [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF], a general linear multistep IMEX method applied to Eq. (33a) can be written as

(1 -∆tc -1 L) x i+1 = K-1 j=0 a j x i-j + ∆tb j N(x i-j ) + ∆tc j Lx i-j , (34)
where c -1 0. It is noteworthy that a K-step IMEX method cannot have order of accuracy greater than K [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF]. The IMEX multistep methods we use for this study have the same order as the number of steps K.

In this work we shall consider multistep methods of order 2, 3 and 4: the popular Crank-Nicolson Adams-Bashforth method of order 2 (CNAB2) already seen in the introduction, and the semi-implicit BDF (SBDF)

schemes given e.g. in [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF] of order 2, 3 and 4. The three SBDF schemes apply a backward differentiation formula to the implicit part, and an extrapolation formula to the explicit part.

In our convergence analysis, ∆t will remain fixed. We shall activate variable time-step for the equilibrated regime of the most turbulent of our reference cases to be reached (section 3.1), and for the stability analysis of section 3.6. The vectors of coefficients a, b and c for the four selected schemes are given in Appendix A.

IMEX multistage methods

IMEX multistage methods rely on quadrature rules to evaluate intermediate stages of x (and z) between discrete times t i and t i+1 . The multistage methods of interest for this work are often referred to as IMEX-RK (for IMEX Runge-Kutta) methods, which indicates that they involve a diagonally implicit Runge-Kutta (DIRK) and an explicit Runge-Kutta (ERK) schemes [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF]. Unlike multistep methods, their stability region can increase slightly with their order.

Let K denote the number of internal stages of an IMEX-RK method. At each substage k ∈ {1, . . . , K} of an IMEX-RK scheme applied to Eq. (33a), one has

(1 -∆ta I kk L)y k = x i + ∆t k-1 j=1 a E k j N(y j ) + ∆t k-1 j=1 a I k j Ly j , (35)
where the coefficients a I k j and a E k j define two matrices for the DIRK and ERK schemes, A I and A E , respectively. The first stage is defined by y 1 = x i . At each stage, the algebraic variables z are updated alongside the differential variables following the strategy detailed in section 2.5.1.

The DIRK and ERK components can be independently summarized using Butcher tableaux [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF]

, § II.1) 0 0 c E 2 a E 21 0 . . . . . . . . . . . . c E K a E K1 • • • a E KK-1 0 b E 1 • • • • • • b E K , 0 0 c I 2 a I 21 a I 22 . . . . . . . . . . . . c I K a I K1 • • • a I KK-1 a I KK b I 1 • • • • • • b I K . ( 36
)
The design of an IMEX-RK method implies a set of constraints on the coefficients

A I , b I , c I , A E , b E , c E , the
number of which depends on the sought order of accuracy. Extra constraints can be added to accommodate discrete algebraic equations [START_REF] Boscarino | On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF]. According to the classification by [START_REF] Boscarino | Error analysis of imex runge-kutta methods derived from differential-algebraic systems[END_REF], the schemes considered in this study belong to the CK type, as the first row of the implicit matrix A I contains zeros. There exists schemes that have non-zero entries in this row, which belong to the so-called type A [START_REF] Boscarino | Error analysis of imex runge-kutta methods derived from differential-algebraic systems[END_REF], the first examples of which were introduced by Pareschi and Russo (2005); we do not investigate such schemes in this study.

In addition, all the schemes considered in this work have c E = c I , which is why c E 1 = c I 1 = a I 11 = 0 above. Also, note that in practice, c E and c I are not required for thermal convection which has no explicit time-dependent forcing and is therefore an autonomous process.

In spite of what Eq. ( 35) suggests, the IMEX-RK methods we analyze in this study do not necessarily comprise the same number of implicit and explicit stages, K I and K E . This is handled by setting the appropriate columns of A I , and entries of b I (or of A E and b E ) to zero. The number of linear inversions per time step, n I , is equal to K I provided a I k1 = 0, ∀k ∈ 1, . . . , K I . DIRK schemes are called "stiffly accurate" if

b I T = e T K A I , (37)
where e K = [0, . . . , 0, 1] T . Following [START_REF] Boscarino | Implicit-explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit[END_REF], we say that a IMEX-RK scheme is globally stiffly accurate (GSA) if b I T = e T K A I and b E T = e T K A E , and c I K = c E K , which implies that the updated state vector is identical to the last internal stage value of the scheme. For non-GSA schemes, the updated differential state vector is assembled as

x i+1 = x i + ∆t K j=1 b E j N(y j ) + ∆t K j=1 b I j Ly j . (38)
If the chosen scheme requires such an assembly then further work is needed to ensure that the updated vorticity meets condition (12) , see section 2.5.4 below.

We follow e.g. [START_REF] Boscarino | Implicit-explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit[END_REF] and identify each IMEX-RK scheme by the initials of the authors (if they are no more than 3), and three numbers (K I , K E , r) denoting, respectively, the number of implicit and explicit stages, and the theoretical order of accuracy. Exceptions to this rule are DBM553 from [START_REF] Vogl | Evaluation of implicit-explicit additive Runge-Kutta integrators for the HOMME-NH dynamical core[END_REF] and BHR553 from [START_REF] Boscarino | On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF] where we kept the initials of the original name; and PC432 which is a second order three stage predictor/corrector scheme constructed using the explicit scheme by [START_REF] Jameson | Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes[END_REF] for its explicit component and a Crank-Nicolson for its implicit counterpart (see Appendix B for details of its Butcher tableaux).

In this work we investigate the properties of 22 IMEX-RK schemes, whose properties are summarized in Table 1.

Treatment of IMEX-RK methods with an assembly stage

IMEX-RK methods which are not stiffly accurate require the assembly stage Eq. ( 38) to be performed. The assembly of the temperature and mean flow is straightforward, as it is a linear combination of nonlinear and linear terms evaluated at the substages. Applying the same linear combination for the vorticity, however, does not guarantee that the boundary conditions ( 12) are properly enforced. Indeed, the assembly stage does not lend itself to the solution method outlined in section 2.5.1 for the vorticity and streamfunction, which allows one to bypass the explicit enforcement of the boundary conditions on vorticity. The last column provides the relevant reference augmented with a section or paragraph number, and possibly the name of the scheme as it appears in the reference. To make sure that the vorticity built at the assembly stage is consistent with the boundary conditions, we follow the strategy outlined by [START_REF] Johnston | Comparison of turbulent thermal convection between conditions of constant temperature and constant flux[END_REF].

Scheme K I K E o n I S. A. S. A. b I = b E storage Reference DIRK ERK ARS222 2 2 2 X 4
We begin by assembling a first guess of the final vorticity y ω,m for each mode m > 0 by means of Eq. ( 38).

Using that intermediate value, we compute y ψ,m by inverting

L ωψ,m y ψ,m = My ω,m , (39)
having modified the first and last lines of L ωψ,m so that ψ m = 0 at s = s i and s = s o .

The knowledge of y ψ,m makes it possible to construct a local interpolant in the vicinity of the two walls, L w , that is constrained by the values of ψ m at the first J + 1, say, Chebyshev-Gauss-Lobatto points, ψ j m (that include the point on the wall) and the extra requirement that ∂ψ m /∂s = 0 on the wall as well. For the inner wall, s = s i , the interpolant reads

L w (s) = J j=0 ψ j m j (s) -(s -s i ) 0 (s)         J j=0 ψ j m j (s i )         , (40)
with a similar expression for the outer wall. In this expression, j is the Lagrange polynomial attached to the j-th point away from the wall, and j is its first derivative. This local interpolant allows us to compute the vorticity on the inner and outer walls,

ω m (s i ) = - d 2 L w ds 2 s=s i = - J j=0 ψ j m j (s i ) -2 j (s i ) 0 (s i ) , (41)
ω m (s o ) = - d 2 L w ds 2 s=s o = - J j=0 ψ j m j (s o ) -2 j (s o ) 0 (s o ) , (42)
where j denotes the second derivative of j . We form a vector of nodal values of vorticity whose interior values are based on y ω,m and whose boundary values are the ones we just computed based on the local interpolant L w . We finally determine y ω,m by applying the inverse of M to this vector of nodal values. In our experiments, we set the value of J to 14.

Fully explicit RK methods

In addition to the IMEX multistep and IMEX-RK multistage techniques detailed above, we found it useful sometimes to consider two well-known fully explicit methods, the explicit Runge-Kutta methods of order 2 and 4, RK2 and RK4 (e.g. Canuto et al., 2006, Eqs. D.2.15 and D.2.17). The solution technique that these methods imply is based on the technique laid out for IMEX-RK methods (no linear solve, except at the assembly stage).

Implementation and validation

The code for solving the problem using the aforementioned pseudospectral methods and time-stepping strategies was written from scratch in the Fortran programming language. The code contains several modules and subroutines where each module has specific dependencies. The fast Fourier and discrete cosine transforms resort to the FFTW3 library [START_REF] Frigo | The design and implementation of fftw3[END_REF]. The matrix equations are solved using standard matrix solvers available in the LAPACK routines dgetrf and dgetrs [START_REF] Anderson | LAPACK Users' Guide[END_REF]. The dgetrf routine is used for computing the LU factorization and the dgetrs routine is used for solving the system using the factored matrix obtained by using the dgetrf routine.

To benchmark the code against peer-reviewed results, we compare it with a reference solution obtained by [START_REF] Alonso | Transition to temporal chaos in an O(2)-symmetric convective system for low Prandtl numbers[END_REF]. They performed their numerical simulations using spectral methods with a fixed radius ratio s i /s o = 0.3 and Prandtl number Pr = 0.025 (which corresponds to liquid Mercury Hg), and the second-order stiffly stable time integrator by [START_REF] Karniadakis | High-order splitting methods for the incompressible Navier-Stokes equations[END_REF]. In table 2 we list the dependency of the equilibrated Nusselt number Nu = Nu o = Nu i to the Rayleigh number reported by [START_REF] Alonso | Transition to temporal chaos in an O(2)-symmetric convective system for low Prandtl numbers[END_REF] and obtained here using the ARS443 IMEX-RK time integrator, together with N s = 32 and N m = 192. Furthermore, for the range of Rayleigh numbers shown in table 2, we observe an oscillation of the solution about the periodic azimuthal direction, which is a characteristic of low Prandtl number fluids. For Ra = 6500, the frequency of oscillation we find is f = 5.15 which exactly matches value published by [START_REF] Alonso | Transition to temporal chaos in an O(2)-symmetric convective system for low Prandtl numbers[END_REF]. Thus we ascertain that the code was benchmarked and ready to be used for the study of various time integration methods.

Results

We begin by a presentation of the 11 cases studied in this work, followed by the analysis of the convergence properties of the time schemes we investigated. We investigate the likely causes of the order reduction observed for some configurations, and finally weigh these findings against a more practical estimate of the computational efficiency.

Cases studied

All cases considered have a radius ratio s i /s o set to 0.35. They are initialized with a temperature perturbation of localized compact support as introduced by [START_REF] Gaspari | Construction of correlation functions in two and three dimensions[END_REF], of width 0.1/ √ 2 and amplitude 10 -4 .

The properties of the cases are summarized in table 3. This table comprises the input control parameters Pr and Ra, the Reynolds number Re (Eq. ( 17)), Nusselt number at the outer boundary (Eq. ( 16)) and the temporal averages of the buoyancy input power (Eq. ( 21)), and heat loss by viscous dissipation (Eq. ( 20)). In addition, we provide in this table the spatial discretization parameters N s and N m introduced in the previous section. Unless otherwise stated, a given case was always run for the same (N s , N m ) pair. That pair was chosen to make spatial discretization error negligible against temporal discretization error. We chose to run 3 cases with a Prandtl number equal to 0.025, which corresponds to liquid metals, 7 cases with Pr equal to 1, which corresponds to a commonly taken value in numerical simulations, and one case with Pr = 40, to have at least one situation in the large Pr limit. The numbering in table 3 was adopted to follow the increase of the Reynolds number. Case 0 is extremely laminar, while case 10 is our most turbulent case with Re > 10 4 . Our goal is to exercise the time schemes over a broad range of regimes.

For this radius ratio, and regardless of the value of the Prandtl number considered (0.025, 1, 40), the most unstable convective mode has a threefold symmetry in the periodic azimuthal direction, that is to say that the value of the critical wavenumber m crit = 3. The corresponding value of the critical Rayleigh number is Ra crit = 1768.

In the range of forcing that we cover, the threefold symmetry is a persisting feature. When increasing the level of turbulence, the energy found in other wavenumbers increases, by virtue of the larger importance taken by turbulent transport of momentum and heat. This is shown in Figure 1, which displays the time averaged kinetic energy spectra of the 11 cases. The kinetic energy in each Fourier mode m is given by

E k (m = 0) = π s o s i u ϕ 2 sds, E k (m > 0) = 2π s o s i |u sm | 2 + |u ϕm | 2 sds.
Note that for case 10 (the most turbulent case) the azimuthal truncation (the value of N m ) chosen enables a 10 6

factor to be achieved between the highest energy level (for m = 3) and the lowest energy level (around m = N m ).

We use a similar criterion to set the truncation in radius, N s . How stiff are these cases numerically? We shall see below that stiffness, as measured by the disparity between linear and nonlinear time scales, remains moderate across the region of parameter space explored by the cases, with a stiffness parameter that varies between 10 -3 and 10 -5 (see section 3.5 and Table 4 below for more details).

We now consider in Fig. 2 a snapshot of the solution obtained for cases 2 and 10. In the latter case, the temperature field (Fig. 2d) shows three major plumes originating from the hot inner boundary, which reflect the maximum energy at m = 3 shown in Fig. 1. Accordingly, the vorticity in Fig. 2c exhibits a large scale m = 3 overturning circulation, with pockets of intense vorticity found in the eyes of the large-scale circulation. Note that in this set-up the plumes are anchored at the inner boundary, and that time-dependency appears mostly in the form of undulations occurring at their tip. In contrast, case 2, which corresponds to a more viscous fluid, and a lower level of forcing, is more laminar; its vorticity is notably concentrated along the edges of the large-scale convective cells (Fig. 2a). The temperature field shown in Fig. 2b appears symmetrical with regard to the top and bottom boundary layers, which are destabilized by similar cold or hot plumes displaying a mushroom head on top of a thin conduit.  2 error (Case 10) 10 -1 10 0 Normalised Δ 10 -2 10 -1 10 0 Normalised Δ Figure 4: L 2 error for (T, u s , ω) normalized by its maximum value (for a given case) as a function of ∆t normalized by its maximum value (for a given case too). Three schemes are considered: SBDF2, SBDF4 and KC664. The darker the color of a symbol, the higher the case number.

ARS222 ARS232 BPR442 PC432 SMR432 ARS233 ARS343 ARS443 BHR553 BPR533 BR343 CB443 CFN343 DBM553 KC443 LZ543 CB664 CFN564 KC664 KC774 LZ764 KC885 CNAB2 SBDF2 SBDF3 SBDF4
are unstable over the range of ∆t investigated here for cases 2 and 10. This is due to the stricter limits imposed on ∆t, which are such that stability coincides with reaching the numerical roundoff error plateau. Case 3 provides a configuration for which fully explicit schemes can be studied; the corresponding convergence curves are provided in Appendix C).

In summary, multistep schemes display the expected convergence rate, and an error level overall higher than IMEX-RK schemes of the same order. To get a better understanding of the error behavior across the 11 cases, we now show in Fig 4 how it evolves for the IMEX-RK KC664 scheme of [START_REF] Kennedy | Additive Runge-Kutta schemes for convection-diffusion-reaction equations[END_REF] and the two multistep schemes SBDF2 and SBDF4. Despite the fact that the error level and admissible time step values vary across the 11 cases, we collapse the information by normalizing the error of a given scheme for a given case by its maximum value, and the timestep ∆t by the value it has when this maximum value is obtained. In addition symbols are represented such that the darker the symbol, the larger the case number. Note also that prior to collapsing the data, we got rid of those unwanted points located on a plateau, such as those present in Fig. 3 for cases 2 and 10.

In the log-log representation of Fig. 4, we first observe that the behavior of SBDF2 and SBDF4 is well captured across the cases by two straight lines, of slopes 2 and 4, respectively, for the three fields of interest, T , u s and ω. This illustrates nicely that they indeed conform to their expected convergence rate over the range of regimes studied. KC664, on the contrary, displays some scatter, lighter symbols (laminar cases) being overall further apart from the 4th order reference defined by SBDF4 than darker symbols (turbulent cases). This is particularly striking for the vorticity field, more moderate for u s , and even less pronounced for T . For cases 9 and 10 (darker symbols) the vorticity appears to transition from 4th order to 2nd order as the value of the normalized time step decreases. This is evidence of order reduction, to an extent that depends strongly on the regime considered.

Order reduction

We now quantify order reduction for all cases and the 22 IMEX-RK schemes considered. To assess the order of convergence of a scheme µ, we consider the two largest values of ∆t used in the convergence analysis and seek a fit for the error of the form e(∆t) ∝ ∆t µ , with µ the sought order, that depends on the scheme and the case, and the field of interest (T , u s and ω). The two largest values of ∆t are admittedly not representative of the asymptotic behavior of a scheme when ∆t → 0.

Several methods display more than one scaling through the range of tested time step sizes (recall Fig. 4 above), which makes it difficult to correctly define the order of convergence. This definition of the order is not perfect, as it is biased towards runs that favor the largest integration length over accuracy, for a given number of time steps performed. The 726 values of µ that we estimated are displayed in Figure 5. To try and synthesize the information, we additionally introduce a χ 2 measure of order reduction, defined by

χ 2 = cases µ m -µ t µ t 2 ,
where the measured order µ m is set to the theoretical order µ t when superconvergence is observed, i.e. when the measured order exceeds the theoretical order. Values of χ 2 are reported in Figure 6. Figure 6 reveals that order reduction, when it affects a scheme, is more severe for ω than it is for u s , which is itself more severe than the order reduction that impacts T , if it impacts T at all.

Overall, we observe in Fig. 5 that IMEX-RK order 2 schemes are immune to order reduction, with SMR432 and ARS232 showing superconvergence, in particular for high Reynolds number cases. Third-order schemes show a variety of behaviors. Two schemes stand out as being particularly impacted by order reduction: CB443 and DBM553. On the contrary, BHR553 is immune to order reduction, and occasionally superconverges, as also reported by e.g. [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF]. Fourth-order schemes are all prone to order reduction, especially CB664. Order reduction manifests itself mostly for our most laminar cases, from 0 to 5, and appears to be stronger in the most laminar cases. This phenomenon is a well-known issue that can arise due to two factors: the discrete algebraic equation (Eq. ( 33)), and the stiffness of the problem (consult Boscarino, 2007, for a thorough theoretical investigation of this issue). As previously noticed by e.g. [START_REF] Kennedy | Additive Runge-Kutta schemes for convection-diffusion-reaction equations[END_REF] in their analysis of order reduction in a convection-diffusion-reaction problem, both differential variables (T , ω here) and algebraic variables (u s here) are impacted.

An attempt to assess the impact of the DAE on order reduction

We try to estimate to which extent order reduction can be ascribed to the DAE by considering the reduced advection-diffusion problem

∂T ∂t = -∇ • (uT ) + 1 Pr ∇ 2 T, (43)
subject to the same boundary conditions for temperature, and the same procedure for spatial discretization as detailed above. To cover the 11 cases investigated, we specify u by extracting the velocity from a random snapshot of the full problem for each case, and take the temperature field from that snapshot as the initial condition. We use this strategy in order to retain the physical properties of the solution to the full problem (in particular the level 0 1 2 3 4 5 6 7 8 9 10 case number Estimated order 0 1 2 3 4 5 6 7 8 9 10 case number 0 1 2 3 4 5 6 7 8 9 10 case number of turbulent transport) while getting rid of its algebraic component. Reference solutions to this reduced problem are produced via its time integration with the SBDF4 scheme using, again, a time step size small enough to allow convergence properties to be determined over two orders of magnitude, for each case. The temperature field T is the only field that remains to evaluate accuracy and convergence properties. Estimated orders are now presented in Figure 7, noting that the findings that we present are not sensitive to the randomly chosen snapshot. Inspection of Figures 7 reveals that order reduction persists, but to a much lesser extent; there is less scatter in Fig. 7 than in the left panel of Fig. 5. In fact order reduction is now restricted to well-identified cases, namely cases 0, 5 and 7, for which the curves of Fig. 7 tend to globally dip. Schemes that exhibit occasional superconvergence (SMR432, BHR553, ARS232) do not superconverge for those cases. Schemes of order 3 and 4 that underperform for the full problem (most notably CB443, DBM553 and CB664) are now on par with other schemes of order 3 and 4.

ARS222 ARS232 BPR442 PC432 SMR432 ARS233 ARS343 ARS443 BHR553 BPR533 BR343 CB443 CFN343 DBM553 KC443 LZ543 CB664 CFN564 KC664 KC774 LZ764 KC885
In summary, using this heuristic method of comparing orders of convergence estimated for a simplified, DAEfree problem against the full problem, we find that a significant fraction of the order reduction that impacts both differential and algebraic variables can be ascribed to the DAE.

An attempt to assess stiffness and its relationship with the observed order reduction

We now try to assess the level of stiffness of the problem we are interested in. As opposed to standard systems of ODEs for which the stiffness is set by means of an input control parameter, such as those analyzed e.g. by [START_REF] Boscarino | On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF], in our case it is the combination of the physical control parameters (Ra,Pr) and the spatial properties of the 2D mesh that sets the level of stiffness, and makes its definition less straightforward.

As seen above, the differential component of the problem at hand reads schematically

dx dt = N(x, z) + Lx, ( 44
)
where x is the differential state vector, and N and L are the nonlinear and linear operators, respectively. To estimate the stiffness we consider the ratio

(t) = τ L τ N (t) , ( 45
)
where τ N and τ L are the smallest time scales associated with the nonlinear and linear terms, respectively. For Boussinesq thermal convection, nonlinearities reflect the transport of momentum and heat by fluid flow, while linearities arise from the diffusion of those same fields. The nonlinear time scale τ N varies along the dynamical trajectory x(t). On the contrary, the linear time scale τ L corresponds to the most negative real eigenvalue of L, and is set once and for all upon prescription of the Prandtl number and the grid properties. A stiff situation occurs when τ L τ N , or 1, which is a strong incentive for an implicit treatment of the linear term Lx.

A first option to estimate τ N and τ L is to follow the logic of e.g. [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF] by writing in dimensionless form

τ N = min grid h s |u s | , h ϕ |u ϕ | , ( 46
)
τ L = Pr 1 + Pr min grid h 2 s , h 2 ϕ , (47)
where min grid is the minimum over the physical s -ϕ grid, and h s and h ϕ are the space-varying grid spacings in the s and ϕ directions, respectively. The nonlinear time scale is estimated based on a local measure of transport in the two directions of space, while the linear time scale assumes that the most negative eigenvalues of the linear operator correspond to an effective diffusivity equal to κ + ν, as suggested by Eq.( 24) of [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF].

We list the nonlinear and linear time scales determined for the 11 cases in the leftmost two columns of Table 4, noting that the value of τ N is, for each case, the average value of τ N (t) found for 5 independent snapshots. The corresponding stiffness parameter, τ , spans a moderate range of two decades. The increase of turbulence and reduction of τ N goes alongside an increase in the resolution that induces a concomitant decrease of τ L . We are in what authors currently refer to as an intermediate range of stiffness [START_REF] Kennedy | Higher-order additive Runge-Kutta schemes for ordinary differential equations[END_REF], that can indeed be detrimental to the order of convergence of some IMEX-RK methods.

Our second estimate of stiffness compares the maximum time-step allowed for stable computation using the ARS343 IMEX-RK method, ∆t imex max , with the maximum time step ∆t expl max allowed if one resorts to the fully explicit RK4 integrator. By considering the ratio

∆t = ∆t expl max ∆t imex max ,
we have an estimate of the stiffness based on an indirect probing of the stability regions of the timeschemes considered. As a matter of fact, the IMEX-RK scheme chosen here is ARS343 because the stability region of its explicit component matches by design that of RK4 (Ascher et al., 1997, §2.7); therefore we should not expect an offset of ∆t by an unwanted factor. Values of ∆t are tabulated alongside values of τ in Table 4. They fall within a factor of 3 within the values of τ , in a non-systematic way.

The third and final option we consider is to investigate directly the eigenvalues of the operators at hand. In the vicinity of a point x (t ), we approximate N by its tangent linear operator N such that

N(x + δx) = N(x ) + N (x )δx. ( 48
) Case τ N τ L τ ∆t λ
λ (red.) 0 9.14 10 -3 2.06 10 -6 2.22 10 -4 1.97 10 -4 2.84 10 -4 2.92 10 -4 1 8.92 10 -4 6.23 10 -7 6.99 10 -4 8.14 10 -4 1.69 10 -3 1.05 10 -3 2 5.99 10 -5 1.40 10 -9 2.34 10 -5 3.66 10 -5 3.50 10 -5 1.40 10 -3 3 1.55 10 -4 1.93 10 -7 1.25 10 -3 2.99 10 -3 4.28 10 -3 2.65 10 -3 4 1.85 10 -5 3.74 10 -8 2.02 10 -3 3.59 10 -3 5.93 10 -3 4.97 10 -3 5 1.07 10 -5 9.42 10 -9 8.78 10 -4 7.94 10 -4 1.36 10 -3 4.27 10 -4 6 4.10 10 -6 1.17 10 -8 2.86 10 -3 6.13 10 -3 9.18 10 -3 7.99 10 -3 7 1.39 10 -6 5.71 10 -10 4.10 10 -4 6.87 10 -4 6.98 10 -4 2.80 10 -4 8 7.42 10 -7 7.20 10 -10 9.71 10 -4 3.04 10 -3 3.54 10 -3 3.47 10 -3 9 2.26 10 -7 2.32 10 -10 1.03 10 -3 2.05 10 -3 1.78 10 -3 1.06 10 -3 10 1.35 10 -7 1.41 10 -10 1.05 10 -3 1.57 10 -3 3.82 10 -3 4.14 10 -3 Under these circumstances, if we define Q(x ) = N (x ) + L, the behavior of the solution in the vicinity of x

obeys dδx dt = Q(x )δx. ( 49
)
For each case, we constructed a two-dimensional, second-order finite-difference approximation of Q(x ), upon the Chebyshev-Fourier grid used in our spatial approximation. We did so in order to obtain a sparse, 2D, operator amenable to eigenvalue analysis for both full and reduced problems. The full problem, though, was approximated using a formulation based on the streamfunction alone (see e.g. Canuto et al., 2007, §1.4), in order to facilitate its implementation. We benchmarked our finite difference approximation by computing the critical parameters for convection, with a critical wavenumber m crit = 3 and Rayleign number Ra crit = 1768 (recall section 3.1). To obtain the eigenvalues of Q(x ), we resorted to the sparse library of the scientific python package Scipy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in python[END_REF], in conjunction with the SLEPc toolbox for python [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF][START_REF] Dalcin | Parallel distributed computing using Python[END_REF].

Cases of modest resolution lend themselves to the full calculation of the spectrum, and the example of case 3 is given in Figure 8. We observe that the distribution of eigenvalues is symmetrical with respect to the real axis, and that eigenvalues are concentrated in the vicinity of the origin, at the exception of a few purely real values that are quite distant, and correspond to the inverse value of the diffusive time scale on the smallest grid spacing, of non-dimensional value π 2 /N 4 s . These negative eigenvalues of large magnitude are due to the linear component of the problem at hand, and the ones responsible for stiffness. We associate the eigenvalues of largest imaginary part with the advective component of Q, and therefore the reciprocal value of τ N .

For cases of larger size (for case 10 the size of the sparse matrix to deal with is 878206 × 878206) we computed only the 100 eigenvalues of largest negative real parts, λ r , and 100 eigenvalues of largest imaginary parts, λ i . Our estimate of stiffness is given by in which () is the imaginary part. For this calculation, we used the same 5 independent snapshots for each case as the ones used to estimate τ . Results are listed in Table 4, and are in agreement, again within a factor of 3 with estimates based on τ and ∆t . The three estimates of point to case 6 as being the least stiff, even though the stiffness parameter does not exceed 10 -2 . Case 6 had a Prandtl number of unity, which implies that the diffusivities of heat and momentum are equal. The stiffest case is case 2, which has Pr = 40. For case 2, a relatively large resolution is needed to resolve the small-scale temperature anomalies within a pretty laminar flow (recall Figures 2a and2b). Table 4 also comprises values of λ computed for the reduced problem (advectiondiffusion equation of temperature with a prescribed flow), that we used previously to try and assess the impact of the DAE on the observed order reductions. The values of λ for the reduced problem are consistent with those found for the complete problem, with one exception. In this simplified system, case 6 remains the least stiff case, and it is now the barely supercritical case 0 that is the stiffest case. Case 2 ceases to be the stiffest case, since its most negative eigenvalues are associated with the diffusion of momentum, not heat. Accordingly, the value of λ we find for the reduced case 2 is precisely multiplied by a factor of Pr = 40, compared with the value found for the full problem. This increase is not seen in the three cases that have Pr below unity (5, 7 and 9), for which the largest negative eigenvalues remain connected with the diffusion of heat. We conclude the analysis of the impact of stiffness by showing the measured order of convergence for the reduced, DAE-free, and full problems against the corresponding values of λ in figure 9, with the hope that this representation will get rid of the jaggedness of figures 5 (left panel) and 7. Within the modest range of λ that our investigations enabled, we observe in Fig. 9, left panel, that for the reduced problem, KC774, CB664 and KC885 are close to meeting their expected order of convergence for λ 10 -3 . SMR432 exhibits an order of convergence larger than 2 for λ 10 -3 as well. ARS343

λ = max (λ i ) max | (λ r )| , (50)
superconverges only for case 10, while BHR553 superconverges for all cases but case 0. This scheme is by design supposed to be immune to order reductions caused by stiffness and the DAE [START_REF] Boscarino | On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF]. For a function of the stiffness parameter λ estimated for each problem for all 11 cases and for some selected time schemes. See text for details. λ sur the full problem, the range of λ covered is a bit broader, and we find it noteworthy that for the stiffest case 2, order is restored. In fact, the curves we obtain in particular for schemes KC774 and KC885 show a reasonable similarity with the curves obtained for those same schemes applied to Kap's problem by their genitors [START_REF] Kennedy | Higher-order additive Runge-Kutta schemes for ordinary differential equations[END_REF].

In summary, we think that the stiffness of the cases we considered in this study covers a moderate, intermediate range of values spanning slightly less than two decades for the reduced problem and two and a half decades for the full Boussinesq convection problem. We find that even if stiffness has an impact on the convergence of some schemes over a limited range, typically 10 -4 λ 10 -3 for the stiffness parameter we determined, it is mostly the DAE that causes the degradation of convergence. The reduction in order affects both differential and algebraic variables, probably by virtue of the coupling induced by the equations of the problem at hand. Yet, schemes of theoretical order 2 are not affected by order reduction. Higher-order time integrators that are by design immune to such problems, such as BHR553 or the IMEX multistage methods, may appear as the schemes of choice.

This statement remains to be weighted against a measure of the stability and computational efficiency of those schemes, which is the topic of the next subsection.

Stability and computational efficiency

The explicit treatment of nonlinearities imposes a restriction on the available time step that is subject to a ν (t) is a reference time series for viscous dissipation obtained with the SBDF4 time scheme. Note that the scale on the y-axis is logarithmic. For a scheme, the maximum admissible value of α cfl is the largest one that precedes the steep increase in the curve.

where the O(1) prefactor α cfl depends on the time integrator and the case considered. In order to determine empirically the maximum admissible value of α cfl , α max cfl , for the 11 × 26 combinations of this study, we followed Gastine (2019) and inspected the timeseries of viscous dissipation D ν (t) and its fluctuations for different values of α cfl . The maximum admissible value of α cfl , α max cfl is determined to 0.02 accuracy by requesting that the timeseries of D ν (t) does not exhibit any flagrant spike, over a time window of width roughly equal to 5 convective turnover times; the time window is case-dependent, but for a given case, it is the same for all time integrators. Figure 10a illustrates this methodology for case 9 and the ARS343 scheme, a configuration for which we find α max cfl = 1.10. This arguably tedious methodology is meant at preserving the accuracy of the solution, and can not lead to the disturbing occurrence of stable yet inaccurate IMEX-RK schemes reported by [START_REF] Grooms | Linearly implicit methods for nonlinear PDEs with linear dispersion and dissipation[END_REF], and mentioned in the introduction. In this regard, our estimate of α max cfl is conservative. We refer readers interested in a more standard assessment of stability and efficiency to appendix Appendix E, where we report accuracy versus runtime for cases 2 and 10.

We now propose an automated way of reaching the same conclusions. We begin by establishing a master curve for D ν (t) over the interval of interest using the SBDF4 time scheme and the smallest α cfl . This master curve is denoted by D r ν (t) where again, the superscript r stands for reference. Given the D ν (t) computed for an integrator and a value of α cfl , we evaluate the time average of |D ν -D r ν | over the window of interest using splines for the numerical integration. As an example, we show |D ν -D r ν | in Figure 10b for case 9 and 6 schemes. We observe a sharp transition in the behavior of this quantity, for relatively low values of α cfl for the two multistep schemes (CNAB2 and SBDF3) and larger values for the IMEX-RK schemes. The largest value of α cfl before the transition matches the α max cfl obtained by visual inspection of the timeseries of D ν (t). The value of α max cfl can be converted into a maximum attainable value of the time step size, ∆t max : we take it to

ARS222 ARS232 BPR442 PC432 SMR432 ARS233 ARS343 ARS443 BHR553 BPR533 BR343 CB443 CFN343 DBM553 KC443 CB664 CFN564 KC664 KC774 LZ764 KC885 CNAB2 SBDF2 SBDF3 SBDF4
Figure 11: Ratio of the maximum ordinate of the stability curve of the explicit component of a time integrator to that of the CNAB2 scheme as a function of the ratio of the maximum timestep size admissible to the maximum time step admissible for CNAB2. The case considered is case 10 (Re = 13320, stiffness parameter ∼ 3 × 10 -3 ). The dashed grey line is the prediction based on Eq. ( 52).

be the average ∆t obtained in the same setups that led to the determination of α max cfl , meaning that we find 26 values of ∆t max (one per scheme) per case. In those cases where nonlinearities dominate, with stiffness parameters of the order of 10 -2 (recall Table 4 above), it is actually possible to make a decent guess of ∆t max based on the boundary of the stability region of the explicit component of the scheme under scrutiny. This stability region is bounded by a curve that mostly lies in the left hand side of the complex plane, that of negative real parts. We anticipate that in moderately stiff situations where nonlinearities play a major part in the dynamics, it is this curve that will control the stability of the implicit explicit scheme, even if it does not correspond to the stability curve of the combined scheme, as studied by e.g. [START_REF] Karniadakis | High-order splitting methods for the incompressible Navier-Stokes equations[END_REF] and [START_REF] Izzo | Highly stable implicit-explicit Runge-Kutta methods[END_REF] for test problems. As discussed above in our analysis of stiffness, transport will effectively probe the eigenvalues of the semi-discrete tangent linear operator with the largest imaginary parts (in absolute value). Let λ m denote the eigenvalue of largest imaginary part. As a rule of thumb we want λ m ∆t to be close to the stability curve. For transport-dominated physics, we thus make the tentative prediction that for any scheme

∆t max (scheme) = f × y(x m ) (scheme) , (51)
where the factor f is a function of the spatial discretization only, x m = (λ m ∆t) and the ordinate y(x m ) on the stability curve can be determined numerically. These curves are provided in Appendix D for completeness. If this equality holds, provided we know ∆t max of a case for one scheme (CNAB2, say) we may expect that

∆t max (scheme) = y(x m )(scheme) y(x m )(CNAB2) ∆t max (CNAB2) . ( 52
)
The relevance of this line of reasoning is shown in Figure 11 where this empirical prediction is compared with the measured value for case 10, our most turbulent case. We find that the overall trend is that the prediction slightly overestimates the actual value, typically by within 10 to 20 %. Figure 11 highlights the fact that the most stable order 3 scheme is DBM553, a result that can be understood by inspection of the stability region of its explicit component given in Figure D.15, middle panel. Of all the third order schemes we considered, DBM553 has the most elongated region of stability in the vicinity of the y-axis of the complex plane. In fact, it was designed for the very purpose of accommodating the constraint on the available time step size arising from the location of eigenvalues of the explicit component of a model being located along the imaginary axis [START_REF] Kinnmark | One step integration methods of third-fourth order accuracy with large hyperbolic stability limits[END_REF]. A closer inspection of the IMEX-RK schemes which share the same stability domain for their explicit component reveals very similar CFL coefficients for case 10 with for instance α max cfl (ARS232) = 0.54 and α max cfl (SMR432) = 0.56; or α max cfl (ARS343) = 0.8 and α max cfl (BR343) = 0.82. This is another indication that the stability domain of the explicit part only provides a decent estimate of the actual stability of an IMEX-RK scheme in the limit of advection-dominated flows. To conclude this paragraph, we note that for cases of more dramatic stiffness, one should probably consider the stability region of the complete IMEX scheme, an endeavor that we did not pursue.

To evaluate the efficiency of a given scheme, we consider the following ratio

eff = α max cfl cost , (53)
where cost refers to the amount of work required to advance the solution by one time step ∆t. In practice cost is the average cpu time measured over 1000 iterations using reproducible runtime conditions (same compute nodes, exclusive access to the compute node, one single OpenMP thread). Note that the linear matrix solves amount for the majority of the walltime, and hence the number of solves per iteration n I in Table 1 provides a decent estimate of the actual relative walltime. We evaluated the efficiency of the 26 schemes considering the 11 cases.

We investigate how the gain one may obtain in terms of a larger α max cfl using an IMEX-RK scheme trades off with the extra operations that are needed. (Again, at this stage, we do not consider the benefit in terms of accuracy.) Given the popularity of the CNAB2 integrator in our community, we normalize the efficiency by the efficiency of CNAB2.

Results are displayed in Figure 12 for cases, 2, 5 and 10 that have Re = 26, 513 and 13320, respectively, and whose stiffness parameter we estimated in section 3.5 to be ∼ 3 × 10 -5 , ∼ 10 -3 and ∼ 3 × 10 -3 , respectively. Relative efficiency is in all cases bounded between 0.5 and 1.50. We observe that the number of IMEX-RK schemes that outperform CNAB2 increases dramatically with the Reynolds number from 5 for Case 2 to 14 for case 10. For the latter, the 14 schemes comprise 3, 7 and 4 schemes of order 2, 3 and 4, respectively. This indicates that the gain in stability in transport-dominated regimes outweigh the increase in the number of operations as the resolution increases. For case 10, we note a relative grouping of multistep schemes around 1, with SBDF3 being slight more efficient than CNAB2, in agreement with its more elongated stability domain close to the imaginary axis (Appendix D). A scheme that appears to be consistently efficient across the 3 cases is ARS343. The DBM553 scheme by [START_REF] Vogl | Evaluation of implicit-explicit additive Runge-Kutta integrators for the HOMME-NH dynamical core[END_REF] and the BR343 scheme by [START_REF] Boscarino | On the uniform accuracy of IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF] transition from a poor efficiency in the laminar case 2 to an excellent one in the turbulent case 10.

Trade-off between accuracy and efficiency

We now weigh efficiency and accuracy for recommendations to be made to the practitioner. The accuracy is characterized by the error one expects for a simulation performed at the maximum available time step size, ∆t max , as defined in the previous section. This error is computed based on the scaling of the error on temperature with ∆t that can be obtained for the various convergence curves introduced in section 3.2 above, by taking ∆t = ∆t max . Therefore, our target practitioner wishes to integrate the solution for the longest time span possible given his/her computing resources with the hope that the error will remain as small as possible. We normalize the error with the error estimated for CNAB2 and show in Figure 13 where the various schemes considered in this work are located in the (relative efficiency, relative error) plane for cases 2, 5 and 10 again.

Schemes located in the top right quadrant of each panel are more accurate and yet more efficient than CNAB2.

Six IMEX-RK schemes are systematically located in this quadrant, one scheme of order 2, SMR432 [START_REF] Spalart | Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions[END_REF], four schemes of order 3, ARS343 [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF], CB443 [START_REF] Cavaglieri | Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff highdimensional ODE systems[END_REF], CFN343 [START_REF] Calvo | Linearly implicit Runge-Kutta methods for advection-reaction-diffusion equations[END_REF], KC443 [START_REF] Kennedy | Additive Runge-Kutta schemes for convection-diffusion-reaction equations[END_REF] and one scheme of order 4, KC664 [START_REF] Kennedy | Additive Runge-Kutta schemes for convection-diffusion-reaction equations[END_REF]. At the operational, dissipation-based limit of stability, it is remarkable to note that ARS343 and KC664 enable a significant improvement of accuracy at a lower cost, regardless of the configuration. In the turbulent regime, where the explicit components of time integrators prevail, CNAB2 is outperformed by no less than 14 schemes, which include 2 schemes of order 2, 8 schemes of order 3 (including SBDF3), and 4 schemes of order 4.

Summary and recommendations

We have applied 22 multistage, 4 multistep implicit-explicit, and 2 fully explicit integrators to the problem of Boussinesq thermal convection in a two-dimensional cylindrical annulus, over a broad range of regimes whose exploration was made possible by the definition of 11 different physical setups. Our spatial discretization rests on a pseudo-spectral collocation method applied to the vorticity streamfunction formulation of the problem. We summarize our findings as follows:

• Over the range of cases considered, IMEX multistep methods exhibit their expected order of convergence.

IMEX-RK methods of second order also show the expected convergence. Some IMEX-RK methods of order 3 and 4 show order reduction, that affect both the algebraic and differential variables of the chosen formulation.

• We have attempted to evaluate the stiffness of the cases using three different strategies that lead to the same conclusion: the small parameter that quantifies stiffness, , covers a moderate and intermediate range of values spanning two and a half decades for the full Boussinesq convection problem, between 3 × 10 -5 and 10 -2 .

• We find that even if stiffness has an impact on the convergence of some schemes over a limited range of , typically 10 -4 10 -3 , it is mostly the discrete algebraic equation that causes the degradation of convergence. The reduction in order affects both differential and algebraic variables, probably by virtue of the coupling induced by the equations of the problem at hand.

• IMEX-RK time integrators that are by design immune to such problems, such as BHR553 [START_REF] Boscarino | On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation[END_REF], exhibit nominal convergence properties, or even better-than-nominal properties in the least stiff (transport-dominated) situations.

• We have defined the efficiency of a scheme by the ratio of the maximum admissible value of the Courant number, α max cfl , divided by the cost of performing one time step. By maximum admissible value we understand a value that generates a smooth timeseries of viscous dissipation within the system. Viscous dissipation is a demanding quantity, whose behavior help us define an acceptable (or not acceptable) solution. With this at hand, we have reported the efficiency of the schemes relative to the popular CNAB2 for 3 cases that we consider representative. We found that the relative efficiency was bounded between 0.5 and 1.5. Also, CNAB2 is not the method of choice when going to transport-dominated cases, as no less than 14 schemes are more efficient than CNAB2 for our most turbulent case.

• This last statement becomes even stronger when relative accuracy is added to the analysis. By relative accuracy here we mean the ratio of the error anticipated for a scheme running at the operational, dissipationbased limit of stability, whose expected error can be estimated based on the convergence analysis, to the same error anticipated for CNAB2 at its limit of stability (and that for each case). For the same three representative cases, we find that 6 schemes combine the assets of being less expensive and more accurate than CNAB2: SMR432 [START_REF] Spalart | Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions[END_REF], CFN343 [START_REF] Calvo | Linearly implicit Runge-Kutta methods for advection-reaction-diffusion equations[END_REF], ARS343 [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF], CB443 [START_REF] Cavaglieri | Low-storage implicit/explicit Runge-Kutta schemes for the simulation of stiff highdimensional ODE systems[END_REF], KC443 and KC664 [START_REF] Kennedy | Additive Runge-Kutta schemes for convection-diffusion-reaction equations[END_REF].

• For the problem of thermal convection in a cylindrical annulus with a spectral discretization in space, it appears that the default integrator should be the third-order ARS343, or possibly KC664 if higher accuracy is sought.

For turbulent cases, or in more general terms, transport-dominated cases, the performance of an implicitexplicit scheme is dictated by its explicit component, as the part treated explicitly here is purely advective. The previous general recommendation can be amended on a case-by-case basis: for instance the third-order scheme DBM553 proposed by [START_REF] Vogl | Evaluation of implicit-explicit additive Runge-Kutta integrators for the HOMME-NH dynamical core[END_REF] may well prove superior to any third-order scheme for turbulent transportdominated 2D problems reaching Re ∼ 10 5 and beyond. In fact, on this path, provided ∼ 1 is reached, RK4 may well prove competitive, even more so if a regridding is performed to alleviate the stringent time step limitations due to the clustering of grid points near the boundaries, as done e.g. by [START_REF] Johnston | Comparison of turbulent thermal convection between conditions of constant temperature and constant flux[END_REF] using the mapping proposed by [START_REF] Kosloff | A modified Chebyshev pseudospectral method with an o(n -1) time step restriction[END_REF].

Before discussing a tentative extrapolation of our results to three-dimensional geometry, we should add the following important caveat: stiffness in our setup is controlled by the most negative real eigenvalues of the linear part that is treated implicitly, and our recommendations may not be suited for those problems where stiffness comes from fast waves, i.e. large imaginary eigenvalues of the linear terms. We may wonder how our findings can carry over to these simulations. In Appendix F, we attempt to convert our 2D tradeoff diagrams into their 3D counterparts, by making the strong assumptions that accuracy and stability remain the same, i.e. that the convergence properties and values of α max cfl are not affected. The sole impact of the change of geometry that factor in the analysis lies in the cost, where we acknowledge that in a three-dimensional pseudo-spectral code of planetary core dynamics, explicit stages are the most expensive steps in a calculation, as there is no efficient equivalent of the fast Fourier transform for the Legendre transform. Taking this into account, non globally stiffly accurate IMEX-RK schemes, that require an extra assembly stage, are penalized compared with stiffly accurate ones.

For the range of Re typical of recent parametric studies, and under the strong assumptions that we made, we predict that the schemes of choice for Re ∼ 1000 in three dimensions should be PC432 (order 2), ARS343 and BPR533 (order 3) and KC664 (order 4). Relative to the classic CNAB2, these schemes should yield a reduction of the time to solution while enabling more accurate solutions. This analysis ignores the memory imprint of the schemes (recall Tab. 1) , as the codes in question are massively parallel and should be immune to this issue, at least for the level of Re considered. If the level of turbulence should increase for 3D simulations, then we anticipate again that it is the robustness and efficiency of the explicit component of IMEX schemes that will matter. Also, depending on the problem studied, hybrid splitting strategies, that would treat explicitly some linear terms with the aim of increasing the efficiency of the calculation, look like an interesting avenue for investigation, especially in turbulent situations. Also worth investigating are the IMEX general linear method (GLM) that have recently come to the fore. They appear to overcome some of the inherent limitations of IMEX-RK schemes, such as order reduction. [START_REF] Zhang | High order implicit-explicit general linear methods with optimized stability regions[END_REF] for instance explored several IMEX strategies to compute 2-D and 3-D simulations of thermal rising bubbles [START_REF] Giraldo | A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases[END_REF] and showed that the IMEX-GLM were immune to order reduction and exhibited the best accuracy/efficiency tradeoff among the tested schemes (see their Fig. 10). We conclude by hoping that our findings will serve as an incentive for the community of stars and planetary fluid interiors modelers to transition towards IMEX-RK integrators, which possess the extra convenient property of being self-restarting (their initialization only requires knowledge of the current state vector), regardless of their order of accuracy.

  function, y ψ,m . Each quadrant in the above equation has size N s × N s . The top-left quadrant originates from the time discretization of the vorticity equation. The right-hand side r.h.s ω contains a mix of linear and nonlinear components, including the contribution of the updated temperature field via the buoyancy term. The constraint (30d) is enforced by means of the second row of the block-matrix system. The boundary conditions are enforced

Figure 1 :

 1 Figure 1: Time averaged kinetic energy versus azimuthal wavenumber m for the 11 configurations considered in this study. Every third mode shown for clarity (m = 0, 3, 6, . . . ) for cases 0, 1, 3, 4 and 6. The scale on both axes is logarithmic. Note that cases 0, 1 and 3 have zero energy in the axisymmetric m = 0 mode.

Figure 3 :

 3 Figure 3: Convergence of the L 2 error for the temperature field (top panels), the cylindrical radial velocity u s (middle panels) and the vorticity ω (bottom panels) for Case 2 (left column) and Case 10 (right column) as a function of the timestep size ∆t. The markers correspond to the class of IMEX, with squares denoting IMEX multistep and circles IMEX-RK multistage schemes. The linestyles highlight the theoretical orderwith dashed lines for second order, dotted lines for third order, solid lines for fourth order and dash-dotted lines for fifth order.

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Measured order of convergence for the 11 cases and 22 IMEX-RK time integrators, based on the error impacting the temperature T (left panel), cylindrical radial velocity u s (middle panel) and vorticity ω (right panel). Thick horizontal lines highlight integer values of 2, 3, 4 and 5. The 22 schemes are listed to the right, and their color reflects the theoretical order of convergence.

Figure 8 :

 8 Figure 8: Eigenvalues λ of the tangent linear operator computed for case 3, for the full problem and the reduced problem restricted to the heat equation considered alone. Left: full spectrum. Right: zoom in the vicinity of the origin, in the region of the complex plane defined by the dashed rectangle in the left panel.

Figure 9 :

 9 Figure 9: Order of convergence for the temperature field T , in the reduced problem set-up (left panel), and for the full problem (right panel), as

Figure 10

 10 Figure 10: a: Time series of viscous dissipation |D ν | for case 9 advanced with the ARS343 IMEX-RK scheme, considering various values of α cfl . The highest admissible value we find in this setup is α cfl = 1.10. b: Time average of |D ν -D rν | as a function of α cfl for case 9 and 6 time integrators. D r ν (t) is a reference time series for viscous dissipation obtained with the SBDF4 time scheme. Note that the scale on the y-axis is logarithmic. For a scheme, the maximum admissible value of α cfl is the largest one that precedes the steep increase in the curve.

Figure 12 :Figure 13 :

 1213 Figure 12: Efficiency of the time integrators relative to the efficiency of CNAB2 for cases 2, 5 and 10 from left to right. Horizontal dashed lines correspond to a value of unity. See text for details.

  Three-dimensional problems concerned with the modeling of planetary interiors in global spherical geometry are nowhere near reaching values of Re as extreme as 10 6 . Recent parametric studies of convection-driven dynamo in a spherical shell geometry by Schwaiger et al. (2019), Gastine et al. (2020) and Tassin et al. (2021) have values of Re in the range 100 -3000, while single-case, rapidly-rotating trophy studies by Schaeffer et al. (2017); Sheyko et al. (2018) (DNS) and Aubert (2019) (LES) report Re ∼ 5000 and Re ∼ 20000 for the DNS and LES, respectively.

Figure

  Figure C.14: Convergence of the L 2 error for the temperature field (left panel), the cylindrical radial velocity u s (middle panel) and the vorticity ω (right panel) for Case 3. In addition to the usual 26 IMEX schemes, this figure also features the fully explicit RK2 and RK4 methods, whose error levels are marked by large triangles which appear in the bottom left corner of each panel.

  Figure E.16: Convergence of the L 2 error for the temperature field (top panels), the cylindrical radial velocity u s (middle panels) and the vorticity ω (bottom panels) for Case 2 (left column) and Case 10 (right column) as a function of total runtime expressed in seconds. The markers correspond to the class of IMEX, with squares denoting IMEX multistep and circles IMEX-RK multistage schemes. The total runtime is the product of the number of iterations times the average walltime, and ignores the initial computation and factorization of the requisite matrices.

  1/2 . The right-hand side of Eqs. 30a-30c comprises nonlinear and linear terms. The nonlinear terms are vectors of size N s denoted by N; N u ϕ , N ω,m , and N T,m respectively contain the values of -u s ω, -∇ m • (uω) m and -∇ m • (uT ) m

Table 1 :

 1 Multistage IMEX-RK methods used in this study. The leftmost "Scheme" column defines the scheme notation, as used in the main text, tables and figures. K I is the number of stages of the diagonally implicit Runge-Kutta (DIRK) component. K E is the number of stages of the explicit Runge-Kutta (ERK) component. The next column contains the expected order of accuracy o of the combined IMEX-RK scheme.n I is the number of linear solves for each time-step, which differs from K I if a I k1 0 ∀k. The star indicates that the scheme involves several matrices because of the changes on the diagonal of the implicit Butcher table (non S-DIRK schemes). S. A. indicates whether the ERK or DIRK part of the method is stiffly accurate, and b I = b E indicates if the DIRK and ERK methods have the same solution weights to compute the assembly. Storage denotes the number of state vectors that need to be stored simultaneously for the time advance of one physical quantity.

Table 2 :

 2 Nusselt number Nu obtained for Pr = 0.025, s i /s o = 0.3 and an increasing Rayleigh number Ra, by[START_REF] Alonso | Transition to temporal chaos in an O(2)-symmetric convective system for low Prandtl numbers[END_REF] (the reference) and with the code developed for this study.

	Ra	Nu -1	Nu -1
		(ref.)	(this study)
	1892 0.005	0.005
	2510 0.163	0.162
	3268 0.383	0.383
	4013 0.544	0.544
	4106 0.562	0.561
	4500 0.617	0.618
	5000 0.679	0.678
	5500 0.733	0.733
	6000 0.783	0.783
	6500 0.827	0.827
	7000 0.871	0.869

Table 3 :

 3 Properties of the 11 convection cases investigated in this study. From left to right: Case number, Prandtl number (input), Rayleigh number (input), Reynolds number (output), Nusselt number at the outer boundary (output), time average buoyancy input power (output), time average heat loss by viscous dissipation (output), and spatial resolution used.

	Case Pr	Ra	Re	Nu o	P	D ν	(N s , N m )
	0	1	2 ×10 3	2.87	1.16 2.03 × 10 3	-2.03 × 10 3	(36, 36)
	1	1	1 ×10 4	18.85	2.51 9.29 × 10 4	-9.29 × 10 4	(48, 48)
	2	40	1 ×10 7	26.00 12.63 4.39 × 10 5	-4.39 × 10 5 (256, 256)
	3	1	1 ×10 5	77.33	4.64 2.22 × 10 6	-2.22 × 10 6	(64, 64)
	4	1	1 ×10 6	279.76	7.70 4.06 × 10 7	-4.06 × 10 7	(96, 128)
	5	0.025 1 ×10 4	513.44	2.07 1.07 × 10 8	-1.07 × 10 8	(64, 192)
	6	1	1 ×10 7	943.12 13.17 7.33 × 10 8	-7.33 × 10 8 (128, 160)
	7	0.025 1 ×10 5	2023.52	3.97 2.89 × 10 9	-2.89 × 10 9 (128, 320)
	8	1	1 ×10 8	3462.47 23.30 1.34 × 10 10 -1.34 × 10 10 (256, 256)
	9	0.025 1 ×10 6	6835.69	6.56 5.36 × 10 10 -5.37 × 10 10 (160, 384)
	10	1	1 ×10 9 13320.12 44.22 2.60 × 10 11 -2.60 × 10 11 (384, 384)

Table 4 :

 4 Three different estimates of the stiffness of the problem at hand, for all cases considered in this study. The rightmost column gives an estimate of stiffness for the problem reduced to the advection-diffusion equation for temperature. Bold face fonts used for largest and smallest values. See text for details.

Table A .

 A 5: Coefficients of the multistep schemes considered in this study when a fixed time step size ∆t is employed.

Present address: Bosch Engineering and Business Solutions, India. e-mail: gopinath.venkatesh2@in.bosch.com

e-mail: gastine@ipgp.fr

Acknowledgments

We thank the two anonymous reviewers whose comments and suggestions helped improve this paper. This work was made possible by support from Indo-French Centre for the Promotion of Advanced Research CE-FIPRA/IFCPAR (project 5307-1). Near the completion of the writing of this manuscript we were saddened to learn the passing of Ms. A. Sathidevi, senior scientific officer at CEFIPRA that had accompanied us over the course of this project, and we would like to dedicate this work to her memory. We also acknowledge fruitful interactions over the long course of this work with Prof. H. Johnston, Dr. P. Livermore, Dr. L. Métivier, Prof. D. Reynolds, Dr. N. Schaeffer and Prof. J.-P. Vilotte. Numerical computations were performed on the S-CAPAD platform, IPGP, France.

Convergence analysis

Our analysis of the convergence of the 26 schemes of interest in this study follows this procedure: for each case, we compute a reference simulation using the 4th order SBDF4 IMEX multistep scheme, using a time-step size ∆t r small enough in order to enable a convergence analysis that spans two orders of magnitude in terms of ∆t.

To equilibrate the solution prior to using ∆t r , we activated the possibility of a variable ∆t for the high-resolution cases. We select a time window [t s , t e ] that typically covers a sizeable fraction of a convective turnover time. The reference state vector at t = t s is taken as the initial condition for the forward integration up to t = t e , performed with each of the 26 schemes. The accuracy of the solution at t = t e is assessed using the L 2 norm. For instance, the error in ω is given by e ω = A [ω(s, ϕ, t e ) -ω r (s, ϕ, t e )] 2 sdsdϕ, where the superscript r corresponds to the reference solution. In the following, unless otherwise stated, we will systematically use this absolute definition of error.

We begin by a global inspection of the error behavior for the two cases we already looked at, cases 2 and cases 10, whose convergence results are shown in Figure 3. As explained above, we tried to assess the convergence properties by having at least two orders of magnitude in the range of ∆t; this is sometimes barely achieved, in particular for IMEX multistep methods SBDF3 and SBDF4 whose stability domain is narrower than IMEX-RK schemes of the same order. Schemes of nominal order 2 systematically display a higher error level than schemes of order 3 and beyond, and never reach the plateau of numerical roundoff error in the range of time step sizes that we considered. Two exceptions are the ARS232 scheme by [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations[END_REF] and SMR432 scheme by [START_REF] Spalart | Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions[END_REF] that exhibit a higher convergence rate; their convergence curves are in fact mixed with those of the schemes of nominal order 3. These schemes aside, we also note that at any given ∆t there can be a factor of 10 difference between the worst (in this sense) order 2 scheme and the best one -PC432 is more accurate by one order of magnitude than SBDF2 or CNAB2. Order 3 schemes find themselves sandwiched between order 2 and order 4 schemes. Their convergence rate is such that the roundoff error plateau is reached for some schemes, for all fields (T , u s and ω) for case 2 and all fields but the vorticity for case 10. BR343, BHR553 and ARS343 appear as the most accurate third-order IMEX-RK schemes, especially in the turbulent case 10. For the latter (ARS343) this makes sense as its explicit component is designed to match the stability properties of the RK4 scheme (Ascher et al., 1997, §2.7); see also Appendix D. For case 2, comparison of the behavior of e u s of IMEX-RK schemes of theoretical order 3 with that of SBDF3 highlights that some do not exhibit third-order accuracy. The 4th-order IMEX-RK schemes display overall similar error levels, below the lowest level attained by third-order schemes, this being marginally true for CFN564. For a given ∆t, if one considers the temperature T (Figure 3, top panel) 4thorder IMEX-RK schemes are more accurate by two orders of magnitude than the SBDF4 multistep scheme. The situation is not so clear when one considers the error in the vorticity. There SBDF4 displays a high convergence rate towards the plateau, whereas 4th order IMEX-RK schemes do not show a clear trend. In fact, the sole scheme that appears to compete with SBDF4 is BHR553. We will return to this later. For now, we complete this preliminary overview by noting that our sole 5th order scheme, KC885, is as expected more accurate than any other scheme considered, with the exception of SBDF4 and BHR553 being more accurate with regard to the vorticity for case 2, over a limited range of ∆t. Note finally that the fully explicit schemes that we have at our disposal (RK2 and RK4)

Appendix A. Multistep schemes

We provide the reader with the vectors of coefficients a, b and c for CNAB2, SBDF2, SBDF3 and SBDF4.

These vectors define a linear multistep method with K steps according to

where c -1 0. Table A.5 enlists these three vectors when a fixed time step size ∆t is used (see e.g. Peyret, 2002, Table 4.4).

Following [START_REF] Wang | Variable step-size implicit-explicit linear multistep methods for time-dependent partial differential equations[END_REF], those vectors can be generalised to the case of variable time step sizes. In the following, we define

and the ratio

For the different IMEX multistep schemes considered here, we then obtain:

where the three constants c 1 , c 2 and c 3 are expressed by

Appendix B. Butcher tableaux of PC432

The PC432 time scheme is assembled using the explicit scheme from [START_REF] Jameson | Numerical solution of the Euler equations by finite volume methods using Runge Kutta time stepping schemes[END_REF] for its explicit component and a Crank-Nicolson scheme for its implicit part. This is a stiffly accurate second order three stage scheme and its Butcher tableaux read We provide in this appendix additional elements to assess the efficiency of the 26 schemes of interest in this study, providing metrics that may be more general than the dissipation-based efficiency introduced in Section 3.6. we can stress that ARS343 appears as a good choice for case 2 and case 10. To obtain higher accuracy with a concomittent moderate increase in the computational cost, SBDF4 (for case 2) and KC664 (for case 10) should be of an explicit state, with respect to the stiffly-accurate schemes. In contrast, the schemes with a lower number of explicit stages, such as BPR533, present an enhanced efficiency compared to the 2-D computations. PC432 stands out as the most efficient scheme for cases 2 and 5, while ARS343, BR343 and DBM553 become more efficient in the least stiff case 10. We anticipate that the overall gain in terms of efficiency compared to the CNAB2 method for 3-D calculations will be smaller than for our 2-D models, bounded to values between 20 and 30%. tions. Focusing our attention on the upper right quadrant, one second order scheme (SMR432), three third order schemes (ARS343, BPR533 and CFN343) and one fourth order scheme (KC664) are always more efficient and more accurate than CNAB2 at the stability limit.

We stress that the reasoning put forward in this Appendix heavily relies on the assumptions that both the stability coefficients and the convergence curves are weakly affected by the change from 2-D to 3-D models. While this is a plausible hypothesis when considering the same physical phenomenon (i.e. non-rotating convection), the incorporation of additional physical effects such as rotation or a magnetic field is likely to significantly impact the convergence and the stability properties of the time integrators.