# Rescaled equations for well-conditioned direct numerical simulations of rapidly rotating convection

Keith Julien, Adrian van Kan, Benjamin Miquel, Edgar Knobloch, Geoffrey Vasil

► To cite this version:

Keith Julien, Adrian van Kan, Benjamin Miquel, Edgar Knobloch, Geoffrey Vasil. Rescaled equations for well-conditioned direct numerical simulations of rapidly rotating convection. Journal of Computational Physics, 2025, 541, pp.114274. 10.1016/j.jcp.2025.114274. hal-05394842

HAL Id: hal-05394842

https://hal.science/hal-05394842v1

Submitted on 3 Dec 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Rescaled Equations for Well-Conditioned Direct Numerical
Simulations of Rapidly Rotating Convection

Keith Julien<fnref n="1" />ᵃ, Adrian van Kan<fnref n="2" />ᵇ, Benjamin Miquelᶜ, Edgar Knoblochᵇ, Geoffrey Vasilᵈ

**Abstract**

Convection is a ubiquitous process driving geophysical/astrophysical fluid flows, which are typi-
cally strongly constrained by planetary rotation on large scales. A celebrated model of such flows,
rapidly rotating Rayleigh-Bénard convection, has been extensively studied in direct numerical simu-
lations (DNS) and laboratory experiments, but the parameter values attainable by state-of-the-art
methods are limited to moderately rapid rotation (Ekman numbers $Ek \gtrsim 10^{-8}$), while realistic
geophysical/astrophysical $Ek$ are significantly smaller. Asymptotically reduced equations of mo-
tion, the nonhydrostatic quasi-geostrophic equations (NHQGE), describing the flow evolution in
the limit $Ek \to 0$, do not apply at finite rotation rates. The geophysical/astrophysical regime of
small but finite $Ek$ therefore remains currently inaccessible. Here, we introduce a new, numerically
advantageous formulation of the Navier-Stokes-Boussinesq equations informed by the scalings valid
for $Ek \to 0$, the *Rescaled Rapidly Rotating incompressible Navier-Stokes Equations* (RRRiNSE).
We solve the RRRiNSE using a spectral quasi-inverse method resulting in a sparse, fast algorithm
to perform efficient DNS in this previously unattainable parameter regime. We validate our results
against the literature across a range of $Ek$, and demonstrate that the algorithmic approaches taken
remain accurate and numerically stable at $Ek$ as low as $10^{-15}$. Like the NHQGE, the RRRiNSE
derive their efficiency from adequate conditioning, eliminating spurious growing modes that oth-
erwise induce numerical instabilities at small $Ek$. We show that the time derivative of the mean
temperature is inconsequential for accurately determining the Nusselt number in the stationary

<footnote>ªDepartment of Applied Mathematics, University of Colorado Boulder, Boulder, Colorado 80309, USA</footnote>

<footnote><sup>b</sup>Department of Physics, University of California at Berkeley, Berkeley, California 94720, USA</footnote>

<footnote>cCNRS, Ecole Centrale de Lyon, INSA Lyon, Univ. Claude Bernard Lyon 1, LMFA, UMR5509, 69130
Ecully, France</footnote>

<footnote>dSchool of Mathematics and the Maxwell Institute for Mathematical Sciences, James Clerk Maxwell Building, University of Edinburgh, Edinburgh, EH9 3FD, UK</footnote>

<footnote n="1">¹Deceased on 14th April 2024.</footnote>

<footnote n="2">²Corresponding Author Email: avankan@berkeley.edu</footnote>
state, significantly reducing the required simulation time, and demonstrate that full DNS using RRRiNSE agree with the NHQGE at very small Ek.

Keywords: Rapidly rotating Rayleigh-Bénard convection, rescaled Navier-Stokes equations, asymptotically reduced equations, quasi-inverse method

2010 MSC: 00-01, 99-00

# 1. Introduction

Buoyant convection in the presence of rotation represents a ubiquitous scenario for geophysical and astrophysical fluid flows that is largely responsible for the turbulent dynamics observed in planetary and stellar interiors [1, 2, 3], and in planetary atmospheres [4, 5] and oceans [6, 7]. The dynamics are highly complex with many influential ingredients such as geometry, compressibility, multiple components, and the presence of magnetic fields. In the absence of such complexities, the quintessential paradigm for investigating rotationally influenced buoyant flows is provided by rotating Rayleigh-Bénard convection (RRBC). A large number of studies has been published on this model system, which is very well suited for detailed experimental, numerical and theoretical studies, including [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], to name but a few. In its most distilled form, the problem consists of a rotating plane layer of fluid confined between two parallel horizontal plates which maintain a destabilizing temperature gradient. However, the interpretation of a layer within this paradigm may be broadened to include confined fluid domains such as cylinders, annuli, as well as spherical interiors and shells, which often arise in geophysical and astrophysical applications.

Five nondimensional parameters of geophysical and astrophysical interest highlight the relative importance of the Coriolis, pressure gradient, buoyancy, and dissipation forces in setting the acceleration of the fluid. These are the bulk Rossby, Euler, buoyancy, Reynolds and Péclet numbers, respectively:

$$ Ro_H = \frac{U}{2\Omega H}, \quad Eu = \frac{P}{\rho_0 U^2}, \quad \Gamma_H = \frac{g\alpha\|\nabla T_b\|H^2}{U^2}, \quad Re_H = \frac{UH}{\nu}, \quad Pe_H = \frac{UH}{\kappa} = PrRe_H, \quad (1) $$

which are comprised of intrinsic, extrinsic and characteristic properties of the fluid. Intrinsic material properties include the coefficient of thermal expansion $\alpha$, the kinematic viscosity $\nu$, and the thermal diffusivity $\kappa$, with $Pr = \nu/\kappa$ denoting the Prandtl number. Extrinsic properties include the magnitude $\Omega$ of the rotation rate, the layer depth $H$, the gravitational acceleration $g$ and the applied temperature gradient $\|\nabla T_b\|$. Characteristic properties include the velocity $U$, pressure $P$,
<table><caption>Table 1: Nondimensional parameter estimates for planetary [21] and satellite interiors [7]. Estimates of the Rossby number are derived from the relation $Ro_H = Re_H Ek$.</caption><thead><tr><th>Celestial body</th><th>Ek</th><th>Pr</th><th>Ro<sub>H</sub></th><th>Re<sub>H</sub></th></tr></thead><tbody><tr><td>Earth's outer core</td><td>10<sup>-15</sup></td><td>0.1</td><td>10<sup>-7</sup></td><td>10<sup>8</sup></td></tr><tr><td>Mercury (core)</td><td>10<sup>-12</sup></td><td>0.1</td><td>10<sup>-4</sup></td><td>10<sup>8</sup></td></tr><tr><td>Jupiter (core)</td><td>10<sup>-19</sup></td><td>0.1</td><td>10<sup>-10</sup></td><td>10<sup>9</sup></td></tr><tr><td>Europa (ocean)</td><td>10<sup>-12</sup></td><td>11.0</td><td>10<sup>-2.5</sup>–10<sup>-1.5</sup></td><td>10<sup>9.5</sup>–10<sup>10.5</sup></td></tr><tr><td>Ganymede (ocean)</td><td>10<sup>-10</sup>–10<sup>-13</sup></td><td>10.0</td><td>10<sup>-3.5</sup>–10<sup>1.5</sup></td><td>10<sup>9.5</sup>–10<sup>11.5</sup></td></tr><tr><td>Saturn (core)</td><td>10<sup>-18</sup></td><td>0.1</td><td>10<sup>-9</sup></td><td>10<sup>9</sup></td></tr><tr><td>Enceladus (ocean)</td><td>10<sup>-10</sup>–10<sup>-11</sup></td><td>13.0</td><td>10<sup>-3.5</sup>–10<sup>-1</sup></td><td>10<sup>7.5</sup>–10<sup>9</sup></td></tr><tr><td>Titan (ocean)</td><td>10<sup>-11</sup>–10<sup>-12</sup></td><td>10.0</td><td>10<sup>-3</sup>–1</td><td>10<sup>9</sup>–10<sup>11</sup></td></tr><tr><td>Neptune (core)</td><td>10<sup>-16</sup></td><td>10.0</td><td>10<sup>-6</sup></td><td>10<sup>10</sup></td></tr><tr><td>Uranus (core)</td><td>10<sup>-16</sup></td><td>10.0</td><td>10<sup>-6</sup></td><td>10<sup>10</sup></td></tr></tbody></table>

and the constant reference density $\rho_0$. The subscript 'H' signifies association with the bulk layer depth. Also of importance is the ratio of the viscous and Coriolis forces that provides an a priori external parameter referred to as the Ekman number

$$Ek = \frac{Ro_H}{Re_H} = \frac{\nu}{2\Omega H^2}. \qquad (2)$$

Turbulent flows are characterized by $Re_H \gg 1$, and when strongly influenced by rotation, by the ordering $Ek \ll Ro_H \lesssim 1$ (cf. Eq. (2)). Table I provides estimates of these nondimensional parameters in important geophysical and astrophysical settings. It can be seen that all such flows are rapidly rotating ($Ek \ll 1$), highly turbulent ($Re_H \gg 1$), and in the majority of situations strongly influenced by rotation ($Ro_H \ll 1$). To first approximation, using the theory of isotropic and statistically stationary turbulence as a benchmark, an order of magnitude estimate of the range of scales between the integral and dissipative scales in terms of the number of degrees of freedom per spatial direction and time is given, respectively, by $Re_H^{3/4} \ge \mathcal{O}(10^6)$ and $Re_H^{1/2} \ge \mathcal{O}(10^4)$ [22, 23]. Probing this region of parameter space is further complicated by an extended temporal range. Specifically, the smallness of the Rossby number $10^{-10} \le Ro_H \le 10^{-1}$ indicates an extreme timescale separation between fast inertial waves, associated with the Coriolis force, propagating on
$O(\Omega^{-1})$ time scales and the motion of eddies evolving on the advective time scale $O(H/U)$. From the combined values of $Re_H$ and $Ro_H$, Table 1 indicates that this extended temporal range may span as much as ten logarithmic decades.

<figure><img src="image_2.png" /><figcaption>Figure 1: Overview of the parameter space of RRBC spanned by the Ekman number $Ek$ and the bulk Reynolds number $Re_H$. Experiments, simulations, and dynamo models populate the parameter space characterized by moderately large Ekman numbers and moderately low Reynolds numbers (shaded) but these are far from their geostrophically relevant values. Adapted from [24], based on Table 1.</figcaption></figure>

From the standpoint of direct numerical simulations (DNS) – which are required to resolve all scales of the motion – these estimates are truly daunting. The current capability of state-of-the-art 3D DNS is $O(10^4)$ degrees of freedom in each spatial direction for periodic boundary conditions [25, 26, 27] indicating an upper threshold of $Re_H = O(10^5)$. In the presence of walls, the current state of the art in 3D DNS is $O(10^3)$ degrees of freedom in each spatial direction [28, 29, 30], further restricting the accessible range of $Re_H$. Hence, directly accessing the geophysical and astrophysical parameter regime is out of reach for the foreseeable future, even with impending advances to exascale supercomputing. Figure 1 captures this void visually. Recent DNS studies place the threshold at $Ek \gtrsim 10^{-8}$, $Re_H \lesssim 10^4$, e.g. [31, 32, 33, 34, 35, 36, 37, 30]. Inclusion of spherical geometry
and the capability for dynamo action in the simulations further restricts the reported range to $Ek \gtrsim 10^{-7}$, $Re_H \lesssim 10^3$ [38, 39, 40, 41], although hyperviscous simulations have been used to emulate lower Ekman numbers [42]. Given these limitations, a popular DNS strategy has been to vary $Re$ over the accessible range and attempt to uncover scaling laws in global quantities such as momentum and heat transport with a goal of extrapolating the results to the geophysical and astrophysical settings of Table 1. However, to be physically meaningful, such an extrapolation must be performed while respecting the strong rotational constraint $Ro_H = Re_H Ek \ll 1$. Inspection of this expression indicates that this gives rise to the challenging and somewhat incompatible requirement that the Ekman number be repeatedly lowered as $Re_H$ increases (see also Fig. 1). This leads to an amplification in the stiffness of the governing equations due to an increased separation between the time scales of inertial waves and advection, as well as between the advective and dissipation time scales (see Section 3). As a result, this requirement imposes severe time-stepping constraints on the majority of numerical algorithms currently in use.

The issue resides, in particular, with the precise implementation of the time-stepping scheme. Specifically, the linear Coriolis force $2\Omega U\hat{z} \times u$, of relative order $Ro_H^{-1}$ compared to inertial forces, is often treated explicitly, e.g. in [16, 33, 43, 44, 45, 46], while the advective timescale associated with the nonlinear advection term, $u \cdot \nabla u$, is invariably treated explicitly and so is known from a prior time step. Algorithmically, this avoids the complexities of implementing a coupled numerical solver for the momentum equations, and by contrast permits the use of a decoupled solver that updates fluid variables sequentially at each time step. However, several recent codes for simulating rapidly rotating convection, including [47, 48, 49], treat the Coriolis force implicitly. This formulation has a number of advantages and we use it below to identify a rescaled RRBC model and numerical algorithms capable of accessing regimes characterized by $Re_H \gg 1$ and $Ek \ll Ro_H \ll 1$.

The remainder of this paper is organized as follows. In Section 2, the incompressible Navier-Stokes equations (iNSE) in the Boussinesq approximation are given in a rotating frame. In Section 3, the detailed spatiotemporal resolution requirements for buoyantly driven, rotationally constrained flows are discussed and the need for implicit time-stepping treatments is highlighted. In Section 4, an asymptotically reduced set of equations, the nonhydrostatic quasi-geostrophic equations (NHQGE), is established as an instrumental guide for deducing a reformulation of the full iNSE. Informed by the asymptotic equations, Section 5 introduces a novel formulation of the iNSE termed the rescaled rapidly rotating incompressible Navier-Stokes equations (RRRiNSE). Section 6 highlights some of
the advantageous numerical properties of RRRiNSE, establishing that the numerical discretization is well conditioned. Section 7 contains a detailed comparison of fully nonlinear DNS using the newly introduced reformulation with established results from the literature, together with an analysis of the mean temperature equation, along with novel DNS results for the Nusselt and Reynolds numbers at unprecedented Ekman numbers (Ek as low as $10^{-15}$ and smaller). Finally, Section 8 concludes with a discussion of the implications of our results for future numerical simulations of rapidly rotating convection. Where necessary, relevant detailed calculations are relegated to Appendices.

## 2. The incompressible Navier-Stokes Equations: iNSE

In the classic paradigm of rotating Rayleigh-Bénard convection in a horizontal plane layer the fluid motion is accurately captured by the Boussinesq approximation that assumes small density fluctuations about a static background state, resulting in the incompressible Navier-Stokes equations (iNSE)

$$ (\partial_t + \mathbf{u} \cdot \nabla) \mathbf{u} = -2\Omega \hat{\mathbf{z}} \times \mathbf{u} - \nabla\pi + g\alpha\vartheta \hat{\mathbf{z}} + \nu \nabla^2 \mathbf{u}, \quad (3a) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (3b) $$

$$ (\partial_t + \mathbf{u} \cdot \nabla) \vartheta + \mathbf{u} \cdot \nabla T_b = \kappa \nabla^2 \vartheta, \quad (3c) $$

where $\mathbf{u}$ represents the convective fluid velocity, $\pi$ is the modified pressure, and $T = T_b(z) + \vartheta(x, t)$, i.e., the temperature is split into a static background profile $T_b(z)$ in the vertical direction and a convective temperature contribution $\vartheta$. The system rotates at a constant frequency $\Omega$ about the vertical direction $\hat{\mathbf{z}}$; the rotational Froude number is assumed to be sufficiently small that the centrifugal force can be neglected.

The equations of motion can be nondimensionalized by a characteristic but as yet undetermined flow velocity scale $U$, the layer depth $H$, and the characteristic temperature gradient $||\nabla T_b||$ giving

$$ (\partial_t + \mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{Ro_H} \hat{\mathbf{z}} \times \mathbf{u} - Eu\nabla\pi + \Gamma_H \vartheta \hat{\mathbf{z}} + \frac{1}{Re_H} \nabla^2 \mathbf{u}, \quad (4a) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (4b) $$

$$ (\partial_t + \mathbf{u} \cdot \nabla) \vartheta + \mathbf{u} \cdot \nabla T_b = \frac{1}{Pe_H} \nabla^2 \vartheta. \quad (4c) $$

The nondimensional parameters are defined in (1). In the next section, we describe the challenges
associated with solving the above set of nondimensional equations numerically in the regime $Ro_H \ll 1$, $Re_H \gg 1$.

### 3. Spatiotemporal resolution requirements for buoyantly driven flow

A review of the rotating convection literature illustrates why the $Ro_H \ll 1$ regime has proven to be so challenging for DNS [24, 37, 30, 50, 51]. Hereafter, for simplicity of exposition, we focus our discussion on the case where the rotation axis is antiparallel with gravity, i.e., the polar regime where $\nabla T_b \equiv \partial_z T_b \hat{z}$. For fixed temperature boundary conditions $\partial_z T_b = (T_b(H) - T_b(0))/H$, while for a fixed heat flux $F$, $\partial_z T_b = -F/\kappa$ instead. Within this regime a dynamical balance exists between the ageostrophic Coriolis, inertial and Archimedean (buoyancy) forces. This so-called CIA balance [52] establishes the rotational free-fall velocity

$$U_{rff} = \frac{g\alpha \|\partial_z T_b\| H}{2\Omega} \equiv Ro_c U_{ff} \quad (5)$$

as an appropriate estimate for the characteristic velocity, an estimate that has been verified both numerically [53, 54] and experimentally [55, 56, 57].

We define the convective Rossby number $Ro_c$ as the ratio of the rotational free-fall velocity to the buoyancy free-fall velocity $U_{ff}$ observed in rotationally unaffected regimes,

$$U_{ff} = \sqrt{g\alpha \|\partial_z T_b\| H^2}. \quad (6)$$

Thus

$$Ro_c \equiv \frac{U_{rff}}{U_{ff}} = \sqrt{\frac{Ra}{Pr}} Ek. \quad (7)$$

Here

$$Ra = \frac{g\alpha \|\partial_z T_b\| H^4}{\nu\kappa} \quad (8)$$

is the thermal Rayleigh number. This definition of $Ro_c$ is physically more precise than the equivalent definition $Ro_c = U_{ff}/(2\Omega H)$, i.e. as the Rossby number $Ro_H$ based on $U_{ff}$.

The convective Rossby number provides an external measure of the rotational constraint based on the imposed thermal Rayleigh number. Since $Ro_c \ll 1$ for rotationally constrained flows, it follows that $U_{rff} \ll U_{ff}$.

With $U = U_{rff}$ as the correct characteristic velocity scale it follows from (1) that

$$\Gamma_H = \frac{1}{Ro_c^2}, \quad Re_H = \frac{Ro_c^2}{Ek}. \quad (9)$$
Moreover, asymptotic linear theory [8] and simulations [58, 59] both indicate that rotating flows are highly anisotropic with

$$ \nabla_{\perp} \sim \frac{H}{\ell} \sim Ek^{-1/3}, \quad \partial_z \sim 1, \qquad (10) $$

and an $O(Ek^{-4/3})$ onset Rayleigh number, indicating that horizontal variations occur on the scale $\ell \ll H$ while vertical variations occur on the scale of the layer depth. Together with the corresponding nonlinear theory [60, 61, 62] these results lead to the introduction of the reduced Rayleigh number

$$ \widetilde{Ra} = RaEk^{4/3}, \qquad (11) $$

with $\widetilde{Ra} = O(1)$ defining the strongly forced ($Ra \gg 1$) but still rotationally constrained ($Ro_c \ll 1$) regime of interest. This regime extends from the convective threshold to highly supercritical Rayleigh numbers subject to the requirement that $\widetilde{Ra}$ is no larger than $\widetilde{Ra} = o(Ek^{-1/3})$. This upper bound represents the constraint required to maintain the local rotational constraint

$$ Ro_{\ell} \equiv \frac{U_{rff}}{2\Omega\ell} = Ro_c^2 \frac{H}{\ell} \sim Ro_c^2 Ek^{-1/3} = o(1). \qquad (12) $$

It is now clear that the external order parameter threaded throughout the rapidly rotating regime is $Ek^{1/3}$. Anticipating the derivation and discussion of the rescaled formulation, this observation suggests the definition of the small parameter:

$$ \epsilon = Ek^{1/3}. \qquad (13) $$

For given (nondimensional) grid resolutions $\Delta x$, $\Delta z \ll 1$ and temporal resolution $\Delta t \ll 1$, the assumption of an explicit time-stepping algorithm leads in Table 2 to the spatio-temporal constraints known as the Courant-Friedrich-Lewy (CFL) criteria required for accurately discretizing the various forces in the iNSE. The most to least restrictive CFL conditions are shown in dimensional (row 2) and nondimensional forms according to the vertical or horizontal diffusion timescale (rows 3 & 4). This ordering holds provided $\Delta x \sim \Delta z < \widetilde{Ra}^{-1}$. In the limit $Ek \to 0$, it is clear from Table 2 that the Coriolis term (column 2) imposes the most restrictive constraint on the time step $\Delta t$ (provided $\Delta x \sim \Delta z \gg \widetilde{Ra}E^{1/3}/Pr$). This suggests that it is numerically advantageous to treat this linear term implicitly with the additional expense of numerically coupling the momentum equations. It is also evident that, compared to nonlinear horizontal advection, an implicit treatment is desirable for the linear horizontal dissipation, the next most prohibitive constraint, $\propto (\Delta x_{\perp})^2$
<table><caption>Table 2: Ordering of the most to least restrictive CFL conditions for the explicit time step $\Delta t_f$ associated with the forcing term $f$ in the incompressible Navier-Stokes Equation (iNSE): $\partial_t u = f \equiv -2\Omega \hat{z} \times u - \nabla p + \nu \nabla_{\|}^2 u - u \cdot \nabla_{\|} u - g\alpha\theta \hat{r} - u \cdot \nabla_{\|} u + \nu \nabla_{\|}^2 u$. The time constraint for the pressure force is identical to that of the Coriolis force. Row 2 gives the dimensional time step estimate. Later rows express the nondimensional estimates based on vertical ($H^2/\nu$) and horizontal ($\ell^2/\nu$) diffusion times. Here $\ell \sim Ek^{1/3}H$ such that $\Delta x_{\perp}^* = \ell \Delta x_{\perp}$, $\Delta z^* = H \Delta z$ where $\Delta x_{\perp} \propto N_{x_{\perp}}^{-1}$ and $\Delta z \propto N_z^{-1}$. From [59, 52], $U_{rff} \sim Ro_c U_{ff}$ and the temperature fluctuation $\theta$, defined in Eq. (16), scales as $\theta \sim Ek^{1/3} \|\partial_z T_b\| H$, where $Ro_c = \sqrt{Ra/PrEk}$. Forces in need of an implicit treatment are presented to the left of the vertical separator $\|$. This holds provided $\Delta x_{\perp} \ll \tilde{R}a^{-1}$, otherwise, no advantages arise from the implicit treatment of horizontal dissipation (column 3) given that its CFL constraint becomes as restrictive as nonlinear horizontal advection (column 4).</caption><thead><tr><th>CFL time step</th><th>rotation</th><th>horiz. diff.</th><th>horiz. adv.</th><th>buoyancy</th><th>vert. adv.</th><th>vert. diff.</th></tr><tr><th>Δt<sub>f</sub></th><th>Δt<sub>Ω</sub></th><th>Δt<sub>ν⊥</sub></th><th>Δt<sub>adv⊥</sub></th><th>Δt<sub>g</sub></th><th>Δt<sub>adv‖</sub></th><th>Δt<sub>ν‖</sub></th></tr></thead><tbody><tr><td>Dimensional</td><td>1/2Ω</td><td>(Δx<sub>⊥</sub><sup>*</sup>)<sup>2</sup>/ν</td><td>Δx<sub>⊥</sub><sup>*</sup>/U<sub>rff</sub></td><td>U<sub>rff</sub>/gαθ</td><td>Δz<sup>*</sup>/U<sub>rff</sub></td><td>(Δz<sup>*</sup>)<sup>2</sup>/ν</td></tr><tr><td>Nondim., H<sup>2</sup>/ν</td><td>Ek</td><td>Ek<sup>2/3</sup>(Δx<sub>⊥</sub>)<sup>2</sup></td><td>Ek<sup>4/3</sup>/Ro<sub>c</sub><sup>2</sup>Δx<sub>⊥</sub></td><td>Ek<sup>2/3</sup></td><td>Ek/Ro<sub>c</sub><sup>2</sup>Δz</td><td>(Δz)<sup>2</sup></td></tr><tr><td>Nondim., l<sup>2</sup>/ν</td><td>Ek<sup>1/3</sup></td><td>(Δx<sub>⊥</sub>)<sup>2</sup></td><td>Ek<sup>2/3</sup>/Ro<sub>c</sub><sup>2</sup>Δx<sub>⊥</sub></td><td>1</td><td>Ek<sup>1/3</sup>/Ro<sub>c</sub><sup>2</sup>Δz</td><td>Ek<sup>-2/3</sup>(Δz)<sup>2</sup></td></tr></tbody></table>

provided $\Delta x_{\perp} = o(\tilde{R}a^{-1})$. If this strategy is adopted then all remaining time-stepping bounds for the linear terms are less severe than the $O(\Delta x_{\perp})$ nonlinear horizontal advection timescale. Thus all remaining linear terms can be treated explicitly without a numerical penalty.

The mechanical conditions at an impenetrable boundary also result in additional resolution constraints in space. Specifically, no-slip boundaries and/or stress-free boundaries that are not perpendicularly aligned to the axis of rotation result in $O(Ek^{1/2})$ Ekman boundary layers. For no-slip boundaries, it has recently been established that this prohibitive constraint can be relaxed by parameterizing its effect on the bulk through the pumping boundary conditions $w = \pm Ek^{1/2} \hat{z} \cdot \nabla \times u / \sqrt{2}$ [45, 63]. This complication does not arise in the present work, since we focus exclusively
on stress-free top and bottom boundaries.

Given the enormous challenges faced by DNS in the $Ro \ll 1$ regime, an attractive alternative is to resort to large-eddy simulations (LES), which resolve only the large turbulent scales, and employ subgrid-scale models for the smaller turbulent scales below a certain threshold scale. This technique has been applied in the context of nonrotating Rayleigh-Bénard convection [64, 65, 66]. However, it must be stressed that, for the highly anisotropic turbulent flows encountered in the geophysical and astrophysical context, LES are still in their infancy and ill-understood due to the complex structure across scales which such flows exhibit. Even when LES can be applied, the results thus obtained still need to be extrapolated to the extreme parameter regimes of geophysical and astrophysical flows. Importantly, LES and subgrid-scale modelling are particularly challenging because there is a notable paucity of validation data in the relevant regimes.

## 4. Asymptotically reduced model as a guide

Attempts to increase the achievable Reynolds number in DNS (or LES) of RRBC while lowering the Ekman number to sustain the low Rossby number environment must result from improving the conditioning of the matrices obtained from numerical discretization. Ultimately this means reducing or removing the discretization dependence on the Rossby or Ekman number. The asymptotic system of equations for RRBC valid in the limit $Ek \to 0$ derived and extensively studied by Julien & coworkers serves as a template for accomplishing this task [58, 59, 62, 67]. Assuming a local plane layer about the North pole, the system leverages $Ro_c \sim Ek^{1/3} \equiv \epsilon \ll 1$ as the small parameter, along with the characteristic anisotropic scalings (10), and the relations (9) and (11). A primary geostrophic balance is obtained together with horizontal incompressibility on horizontal spatial scales, namely,

$$ \hat{\mathbf{z}} \times \mathbf{u} \approx -\nabla\pi, \quad \nabla_{\perp} \cdot \mathbf{u}_{\perp} \approx 0. \qquad (14) $$

It follows that $\mathbf{u}_{\perp} = \hat{\mathbf{z}} \cdot \nabla_{\perp} \times \pi \equiv \nabla^{\perp}\pi$ where $\nabla^{\perp} = (-\partial_y, \partial_x)$. Moreover, it is found that the modified pressure $\pi = \Psi$ serves as the geostrophic streamfunction with $\mathbf{u} = (\nabla^{\perp}\Psi, w)$. When observed on the characteristic anisotropic spatial scales $\ell$ and $H$, and velocity scale $\nu/\ell$, the reduced
system of equations (the Non-Hydrostatic Quasi-Geostrophic Equations [58]) is given by

$$ \partial_t \nabla_{\perp}^2 \Psi + J [\Psi, \nabla_{\perp}^2 \Psi] - \partial_Z w = \nabla_{\perp}^4 \Psi, \qquad (15a) $$

$$ \partial_t w + J[\Psi, w] + \partial_Z \Psi = \frac{\tilde{R}a}{Pr} \theta + \nabla_{\perp}^2 w, \qquad (15b) $$

$$ \partial_t \theta + J[\Psi, \theta] + w(\partial_Z \bar{\Theta} - 1) = \frac{1}{Pr} \nabla_{\perp}^2 \theta, \qquad (15c) $$

$$ \epsilon^{-2} \partial_t \bar{\Theta} + \partial_Z (\overline{w\theta}) = \frac{1}{Pr} \partial_{ZZ} \bar{\Theta} \qquad (15d) $$

with

$$ w = \bar{\Theta} = 0, \quad \text{on} \quad Z = 0, 1. \qquad (15e) $$

Here $J[\Psi, f] = u_{\perp} \cdot \nabla_{\perp} f, \hat{z} \cdot \omega = \nabla_{\perp}^2 \Psi$ is the vertical vorticity and corrections at $O(Ek^{1/3})$ have been dropped; $Z$ is the $O(1)$ vertical scale. The temperature field is decomposed into a mean (horizontally-averaged) and $O(\epsilon)$ fluctuating component, i.e.,

$$ \vartheta \equiv \bar{\Theta} + \epsilon\theta. \qquad (16) $$

It follows that the Ekman number dependence remains only in the evolution of the mean temperature $\bar{\Theta}$ which can be seen to evolve on a much slower timescale $T = \epsilon^2 t$ (ratio of vertical viscous diffusion time to horizontal diffusion time) compared to the convective dynamics. Importantly, it has been established in [58, 59, 67] that this term can be omitted provided: (i) $\bar{\Theta}$ evolves to a statistically stationary state $\partial_t \bar{\Theta} \approx 0$, and (ii) the fluid domain is sufficiently large that numerous convective cells or plumes contribute to the horizontal spatial averaging.

The result serves as an accurate representation for $\bar{\Theta}$ with $O(Ek^{2/3})$ error (see Appendix D). Assuming an implicit treatment in time for all linear terms, it can now be seen that the most restrictive condition is the Ekman number-independent CFL constraint on the horizontal advection terms, i.e. $\Delta t = \Delta x/\|U_{RMS}\|$, which is consistent with our discussion in Table 2.

The NHQGE (15a)-(15e) have been instrumental in probing and identifying the properties of turbulent Rayleigh-Bénard convection in the rapidly rotating regime, from the identification of regimes of distinct flow morphologies [59], to the understanding of a novel inverse energy cascade in three dimensions [53, 68], through to uncovering the dissipation-free momentum and heat transport scaling laws [53, 54, 69]. However, by design, the NHQGE are constructed to be valid in the asymptotic regime $Ek \ll Ro_H \to 0$, and a complete understanding of its robustness to finite Ek remains
an open question [70]. Bridging the intermediate void in parameter space between the regimes obtained in current laboratory experiments and DNS and actual geophysical and astrophysical settings as highlighted in Fig. 1 is a key scientific objective.

## 5. Rescaled incompressible DNS

Based on the discussion in the previous section of the asymptotically reduced governing equations, we can now reformulate the iNSE in a rescaled form, which is advantageous for simulating very low Ekman and Rossby numbers. For this we will follow the template that produced the reduced system (15). We begin by introducing the anisotropic characteristic length scales: we nondimensionalize vertical lengths by $H$ and horizontal lengths by $l = \epsilon H$, where $\epsilon = Ek^{1/3}$ as before. We also adopt the velocity scale $U = U_ν = v/l$ which differs from a rotational free-fall velocity scale according to $U_ν = (Ra/Pr)U_{reff}$. This implies that

$$Ro_H \equiv \frac{U}{(2\Omega H)} = Ek^{2/3} \equiv \epsilon^2, \quad \Gamma_H = \frac{1}{\epsilon^2} \frac{\widetilde{Ra}}{Pr}, \quad Re_H = \frac{1}{\epsilon}, \quad Pe_H = PrRe_H = \frac{Pr}{\epsilon}, \quad (17)$$

$$Eu \equiv \frac{P}{\rho_0 U^2} = \frac{2\Omega U l^3}{\nu^2} = \frac{1}{Ek^{1/3}} = \frac{1}{\epsilon}. \quad (18)$$

Note that due to the anisotropic rescaling one finds

$$\nabla_{\perp} \mapsto \frac{1}{\epsilon} \nabla_{\perp}, \quad \hat{\mathbf{z}} \cdot \nabla = \partial_Z, \quad \partial_t \mapsto \frac{1}{\epsilon} \partial_t. \quad (19)$$

As in the derivation of the reduced equations, we decompose the temperature deviation from the linear conductive background state according to $\vartheta = \bar{\Theta}(Z, t) + \epsilon\theta(\mathbf{x}, t)$. Finally, we define the ageostrophic velocities

$$U = \frac{1}{\epsilon}(u + \partial_y \pi), \quad V = \frac{1}{\epsilon}(v - \partial_x \pi) \iff \mathbf{U}_{\perp} = \frac{1}{\epsilon}(\mathbf{u}_{\perp} - \nabla^{\perp} \pi). \quad (20a)$$
Recall $\nabla^{\perp} = (-\partial_y, \partial_x)$ and $\nabla_{\perp} = (\partial_x, \partial_y)$. The iNSE then take the form

$$ \partial_t u + \mathcal{N}_\epsilon u - V = \tilde{\nabla}_\epsilon^2 u. \qquad (20b) $$

$$ \partial_t v + \mathcal{N}_\epsilon v + U = \tilde{\nabla}_\epsilon^2 v, \qquad (20c) $$

$$ \partial_t w + \mathcal{N}_\epsilon w + \partial_z \pi = \tilde{\nabla}_\epsilon^2 w + \frac{\tilde{R}a}{Pr} \theta, \qquad (20d) $$

$$ \partial_x U + \partial_y V + \partial_z w = 0, \qquad (20e) $$

$$ \partial_t \theta + \mathcal{N}_\epsilon \theta + (\partial_z \bar{\Theta} - 1)w = \frac{1}{Pr} \tilde{\nabla}_\epsilon^2 \theta, \qquad (20f) $$

$$ \epsilon^{-2} \partial_t \bar{\Theta} + \partial_z \bar{w} \theta = \frac{1}{Pr} \partial_z^2 \bar{\Theta}, \qquad (20g) $$

where

$$ \tilde{\nabla}_\epsilon^2 = \nabla_\perp^2 + \epsilon^2 \partial_z^2, \qquad \nabla_\perp^2 = \partial_x^2 + \partial_y^2, \qquad \mathcal{N}_\epsilon = u\partial_x + v\partial_y + \epsilon w \partial_z. \qquad (21) $$

We consider impenetrable stress-free, fixed-temperature boundary conditions, i.e.,

$$ w = \partial_z u = \partial_z v = 0 \text{ at } Z = 0, 1, \qquad (22) $$

and

$$ \theta = \bar{\Theta} = 0 \text{ at } Z = 0, 1. \qquad (23) $$

The equation set (20a)-(23) is an equivalent reformulation of the iNSE obtained by rescaling terms (without omitting any terms in the process) in accordance with the asymptotic theory, specifically utilizing the distinguished limits (17)-(18) described in the previous section. We refer to these equations as the Rescaled Rapidly Rotating incompressible Navier-Stokes Equations (RRRiNSE).

We complete our exposition of these equations by noting that the formal limit $\epsilon \to 0$ of these equations leads directly to the asymptotic reduced equations (15a)-(15d) describing quasi-geostrophic Rayleigh-Bénard convection. This follows on noting that

$$ \lim_{\epsilon \to 0} \tilde{\nabla}_\epsilon^2 = \nabla_\perp^2; \qquad (24) $$

moreover, introducing the streamfunction $\Psi$ and velocity potential $\chi$ decomposition of the horizontal velocity field,

$$ \mathbf{u}_\perp = \nabla^\perp \Psi + \epsilon \nabla_\perp \chi, \qquad (25) $$

where $\nabla_\perp = (\partial_x, \partial_y, 0)^T$ and $\nabla^\perp = (-\partial_y, \partial_x, 0)^T$, we see that

$$ \lim_{\epsilon \to 0} \pi = \Psi; \quad \lim_{\epsilon \to 0} \mathbf{U}_\perp = \nabla_\perp \chi; \quad \lim_{\epsilon \to 0} \mathbf{u}_\perp = \nabla^\perp \Psi. \qquad (26) $$
Thus, as in the asymptotic equations, the dependence of the pressure on the velocity changes from quadratic to linear and a leading order geostrophic velocity field is recovered. Three-dimensional incompressibility is maintained through the ageostrophic velocity, $\nabla_{\perp} \cdot U_{\perp} = \nabla_{\perp}^2 \chi = -\partial_z w$. In Section 7 we demonstrate empirically that in sufficiently large domains $\epsilon^{-2}\partial_t \bar{\Theta} \approx 0$ in the statistically steady state, including at very small values of Ek.

## 6. Conditioning properties of rescaled equations: RRRiNSE

The advantage of the RRRiNSE formulation can be displayed through the properties of its spatio-temporal discretization. The findings of Table 2 suggest an implicit-explicit time discretization scheme for the governing equations of the form

$$ (\partial_t \mathcal{M} - \mathcal{L}_I) \mathcal{X}^{(n+1)} = \mathcal{L}_E \mathcal{X}^{(n)} + \mathcal{N}(\mathcal{X}^{(n)}, \mathcal{X}^{(n)}) \quad (27) $$

with implicit and explicit vectors of state $\mathcal{X}^{(n+1)} \equiv (\mathbf{u}^{(n+1)}, U_{\perp}^{(n+1)}, \pi^{(n+1)}, \theta^{(n+1)})^T$ and similarly for $\mathcal{X}^{(n)}$. The exact list of the variables that enter $\mathcal{X}$ and the expressions for differential operators $\mathcal{M}, \mathcal{L}_I, \mathcal{L}_E$ and $\mathcal{N}$ all depend on the adopted formulation. The specific details for various forms of RRRiNSE and the asymptotic model NHQGE are relegated to Appendix C.

Numerically, given that $\mathcal{L}_I$ is a non-diagonal operator, this requires the utilization of a coupled solver at each time step. The non-diagonal component is associated with the Coriolis force that would impose an $\mathcal{O}(\epsilon^{-1})$ explicit time-stepping constraint. The linear operator $\mathcal{L}_E$ represents vertical diffusion that, consistent with Table 2, imposes a nonrestrictive $\mathcal{O}(\epsilon^{-2}(\Delta z)^2)$ explicit time-stepping constraint. However, we note an apparent conflict in pursuing an explicit treatment of vertical diffusion, a term via which mechanical boundary conditions at the top and bottom are imposed, be it grid-based or via a basis of special functions. Two possibilities arise as resolutions to this predicament: (i) a reversion to an implicit treatment of $\mathcal{L}_E$, or (ii) the construction of Galerkin basis functions that automatically impose the mechanical boundary constraint. Implementation of the latter is insensitive to implicit/explicit treatments while the former requires a near-boundary resolution

$$ \epsilon^{-2} (\Delta z_b)^2 \lesssim \Delta t_{adv_{\perp}} \propto \frac{\Delta x}{\|U_{\perp}\|_{\infty}} \implies \Delta z_b \lesssim \epsilon (\Delta t_{adv_{\perp}})^{1/2} \quad (28) $$

to ensure that the numerical scheme is aware of the boundary constraint. This is more prohibitive than the $\mathcal{O}(\epsilon^{3/2})$ resolution constraint of an Ekman boundary layer in the limit $\epsilon \to 0$. We therefore
pursue a Galerkin function approach and, in addition, adopt an implicit treatment of $\mathcal{L}_E$, unless specified otherwise. Parameterized Ekman pumping boundary conditions [45], while not considered here, go beyond the Galerkin basis approach but can be handled via the use of the tau method [71] implemented in widely available solvers such as Dedalus [48].

Equations (27) are solved with the numerical code Coral [49], a flexible platform for solving systems of PDEs with spectral accuracy, i.e., with exponential error convergence. All fluid variables are discretized with a Fourier mode expansion in the horizontal and a Chebyshev-Galerkin polynomial expansion in the vertical direction, i.e.,

$$ \mathbf{v} = \sum \hat{\mathbf{v}}_{\mathbf{k}_{\perp},j}(t) e^{i\mathbf{k}_{\perp} \cdot \mathbf{x}_{\perp}} \Phi_j(Z). \quad (29) $$

Here the Chebyshev-Galerkin basis functions $\Phi_j(Z)$ can be constructed to satisfy Dirichlet ($\Phi_j = 0$), Neumann ($\partial_Z\Phi_j = 0$) or stress-free conditions ($\Phi_j = \partial_{ZZ}\Phi_j = 0$) or mixtures thereof on the boundaries; here, we consider the latter case only. The code temporally evolves the spectral coefficients of the modes $\hat{\mathbf{v}}_{\mathbf{k}_{\perp},j}(t)$ in spectral space, here via the third order-four stage implicit-explicit Runge-Kutta time-stepping scheme RK443 [72]. For constant-coefficient differential equations, as considered here, Coral adopts the quasi-inverse method presented in [73], based on an integral formulation of the problem, applied to Chebyshev-Galerkin bases obtained by basis recombination [73, 48, 49]. This procedure, detailed in Appendix B, is implemented in Coral and results in a sparse banded structure for the coupling matrices $\mathcal{M}$ and $\mathcal{L}_I$ in (27).

We note that, given the dependence of the time step $\Delta t$ on $Ek$ through its presence in $\mathcal{L}_I$, the ability to take the limit $Ek \to 0$ is ultimately bounded by the accuracy of the time integration due to the specific time-discretization error associated with the scheme and round-off errors. Such errors are ultimately related to the condition number of the matrix $\mathbf{A}$ in the linear algebraic system $\mathbf{A}\mathbf{v} = \mathbf{b}$ that results from the spatial-temporal discretization of (27). Here $\mathbf{A}$ is the discretization of $\partial_t\mathcal{M} - \mathcal{L}_I$ and $\mathbf{b}$ is the explicit right-hand side of (27). Alternatively, the sensitivity of an implicit time-stepping scheme can be explored through the eigenspectrum of the generalized eigenproblem $\lambda\mathcal{M} - \mathcal{L}_I$ deduced from the discretization of $(\partial_t\mathcal{M} - \mathcal{L}_I)\mathbf{v} = 0$.

Figure 2 shows this eigenspectrum, obtained with the RRRiNSE formulation, in the complex plane for four different Ekman numbers $Ek = 10^{-6}, 10^{-9}, 10^{-12}, 10^{-15}$ at $\tilde{Ra} = 0$ (top row) and $\tilde{Ra} = 5$ (bottom row), with the numerical results indicated by blue crosses (for the mixed velocity-vorticity formulation (A.1) with $N_Z = 256$ Chebyshev modes in the vertical), magenta circles and
<figure><img src="image_3.png" /><figcaption>Figure 2: Linear spectra from the rescaled equations (20) and (A.1) in the complex plane obtained numerically using the quasi-inverse method with Chebyshev-Galerkin basis (see appendix B). The top row illustrates the case of pure inertial waves ($\tilde{Ra} = 0$) while thermal stratification below the convective onset ($\tilde{Ra} = 5$) is included in the bottom row. All panels use $\tilde{k}_{\perp} = 1.3$ and $Pr = 1$. Numerical solutions are obtained for both the primitive variable formulation (20) with $N_Z = 256$ (magenta circles) and $N_Z = 512$ (green squares), and the mixed velocity-vorticity formulation (A.1) with $N_Z = 256$ (blue crosses). For reference, the analytical dispersion relation is represented with black dots appearing as a continuous black line. In both cases the computation of the numerical spectra remains stable as $Ek$ reaches values as low as $10^{-15}$. The mixed velocity-vorticity formulation leads to a remarkably accurate numerical spectrum, at the cost of larger memory usage.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Linear spectra from the unscaled equations (4) in the complex plane obtained numerically using the quasi-inverse method with Chebyshev-Galerkin basis (see appendix B). The top row illustrates the case of pure inertial waves ($\tilde{Ra} = 0$) while thermal stratification below the convective onset ($\tilde{Ra} = 5$) is included in the bottom row. All panels use $\tilde{k}_{\perp} = 1.3$ and $Pr = 1$. Numerical solutions obtained with $N_Z = 256$ (magenta circles) and $N_Z = 512$ (green squares) are compared against the analytical dispersion relation (black dots appearing as a continuous black line). In both cases spurious unstable modes proliferate with decreasing $Ek$.</figcaption></figure>
green squares (for the primitive-variable form (20) with $N_Z = 256$ and $N_Z = 512$ Chebyshev modes, respectively), in excellent agreement with the analytical result shown by the black line. This result should be compared with Fig. 3, which displays the eigenspectrum obtained with the unscaled Boussinesq equations at Ekman numbers $Ek = 10^{-7}, 10^{-9}, 10^{-10}, 10^{-11}$, again at $\tilde{Ra} = 0$ (top row) and $\tilde{Ra} = 5$ (bottom row), with magenta circles ($N_Z = 256$) and green squares ($N_Z = 512$) indicating the numerical data for comparison with the analytical result shown by the black line. In contrast with Fig. 2, Fig. 3 shows that the accuracy of the numerical spectra deteriorates significantly as the Ekman number decreases below $Ek \lesssim 10^{-7}$. Particularly damaging to time-stepping the solution are the spurious, linearly unstable modes (i.e., modes with $(\Re(s) > 0)$) visible in Fig. 3 for $\tilde{Ra} = 0$ when $Ek = 10^{-11}$ (top row) and even more spectacularly for the thermally forced case $\tilde{Ra} = 5$ when $Ek = 10^{-9}$ (bottom row), i.e., at a substantially slower rotation rate than in the purely hydrodynamical case. These results indicate why traditional DNS has proved unable to reach $Ek = 10^{-9}$.

The behavior of the spectra associated with the standard and the rescaled formulations is quantified and summarized in Fig. 4 which presents the condition number of the operator $L_I$ as rotation increases. High values for the condition number of matrices are commonly associated with unstable numerical computations. With the standard formulation (4), the conditioning of the discretized operator degrades rapidly as the Ekman number decreases following an approximate $Ek^{-3/2}$ scaling law at $\tilde{Ra} = 0$, and an even steeper scaling close to $Ek^{-2}$ at $\tilde{Ra} = 5$, as $Ek \to 0$. For comparison, the condition number associated to the RRRiNSE obeys a more moderate $Ek^{-1/2}$ scaling for the primitive variable form (20), and only a somewhat steeper scaling for the mixed velocity-vorticity form (A.1). The relative values of the condition number also speak clearly in favor of the rescaled formulation: the condition number computed with the RRRiNSE (in either the primitive variable or mixed velocity-vorticity forms) for geophysically relevant rotation strengths ($Ek = 10^{-16}$) appears to be smaller than its counterpart computed with the standard formulation even at modest rotation rate ($Ek \approx 10^{-6}$).

These results reflect the well-conditioned nature and the practicability of the RRRiNSE equations in the small $Ek$ limit, which the unscaled equations do not possess. This fact provides a strong motivation for using the RRRiNSE system to perform accurate DNS in the limit of small $Ek$.
<figure><img src="image_5.png" /><figcaption>Figure 4: Condition number of the operator $L_I$ computed for both the standard formulation [equation (4), open symbols] and the rescaled formulations [equations (20) and (A.1), filled symbols]. Both panels use $\tilde{k}_{\perp} = 1.3$ and $Pr = 1$. (a) Pure inertial waves for $\tilde{Ra} = 0$. (b) Thermally stratified case $\tilde{Ra} = 5$. The dash-dotted, dashed, and dotted lines are eye-guides with slope 1/2, 3/2 and 2, respectively.</figcaption></figure>

## 7. Results

To further validate the RRRiNSE formulation, going beyond the improved conditioning properties of RRRiNSE presented in Section 6, we perform extensive direct numerical simulations of RRBC using Coral. All runs described below were performed with $Pr = 1$ and a rescaled (nondimensional) domain with dimensions $10\ell_c \times 10\ell_c \times 1$ was used throughout, with critical convective onset length scale $\ell_c \approx 4.82$, unless specified otherwise. An explicit treatment is used for all advection RRRiNSE terms $N_i$, $i = u, v, w, \theta, \Theta$ and also for the mean temperature advection term $(\partial_z \bar{\Theta} - 1)w$ in the $\theta$ equation. The CFL condition is imposed based on the horizontal velocity components. Stress-free boundary conditions are adopted in all runs.

First, in Sections 7.1 and 7.2 we discuss different numerical schemes that can be used for solving the RRRiNSE. These differ in the treatment of the mean temperature equation and vertical derivatives in the diffusion terms. Our results focus on the global heat and momentum transport as defined by the nondimensional Nusselt and Reynolds numbers, namely,

$$Nu = 1 + Pr \langle w\theta \rangle_{Z,t} \equiv 1 + \partial_z \bar{\Theta}|_{0,1}, \quad Re_w = \langle \langle w^2 \rangle_Z^{1/2} \rangle_t, \quad (30)$$

where $\langle \cdot \rangle_{Z,t}$ denotes the combined average in the vertical and in time, while $\langle \cdot \rangle_Z$ and $\langle \cdot \rangle_t$ denote averaging in the vertical or in time separately. Along with the Nusselt number $Nu$, the quantity $Re_w$
saturates significantly earlier than the horizontal velocity components that are strongly impacted
by an inverse kinetic energy cascade [53, 54]. In Section 7.3, we compare the Nusselt number
obtained in our simulation to published results in the literature. Next, in Section 7.4, we present
the Nusselt and Reynolds numbers for different values of Ek and $\tilde{Ra}$, verifying the convergence of
the RRRiNSE to the asymptotically reduced equations presented in Section 4. Finally, we provide
visualizations of our RRRiNSE simulation results at a very low Ekman number, $Ek = 10^{-15}$, well
within the geostrophic turbulence regime.

7.1. Slaving: the role of $\epsilon^{-2}\partial_t\bar{\Theta}$ in the mean temperature equation

In the statistically steady state, quantities of interest, such as the Nusselt and Reynolds numbers
are typically given as space-time averages. It is natural to expect that this averaging improves as the
domain size is increased. This motivates the hypothesis, discussed in Appendix D, that the term
$\epsilon^{-2}\partial_t\bar{\Theta}$ in the equation may become subdominant in calculating, for instance, the Nusselt number,
provided the domain size is sufficiently large. The strategy of omitting the temporal variation of the
mean temperature, which we will refer to as the *slaving strategy*, has the significant advantage of
accelerating the convergence to the steady state at small $\epsilon$ by orders of magnitude, due to the fact
that $\epsilon^{-2}\partial_t = \partial_T$, where $T = \epsilon^2 t$, is a derivative with respect to a slow time variable. The slaving
approach has already been used successfully for the reduced equations of Section 4 in a number of
works, including [58, 67, 69, 68]. Here, we begin with a detailed verification of the slaving approach
for the full Boussinesq system in the RRRiNSE formulation.

Figures 5(a)-(c) show the time series of the Nusselt number for long simulations in domains
of three different sizes, $5\ell_c \times 5\ell_c \times 1$, $10\ell_c \times 10\ell_c \times 1$, and $15\ell_c \times 15\ell_c \times 1$ at $\tilde{Ra} = 40$, $Pr = 1$
and $Ek = 10^{-9}$. With this set of parameters, an inverse cascade of energy is observed, which
leads to the accumulation of energy at large scales and the formation of a large-scale vortex dipole
(LSV), also observed in the nonhydrostatic quasi-geostrophic equations (15), cf. [69, 53, 68]. Each
of the three simulation sets consists of two parts: first, each set is initialized with small-amplitude
initial conditions and integrated for a long time with the slaving approach, until the LSV has
saturated. Then, at the time indicated in Figs. 5(a)-(c) by an arrow, the time derivative of the mean
temperature is again included, restoring the full RRRiNSE equations, and the run is continued. In
small domains there is a notable discrepancy between the slaving strategy and the solution of the
full RRRiNSE equations, but this discrepancy decreases as the domain size increases (see dashed
<figure><img src="image_6.png" /><figcaption>Figure 5: Panels (a)-(c): time series of $Nu - 1$ at $\tilde{Ra} = 40$, $Pr = 1$, $Ek = 10^{-9}$ for varying horizontal domain sizes. (a) $L = 5\ell_c$ (resolution $64^2 \times 128$ modes), (b) $L = 10\ell_c$ (resolution $128^3$ modes), (c) $L = 15\ell_c$ (resolution $192^3$ modes), with the rescaled most unstable length scale $\ell_c \approx 4.82$. Each simulation consists of two parts. In the first part, which extends from $t = 0$ to the time $t = t_c$ indicated by a black arrow, the time derivative of $\bar{\Theta}$ in Eq. (A.1i) is omitted, leading to a slaving relation between $\bar{\Theta}$ and the heat flux. At $t = t_c$ (at which the large scale vortex condensate reaches saturation), the time derivative of $\bar{\Theta}$ is restored. Horizontal dashed lines indicate time averages over each of the two segments of the simulation, while the blue and orange shaded areas indicate the observed standard deviation. The agreement between the results of the two numerical schemes improves noticeably as the domain size increases. Fluctuations are seen to be larger in smaller domains since the volume average contains fewer points. Panels (d)-(f): histograms of $Nu - 1$ computed from each time series with and without the mean temperature time derivative, with dashed lines indicating a Gaussian fit. The histograms illustrate the convergence of the mean between the two schemes with increasing domain size. The histograms also show that fluctuations in the presence of the $\epsilon^{-2}\partial_t\bar{\Theta}$ term are of larger amplitude, a fact consistent with the time series in the top row.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Time series of $Nu - 1$ for various horizontal domain sizes at $Ek = 10^{-9}$. (a) $L = 5\ell_c$ (resolution $64^2 \times 128$). (b) $L = 10\ell_c$ (resolution $128^3$). (c) $L = 15\ell_c$ (resolution $192^3$). First part of each simulation is performed without the $\partial_t \bar{\Theta}$ term in the mean temperature equation. In the second part of each run, the term $\epsilon^{-2}\partial_t \bar{\Theta}$ is replaced by $\partial_t \bar{\Theta}$. The mean and variance of each of the two segments are found to be very close.</figcaption></figure>

lines). This is accompanied by a decrease in the statistical fluctuations about the mean Nusselt number, as expected given the improved horizontal averaging in larger domains. The histograms in Figs. 5(d)-(f) correspond to each of the two parts of the time series above and illustrate both of these trends: the averages of the two histograms approach each other as the domain size is increased, and the variance decreases, being somewhat larger for the full equations than with the slaving strategy. Thus, for large enough domains (horizontal domain size $L \gtrsim 10\ell_c$ in terms of the critical length scale $\ell_c$), the slaving scheme yields approximately the same answers for mean quantities as the unaltered equations. On the other hand, differences remain in the fluctuations about that mean, owing to the additional slow time scale arising from the time derivative of the mean temperature, eliminated in the slaving strategy. We also note that in all cases, the peak of the histogram is close to a Gaussian, while in the presence of $\epsilon^{-2}\partial_t \bar{\Theta}$ deviations from that shape are seen, most strikingly in panel (e) in the intermediate domain, leading to a certain degree of skewness. This is not observed to the same degree in the smaller or larger domain, and remains to be better understood in future investigations.

Figure 6 shows long time series similar to those in Figs. 5(a)-(c), again at $Ek = 10^{-9}$. The first, earlier part of each simulation (orange curves in Fig. 6) is identical between the two figures and was computed using the slaving approach. In contrast, the second part (in blue) was computed with $\epsilon^{-2}\partial_t \bar{\Theta}$ replaced by $\partial_t \bar{\Theta}$ (without the $\epsilon^{-2}$ prefactor). The Nusselt number evolution in the two segments is found to be close to indistinguishable. This indicates that, in the statistically stationary state, $\partial_t \bar{\Theta}$ is small compared to the remaining terms in the mean temperature equation.
In summary, the slaving strategy is a highly attractive scheme for accelerating transient dynamics
in the approach to a statistically steady state, and it is often preferable to adopt this strategy for a
sizeable efficiency gain. We therefore adopt the slaving strategy in all the runs described below. In
Appendix E, we examine the impact of rescaling and slaving on the stability of strongly nonlinear
simulations, showing that the slaving strategy yields a significant improvement in stability.

7.2. Implicit and explicit vertical derivatives in the diffusion terms

A particular feature of the RRRiNSE equations is the appearance of vertical derivatives with
a prefactor of $\epsilon = Ek^{1/3}$. This indicates that, for sufficiently small $\epsilon$, terms involving vertical
derivatives become subdominant and their numerical treatment becomes irrelevant in the Courant-
Friedrichs-Lewy (CFL) constraint on the integration of the equations of motion.

To test whether this intuition is correct, we perform two sets of runs at $\tilde{Ra} = 60, Pr = 1$, varying $Ek$ between $10^{-1}$ and $10^{-12}$, as listed in Table 3. In the first set of runs, all diffusive terms are treated implicitly (as in all other runs described in later sections), while in the second set the diffusive terms in Eqs. (A.1e)-(A.1i) involving a vertical derivative are treated explicitly (the CFL condition is still only applied based on the horizontal velocity field). This is expected to produce no significant difference in the simulation outcome, provided $Ek$ is sufficiently small.

The data provided in Table 3 and its visualization in Fig. 7 show that the runs with explicit and
implicit vertical diffusion schemes produce Nusselt and Reynolds numbers which are compatible
with each other within the margin of error (computed as the standard deviation in steady state),
provided that $Ek = 1/\sqrt{Ta} \lesssim 10^{-6}$. For $Ek \gtrsim 10^{-6}$, i.e. $\epsilon \gtrsim 0.01$, the simulation becomes unstable
with a CFL prefactor of 0.2 and a time step tolerance of 0.3, leading to an unphysical blow-up.

In short, diffusion terms involving vertical derivatives become irrelevant in the limit of small $Ek$. This is consistent with the discussion in Section 4, which highlighted that, as $Ek \to 0$, the RRRiNSE formulation directly converges to the asymptotically reduced equations, which do not contain vertical diffusion terms except in the mean temperature equation.

7.3. *Comparison with published Nusselt numbers*

To further ascertain the validity of the RRRiNSE formulation, we reproduce results from the
literature. In [34], Kunnen and co-workers provide Nusselt numbers obtained from direct numerical
simulations of the same set-up as ours, using a second-order energy-conserving, finite-difference
<table><caption>Table 3: Overview of runs with $\tilde{Ra} = 60$, $Pr = 1$ and implicit or explicit vertical diffusion schemes (in a rescaled domain of size $10\ell_c \times 10\ell_c \times 1$, where $\ell_c = 4.82$) for Ekman numbers between $10^{-12}$ and $10^{-1}$. Values of $Nu-1$ and $Re_w$, defined in Eq. (30), refer to the average in the quasi-steady state during the early phase of the nonlinear evolution where a large-scale vortex condensate slowly grows in amplitude if an inverse cascade is present. Uncertainties represent the standard deviation of the time series.</caption><thead><tr><th>∂<sub>Z</sub> (diffusion)</th><th>Ek</th><th>Ra</th><th>Ra</th><th>N<sub>x</sub> × N<sub>y</sub> × N<sub>Z</sub></th><th>Stability</th><th>(Nu − 1) ± ΔNu</th><th>Re<sub>w</sub> ± ΔRe<sub>w</sub></th></tr></thead><tbody><tr><td>Implicit</td><td>1.0 × 10<sup>−1</sup></td><td>60</td><td>1.29 × 10<sup>3</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>0.6 ± 0.1</td><td>1.8 ± 0.1</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−2</sup></td><td>60</td><td>2.78 × 10<sup>4</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>3.2 ± 0.2</td><td>5.5 ± 0.1</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−3</sup></td><td>60</td><td>6.0 × 10<sup>5</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>8.3 ± 0.3</td><td>9.1 ± 0.2</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−4</sup></td><td>60</td><td>1.29 × 10<sup>7</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>20.7 ± 0.6</td><td>17.7 ± 0.3</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−5</sup></td><td>60</td><td>2.78 × 10<sup>8</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>26.1 ± 1.0</td><td>18.3 ± 0.4</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−6</sup></td><td>60</td><td>6.00 × 10<sup>9</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>25.4 ± 1.4</td><td>18.8 ± 0.7</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−7</sup></td><td>60</td><td>1.29 × 10<sup>11</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>21.9 ± 1.4</td><td>17.4 ± 0.5</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−8</sup></td><td>60</td><td>2.78 × 10<sup>12</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.6 ± 1.0</td><td>16.6 ± 0.5</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−9</sup></td><td>60</td><td>6.00 × 10<sup>13</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.6 ± 0.9</td><td>16.6 ± 0.5</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−10</sup></td><td>60</td><td>1.29 × 10<sup>15</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.4 ± 1.0</td><td>16.7 ± 0.5</td></tr><tr><td>Implicit</td><td>1.0 × 10<sup>−12</sup></td><td>60</td><td>6.00 × 10<sup>17</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.5 ± 0.9</td><td>16.5 ± 0.5</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−1</sup></td><td>60</td><td>1.29 × 10<sup>3</sup></td><td>128 × 128 × 128</td><td>Unstable</td><td>-</td><td>-</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−2</sup></td><td>60</td><td>2.78 × 10<sup>4</sup></td><td>128 × 128 × 128</td><td>Unstable</td><td>-</td><td>-</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−3</sup></td><td>60</td><td>6.00 × 10<sup>5</sup></td><td>128 × 128 × 128</td><td>Unstable</td><td>-</td><td>-</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−4</sup></td><td>60</td><td>1.29 × 10<sup>7</sup></td><td>128 × 128 × 128</td><td>Unstable</td><td>-</td><td>-</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−5</sup></td><td>60</td><td>2.78 × 10<sup>8</sup></td><td>128 × 128 × 128</td><td>Unstable</td><td>-</td><td>-</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−6</sup></td><td>60</td><td>6.00 × 10<sup>9</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>25.4 ± 1.3</td><td>18.6 ± 0.6</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−7</sup></td><td>60</td><td>1.29 × 10<sup>11</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>21.7 ± 1.2</td><td>17.5 ± 0.6</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−8</sup></td><td>60</td><td>2.78 × 10<sup>12</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.7 ± 1.0</td><td>16.9 ± 0.5</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−9</sup></td><td>60</td><td>6.00 × 10<sup>13</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.6 ± 1.0</td><td>16.8 ± 0.5</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−10</sup></td><td>60</td><td>1.29 × 10<sup>15</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.4 ± 0.9</td><td>16.8 ± 0.5</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−11</sup></td><td>60</td><td>2.78 × 10<sup>16</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.2 ± 0.9</td><td>16.6 ± 0.4</td></tr><tr><td>Explicit</td><td>1.0 × 10<sup>−12</sup></td><td>60</td><td>6.00 × 10<sup>17</sup></td><td>128 × 128 × 128</td><td>Stable</td><td>19.3 ± 1.0</td><td>16.7 ± 0.5</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Figure 7: Visualization of the data listed in Table 3, comparing the results of implicit and explicit vertical diffusion schemes at different Taylor numbers $Ta = Ek^{-2}$. At $Ta \ge 10^{12}$, implicit and explicit diffusion schemes are both stable and yield very similar results which are compatible within the margin of error (given by the standard deviation).</figcaption></figure>

code with fractional time-stepping. We stress that, while the authors of [34] discuss the transition to geostrophic turbulence, they are only able to reach relatively modest Ekman numbers, $Ek \gtrsim 1.34 \times 10^{-7}$, which is large relative to the values of $Ek$ that can be reached using RRRiNSE, as shown in Sections 7.2 and 7.4. Nonetheless, the results presented in [34] provide a valuable benchmark as one bookend at moderate $Ek$.

We perform runs at the same parameters as those given in [34]. In Table 4, set A, we list those runs along with the Nusselt and Reynolds numbers obtained from our simulations. We choose up to $384 \times 384$ dealiased Fourier modes in the horizontal directions and up to 256 dealiased Chebyshev modes in the vertical derivation, while Kunnen et al. [34] consider up to $512 \times 512 \times 1024$ spatial grid points. We stress that these numbers cannot be easily compared between a spectral code such as Coral and finite-difference codes such as that employed in [34]. However, the differences reside in the exponential vs algebraic error convergence properties of the two algorithmic approaches.

Despite the different codes and resolution requirements, Fig. 8 shows that the Nusselt numbers obtained using RRRiNSE and those of Kunnen et al. agree well within the margin of error (the standard deviation of the Nusselt number time series). This provides a first bookend at relatively large Ekman numbers, where the RRRiNSE formulation correctly reproduces known results.
<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison of Nusselt numbers from [34] with the RRRiNSE results. The simulations shown correspond to runs A1–A12 in Table 4. Error bars for RRRiNSE data represent the observed standard deviation.</figcaption></figure>

7.4. Convergence to the geostrophic branch

As discussed in Section 4, it is expected that rotating convective flows converge to the well-
studied, asymptotically reduced equations in the limit of small Ekman numbers. However, to date,
it has not been possible to achieve sufficiently small Ekman numbers in direct numerical simulations
of the full Boussinesq equations to observe this convergence.

Owing to the improved conditioning of the RRRiNSE formulation in the small *Ek* limit, it is shown below that it becomes possible, for the first time, to observe this convergence. We perform four sets of simulations of the RRRiNSE at *Ra* = 40, 60, 80, 120, *Pr* = 1 and *Ek* varying between 10<sup>-1</sup> and 10<sup>-15</sup> (corresponding to *Ta* between 10<sup>2</sup> and 10<sup>30</sup>), summarized in Tables 4, 5 and 6. We note that the RRRiNSE formulation remained numerically stable even at *Ek* = 10<sup>-24</sup>, and yielded approximately the same Nusselt and Reynolds numbers as the case *Ek* = 10<sup>-15</sup>, but in order to avoid potential issues due to machine precision, these results are not shown here. In addition, we perform simulations of the asymptotically reduced equations described in Section 4 with *Pr* = 1 and *Ra* = 40, 60, 80, 120 and compare the observed Nusselt and Reynolds numbers with the RRRiNSE results.

Figures 9 and 10 show, respectively, the Nusselt and Reynolds numbers (based on vertical velocity and domain height), obtained at each $\tilde{Ra}$ as a function of $Ta$ (symbols), with error bars
<figure><img src="image_10.png" /><figcaption>Figure 9: Nondimensional heat flux measured by the Nusselt number $Nu$, defined in Eq. (30), compensated by the turbulent scaling law $\tilde{Ra}^{3/2}$, cf. [59], versus the Taylor number $Ta = Ek^{-2}$ for fixed $\tilde{Ra} \equiv RaEk^{4/3}$ (blue $\tilde{Ra} = 40$, orange $\tilde{Ra} = 60$, red $\tilde{Ra} = 80$, green $\tilde{Ra} = 120$), computed with implicit vertical diffusion. Dashed lines show the average in steady state predicted by the reduced equations (15a)-(15d) with the shaded region indicating one standard deviation about the mean. At low Ekman numbers, the RRRiNSE predictions converge to the reduced equations.</figcaption></figure>

indicating the observed standard deviation of the time series. The Nusselt number is rescaled by $\tilde{Ra}^{3/2}$, the turbulent scaling law [52], leading to an approximate data collapse between different $\tilde{Ra}$ at large $Ta$ (small $Ek$), in agreement with the asymptotically reduced equations [69]. The Reynolds numbers are rescaled by $\tilde{Ra}$, leading to a less satisfactory collapse, which is known to be related to the presence of the inverse energy cascade [53, 54].

In addition, the dashed lines in Figs. 9 and 10 indicate the results obtained using the asymptotically reduced equations, with the shaded area showing the standard deviation. The RRRiNSE results are seen to converge to the values observed in the asymptotically reduced equations above a certain threshold in the Taylor number $Ta$, within the error margins given by the standard deviation of the time series from the reduced equations. The threshold $Ta$ required for this convergence appears to increase with $\tilde{Ra}$, but a more detailed investigation will be required in the future to quantitatively investigate this behavior. The observed convergence of the RRRiNSE to the reduced equations provides an additional bookend validating the accuracy of the RRRiNSE formulation against an established body of work in the limit $Ek \to 0$ ($Ta \to \infty$). We also verified that the
<table><caption>Table 4: List of simulations described in this work (see also Tables 5 and 6). All simulations are done with $Pr = 1$ in a rescaled domain of size $10\ell_c \times 10\ell_c \times 1$, where $\ell_c \approx 4.82$. Simulations A1 through A12 have $Ra$, $Ek$ and $Pr$ identical to those in [34]. The resolution is specified by the numbers $N_x$, $N_y$ of Fourier modes in the horizontal directions and the number $N_Z$ of Chebyshev modes in the vertical. The values of $Nu$ and $Re_w$, defined in Eq. (30), refer to the average computed in the early, quasi-steady, nonlinear stage of the evolution, in the absence of an LSV. Uncertainties represent the standard deviation of the time series.</caption><thead><tr><th>Simulation nr.</th><th>Ek</th><th>R̃a</th><th>Ra</th><th>N<sub>x</sub> × N<sub>y</sub> × N<sub>Z</sub></th><th>(Nu − 1) ± ΔNu</th><th>Re<sub>w</sub> ± ΔRe<sub>w</sub></th></tr></thead><tbody><tr><td>A1</td><td>4.0 × 10<sup>−7</sup></td><td>29.5</td><td>1 × 10<sup>10</sup></td><td>128 × 128 × 128</td><td>8.1 ± 0.5</td><td>7.6 ± 0.2</td></tr><tr><td>A2</td><td>6.0 × 10<sup>−7</sup></td><td>50.6</td><td>1 × 10<sup>10</sup></td><td>128 × 128 × 128</td><td>18.9 ± 1.1</td><td>15.3 ± 0.5</td></tr><tr><td>A3</td><td>9.0 × 10<sup>−7</sup></td><td>86.9</td><td>1 × 10<sup>10</sup></td><td>128 × 128 × 128</td><td>44.5 ± 2.6</td><td>29.2 ± 1.3</td></tr><tr><td>A4</td><td>1.2 × 10<sup>−6</sup></td><td>127.5</td><td>1 × 10<sup>10</sup></td><td>256 × 256 × 256</td><td>67.2 ± 3.1</td><td>38.5 ± 1.6</td></tr><tr><td>A5</td><td>1.5 × 10<sup>−6</sup></td><td>171.5</td><td>1 × 10<sup>10</sup></td><td>256 × 256 × 256</td><td>91.6 ± 4.7</td><td>53.1 ± 2.4</td></tr><tr><td>A6</td><td>2.0 × 10<sup>−6</sup></td><td>252.2</td><td>1 × 10<sup>10</sup></td><td>384 × 384 × 256</td><td>110.2 ± 5</td><td>59.6 ± 2.5</td></tr><tr><td>A7</td><td>1.34 × 10<sup>−7</sup></td><td>34.3</td><td>5 × 10<sup>10</sup></td><td>128 × 128 × 128</td><td>9.5 ± 0.6</td><td>9.1 ± 0.3</td></tr><tr><td>A8</td><td>1.79 × 10<sup>−7</sup></td><td>50.4</td><td>5 × 10<sup>10</sup></td><td>128 × 128 × 128</td><td>18.2 ± 1.1</td><td>14.8 ± 0.4</td></tr><tr><td>A9</td><td>2.95 × 10<sup>−7</sup></td><td>98.3</td><td>5 × 10<sup>10</sup></td><td>128 × 128 × 128</td><td>52.6 ± 3.1</td><td>33.0 ± 1.6</td></tr><tr><td>A10</td><td>4.02 × 10<sup>−7</sup></td><td>148.3</td><td>5 × 10<sup>10</sup></td><td>256 × 256 × 256</td><td>88.7 ± 5.2</td><td>47.9 ± 1.5</td></tr><tr><td>A11</td><td>4.92 × 10<sup>−7</sup></td><td>194.6</td><td>5 × 10<sup>10</sup></td><td>256 × 256 × 256</td><td>116.5 ± 7.2</td><td>61.4 ± 3.2</td></tr><tr><td>A12</td><td>6.71 × 10<sup>−7</sup></td><td>293.6</td><td>5 × 10<sup>10</sup></td><td>384 × 384 × 256</td><td>148.4 ± 10.2</td><td>78.0 ± 4.2</td></tr><tr><td>B1</td><td>1.0 × 10<sup>−1</sup></td><td>40</td><td>8.62 × 10<sup>2</sup></td><td>128 × 128 × 128</td><td>0.029 ± 0.007</td><td>0.34 ± 0.02</td></tr><tr><td>B2</td><td>1.0 × 10<sup>−2</sup></td><td>40</td><td>1.86 × 10<sup>4</sup></td><td>128 × 128 × 128</td><td>0.287 ± 0.014</td><td>1.48 ± 0.03</td></tr><tr><td>B3</td><td>1.0 × 10<sup>−3</sup></td><td>40</td><td>4.0 × 10<sup>5</sup></td><td>128 × 128 × 128</td><td>6.26 ± 0.17</td><td>7.0 ± 0.1</td></tr><tr><td>B4</td><td>1.0 × 10<sup>−4</sup></td><td>40</td><td>8.62 × 10<sup>6</sup></td><td>128 × 128 × 128</td><td>12.0 ± 0.4</td><td>10.3 ± 0.2</td></tr><tr><td>B5</td><td>3.0 × 10<sup>−5</sup></td><td>40</td><td>4.29 × 10<sup>7</sup></td><td>128 × 128 × 128</td><td>14.7 ± 0.4</td><td>11.4 ± 0.3</td></tr><tr><td>B6</td><td>1.0 × 10<sup>−5</sup></td><td>40</td><td>1.86 × 10<sup>8</sup></td><td>128 × 128 × 128</td><td>15.8 ± 0.7</td><td>12.2 ± 0.3</td></tr><tr><td>B7</td><td>1.0 × 10<sup>−6</sup></td><td>40</td><td>4.0 × 10<sup>9</sup></td><td>128 × 128 × 128</td><td>14.8 ± 1.0</td><td>12.0 ± 0.4</td></tr><tr><td>B8</td><td>1.0 × 10<sup>−7</sup></td><td>40</td><td>8.62 × 10<sup>10</sup></td><td>128 × 128 × 128</td><td>12.2 ± 0.7</td><td>10.9 ± 0.3</td></tr><tr><td>B9</td><td>1.0 × 10<sup>−8</sup></td><td>40</td><td>1.86 × 10<sup>12</sup></td><td>128 × 128 × 128</td><td>11.7 ± 0.5</td><td>10.6 ± 0.2</td></tr><tr><td>B10</td><td>1.0 × 10<sup>−9</sup></td><td>40</td><td>4.0 × 10<sup>13</sup></td><td>128 × 128 × 128</td><td>11.5 ± 0.6</td><td>10.7 ± 0.3</td></tr><tr><td>B11</td><td>1.0 × 10<sup>−10</sup></td><td>40</td><td>8.62 × 10<sup>14</sup></td><td>128 × 128 × 128</td><td>11.65 ± 0.65</td><td>10.7 ± 0.3</td></tr><tr><td>B12</td><td>1.0 × 10<sup>−11</sup></td><td>40</td><td>1.86 × 10<sup>16</sup></td><td>128 × 128 × 128</td><td>11.4 ± 0.6</td><td>10.7 ± 0.3</td></tr><tr><td>B13</td><td>1.0 × 10<sup>−12</sup></td><td>40</td><td>4.0 × 10<sup>17</sup></td><td>128 × 128 × 128</td><td>11.4 ± 0.6</td><td>10.6 ± 0.3</td></tr><tr><td>B14</td><td>1.0 × 10<sup>−15</sup></td><td>40</td><td>4.0 × 10<sup>21</sup></td><td>128 × 128 × 128</td><td>11.4 ± 0.6</td><td>10.6 ± 0.3</td></tr><tr><td>B15</td><td>1.0 × 10<sup>−15</sup></td><td>40</td><td>4.0 × 10<sup>21</sup></td><td>256 × 256 × 256</td><td>11.4 ± 0.5</td><td>10.7 ± 0.3</td></tr></tbody></table>
<table><caption>Table 5: List of simulations described in this work (continued 1); see also Tables 4 and 6. Simulations C1-C3, C6-13 are identical to the runs with an implicit vertical diffusion scheme listed in Table 3.</caption><thead><tr><th>Simulation nr.</th><th>Ek</th><th>Ra</th><th>Ra</th><th>N<sub>x</sub> × N<sub>y</sub> × N<sub>Z</sub></th><th>(Nu − 1) ± ΔNu</th><th>Re<sub>w</sub> ± ΔRe<sub>w</sub></th></tr></thead><tbody><tr><td>C1</td><td>1.0 × 10<sup>−1</sup></td><td>60</td><td>1.29 × 10<sup>3</sup></td><td>128 × 128 × 128</td><td>0.6 ± 0.1</td><td>1.8 ± 0.1</td></tr><tr><td>C2</td><td>1.0 × 10<sup>−2</sup></td><td>60</td><td>2.78 × 10<sup>4</sup></td><td>128 × 128 × 128</td><td>3.2 ± 0.2</td><td>5.5 ± 0.1</td></tr><tr><td>C3</td><td>1.0 × 10<sup>−3</sup></td><td>60</td><td>6.0 × 10<sup>5</sup></td><td>128 × 128 × 128</td><td>8.3 ± 0.3</td><td>9.1 ± 0.2</td></tr><tr><td>C4</td><td>3.0 × 10<sup>−4</sup></td><td>60</td><td>2.99 × 10<sup>6</sup></td><td>128 × 128 × 128</td><td>12.5 ± 0.4</td><td>12.1 ± 0.2</td></tr><tr><td>C5</td><td>1.75 × 10<sup>−4</sup></td><td>60</td><td>6.13 × 10<sup>6</sup></td><td>128 × 128 × 128</td><td>14.7 ± 0.4</td><td>13.3 ± 0.3</td></tr><tr><td>C6</td><td>1.0 × 10<sup>−4</sup></td><td>60</td><td>1.29 × 10<sup>7</sup></td><td>128 × 128 × 128</td><td>20.7 ± 0.6</td><td>17.7 ± 0.3</td></tr><tr><td>C7</td><td>1.0 × 10<sup>−5</sup></td><td>60</td><td>2.78 × 10<sup>8</sup></td><td>128 × 128 × 128</td><td>26.1 ± 1.0</td><td>18.3 ± 0.4</td></tr><tr><td>C8</td><td>1.0 × 10<sup>−6</sup></td><td>60</td><td>6.0 × 10<sup>9</sup></td><td>128 × 128 × 128</td><td>25.4 ± 1.4</td><td>18.8 ± 0.7</td></tr><tr><td>C9</td><td>1.0 × 10<sup>−7</sup></td><td>60</td><td>1.29 × 10<sup>11</sup></td><td>128 × 128 × 128</td><td>21.9 ± 1.4</td><td>17.4 ± 0.5</td></tr><tr><td>C10</td><td>1.0 × 10<sup>−8</sup></td><td>60</td><td>2.78 × 10<sup>12</sup></td><td>128 × 128 × 128</td><td>19.6 ± 1.0</td><td>16.6 ± 0.6</td></tr><tr><td>C11</td><td>1.0 × 10<sup>−9</sup></td><td>60</td><td>6.00 × 10<sup>13</sup></td><td>128 × 128 × 128</td><td>19.6 ± 0.9</td><td>16.5 ± 0.4</td></tr><tr><td>C12</td><td>1.0 × 10<sup>−10</sup></td><td>60</td><td>1.29 × 10<sup>15</sup></td><td>128 × 128 × 128</td><td>19.4 ± 1.0</td><td>16.6 ± 0.5</td></tr><tr><td>C13</td><td>1.0 × 10<sup>−11</sup></td><td>60</td><td>2.78 × 10<sup>16</sup></td><td>128 × 128 × 128</td><td>19.7 ± 1.0</td><td>16.5 ± 0.4</td></tr><tr><td>C14</td><td>1.0 × 10<sup>−12</sup></td><td>60</td><td>6.00 × 10<sup>17</sup></td><td>128 × 128 × 128</td><td>19.5 ± 1.0</td><td>16.5 ± 0.4</td></tr><tr><td>C15</td><td>1.0 × 10<sup>−15</sup></td><td>60</td><td>6.00 × 10<sup>21</sup></td><td>128 × 128 × 128</td><td>19.4 ± 1.0</td><td>16.7 ± 0.5</td></tr><tr><td>D1</td><td>1.0 × 10<sup>−1</sup></td><td>80</td><td>1.72 × 10<sup>3</sup></td><td>192 × 192 × 192</td><td>0.95 ± 0.05</td><td>2.4 ± 0.1</td></tr><tr><td>D2</td><td>1.0 × 10<sup>−2</sup></td><td>80</td><td>3.71 × 10<sup>4</sup></td><td>192 × 192 × 192</td><td>3.8 ± 0.1</td><td>6.8 ± 0.1</td></tr><tr><td>D3</td><td>1.0 × 10<sup>−3</sup></td><td>80</td><td>8.00 × 10<sup>5</sup></td><td>192 × 192 × 192</td><td>9.7 ± 0.2</td><td>11.6 ± 0.2</td></tr><tr><td>D4</td><td>1.0 × 10<sup>−4</sup></td><td>80</td><td>1.72 × 10<sup>7</sup></td><td>192 × 192 × 192</td><td>20.7 ± 0.6</td><td>17.7 ± 0.3</td></tr><tr><td>D5</td><td>1.0 × 10<sup>−5</sup></td><td>80</td><td>3.71 × 10<sup>8</sup></td><td>192 × 192 × 192</td><td>33.8 ± 0.9</td><td>25.7 ± 0.9</td></tr><tr><td>D6</td><td>1.0 × 10<sup>−6</sup></td><td>80</td><td>8.00 × 10<sup>9</sup></td><td>192 × 192 × 192</td><td>38.9 ± 0.9</td><td>24.4 ± 0.7</td></tr><tr><td>D7</td><td>1.0 × 10<sup>−7</sup></td><td>80</td><td>1.72 × 10<sup>11</sup></td><td>192 × 192 × 192</td><td>34.4 ± 1.4</td><td>24.4 ± 0.7</td></tr><tr><td>D8</td><td>1.0 × 10<sup>−8</sup></td><td>80</td><td>3.71 × 10<sup>12</sup></td><td>192 × 192 × 192</td><td>30.8 ± 1.1</td><td>22.8 ± 0.8</td></tr><tr><td>D9</td><td>1.0 × 10<sup>−9</sup></td><td>80</td><td>8.00 × 10<sup>13</sup></td><td>192 × 192 × 192</td><td>28.9 ± 1.2</td><td>22.3 ± 0.6</td></tr><tr><td>D10</td><td>1.0 × 10<sup>−10</sup></td><td>80</td><td>1.72 × 10<sup>15</sup></td><td>192 × 192 × 192</td><td>29.1 ± 1.6</td><td>22.6 ± 0.7</td></tr><tr><td>D11</td><td>1.0 × 10<sup>−11</sup></td><td>80</td><td>3.71 × 10<sup>16</sup></td><td>192 × 192 × 192</td><td>28.5 ± 1.4</td><td>22.4 ± 0.8</td></tr><tr><td>D12</td><td>1.0 × 10<sup>−12</sup></td><td>80</td><td>8.00 × 10<sup>17</sup></td><td>192 × 192 × 192</td><td>29.1 ± 1.2</td><td>22.8 ± 0.7</td></tr><tr><td>D13</td><td>1.0 × 10<sup>−15</sup></td><td>80</td><td>8.00 × 10<sup>21</sup></td><td>192 × 192 × 192</td><td>28.5 ± 1.4</td><td>22.7 ± 0.7</td></tr></tbody></table>
<table><caption>Table 6: List of simulations described in this work (continued 2); see also Tables 4 and 5.</caption><thead><tr><th>Simulation nr.</th><th>Ek</th><th>Ra</th><th>N<sub>x</sub> × N<sub>y</sub> × N<sub>Z</sub></th><th>(Nu − 1) ± ΔNu</th><th>Re<sub>w</sub> ± ΔRe<sub>w</sub></th></tr></thead><tbody><tr><td>E1</td><td>1.0 × 10<sup>−1</sup></td><td>120</td><td>2.59 × 10<sup>3</sup></td><td>128 × 128 × 128</td><td>3.4 ± 0.1</td></tr><tr><td>E2</td><td>1.0 × 10<sup>−2</sup></td><td>120</td><td>5.57 × 10<sup>4</sup></td><td>128 × 128 × 128</td><td>8.7 ± 0.1</td></tr><tr><td>E3</td><td>1.0 × 10<sup>−3</sup></td><td>120</td><td>1.2 × 10<sup>6</sup></td><td>128 × 128 × 128</td><td>14.8 ± 0.2</td></tr><tr><td>E4</td><td>1.0 × 10<sup>−4</sup></td><td>120</td><td>2.59 × 10<sup>7</sup></td><td>128 × 128 × 128</td><td>22.7 ± 0.4</td></tr><tr><td>E5</td><td>1.0 × 10<sup>−5</sup></td><td>120</td><td>5.57 × 10<sup>8</sup></td><td>256 × 256 × 256</td><td>30.2 ± 0.9</td></tr><tr><td>E6</td><td>1.0 × 10<sup>−6</sup></td><td>120</td><td>1.2 × 10<sup>10</sup></td><td>256 × 256 × 256</td><td>37.9 ± 1.6</td></tr><tr><td>E7</td><td>1.0 × 10<sup>−7</sup></td><td>120</td><td>2.59 × 10<sup>11</sup></td><td>256 × 256 × 256</td><td>39.6 ± 1.1</td></tr><tr><td>E8</td><td>1.0 × 10<sup>−8</sup></td><td>120</td><td>5.57 × 10<sup>12</sup></td><td>256 × 256 × 256</td><td>40.5 ± 2.1</td></tr><tr><td>E9</td><td>1.0 × 10<sup>−9</sup></td><td>120</td><td>1.2 × 10<sup>14</sup></td><td>256 × 256 × 256</td><td>40.1 ± 2.2</td></tr><tr><td>E10</td><td>1.0 × 10<sup>−10</sup></td><td>120</td><td>2.59 × 10<sup>15</sup></td><td>256 × 256 × 256</td><td>39.7 ± 1.5</td></tr><tr><td>E11</td><td>1.0 × 10<sup>−11</sup></td><td>120</td><td>2.59 × 10<sup>15</sup></td><td>256 × 256 × 256</td><td>38.7 ± 1.0</td></tr><tr><td>E12</td><td>1.0 × 10<sup>−12</sup></td><td>120</td><td>1.20 × 10<sup>18</sup></td><td>256 × 256 × 256</td><td>40.3 ± 1.1</td></tr><tr><td>E13</td><td>1.0 × 10<sup>−15</sup></td><td>120</td><td>1.20 × 10<sup>22</sup></td><td>256 × 256 × 256</td><td>40.4 ± 1.0</td></tr><tr><td>F1</td><td>Ek ≪ 1</td><td>40</td><td>Ra ≫ 1</td><td>128 × 128 × 128</td><td>10.6 ± 0.3</td></tr><tr><td>F2</td><td>Ek ≪ 1</td><td>60</td><td>Ra ≫ 1</td><td>128 × 128 × 128</td><td>16.7 ± 0.5</td></tr><tr><td>F3</td><td>Ek ≪ 1</td><td>80</td><td>Ra ≫ 1</td><td>128 × 128 × 128</td><td>23.1 ± 0.8</td></tr><tr><td>F4</td><td>Ek ≪ 1</td><td>120</td><td>Ra ≫ 1</td><td>256 × 256 × 256</td><td>40.9 ± 2.1</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 10: Reynolds number $Re_w$, defined in Eq. (30), versus $Ta = Ek^{-2}$ at fixed $\tilde{Ra} \equiv RaEk^{4/3} = 60$. The dashed lines show the corresponding mean value obtained from the asymptotically reduced equations, with the shaded area indicating one standard deviation above and below that value.</figcaption></figure>

flow statistics obtained in the low $Ek$ regime are independent of the time step using runs with CFL prefactor 0.4 or 0.1 instead of 0.2 (which was used in all other runs) using $Ek = 10^{-15}$, $\tilde{Ra} = 80$, both of which gave the same Nusselt and Reynolds numbers within one standard deviation (not shown).

The above simulations were performed with between 128 and 384 dealiased Fourier modes in the $x$ and $y$ directions and between 128 and 256 dealiased Chebyshev modes in the vertical direction. It was verified for each simulation that the thermal boundary layer (defined in terms of the root-mean-square temperature fluctuation, cf. [59]) was well resolved, with at least 10 grid points.

Beyond the convergence of the Nusselt and Reynolds numbers to the values predicted by the asymptotic equations, an interesting pattern emerges. Both $Nu$ and $Re$ are small when $Ta$ is small (weak rotation). As $Ta$ increases, $Nu$ and $Re$ increase as well and exhibit an overshoot before converging to the asymptotic value. The amplitude of the overshoot is seen to decrease with increasing $\tilde{Ra}$. Similar results were very recently reported in [74] for high $Pr$ rotating Rayleigh-Bénard convection with no-slip boundary conditions, but the physical origin of these features remains to be explained. An enhancement of the Nusselt number with increasing rotation rate has also been observed for $Pr = 4.38$ and 6.4 [75]. The RRRiNSE formulation allows us, for the first time, to
observe the full range of $Ek$ from order one values down to the asymptotic regime within a single code, opening the door to detailed numerical studies of the classical problem of rapidly rotating convection, which has long posed a major challenge to the fluid dynamics community.

### 7.5. Visualizations

Here, we provide some visualizations of the various fields from run D12, at the low Ekman number of $Ek = 10^{-15}$ with $Pr = 1$ and $\tilde{Ra} = 80$. This run is in the geostrophic turbulence regime and the visualizations are produced in the statistically steady state, where a pronounced large-scale vortex (LSV) is present. The software Vapor [76] was used to generate the visualizations. We indicate positive values by orange and red contours and negative values by light and dark blue contours.

The left column of Fig. 11 shows the pressure field, where the large-scale columnar vortices are clearly visible. The right column of Fig. 11 shows the $y$-component $v$ of the velocity, which displays smaller-scale features than the pressure field. Figure 12 shows the vertical vorticity $\omega_z$ (left column) and the vertical component $w$ of the velocity (right column). In the $u$ field, a clear trace of the large-scale vortex column is visible, while this is less obvious in the $w$ snapshot. Finally, Fig. 13 shows the temperature perturbation field $\theta$ (left column) and the ageostrophic $x$ velocity $U$. Both fields show small-scale structures without any visible trace of the large-scale vortex.

## 8. Conclusions

In this work, we introduced the *Rescaled Rapidly Rotating incompressible Navier-Stokes Equations* referred to as RRRiNSE – a new formulation of the Navier-Stokes equations in the Boussinesq approximation describing rotating Rayleigh-Bénard convection, informed by the scalings valid in the asymptotic limit $Ek \to 0$. We solved these equations for stress-free boundary conditions using the quasi-inverse method to perform efficient DNS in a previously unattainable parameter regime of extremely small Ekman numbers $Ek$ relevant to geophysical and astrophysical fluid flows. We showed that the reduced equations of motion derive their increased efficiency from being well conditioned, thereby eliminating spurious growing modes that otherwise lead to numerical instabilities at small $Ek$. We have validated our simulation results against published results in the literature, and showed that the vertical diffusion terms can be treated either implicitly or explicitly for small $Ek$ owing to their smallness. We demonstrated for the first time that the full DNS of the RRRiNSE
<figure><img src="image_12.png" /><figcaption>Figure 11: Snapshots of the pressure (left column) and y-component $v$ of the velocity (right column) from run D12 (with $Ek = 10^{-15}$, $Ra = 80$, $Pr = 1$) in the steady state, where a saturated LSV is present. The axes in all panels are indicated by black arrows. Top row: side view. Bottom row: same data as in the top row (viewed from top). Blue color indicates negative values while orange and red colors indicate positive values.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Snapshots of the vertical vorticity $\omega_z$ (left column) and the vertical velocity $w$ from run D12 (with $Ek = 10^{-15}$, $Ra = 80$, $Pr = 1$) in the steady state, where a saturated LSV is present. The orientation in all panels is identical to Fig. 11. Top row: side view of $\omega_z$ (left) and $w$ (right). Bottom row: same data as in top row (viewed from top). Blue color indicates negative values while orange and red colors indicate positive values.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Snapshots of the temperature fluctuation $θ$ (left column) and the ageostrophic velocity component $U$ (right column) from run D12 (with $Ek = 10^{-15}$, $Ra = 80$, $Pr = 1$) in the steady state, where a saturated LSV is present. The orientation in all panels is identical to Fig. 11. Top row: side view. Bottom row: same data as in the top row, viewed from top. Blue color indicates negative values while orange and red colors indicate positive values.</figcaption></figure>
converge to the asymptotically reduced equations for small $Ek$, and showed that the time derivative in the mean temperature is inconsequential for the accurate determination of the average Nusselt number in the statistically stationary state, thus allowing a reduction by orders of magnitude in the simulation time required.

The results presented here provide an important advance in the numerical treatment of rotating convection in the rapid rotation regime, which will make it possible to explore for the first time the previously unattainable parameter regime of small but finite $Ek$man and $Rossby$ numbers. The first study applying the approach described here to this regime, exploring the flow physics found there, is given in [77]. Future studies will address the physics of the transition to the asymptotic parameter regime and the properties of optimal heat transport in rapidly rotating convection and the associated Reynolds numbers. Another direction for future investigation concerns the possibility of misalignment between the rotation axis and gravity, which has previously been studied in the context of the tilted $f$-plane using both the asymptotically reduced equations of motion [60, 78] and the full rotating Boussinesq equations [79, 80, 81].

## Acknowledgements

K. J. passed away before this manuscript was finalized. We have attempted to present the results of our collaboration in accordance with his high standards. Any errors or misinterpretations remain our own. This work was supported by the National Science Foundation (Grants DMS-2009319 and DMS-2308338 (KJ), Grants DMS-2009563 and DMS-2308337 (EK)), by the German Research Foundation (DFG Projektnummer: 522026592, AvK) and by Agence Nationale de le Recherche (Grant ANR-23-CE30-0016-01, BM). This research used the Savio computational cluster resource provided by the Berkeley Research Computing program at the University of California Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor for Research, and Chief Information Officer). This research also utilized the Alpine high performance computing resource at the University of Colorado Boulder. Alpine is jointly funded by the University of Colorado Boulder, the University of Colorado Anschutz, and Colorado State University. Data storage for this project was supported by the University of Colorado Boulder PetaLibrary. This project also made use of computational resources from TGCC provided by GENCI (2024-A0162A10803).
Declaration of interests

The authors declare no conflict of interests.

**Appendix A. Mixed vorticity-velocity formulation of the RRRiNSE**

The primitive variable formulation of the RRRiNSE in terms of $\mathbf{u} = (u, v, w)$, $U_{\perp} = (U, V)$, $\pi$, and $\bar{\Theta}, \theta$ given by equations (20) in the main text is of 11$^{th}$ order in $Z$. Specifically, the continuity equation requires, e.g., the imposition of an 11$^{th}$ auxiliary boundary condition applied to the pressure function. Instead of pursuing this option, we numerically solve the following modified set of equations for the variables $\mathbf{u} = (u, v, w)$, $\omega = (\omega_x, \omega_y, \omega_z)$, $U_{\perp} = (U, V)$, $\pi$ and $\bar{\Theta}, \theta$:

$$U = \frac{1}{\epsilon}(u + \partial_y \pi), \qquad V = \frac{1}{\epsilon}(v - \partial_x \pi), \qquad \omega_z = \partial_x v - \partial_y u, \qquad (A.1a)$$

$$\epsilon \partial_Z v - \partial_y w + \omega_x = 0, \qquad (A.1b)$$

$$\epsilon \partial_Z u - \omega_y - \partial_x w = 0, \qquad (A.1c)$$

$$\partial_x U + \partial_y V + \partial_z w = 0, \qquad (A.1d)$$

$$\partial_t u - V - \mathcal{D}_u = - \mathcal{N}_u, \qquad (A.1e)$$

$$\partial_t v + U - \mathcal{D}_v = - \mathcal{N}_v, \qquad (A.1f)$$

$$\partial_t w + \partial_z \pi - \frac{\widetilde{Ra}}{Pr} \theta - \mathcal{D}_w = - \mathcal{N}_w, \qquad (A.1g)$$

$$\partial_t \theta + (\partial_z \bar{\Theta} - 1)w - \mathcal{D}_\theta = - \mathcal{N}_\theta, \qquad (A.1h)$$

$$\epsilon^{-2} \partial_t \bar{\Theta} - \mathcal{D}_\Theta = - \mathcal{N}_\Theta, \qquad (A.1i)$$

where the linear diffusion and nonlinear advection terms are given by

$$-\mathcal{D}_u = \partial_y \omega_z - \epsilon \partial_z \omega_y, \qquad -\mathcal{N}_u = \omega_z v - \omega_y w,$$

$$-\mathcal{D}_v = \epsilon \partial_z \omega_x - \partial_x \omega_z, \qquad -\mathcal{N}_v = \omega_x w - \omega_z u,$$

$$-\mathcal{D}_w = \partial_x \omega_y - \partial_y \omega_x, \qquad -\mathcal{N}_w = \omega_y u - \omega_x v,$$

$$-\mathcal{D}_\theta = - \frac{1}{Pr}(\partial_x^2 + \partial_y^2 + \epsilon^2 \partial_z^2), \qquad -\mathcal{D}_\theta = - \partial_x(u\theta) - \partial_y(v\theta) - \epsilon \partial_z(w\theta),$$

$$-\mathcal{D}_\Theta = - \frac{1}{Pr} \partial_z^2, \qquad -\mathcal{N}_\Theta = - \partial_z(\overline{w\theta}).$$
The above equations are of 10th order in Z and do not require an auxiliary pressure boundary condition. We apply impenetrable, stress-free, fixed-temperature boundary conditions at the top and bottom which provide the 10 conditions

$$w = \omega_x = \omega_y = \theta = \bar{\Theta} = 0 \text{ at } Z = 0, 1. \qquad (A.2)$$

An immediate consequence of this formulation is the fact that $\partial_Z p = 0$ on $Z = 0, 1$. We note that the CFL constraints on the linear and nonlinear terms are identical to those presented in Table 2.

## Appendix B. The quasi-inverse method with Chebyshev-Galerkin bases

We illustrate how the direction $z$ is treated in Coral using the quasi-inverse method and Galerkin bases. This technique can readily be applied to coupled sets of equations of arbitrary order. For brevity and clarity, however, we consider the simple case of the second-order, scalar heat equation:

$$\partial_t \phi - (\partial_{zz} - k_{\perp}^2) \phi = b(z) \qquad (B.1)$$

on the interval $z \in [-1, 1]$, where the right-hand side $b$ contains explicit contributions (e.g. source terms or advection).

We suppose this second-order equation is supplemented with two linear and homogeneous boundary conditions, a common case in fluid mechanics. By computing linear combinations of $N$ Chebyshev polynomials $(T_n)_{0 \le n < N}$, one defines a Galerkin family of function $(\Phi_m)_{2,N}$:

$$\Phi_m(z) = \sum_{0 \le n < N} S_{mn} T_n(z), \qquad (B.2)$$

each of which obeys the linear, homogeneous boundary conditions. It is crucial to note here that, as a result of enforcing these two boundary conditions, the Galerkin basis has been reduced as compared with the initial Chebyshev basis and now contains only $N - 2$ polynomials. Next, we expand the variable $\phi$ in this Galerkin basis:

$$\phi(z, t) = \sum_{2 \le m < N} \tilde{\phi}_m(t) \Phi_m(z). \qquad (B.3)$$

The standard discretization of Eq. (B.1) consists in using expansion (B.3) and projecting on Chebyshev polynomials. Owing to the presence of derivatives $\partial_z$, this is conducive to dense (triangular) and, perhaps more importantly, ill-conditioned matrices [82]. The spirit of the quasi-inverse method
consists in integrating the differential equation repeatedly, until the reformulated problem is clear of derivatives. In our case, we integrate Eq. (B.1) with respect to z twice:

$$ \partial_t \iint \phi - \left(1 - \iint k_{\perp}^2\right) \phi = \iint b(z) + a_0 + a_1 z, \quad (B.4) $$

where $a_0$ and $a_1$ are two arbitrary integration constants. Fortunately, these unknown constants appear in (and pollute) the $T_0(z)$ (constant) and $T_1(z)$ (linear) projections only. By projecting Eq. (B.4) on the unpolluted $N-2$ higher Chebyshev polynomials, one obtains an algebraic system for the $N-2$ unknown Galerkin coefficients $\tilde{\phi}_m$. Denoting the natural scalar product for Chebyshev polynomials with brackets, $\langle...\rangle$, we have for $2 \le m, p < N$:

$$ \sum_{0 \le n < N} S_{mn} \left\langle T_p(z), \iint T_n(z) \right\rangle \partial_t \tilde{\phi}_m - S_{mp} \tilde{\phi}_m + k_{\perp}^2 \sum_{0 \le n < N} S_{mn} \left\langle T_p(z), \iint T_n(z) \right\rangle \tilde{\phi}_m = \left\langle T_p(z), \iint b(z) \right\rangle, \quad (B.5) $$

where we have used the orthonormality condition $\langle T_p(z), T_n(z) \rangle = \delta_{pn}$. Crucially, the matrix representing the double integration,

$$ \left\langle T_p(z), \iint T_n(z) \right\rangle, \quad (B.6) $$

is penta-diagonal and, more importantly, well-conditioned.

Finally, a discussion is in order considering the Galerkin stencil $S_{mn}$. Some care must be taken when defining the Galerkin basis, among all the possibilities. Considering the simple case of Dirichlet boundary conditions on both boundaries, one may be tempted by the following simple recombination:

$$ \Phi_{2p}(z) = T_{2p}(z) - T_0(z) \quad \text{and} \quad \Phi_{2p+1}(z) = T_{2p+1}(z) - T_1(z). \quad (B.7) $$

However, a dense discretization would result and therefore this stencil should be avoided. Instead, one should use

$$ \Phi_m(z) = T_m(z) - T_{m-2}(z), \quad (B.8) $$

which is a well-conditioned and banded stencil. Thus, all coupling matrices appearing in Eq. (B.5) are also banded and the system can be efficiently marched in time implicitly, e.g. with a Runge-Kutta scheme.
We emphasize that this procedure, exemplified on a simple scalar equation, can be generalized to systems of coupled PDEs without noticeable difficulty (but at the cost of increased book-keeping), as long as linear and homogeneous boundary conditions are imposed.

## Appendix C. Implicit-explicit time discretizations and the quasi-inverse method

In this appendix, we summarize the specific formulation of implicit-explicit time discretization for both the NHQG model (15) and the rescaled equations (20) and (A.1). In all generality, these governing equations are represented by a system of the form:

$$ (\partial_t \mathcal{M} - \mathcal{L}_I) \mathbf{v}^{(n+1)} = \mathcal{L}_E \mathbf{v}^{(n)} + \mathcal{N}(\mathbf{v}^{(n)}, \mathbf{v}^{(n)}) + \mathcal{F}_{\theta}^{(n)}. \quad (C.1) $$

where the superscripts (n + 1) and (n) denote implicit (unknown) and explicit (known) variables from a prior time step, and $\mathbf{v}$ represents the dependent variables associated with the hydrodynamics problem only. Generalization to the full problem including thermal effects (27) is straightforward and omitted for brevity. We summarize below the expressions for the various operators appearing in this equation for the different equations studied here.

### A. NHQG-RRBC: $\mathbf{v} = (\Psi, w)^T$ and differential operators for equations (15a,b)

$$ \mathcal{M} = \begin{pmatrix} \nabla_{\perp}^2 & 0 \\ 0 & 1 \end{pmatrix}, \quad \mathcal{L}_I = \begin{pmatrix} \nabla_{\perp}^4 & \partial_Z \\ -\partial_Z & \nabla_{\perp}^2 \end{pmatrix}, \quad \mathcal{L}_E = \mathbf{0}_2, \quad \mathcal{N} = \begin{pmatrix} J[\Psi, \cdot] & 0 \\ 0 & J[\Psi, \cdot] \end{pmatrix}. \quad (C.2) $$

### B. RiNSE: Primitive variable formulation: $\mathbf{v}^{(j)} = (\mathbf{u}, \mathbf{U}_{\perp}, \pi)$ for equations (20a-d) and

$$ \mathcal{M} = \begin{pmatrix} \mathcal{I}_3 & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{pmatrix}, \quad \mathcal{L}_I = \begin{pmatrix} -\nabla_{\perp}^2 \mathcal{I}_3 & \begin{pmatrix} \mathcal{J}_2 & 0 \\ 0 & 0 \end{pmatrix} \\ \begin{pmatrix} \mathcal{J}_2 & 0 \\ 0 & 0 \end{pmatrix} & \begin{pmatrix} -\epsilon \mathcal{J}_2 & \nabla_{\perp}^T \\ \nabla_{\perp} & 0 \end{pmatrix} \end{pmatrix}, \quad (C.3) $$

$$ \mathcal{L}_E = \epsilon^2 \begin{pmatrix} \mathcal{I}_3 \partial_Z^2 & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{pmatrix}, \quad \mathcal{N} = - \begin{pmatrix} (\mathbf{u}_{\perp} \cdot \nabla_{\perp} + \epsilon w \partial_Z) \mathcal{I}_3 & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{pmatrix}. \quad (C.4) $$

Here $\mathcal{I}_3$ is the order three identity matrix and

$$ \mathcal{J}_2 = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}. \quad (C.5) $$
<figure><img src="image_15.png" /><figcaption>Figure C.14: Spy plots of the mass $\mathcal{M}$ (left) and stiffness $\mathcal{L}_I$ (center) matrices, both in banded format. The fine structure of the stiffness matrix is illustrated in the close-up in the right panel.</figcaption></figure>

Figure C.14 demonstrates the sparsity of the quasi-inverse approach via spyplots for the mass and stiffness matrices $\mathcal{M}$ and $\mathcal{L}_I$, respectively.

C. RiNSE: Mixed vorticity-velocity formulation: $\mathbf{v}^{(j)} = (\mathbf{u} | \mathbf{U}_{\perp}, \pi | \omega)$ (see Appendix A) and

$$ \mathcal{M} = \begin{pmatrix} \mathcal{I}_3 & \mathbf{0}_3 & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 & \mathbf{0}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 & \mathbf{0}_3 \end{pmatrix}, \quad -\mathcal{L}_E = \epsilon \begin{pmatrix} \mathbf{0}_3 & \mathbf{0}_3 & \begin{pmatrix} \partial_Z \mathcal{J}_2 & 0 \\ 0 & 0 \end{pmatrix} \\ \mathbf{0}_3 & \mathbf{0}_3 & \mathbf{0}_3 \\ \begin{pmatrix} \partial_Z \mathcal{J}_2 & 0 \\ 0 & 0 \end{pmatrix} & \mathbf{0}_3 \end{pmatrix}, \quad (C.6) $$

$$ -\mathcal{L}_I = \begin{pmatrix} \mathbf{0}_3 & \begin{pmatrix} \mathcal{J}_2 & 0 \\ 0 & \partial_Z \end{pmatrix} & \begin{pmatrix} 0 & -\nabla^{\perp T} \\ \nabla^{\perp} & 0 \end{pmatrix} \\ \begin{pmatrix} \mathcal{J}_2 & 0 \\ 0 & \partial_Z \end{pmatrix} & \begin{pmatrix} -\epsilon \mathcal{J}_2 & \nabla_{\perp}^T \\ \nabla_{\perp} & 0 \end{pmatrix} & \mathbf{0}_3 \\ \begin{pmatrix} 0 & -\nabla^{\perp T} \\ \nabla^{\perp} & 0 \end{pmatrix} & \mathbf{0}_3 & \mathcal{I}_3 \end{pmatrix}, \quad (C.7) $$

and $\mathcal{N} = \text{Diag}[\{\mathcal{N}_u, \mathcal{N}_v, \mathcal{N}_w\}, \mathbf{0}_3, \mathbf{0}_3]$.

From $\mathcal{L}_I$ we see that the system (C.1) is in each case of second order in $Z$ requiring impenetrable boundary conditions $w = 0$ at $Z = 0, 1$. For cases B and C stress-free boundary conditions
($\partial_Z \mathbf{u}_\perp = 0$) are enforced via an appropriate Chebyshev-Galerkin basis $\Phi_j(Z)$ for each variable (no-slip boundary conditions, $\mathbf{u}_\perp = 0$, can also be considered in this approach, but are not analyzed here). The above holds regardless of whether $\mathcal{L}_E$ is treated explicitly or implicitly. For case B only, an implicit treatment of $\mathcal{L}_E$ increases the order of the system to seven and an additional auxillary boundary condition on the pressure is required. This does not occur for Case C which is preferred.

## Appendix D. Analysis of the mean temperature equation

The mean temperature equation

$$ \epsilon^{-2} \partial_t \bar{\Theta} + \partial_Z \left( \overline{w\theta} - \frac{1}{Pr} \partial_Z \bar{\Theta} \right) = 0 \qquad (D.1) $$

in a statistically stationary state implies

$$ Nu_t - 1 = Pr \overline{w\theta}^t - \overline{\partial_Z \bar{\Theta}}^t \quad \Rightarrow \quad Nu_t - 1 = Pr \left\langle \overline{w\theta}^t \right\rangle_Z . \qquad (D.2) $$

With this interpretation $Nu_t$ is strictly a constant. It follows

$$ \epsilon^{-2} \partial_t \bar{\Theta} + \partial_Z \left( \overline{w\theta} - \overline{\overline{w\theta}^t} - \frac{1}{Pr} \left( \partial_Z \bar{\Theta} - \overline{\partial_Z \bar{\Theta}^t} \right) \right) = 0 \qquad (D.3) $$

and given $\epsilon^2 \bar{\Theta}_{2+} \equiv \sum_{j \ge 2} \epsilon^j \bar{\Theta}_j$, $\partial_t \bar{\Theta}_{(0,1)} = 0$ such that $\partial_Z \bar{\Theta}_{(0,1)} = \partial_Z \overline{\bar{\Theta}}_{(0,1)}^t$. This implies

$$ \partial_t \bar{\Theta}_{2+} + \partial_Z \left( \overline{w\theta} - \overline{\overline{w\theta}^t} - \frac{\epsilon^2}{Pr} \left( \partial_Z \bar{\Theta}_{2+} - \overline{\partial_Z \bar{\Theta}_{2+}^t} \right) \right) = 0. \qquad (D.4) $$

To leading order

$$ \partial_t \bar{\Theta}_2 + \partial_Z \left( \overline{w\theta} - \overline{\overline{w\theta}^t} \right) \approx 0 \qquad (D.5) $$

indicating that $O(1)$ fluctuations in the heat transport about the mean $Nu_t$ are accounted for by mean temporal variations in the mean temperature at $O(\epsilon^2)$, i.e., $\bar{\Theta}_2$. However, for numerical efficiency it is found that the temporal fluctuations of $\epsilon^{-2} \partial_t \bar{\Theta}$ can be neglected. Hence,

$$ \partial_Z \left( \overline{w\theta} - \frac{1}{Pr} \partial_Z \bar{\Theta} \right) = 0, \qquad (D.6) $$

resulting in the time-dependent Nusselt number

$$ Nu(t) - 1 = Pr \overline{w\theta} - \partial_Z \bar{\Theta}. \qquad (D.7) $$

This implies

$$ Nu(t) - 1 = Pr \langle \overline{w\theta} \rangle_Z \quad \Rightarrow \quad \overline{Nu(t)}^t - 1 = Pr \langle \overline{w\theta} \rangle_Z^t. \qquad (D.8) $$
The difference in averaged Nusselt numbers from the two methods is given by

$$ \overline{Nu(t)}^t - Nu_t = Pr \left( \left\langle \overline{w\theta} \right\rangle_Z^t - \left\langle \overline{w\theta} \right\rangle_Z^t \right) = \left( \partial_Z \overline{\Theta} - \partial_Z \overline{\Theta}^t \right) \Big|_{0,1} . \quad (D.9) $$

If the operations of depth- and time-averaging commute then assuming equivalence in the thermal and velocity mean statistics implies $\overline{Nu(t)}^t = Nu_t$.

Moreover, if $\epsilon \mapsto \epsilon^*$ such that $\epsilon^* > \epsilon$ then it follows

$$ Nu_t^* - Nu_t = Pr \left( \left\langle \overline{w\theta}^{*t} \right\rangle_Z - \left\langle \overline{w\theta} \right\rangle_Z^t \right) = \left( \partial_Z \overline{\Theta}^{*t} - \partial_Z \overline{\Theta}^t \right) \Big|_{0,1} . \quad (D.10) $$

If the period of time-averaging is sufficiently long then assuming equivalence in the thermal and velocity mean statistics this implies $Nu_t^* = Nu_t$.

It is found that the magnitudes $|\overline{Nu(t)}^t - Nu_t|$ or $|Nu_t^* - Nu_t|$ depend on the horizontal domain size upon which area-averaging is performed. Increasingly larger domains contain greater statistical sampling advantageous to the aforementioned commutation that results in asymptotic error convergence.

## Appendix E. Comparison of iNSE and RRRiNSE with and without slaving

Here we present results on different variations of the numerical approach presented in the main text. Specifically, we compare simulations of the iNSE and of the RRRiNSE with and without slaving of the mean temperature field. For the tests described below, we fix $\tilde{Ra} = 40$, while setting $N_x = N_y = 128$, $L_x = L_y = 10\ell_c$, and we employ the second order ARS-222 time-stepping scheme. We set the Ekman number to low values, either $Ek = 10^{-9}$ or $Ek = 10^{-12}$. We have tested different values for the vertical resolution $N_z$, since an increased $N_z$ is seen in Fig. 1 to lead to less accurate spectra in the unscaled equations.

The results of these tests are summarized in Fig. E.15 showing the time series of the kinetic energy in horizontal (top row) and vertical (bottom row) motions from different simulations. First, we note that the RRRiNSE are found to be unconditionally stable in all cases. In contrast, the standard iNSE are unstable when a slow adjustment of the vertical profile is permitted. This is seen to be the case when small amplitude, random perturbations are used for initial conditions at both $Ek = 10^{-9}$ (crimson line in the left column of Fig. E.15) and $Ek = 10^{-12}$ (purple line in the right column of Fig. E.15). In addition to using noise as initial condition, we have also considered
an equilibrated RRRiNSE solution, which was suitably rescaled and used as an initial condition for continued time-marching of the solution with the unscaled iNSE. This procedure yielded the following observations: including the slow mean temperature time derivative (e.g. purple line on the left panels) is associated with eventual blow-up. Rather surprisingly, despite the presence of spurious modes polluting the numerical spectra, we have not observed a numerical blow-up when starting from an equilibrated solution and simultaneously omitting the slow mean temperature time derivative (e.g. salmon lines on both panels). When all terms in the mean temperature equation are retained, including the slow time derivative of the mean temperature, the unscaled equations quickly diverge (see the red and purple lines in the left panel). When the mean temperature derivative is set to zero, the unscaled equations of motion remain surprisingly stable for the times investigated, despite the presence of spurious growth rate values in the spectra obtained when solving the numerically ill-conditioned generalized eigenproblem.

References

[1] C. A. Jones, Planetary magnetic fields and fluid dynamos, Annu. Rev. Fluid Mech. 43 (2011) 583–614.

[2] P. H. Roberts, E. M. King, On the genesis of the Earth’s magnetism, Rep. Prog. Phys. 76 (9) (2013) 096801.

[3] M. S. Miesch, The coupling of solar convection and rotation, Sol. Phys. 192 (2000) 59–89.

[4] A. P. Ingersoll, Atmospheric dynamics of the outer planets, Science 248 (4953) (1990) 308–315.

[5] K. A. Emanuel, et al., Atmospheric Convection, Oxford University Press, 1994.

[6] J. Marshall, F. Schott, Open-ocean convection: Observations, theory, and models, Rev. Geo-phys. 37 (1) (1999) 1–64.

[7] K. M. Soderlund, Ocean dynamics of outer solar system satellites, Geophys. Res. Lett. 46 (15) (2019) 8700–8710.

[8] S. Chandrasekhar, The instability of a layer of fluid heated below and subject to Coriolis forces, Proc. R. Soc. A 217 (1130) (1953) 306–327.

[9] Y. Nakagawa, P. Frenzen, A theoretical and experimental study of cellular convection in rotating fluids, Tellus 7 (1955) 2–21.
<figure><img src="image_16.png" /><figcaption>Figure E.15: Comparison of the RRRiNSE with the unscaled equations including (no approx. in legend) or omitting ($\partial_t\Theta = 0$ in legend) the slow mean temperature time derivative, showing time series of vertical kinetic energy (top row) and horizontal kinetic energy (bottom row). All quantities are plotted in rescaled units, regardless of the dimensionless formulation adopted for the governing equations. The horizontal kinetic energy grows as a consequence of inverse energy cascade.</figcaption></figure>
[10] G. Veronis, Cellular convection with finite amplitude in a rotating fluid, J. Fluid Mech. 5 (1959) 401-435.

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

[12] J. M. Pfotenhauer, P. G. J. Lucas, R. J. Donnelly, Stability and heat transfer of rotating cryogens. Part 2. Effects of rotation on heat-transfer properties of convection in liquid He, J. Fluid Mech. 145 (1984) 239-252.

[13] B. M. Boubnov, G. S. Golitsyn, Experimental study of convective structures in rotating fluids, J. Fluid Mech. 167 (1986) 503-531.

[14] F. Zhong, R. Ecke, V. Steinberg, Asymmetric modes and the transition to vortex structures in rotating Rayleigh-Bénard convection, Phys. Rev. Lett. 67 (1991) 2473-2476.

[15] B. M. Boubnov, G. S. Golitsyn, Convection in Rotating Fluids, Vol. 29, Springer Science & Business Media, 2012.

[16] K. Julien, S. Legg, J. McWilliams, J. Werne, Rapidly rotating turbulent Rayleigh-Bénard convection, J. Fluid Mech. 322 (1996) 243-273.

[17] E. Knobloch, Rotating convection: recent developments, Int. J. Eng. Sci. 36 (12-14) (1998) 1421-1450.

[18] J. E. Hart, S. Kittelman, D. R. Ohlsen, Mean flow precession and temperature probability density functions in turbulent rotating convection, Phys. Fluids 14 (2002) 955-962.

[19] P. Vorobieff, R. E. Ecke, Turbulent rotating convection: an experimental study, J. Fluid Mech. 458 (2002) 191-218.

[20] R. J. A. M. Stevens, H. J. H. Clercx, D. Lohse, Heat transport and flow structure in rotating Rayleigh-Bénard convection, Eur. J. Mech. B Fluids 40 (2013) 41-49.

[21] G. Schubert, K. M. Soderlund, Planetary magnetic fields: Observations and models, Phys. Earth Planet. Inter. 187 (2011) 92-108.

[22] S. B. Pope, Turbulent Flows, Cambridge University Press, Cambridge, 2000.
[23] U. Frisch, Turbulence: the Legacy of AN Kolmogorov, Cambridge University Press, 1995.

[24] J. M. Aurnou, M. A. Calkins, J. S. Cheng, K. Julien, E. M. King, D. Nieves, K. M. Soderlund, S. Stellmach, Rotating convective turbulence in Earth and planetary cores, Phys. Earth Planet. Inter. 246 (2015) 52–71.

[25] Y. Kaneda, T. Ishihara, High-resolution direct numerical simulation of turbulence, J. Turbulence 7 (20) (2006) 1–17.

[26] P. K. Yeung, X. M. Zhai, K. R. Sreenivasan, Extreme events in computational turbulence, Proc. Natl. Acad. Sci. U.S.A. 112 (41) (2015) 12633–12638.

[27] P. K. Yeung, K. Ravikumar, S. Nichols, R. Uma-Vaideswaran, GPU-enabled extreme-scale turbulence simulations: Fourier pseudo-spectral algorithms at the exascale using OpenMP offloading, Comput. Phys. Commun. 306 (2025) 109364.

[28] A. Blass, P. Tabak, R. Verzicco, R. J. A. M. Stevens, D. Lohse, The effect of Prandtl number on turbulent sheared thermal convection, J. Fluid Mech. 910 (2021) A37.

[29] C. J. Howland, G. S. Yerragolam, R. Verzicco, D. Lohse, Turbulent mixed convection in vertical and horizontal channels, J. Fluid Mech. 998 (2024) A48.

[30] J. Song, O. Shishkina, X. Zhu, Direct numerical simulations of rapidly rotating Rayleigh-Bénard convection with Rayleigh number up to $5 \times 10^{13}$, J. Fluid Mech. 989 (2024) A3.

[31] E. M. King, S. Stellmach, J. M. Aurnou, Heat transfer by rapidly rotating Rayleigh-Bénard convection, J. Fluid Mech. 691 (2012) 568–582.

[32] B. Favier, L. J. Silvers, M. R. E. Proctor, Inverse cascade and symmetry breaking in rapidly rotating Boussinesq convection, Phys. Fluids 26 (2014) 096605.

[33] C. Guervilly, D. W. Hughes, C. A. Jones, Large-scale vortices in rapidly rotating Rayleigh-Bénard convection, J. Fluid Mech. 758 (2014) 407–435.

[34] R. P. J. Kunnen, R. Ostilla-Mónico, E. P. Van Der Poel, R. Verzicco, D. Lohse, Transition to geostrophic convection: the role of the boundary conditions, J. Fluid Mech. 799 (2016) 413–432.
[35] B. Favier, E. Knobloch, Robust wall states in rapidly rotating Rayleigh-Bénard convection, J. Fluid Mech. 895 (2020) R1.

[36] A. J. A. Guzmán, M. Madonia, J. S. Cheng, R. Ostilla-Mónico, H. J. H. Clercx, R. P. J. Kunnen, Force balance in rapidly rotating Rayleigh-Bénard convection, J. Fluid Mech. 928 (2021) A16.

[37] J. Song, O. Shishkina, X. Zhu, Scaling regimes in rapidly rotating thermal convection at extreme Rayleigh numbers, J. Fluid Mech. 984 (2024) A45.

[38] R. G. Cooper, P. J. Bushby, C. Guervilly, Subcritical dynamos in rapidly rotating planar convection, Phys. Rev. Fluids 5 (2020) 113702.

[39] S. J. Mason, C. Guervilly, G. R. Sarson, Magnetoconvection in a rotating spherical shell in the presence of a uniform axial magnetic field, Geophys. Astrophys. Fluid Dyn. 116 (2022) 458-498.

[40] P. Kolhey, S. Stellmach, D. Heyner, Influence of boundary conditions on rapidly rotating convection and its dynamo action in a plane fluid layer, Phys. Rev. Fluids 7 (2022) 043502.

[41] D. Majumder, B. Sreenivasan, G. Maurya, Self-similarity of the dipole-multipole transition in rapidly rotating dynamos, J. Fluid Mech. 980 (2024) A30.

[42] J. Aubert, T. Gastine, A. Fournier, Spherical convective dynamos in the rapidly rotating asymptotic regime, J. Fluid Mech. 813 (2017) 558-593.

[43] R. Verzicco, P. Orlandi, A finite-difference scheme for three-dimensional incompressible flows in cylindrical coordinates, J. Comp. Phys. 123 (1996) 402-414.

[44] E. M. King, K. M. Soderlund, U. R. Christensen, J. Wicht, J. M. Aurnou, Convective heat transfer in planetary dynamo models, Geochem, Geophys. Geosys. 11 (2010) Q06016.

[45] M. Plumley, K. Julien, P. Marti, S. Stellmach, The effects of Ekman pumping on quasi-geostrophic Rayleigh-Bénard convection, J. Fluid Mech. 803 (2016) 51-71.

[46] X. Zhu, E. Phillips, V. Spandan, J. Donners, G. Ruetsch, J. Romero, R. Ostilla-Mónico, Y. Yang, D. Lohse, R. Verzicco, et al., AFiD-GPU: a versatile Navier-Stokes solver for wall-bounded turbulent flows on GPU clusters, Comput. Phys. Commun. 229 (2018) 199-210.
[47] P. Marti, M. A. Calkins, K. Julien, A computationally efficient spectral method for modeling core dynamics, Geochem. Geophys. Geosys. 17 (8) (2016) 3031–3053.

[48] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Res. 2 (2020) 023068.

[49] B. Miquel, Coral: A parallel spectral solver for fluid dynamics and partial differential equations, J. Open Source Softw. 66 (2021) 2978. doi:10.21105/joss.02978.

[50] R. P. J. Kunnen, The geostrophic regime of rapidly rotating turbulent convection, J. Turbulence 22 (2021) 267–296.

[51] R. E. Ecke, O. Shishkina, Turbulent rotating Rayleigh–Bénard convection, Annu. Rev. Fluid Mech. 55 (2023) 603–638.

[52] J. M. Aurnou, S. Horn, K. Julien, Connections between nonrotating, slowly rotating, and rapidly rotating turbulent convection transport scalings, Phys. Rev. Res. 2 (2020) 043115.

[53] S. Maffei, M. J. Krouss, K. Julien, M. A. Calkins, On the inverse cascade and flow speed scaling behaviour in rapidly rotating Rayleigh–Bénard convection, J. Fluid Mech. 913 (2021) A18.

[54] T. G. Oliver, A. S. Jacobi, K. Julien, M. A. Calkins, Small scale quasigeostrophic convective turbulence at large Rayleigh number, Phys. Rev. Fluids 8 (2023) 093502.

[55] E. K. Hawkins, J. S. Cheng, J. A. Abbate, T. Pilegard, S. Stellmach, K. Julien, J. M. Aurnou, Laboratory models of planetary core-style convective turbulence, Fluids 8 (4) (2023) 106.

[56] M. Madonia, A. J. Aguirre Guzmán, H. J. Clercx, R. P. Kunnen, Reynolds number scaling and energy spectra in geostrophic convection, J. Fluid Mech. 962 (2023) A36.

[57] J. A. Abbate, J. M. Aurnou, Rotating convective turbulence in moderate to high Prandtl number fluids, Geophys. Astrophys. Fluid Dyn. 117 (2023) 397–436.

[58] M. Sprague, K. Julien, E. Knobloch, J. Werne, Numerical simulation of an asymptotically reduced system for rotationally constrained convection, J. Fluid Mech. 551 (2006) 141–174.

[59] K. Julien, A. M. Rubio, I. Grooms, E. Knobloch, Statistical and physical balances in low Rossby number Rayleigh–Bénard convection, Geophys. Astrophys. Fluid Dyn. 106 (2012) 392–428.
[60] K. Julien, E. Knobloch, Strongly nonlinear convection cells in a rapidly rotating fluid layer: the tilted f-plane, J. Fluid Mech. 360 (1998) 141–178.

[61] K. Julien, E. Knobloch, Fully nonlinear three-dimensional convection in a rapidly rotating layer, Phys. Fluids 11 (1999) 1469–1483.

[62] K. Julien, E. Knobloch, Reduced models for fluid flows with strong constraints, J. Math. Phys. 48 (2007) 065405.

[63] K. Julien, J. M. Aurnou, M. A. Calkins, E. Knobloch, P. Marti, S. Stellmach, G. M. Vasil, A nonlinear model for rotationally constrained convection with Ekman pumping, J. Fluid Mech. 798 (2016) 50–87.

[64] S.-H. Peng, K. Hanjalic, L. Davidson, Large-eddy simulation and deduced scaling analysis of Rayleigh–Bénard convection up to Ra=10⁹, J. Turbulence 7 (2006) N66.

[65] N. Foroozani, J. J. Niemela, V. Armenio, K. R. Sreenivasan, Large eddy simulation of turbulent Rayleigh–Bénard convection in a cubic cell, in: Direct and Large-Eddy Simulation X, Springer, 2018, pp. 559–565.

[66] R. Samuel, R. Samtaney, M. K. Verma, Large-eddy simulation of Rayleigh–Bénard convection at extreme Rayleigh numbers, Phys. Fluids 34 (2022) 075133.

[67] K. Julien, E. Knobloch, J. Werne, A new class of equations for rotationally constrained flows, Theor. Comput. Fluid Dyn. 11 (1998) 251–261.

[68] A. M. Rubio, K. Julien, E. Knobloch, J. B. Weiss, Upscale energy transfer in three-dimensional rapidly rotating turbulent convection, Phys. Rev. Lett. 112 (2014) 144501.

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

[70] S. Stellmach, M. Lischper, K. Julien, G. Vasil, J. S. Cheng, A. Ribeiro, E. M. King, J. M. Aurnou, Approaching the asymptotic regime of rapidly rotating convection: Boundary layers versus interior dynamics, Phys. Rev. Lett. 113 (2014) 254501.

[71] E. L. Ortiz, The tau method, SIAM J. Numer. Anal. 6 (3) (1969) 480–492.
[72] U. M. Ascher, S. J. Ruuth, R. J. Spiteri, Implicit-explicit Runge-Kutta methods for time-
dependent partial differential equations, Appl. Numer. Math. 25 (2) (1997) 151–167, Special
Issue on Time Integration.

[73] K. Julien, M. Watson, Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods, J. Comp. Phys. 228 (2009) 1480–1503.

[74] M. Anas, P. Joshi, Critical Prandtl number for heat transfer enhancement in rotating convection, Phys. Rev. Lett. 132 (2024) 034001.

[75] R. Hartmann, G. S. Yerragolam, R. Verzicco, D. Lohse, R. J. A. M. Stevens, Optimal heat transport in rotating Rayleigh-Bénard convection at large Rayleigh numbers, Phys. Rev. Fluids 8 (2023) 083501.

[76] S. Li, S. Jaroszynski, S. Pearse, L. Orf, J. Clyne, VAPOR: A visualization package tailored to analyze simulation data in Earth System Science, Atmosphere 10 (9) (2019) 488.

[77] A. van Kan, K. Julien, B. Miquel, E. Knobloch, Bridging the Rossby number gap in rapidly rotating thermal convection, J. Fluid Mech., to appear (2025).

[78] K. Julien, A. Ellison, B. Miquel, M. Calkins, E. Knobloch, Quasi-geostrophic convection on the tilted $f$-plane., to be submitted to J. Fluid Mech. (2025).

[79] L. Novi, J. von Hardenberg, D. W. Hughes, A. Provenzale, E. A. Spiegel, Rapidly rotating Rayleigh-Bénard convection with a tilted axis, Phys. Rev. E 99 (2019) 053116.

[80] L. K. Currie, A. J. Barker, Y. Lithwick, M. K. Browning, Convection with misaligned gravity and rotation: simulations and rotating mixing length theory, Mon. Not. R. Astron. Soc. 493 (2020) 5233–5256.

[81] T. Cai, Penetrative convection for rotating Boussinesq flow in tilted $f$-planes, Astrophys. J. 898 (1) (2020) 22.

[82] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Courier Corporation, 2001.