# Evolution of weak, homogeneous turbulence subject to rotation and stratification: the weakly dispersive case

J. F. Scott†

Ecole Centrale de Lyon, CNRS, Université Claude Bernard Lyon 1, INSA Lyon, Laboratoire de Mécanique des Fluides et d'Acoustique (UMR5509), 36 avenue Guy de Collongue, 69134 Ecully CEDEX, France

## Abstract

This article follows on from Scott and Cambon (2024) and Scott (2025). Like those articles, it concerns weak, decaying homogeneous turbulence in a rotating, stably stratified fluid with constant Brunt-Vaisala frequency, $N$. The difference is that here we consider the case in which $\beta = 2\Omega/N$ is close to 1, where $\Omega$ is the rotation rate. Because this renders inertial-gravity waves only weakly dispersive, wave-turbulence theory, which played a prominent role in the earlier studies, no longer applies. Indeed, wave-turbulence analysis does not appear here. Nonetheless, much of the analytical framework, based on modal decomposition, carries over, as do many of the conclusions. The flow is expressed as a sum of wave and non-propagating (NP) modes and their weak-turbulence mode-amplitude evolution equations are derived for small $\beta - 1$. The NP component is found to evolve independently of the wave one, following an amplitude equation which is precisely that of the previous studies in the limit $\beta \to 1$. The NP component induces coupling between wave modes and, without it, the wave component has purely linear decay. The mode-amplitude equations are integrated numerically using a scheme similar to classical direct numerical simulation (DNS) and results given. We find an inverse energy cascade of the NP component, whereas the presence of that component induces a forward cascade, hence significant dissipation, of the wave component. Detailed results are given for the energy, energy spectra and energy fluxes of the two components.

## 1. Introduction

As in Scott and Cambon (2024) and Scott (2025) (henceforth referred to as [I] and [II]), this paper concerns weak (small Rossby or Froude number) unbounded, homogeneous turbulence subject to rotation and stable stratification with constant Brunt-Vaisala frequency, $N$, and a constant and vertical rotation vector, $\Omega$, whose norm gives the rotation rate, $\Omega = |\Omega|$. The approach used is a spectral one, i.e. based on Fourier transforms, and a general overview of homogeneous turbulence with this viewpoint in mind can be found in the book by Sagaut and Cambon (2018). The difference compared with [I] and [II] is that we focus on the case $\beta = 2\Omega/N$ close to 1, for which the weak-turbulence analysis used in the previous articles is no longer applicable.

<footnote>† julian.scott@ec-lyon.fr</footnote>
Different approaches have been used to address the problem of homogeneous turbulence with rotation and/or stratification. These include the Eddy-Damped, Quasi-Normal, Markovian (EDQNM) closure (see Cambon and Jacquin (1989), Godeferd and Cambon (1994), Cambon, Mansour and Godeferd (1997), Godeferd and Staquet (2003) for studies of turbulence with both rotation and stratification). Another is direct numerical simulation (DNS), which approximates the flow as periodic using truncated Fourier series (see e.g. Orszag and Patterson (1972), Canuto et al. (1988)). DNS studies have been undertaken of purely stratified turbulence (see e.g. Waite and Bartello (2004 and 2006a)) and of the purely rotating case (see e.g. Smith and Waleffe (1999), Bourouiba and Bartello (2007)). Of course, pure rotation or stratification does not include the case of small $\beta-1$ considered here. Applications of DNS to stratified, rotating turbulence can be found in Coleman, Ferziger and Spalart (1992), Smith and Waleffe (2002), Liechtenstein, Godeferd and Cambon (2005), Waite and Bartello (2006b), Thomas and Daniel (2020 and 2021) and Thomas (2023). Several of these papers discuss wave-NP mode interactions and energy transfers, these being central themes of the present work. This being said, none of these studies are specifically aimed at the case of weak turbulence with small $\beta-1$ studied here.

As discussed in the introduction of [I], in the case of weak turbulence, inertial-gravity wave oscillations lead to problems for classical DNS. Firstly, the time step must be sufficiently small as to resolve the oscillations, while the total time for evolution needs to be long enough that nonlinearity is significant. Since nonlinearity requires many wave periods for weak turbulence, the required number of time steps is considerably more that in the absence of waves. Secondly, DNS involves discretisation in spectral space. Wave oscillations in that space become increasingly rapid as time increases and the spectral resolution should be sufficiently fine that it resolves the oscillations at the times for which nonlinearity intervenes. These constraints mean that accuracy of classical DNS requires more and more computer resources as the turbulence weakens, hence the interest of studies, such as the present one, which deliberately focus on the weak-turbulence limit.

With weak turbulence in mind, the analysis of [I] expressed the flow as a sum over modes, which are solutions of the governing equations without nonlinearity or visco-diffusion (see Bartello (1995) and [I]). Modes have the form $\exp(i(k \cdot x - s\omega(k)t))$, where $\omega(k)$ is the dispersion relation of inertial-gravity waves. They are indexed by the wave vector $k$ and the integer $s=0, \pm 1$. Those with $s=\pm 1$ are inertial-gravity waves, while $s=0$ modes have zero frequency, hence zero group velocity when regarded as waves. Thus, the latter modes are non-propagating, here abbreviated to NP. They are often called vortical or geostrophic modes in the geophysical literature, where vortical refers to potential vorticity (see e.g. Pedlosky (1987)). In any case, there are two components of the flow, wave and NP, whose physics is rather different in the weak-turbulence limit.

Mode-amplitude equations were derived in [I]. These equations describe the temporal evolution of the mode amplitudes once nonlinearity and visco-diffusion are allowed for. A spectral matrix, whose diagonal elements give the energy density of the different modes in
spectral space was defined, the total energy being a sum of wave and NP contributions. This
analytical framework forms the basis of the present study and is summarised in section 2.

Assuming weak turbulence and $\beta = 2\Omega/N$ not too close to 1, it was shown in [I], section 3, that the NP component decouples from the wave one at leading order and a modified version of DNS was developed for that component. On the other hand, the wave component was described by wave-turbulence analysis, which has a long and ongoing history (see e.g. Benney and Saffman (1966), Benney and Newell (1969), Zhakarov, Lvov and Falkovich (1992), Newell and Rumpf (2011), Nazarenko (2011) and Galtier (2022)). To close the wave-turbulence equations, it was assumed in [I] that the amplitude of the NP component was small compared with the wave one. This restriction was lifted in [II], which allowed for comparable amplitudes of the two components. In that case, the NP component was found to induce coupling between wave modes, but itself still evolved as if the wave component did not exist.

However, both [I] and [II] assume that $\beta-1$ is not too small, otherwise the wave modes
are only weakly dispersive and wave-turbulence analysis, which assumes dispersive waves,
no longer applies. The case in which $\beta-1$ is small or zero is the subject of this article and
wave-turbulence analysis makes no appearance. Instead, we derive weak-turbulence
equations for the mode amplitudes and examine their consequences, at first analytically, then
using numerical integration.

As in [II], we find that, owing to vanishing of a particular NP-wave coupling coefficient
for $\beta=1$, the NP component is decoupled from the wave one at leading order, though the
latter is affected by the former. We should however note that, being based on leading-order
analysis, the theory describes what happens up to times of $O(\epsilon^{-1})$ times the wave period,
where $\epsilon$ is a small parameter characterising the amplitude of turbulence (the Rossby or
Froude number). $\epsilon^{-1}$ is the time scale for nonlinear evolution, but it may be that, as suggested
by Thomas (2023), NP decoupling does not persist to larger times.

The approach used here is, in the end, numerical and similar to classical DNS, though
weak-turbulence analysis is used to suppress the oscillations due to inertial-gravity waves. As
a result, it avoids the problems of direct application of classical DNS referred to above.

This paper is organised as follows. In section 2, the analytical framework using modal
projection is summarised and weak-turbulence approximations of the mode-amplitude
equations for small $\beta-1$ derived. Decoupling of the NP from the wave component is shown
in appendix A, while energy conservation by nonlinearity is the subject of appendix B.
Finally, appendix C describes the numerical procedure used to integrate the mode-amplitude
equations and results are given in section 3.

## 2. Analytical treatment

Consider homogeneous turbulence in a rotating, stably stratified fluid having constant Brunt-
Vaisala frequency, $N$, and constant rotation vector $\Omega$, which is supposed parallel to gravity
and has unit vector $\mathbf{e} = \mathbf{\Omega}/\mathbf{\Omega}$, where $\mathbf{\Omega} = |\mathbf{\Omega}|$. Henceforth, spatial coordinates, time and velocity are nondimensionalised using $L$, $(N^2 + 4\mathbf{\Omega}^2)^{-1/2}$ and $L(N^2 + 4\mathbf{\Omega}^2)^{1/2}$, where $L$ is a length scale characterising the large scales of the initial turbulence.

Using a rotating frame of reference and Cartesian coordinates, as well as the summation convention, the nondimensional Boussinesq equations of motion are

$$ \frac{\partial u_i}{\partial t} + \hat{\Omega} \varepsilon_{ijk} e_j u_k + \hat{N} \eta e_i + \frac{\partial p}{\partial x_i} = - \frac{\partial}{\partial x_j} (u_i u_j) + D_u \frac{\partial^2 u_i}{\partial x_j \partial x_j}, \quad (2.1) $$

$$ \frac{\partial u_i}{\partial x_i} = 0, \quad (2.2) $$

$$ \frac{\partial \eta}{\partial t} - \hat{N} e_i u_i = - \frac{\partial}{\partial x_j} (\eta u_j) + D_\eta \frac{\partial^2 \eta}{\partial x_j \partial x_j}, \quad (2.3) $$

where $\hat{\Omega} = 2\Omega(N^2 + 4\Omega^2)^{-1/2}$, $\hat{N} = N(N^2 + 4\Omega^2)^{-1/2}$, $u_i$, $p$ and $\eta$ denote velocity, pressure and the buoyancy variable, while $\varepsilon_{ijk}$ is the alternating tensor. Furthermore, $D_u = \nu(N^2 + 4\Omega^2)^{-1/2} L^{-2}$ and $D_\eta = \kappa(N^2 + 4\Omega^2)^{-1/2} L^{-2}$, where $\nu$ is the kinematic viscosity and $\kappa$ the diffusivity associated with the buoyancy variable $\eta$. Note that $\hat{N} = (\beta^2 + 1)^{-1/2}$ and $\hat{\Omega} = \beta(\beta^2 + 1)^{-1/2}$, where $\beta = \hat{\Omega}/\hat{N} = 2\Omega/N$. Thus, the nondimensional governing equations, (2.1)-(2.3), only depend on the parameters $\beta$, $D_u$ and $D_\eta$. For simplicity's sake, $\hat{\Omega}$, $\hat{N}$ are denoted $\Omega$, $N$ in what follows. Because we will only be working with nondimensional quantities, this should not lead to confusion. Note that the right-hand sides of (2.1) and (2.3) express nonlinearity, viscosity and diffusion. The visco-diffusive terms dissipate energy.

As described in [I], in the absence of nonlinearity and visco-diffusion, (2.1)-(2.3) have particular solutions, referred to as modes, which are indexed by an integer $s = 0, \pm 1$ and a wave vector $\mathbf{k}$. Modes have the velocity and buoyancy variable

$$ u_i = v_i^{(s)}(\mathbf{k}) \exp\left(i(\mathbf{k} \cdot \mathbf{x} - s\omega(\mathbf{k})t)\right), \quad (2.4) $$

$$ \eta = \eta^{(s)}(\mathbf{k}) \exp\left(i(\mathbf{k} \cdot \mathbf{x} - s\omega(\mathbf{k})t)\right), \quad (2.5) $$

where

$$ \omega(\mathbf{k}) = (N^2 \sin^2 \theta_k + \Omega^2 \cos^2 \theta_k)^{1/2}, \quad (2.6) $$
$\theta_k$ is the angle between $k$ and the rotation vector and $v_i^{(s)}$, $\eta^{(s)}$ are given by (I.2.12)-(I.2.15) (where, here and henceforth, (I.x.y) refers to equation (x.y) of [I]).

Although modes do not allow for nonlinearity or visco-diffusion, they form a complete set for $u_i(\mathbf{x})$, $\eta(\mathbf{x})$ satisfying (2.2). Thus, the flow can be expressed at any given time as

$$u_i = \sum_{s=0,\pm 1} u_{i,s}, \quad \eta = \sum_{s=0,\pm 1} \eta_s, \qquad (2.7)$$

where

$$u_{i,s} = \int a_s(\mathbf{k}, t) v_i^{(s)}(\mathbf{k}) \exp(i(\mathbf{k} \cdot \mathbf{x} - s\omega(\mathbf{k})t)) d^3\mathbf{k}, \qquad (2.8)$$

$$\eta_s = \int a_s(\mathbf{k}, t) \eta^{(s)}(\mathbf{k}) \exp(i(\mathbf{k} \cdot \mathbf{x} - s\omega(\mathbf{k})t)) d^3\mathbf{k}, \qquad (2.9)$$

and $a_s$ are the mode amplitudes at time $t$. The $a_s$ evolve in time due to nonlinearity and visco-diffusion.

(2.7)-(2.9) express the flow as a linear combination of modes indexed by $s = 0, \pm 1$ and the wave vector $\mathbf{k}$, each mode being a solution of the governing equations without visco-diffusion and nonlinearity of the form $\exp(i(\mathbf{k} \cdot \mathbf{x} - s\omega(\mathbf{k})t))$. The modes with $s = \pm 1$ are inertial-gravity waves, whose dispersion relation is (2.6). The nondimensionalisation of time makes the period of oscillation of such waves of order 1. The mode $s = 0$ has zero frequency, hence zero group velocity. The corresponding component of the flow is non-propagating, hence its designation by NP. As in [I], since $u_i$ and $\eta$ are real, $a_s(-\mathbf{k}) = a_{-s}^*(\mathbf{k})$, where * denotes complex conjugation. Using the identities $v_i^{(s)}(-\mathbf{k}) = v_i^{(-s)*}(\mathbf{k})$, $\eta^{(s)}(-\mathbf{k}) = \eta^{(-s)*}(\mathbf{k})$ and $\omega(-\mathbf{k}) = \omega(\mathbf{k})$, (2.8) and (2.9) imply $u_{i,s} = u_{i,-s}^*$, $\eta_s = \eta_{-s}^*$. Thus, the NP component, defined by $u_{i,0}$ and $\eta_0$, is real, whereas $s = +1$ and $s = -1$ provide conjugate contributions to (2.7), whose sum is the wave component and is also real.

Since we are dealing with turbulent flow, a statistical description is to be expected. In particular averaging plays an important role. The mean flow is supposed zero, hence $\overline{a_s(\mathbf{k})} = 0$, where the overbar denotes an ensemble average. Given statistical homogeneity, the second-order moments have the form

$$\overline{a_s^*(\mathbf{k}) a_s(\mathbf{k}')} = A_{ss'}(\mathbf{k}) \delta(\mathbf{k} - \mathbf{k}'), \qquad (2.10)$$
where $A_{ss'}$ is the spectral matrix, which is Hermitian and positive semi-definite. Given $a_s(-\mathbf{k}) = a_{-s}(\mathbf{k})$, $A_{ss'}(-\mathbf{k}) = A_{-s',-s}(\mathbf{k})$. The average, nondimensional energy density in physical space is the sum of kinetic and potential energies:

$$ \underbrace{\frac{1}{2} \overline{u_i u_i}}_{\text{Kinetic}} + \underbrace{\frac{1}{2} \overline{\eta^2}}_{\text{Potential}}. \qquad (2.11) $$

As shown in [I], (2.11) can be expressed as

$$ \frac{1}{2} (\overline{u_i u_i} + \overline{\eta^2}) = \int e(\mathbf{k}) d^3\mathbf{k}, \qquad (2.12) $$

where $e(\mathbf{k})$ is the energy density in spectral space and is given by

$$ e = \frac{1}{2} \sum_{s=0,\pm 1} A_{ss}. \qquad (2.13) $$

The diagonal elements of the spectral matrix, which appear in (2.13), are real, non-negative and express the energy densities of the different modes in spectral space. Given $A_{ss'}(-\mathbf{k}) = A_{-s',-s}(\mathbf{k})$, they are such that $A_{00}(-\mathbf{k}) = A_{00}(\mathbf{k})$ and $A_{-1,-1}(-\mathbf{k}) = A_{11}(\mathbf{k})$. On the other hand, the off-diagonal elements of $A_{ss'}$ can be complex and represent correlations between modes. The energy densities of the wave and NP components are

$$ e_W = \frac{1}{2} \sum_{s=\pm 1} A_{ss}, \qquad e_{NP} = \frac{1}{2} A_{00}, \qquad (2.14) $$

whose sum provides $e$. The energy densities are such that $e_W(-\mathbf{k}) = e_W(\mathbf{k})$, $e_{NP}(-\mathbf{k}) = e_{NP}(\mathbf{k})$ and $e(-\mathbf{k}) = e(\mathbf{k})$.

The mode amplitudes evolve according to (I.2.26), i.e.

$$ \frac{\partial a_s}{\partial t} = -ik_j \left( v_i^{(s)*} \widetilde{u}_i u_j + \eta^{(s)*} \widetilde{\eta} u_j \right) \exp(is\omega t) - \sum_{\hat{s}=0,\pm 1} D_{ss'} \exp(i(s-\hat{s})\omega t) a_{\hat{s}}, \quad (2.15) $$

where

$$ \tilde{f}(\mathbf{k}) = \frac{1}{8\pi^3} \int f(\mathbf{x}) \exp(-i\mathbf{k}\cdot\mathbf{x}) d^3\mathbf{x} \qquad (2.16) $$

denotes spatial Fourier transformation and $D_{ss'}$ is the Hermitian, positive-definite matrix defined by (I.2.20). The first term on the right-hand side of (2.15) represents nonlinearity, whereas the second expresses visco-diffusive effects, which are dissipative.
Assuming, as we do from here on, that the turbulence is weak (small Rossby or Froude number), the nonlinear term in (2.15) is small and many wave periods (large $t$) are required for it to be effective. The visco-diffusive term must then also be small (i.e. small $D_{\tilde{s}\tilde{s}}$), otherwise dissipation will kill the turbulence before nonlinearity intervenes. Terms in the final sum of (2.15) are then small and oscillatory when $\hat{s} \neq s$ and thus have little effect on the long-time evolution considered here. Neglecting such terms, (2.15) becomes

$$ \frac{\partial a_s}{\partial t} + D_{\tilde{s}\tilde{s}} a_s = -ik_j \left( v_i^{(s)*} \tilde{u}_i \tilde{u}_j + \eta^{(s)*} \tilde{\eta} u_j \right) \exp(is\omega t). \quad (2.17) $$

Since the matrix $D_{\tilde{s}\tilde{s}}$ is Hermitian and positive-definite, the diagonal elements appearing in (2.17) are real and positive. They can be written as $D_{00} = D_{NP}$, $D_{11} = D_{-1,-1} = D_W$, where $D_{NP}$ and $D_W$ are the NP- and wave-mode damping factors. It remains to consider the nonlinear term, i.e. the right-hand side of (2.17).

From here on, we suppose small $\beta - 1$ and derive the leading-order amplitude equations. Given $N = (\beta^2 + 1)^{-1/2}$ and $\Omega = \beta(\beta^2 + 1)^{-1/2}$, $\beta = 1$ leads to $\omega(\mathbf{k}) = 2^{-1/2}$ according to (2.6). This value, independent of $\mathbf{k}$, means that the wave modes are nondispersive when $\beta = 1$, hence wave-turbulence theory does not apply. For small $\beta - 1$, (2.6) implies

$$ \omega(\mathbf{k}) \sim 2^{-1/2} \left( 1 + \frac{1}{2}(\beta - 1)\cos 2\theta_{\mathbf{k}} \right), \quad (2.18) $$

which gives the constant value $\omega = 2^{-1/2}$, noted above, when $\beta = 1$ and shows the small departures from that value when $|\beta - 1|$ is small. Using (2.18) in (2.8) and (2.9) gives

$$ u_{i,s} = w_{i,s} \exp(-ist/2^{1/2}), \quad \eta_s = \zeta_s \exp(-ist/2^{1/2}) \quad (2.19) $$

where

$$ w_{i,s} = \int c_s(\mathbf{k}) v_i^{(s)}(\mathbf{k}) \exp(i\mathbf{k}\cdot\mathbf{x}) d^3\mathbf{k}, \quad (2.20) $$

$$ \zeta_s = \int c_s(\mathbf{k}) \eta^{(s)}(\mathbf{k}) \exp(i\mathbf{k}\cdot\mathbf{x}) d^3\mathbf{k}, \quad (2.21) $$

$$ c_s(\mathbf{k}) = a_s(\mathbf{k}) \exp(-i2^{-3/2} s (\beta - 1) t \cos 2\theta_{\mathbf{k}}). \quad (2.22) $$

According to (2.17), given small nonlinearity and visco-diffusion, $a_s$ is a slowly varying function of time. Combining this with small $\beta - 1$, $c_s$ is also slowly varying, hence so are $w_{i,s}$ and $\zeta_s$ according to (2.20) and (2.21). Note that $w_{i,-s} = w_{i,s}^*$, $\zeta_{-s} = \zeta_s^*$ because $u_{i,s} = u_{i,-s}^*$
and $\eta_s = \eta_{-s}^*$. Note also that small $\beta-1$ allows us to approximate $D_{ss}$, $v_i^{(s)}$ and $\eta^{(s)}$ in (2.17), (2.20) and (2.21) by their values for $\beta=1$ (values obtained using (I.2.12)-(I.2.15) and (I.2.20) with $N=\Omega=2^{-1/2}$). From here on, $D_{ss}$, $v_i^{(s)}$ and $\eta^{(s)}$ should be understood as having these values. The resulting $v_i^{(s)}$ and $\eta^{(s)}$ are given in equations (A.13) and (A.14). This approximation, which is valid at leading order, has important consequences. In particular, it leads to decoupling of the NP component from the wave one. Note that $c_s(-\mathbf{k}) = c_{-s}(\mathbf{k})$ follows from (2.22) and $a_s(-\mathbf{k}) = a_{-s}(\mathbf{k})$.

Using (2.6), (2.18) and (2.19),

$$ -ik_j \left( v_i^{(s)*} \widetilde{u_i u_j} + \eta^{(s)*} \widetilde{\eta u_j} \right) \exp(is\omega t) = -ik_j \exp(i2^{-3/2}s(\beta-1)t \cos 2\theta_k) \\ \sum_{s',s'=0,\pm 1} \exp(i2^{-1/2}(s-s'-s'')t) \left( v_i^{(s)*} \widetilde{w_{i,s'} w_{j,s''}} + \eta^{(s)*} \widetilde{\zeta_{s'} w_{j,s''}} \right) . \quad (2.23) $$

Given that $w_{i,s}$ and $\zeta_s$ are slowly varying, the sum consists of oscillatory ($s-s'-s'' \neq 0$) and non-oscillatory ($s-s'-s'' = 0$) terms. Since nonlinearity is small, it takes a long time to act and we thus focus on large-$t$ evolution. Whereas the oscillatory terms in (2.23) only induce small oscillations of $a_s$ when (2.17) is integrated, the non-oscillatory ones can lead to significant cumulative evolution. Thus, we neglect the oscillatory terms, leaving only the non-oscillatory ones, $s-s'-s'' = 0$, hence (2.17), (2.22) and (2.23) yield

$$ \frac{\partial c_s}{\partial t} + (D_{ss} + i2^{-3/2}s(\beta-1)\cos 2\theta_k)c_s = -ik_j \sum_{s',s'=0,\pm 1} \delta_{s-s'-s'} \left( v_i^{(s)*} \widetilde{w_{i,s'} w_{j,s''}} + \eta^{(s)*} \widetilde{\zeta_{s'} w_{j,s''}} \right), \quad (2.24) $$

where $\delta_0 = 1$ and $\delta_\sigma = 0$ for $\sigma \neq 0$.

When $s=0$, it is shown in appendix A that (2.24) gives

$$ \frac{\partial c_0}{\partial t} + D_{NP}c_0 = -ik_j \left( v_i^{(0)*} \widetilde{w_{i,0} w_{j,0}} + \eta^{(0)*} \widetilde{\zeta_0 w_{j,0}} \right), \quad (2.25) $$

hence, as in [I] and [II], the NP component evolves independently of the wave one. This is due to vanishing of the NP-wave coupling coefficient $M_{0s',-s'}$ ($s'=\pm 1$) for $\beta=1$, where $M_{ss's'}$ is given by (A.11). The resulting decoupling of the NP component from the wave one is one of the main results of this paper, though it is far from obvious a priori. The result indicates that the treatment of the NP component in [I] and [II] continues to apply for small $\beta-1$. It should however be recalled that $D_{ss}$, $v_i^{(s)}$ and $\eta^{(s)}$ here take their values for $\beta=1$. Thus, the NP component in the present work is the $\beta \to 1$ limit of that of [I] and [II]. Note
that decoupling of the NP component has been found to hold numerically by Thomas and Daniel (2021) for particular cases in which $\beta$ is not close to 1 and the wave energy is not too large compared with the NP energy (i.e. the regime covered by [I] and [II]). As far as we are aware, that it also holds for small $\beta-1$ is a new result.

As regards the wave component, (2.24) yields

$$ \frac{\partial c_s}{\partial t} + (D_W + i2^{-3/2} s(\beta-1)\cos 2\theta_k) c_s = -ik_j \left( v_i^{(s)*} \overline{w_{i,0} w_{j,s}} + w_{i,s} w_{j,0} + \eta^{(s)*} \overline{\zeta_0 w_{j,s}} + \zeta_s w_{j,0} \right) \quad s = \pm 1. \quad (2.26) $$

Equations (2.20), (2.21), (2.25) and (2.26) form the evolution equations of $c_s$, which are the basis of this study. Equations (2.20), (2.21) and (2.25) show that the NP component is decoupled from the wave one at leading order, whereas (2.20), (2.21) and (2.26) indicate that the NP component induces coupling between wave modes. These conclusions are similar to those of [II], though no use of wave-turbulence analysis has been made because, as discussed earlier, it does not apply in the present case.

Let $\varepsilon$ be a small, positive parameter measuring the amplitude of $u_i$ and $\eta$. It is then natural to rescale variables to be of order 1 by defining $\hat{u}_i = u_i / \varepsilon$ and $\hat{\eta} = \eta / \varepsilon$. As a result, the rescalings $\hat{c}_s = c_s / \varepsilon$, $\hat{w}_{i,s} = w_{i,s} / \varepsilon$, $\hat{\zeta}_s = \zeta_s / \varepsilon$ and $\hat{t} = \varepsilon t$ are introduced, so (2.20), (2.21), (2.25) and (2.26) become

$$ \hat{w}_{i,s} = \int \hat{c}_s(\mathbf{k}) v_i^{(s)}(\mathbf{k}) \exp(i\mathbf{k}\cdot\mathbf{x}) d^3\mathbf{k}, \quad (2.27) $$

$$ \hat{\zeta}_s = \int \hat{c}_s(\mathbf{k}) \eta^{(s)}(\mathbf{k}) \exp(i\mathbf{k}\cdot\mathbf{x}) d^3\mathbf{k}, \quad (2.28) $$

$$ \frac{\partial \hat{c}_0}{\partial \hat{t}} + \hat{D}_{NP} \hat{c}_0 = -ik_j \left( v_i^{(0)*} \overline{\hat{w}_{i,0} \hat{w}_{j,0}} + \eta^{(0)*} \overline{\hat{\zeta}_0 \hat{w}_{j,0}} \right) \quad (2.29) $$

and

$$ \frac{\partial \hat{c}_s}{\partial \hat{t}} + (\hat{D}_W + i2^{-3/2} s \hat{\beta} \cos 2\theta_k) \hat{c}_s = -ik_j \left( v_i^{(s)*} \overline{\hat{w}_{i,0} \hat{w}_{j,s}} + \hat{w}_{i,s} \hat{w}_{j,0} + \eta^{(s)*} \overline{\hat{\zeta}_0 \hat{w}_{j,s}} + \hat{\zeta}_s \hat{w}_{j,0} \right) \quad s = \pm 1, \quad (2.30) $$

where $\hat{D}_{NP} = D_{NP} / \varepsilon$, $\hat{D}_W = D_W / \varepsilon$ and $\hat{\beta} = (\beta - 1) / \varepsilon$. Equations (2.27)-(2.30) describe the time evolution of the coefficients $\hat{c}_s$ and are independent of $\varepsilon$. Given the definition of $\hat{t}$, it is apparent that the time scale for nonlinear evolution is of order $\varepsilon^{-1}$, a nondimensional scaling which corresponds with the classical large-eddy turnover time ($L/u'$, $u'$ being a measure of
the turbulent velocities) for evolution of turbulence without rotation or stratification. This is short compared with the time scale, $O(\epsilon^{-2})$, which is typical of wave-turbulence theory and characterised the evolution of the wave component in [I] and [II]. According to (2.27)-(2.30) the NP component is independent of $\hat{\beta}$, whereas, for the wave one, $\hat{\beta}=(\beta-1)/\epsilon$ is the parameter which determines closeness or otherwise to $\beta=1$. We expect the wave-turbulence formulation of [II] to apply when $|\hat{\beta}| \to \infty$.

It is shown in appendix B that nonlinearity conserves the total energies of both the NP and wave components, resulting in

$$ \frac{d}{d\hat{t}} \int \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k} = -2 \int \hat{D}_{NP}(\mathbf{k}) \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k}, \qquad (2.31) $$

$$ \frac{d}{d\hat{t}} \int \hat{e}_W(\mathbf{k}) d^3\mathbf{k} = -2 \int \hat{D}_W(\mathbf{k}) \hat{e}_W(\mathbf{k}) d^3\mathbf{k}, \qquad (2.32) $$

where $\hat{e}_{NP} = e_{NP} / \epsilon^2$ and $\hat{e}_W = e_W / \epsilon^2$ are scaled energy densities. These equations indicate that the energy of the NP and wave components each decay due to their individual visco-diffusive dissipation. The sum of (2.31) and (2.32) describes the decay of the total energy.

Of course, it should not be concluded from the above energy conservation principles that nonlinearity plays no role in dissipation. Rather, when large-scale dissipation is small (i.e. small $\hat{D}_{ss}(\mathbf{k})$ for $k = O(1)$, where $k = |\mathbf{k}|$ is the wavenumber), as is usual in turbulent flows, nonlinearity may induce energy transfer from larger to smaller scales, energy which is finally dissipated by linear mechanisms at the smallest scales. This scenario, involving transfer from large to small scales, is often referred to as a forward energy cascade, the expression “inverse cascade” being used to describe energy transfer from small to large scales. In [I], it was found (numerically) that the NP component has only small dissipation, which suggests that, despite nonlinear spectral transfer, there is no forward energy cascade for that component. Since the description of the NP component is here unchanged, the same conclusion applies in the present study. On the other hand, a forward cascade of wave energy was found in [II] and we expect a similar result here.

Note that if one of the two components, NP or wave, is initially absent, it remains so according to (2.27)-(2.30). Thus, pure-NP and pure-wave flows are possible according to the present theory, as indeed they were in [I] and [II]. In the pure-wave case, the nonlinear term in (2.30) is zero, hence there is no nonlinear transfer and the initial $\hat{e}_W$ simply decays under linear dissipation. Nonlinear effects on wave modes require the NP component, which disallows an energy cascade in the pure-wave case.

For simplicity's sake, it is henceforth assumed that, in addition to being homogeneous, the initial turbulence has two further statistical symmetries: axisymmetry and symmetry under reflection in any plane perpendicular to the axis of gravity/rotation. Both are preserved by
evolution according to the governing equations and so apply at all times. These symmetries lead to $A_{ss} = A_{ss}(k, \theta_k)$, $A_{ss}(k, \pi - \theta_k) = A_{ss}(k, \theta_k)$ and $A_{-1,-1} = A_{11} = e_W$.

The initial ($\hat{t} = 0$) conditions are chosen according to

$$ \hat{A}_{00}(\mathbf{k}) = k^2 \exp(-k^2), \qquad (2.33) $$

$$ \hat{A}_{11}(\mathbf{k}) = \hat{A}_{-1,-1}(\mathbf{k}) = \alpha k^2 \exp(-k^2), \qquad (2.34) $$

where $\hat{A}_{ss} = A_{ss} / \varepsilon^2$ and $\alpha \ge 0$ is a parameter which determines the relative importance of the wave and NP components. This choice respects the statistical axisymmetry and reflection symmetry discussed above. For simplicity's sake, the form, $k^2 \exp(-k^2)$, used here is the same for both components. This form is isotropic and places the initial energy firmly in the large scales. Time evolution is expected to produce anisotropy and to transfer energy towards smaller scales, leading to inertial and dissipative ranges, and perhaps an energy cascade. Note that the spherically averaged spectra resulting from (2.33) and (2.34) have the form $k^4 \exp(-k^2)$, a form often used in theoretical studies of turbulence. Integrating (2.33) and (2.34) over $\mathbf{k}$, equation (2.14) leads to

$$ \int \hat{e}_W(\mathbf{k}) d^3\mathbf{k} = \frac{3}{2} \alpha \pi^{3/2}, \quad \int \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k} = \frac{3}{4} \pi^{3/2} \qquad (2.35) $$

for the total initial scaled energies of the two components. Using $\hat{e}_{NP} = e_{NP} / \varepsilon^2$,

$$ \varepsilon = \left( \frac{4}{3\pi^{1/2}} \int e_{NP}(\mathbf{k}) d^3\mathbf{k} \right)^{1/2} \qquad (2.36) $$

gives $\varepsilon$ in terms of the initial NP energy. This fixes the precise value of $\varepsilon$, up to now an order-of-magnitude parameter.

The full problem, which is integrated numerically, consists of (2.27)-(2.30) with initial conditions for $\hat{c}_s$ which respect from (2.33) and (2.34). The problem, and hence the results of the next section, is independent of $\varepsilon$. Given (2.27)-(2.29) and (2.33), the NP component is independent of $\alpha$ and $\hat{\beta}$. Once that component has been determined, (2.27), (2.28) and (2.30) form a linear system for evolution of $\hat{c}_s$ ($s = \pm 1$). (2.34) implies that the initial $\hat{c}_s$ scale like $\alpha^{1/2}$ and it follows from linearity that $\hat{c}_s(\hat{t}) / \alpha^{1/2}$ is independent of $\alpha$, hence the same is true of $\hat{A}_{ss}(\hat{t}) / \alpha$. Thus, the NP calculation has no physical parameters, whereas, when scaled appropriately, the wave one depends only on $\hat{\beta}$.
### 3. Numerical results

As usual in studies of turbulence, dissipation, represented by $\hat{D}$ in (2.29) and (2.30), is taken to be small for the large scales. In fact, it is only included to mop up energy before it can be transferred to higher $k$ than resolved by the numerics, thus avoiding numerical problems such as instabilities or a pile up of energy near a large-$k$ numerical cut-off. As in the DNS of [I], the numerical treatment of (2.27)-(2.30) uses a hyperviscosity in which $\hat{D}_{ss}$ is replaced by $dk^6$, with small $d$ and, for simplicity's sake, $d$ has the same value for the wave and NP components. As usual (see e.g. Haugen and Brandenburg (2004)), this is done with small large-scale dissipation in mind and in the hope of extending the expected inertial range to higher wavenumbers, without significantly affecting the large and inertial-range scales.

Given $\hat{c}_s$ at some instant of time, (2.27) and (2.28) yield $\hat{w}_{i,s}$ and $\hat{\zeta}_s$ in physical space, where the products appearing in (2.29) and (2.30) are taken. Fourier transformation then provides the right-hand sides of (2.29) and (2.30) and hence the time derivative of $\hat{c}_s$. It will be observed that this procedure is similar to classical DNS and indeed the numerical scheme used to integrate (2.27)-(2.30), and described in appendix C, is based on that approach. The flow is approximated as periodic, resulting in discretisation in spectral space and truncated Fourier series in place of the Fourier integrals of (2.27) and (2.28). The $\hat{c}_s$ are initialised using (2.33) and (2.34) and random phases. Time stepping of (2.29) and (2.30) gives the flow at later times, where the various measures of energy described below are calculated. Numerical parameters include the hyperviscous coefficient, $d$. Others appear in appendix C, where the values of all the numerical parameters used to obtain the results are given at the end.

In an ideal world the numerical calculations would be used to obtain ensemble averages. However, this would require very many runs, which, given the computational cost, is not attempted. Instead, as usual for DNS, the results given below represent just a single run for each value of $\hat{\beta}$. The initial conditions, and hence the subsequent flow, vary randomly from run to run (as described in appendix C, the initial $\hat{c}_s$ involve random phases). Thus, random fluctuations about the average are to be expected. The numerical scheme also introduces discretisation in spectral space.

In this section, different measures of wave and NP energy are used to present the numerical results. Each of these quantities illustrate different aspects of the results. One consists of the spectra, $\hat{e}_W$ and $\hat{e}_{NP}$, which represent the energy density in spectral space. Another is the total energies

$$ \hat{E}_W^{tot} = \int \hat{e}_W(\mathbf{k}) d^3\mathbf{k}, \qquad \hat{E}_{NP}^{tot} = \int \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k}, \qquad (3.1) $$

which decay according to (2.31) and (2.32), while a third is the spherically averaged spectra
$$ \hat{E}_W(k) = 2\pi \int_0^\pi k^2 \hat{e}_W(k, \theta_k) \sin \theta_k d\theta_k, \qquad (3.2) $$

$$ \hat{E}_{NP}(k) = 2\pi \int_{0}^{\pi} k^2 \hat{e}_{NP}(k, \theta_k) \sin \theta_k d\theta_k, \qquad (3.3) $$

which describe the distribution of energy over $k = |k|$ and give the total energies of the components when integrated over $k$.

Spectral energy fluxes are defined as follows. The wave and NP energies inside the sphere of radius $k$ and centred on the origin evolve according to

$$ \frac{d}{d\hat{t}} \int_{|k'|<k} \hat{e}_W(\mathbf{k}') d^3\mathbf{k}' = -\Phi_W(k) - 2\int_{|k'|<k} \hat{D}_W(\mathbf{k}') \hat{e}_W(\mathbf{k}') d^3\mathbf{k}', \qquad (3.4) $$

$$ \frac{d}{d\hat{t}} \int_{|k'|<k} \hat{e}_{NP}(\mathbf{k}') d^3\mathbf{k}' = -\Phi_{NP}(k) - 2\int_{|k'|<k} \hat{D}_{NP}(\mathbf{k}') \hat{e}_{NP}(\mathbf{k}') d^3\mathbf{k}', \qquad (3.5) $$

where $\Phi_W$ and $\Phi_{NP}$ are the spectral fluxes due to nonlinearity and the final terms represent dissipation inside the given spectral sphere. It may be noted that, according to (3.4) and (3.5), the fluxes are zero when $k = 0$. Taking the limit $k \to \infty$ of (3.4) and (3.5) and using (2.31) and (2.32), $\Phi_W(\infty) = \Phi_{NP}(\infty) = 0$ represents conservation of energy by nonlinearity.

Beginning with the NP component, there is no dependence on $\alpha$ or $\hat{\beta}$. Figure 1 shows the time evolution of the spherically averaged spectrum. The expected initial transfer of energy from large to small scales in figure 1a leads to the establishment of inertial and dissipative ranges. The result is a clear inertial range in which a power law, $k^{-4}$, appears to apply. During the later time range of figure 1b, the exponent of this law gradually evolves from -4 to -3, the value reported in [I], while the spectral peak moves towards lower $k$, corresponding to growth in size of the large scales and suggesting an inverse energy cascade. $\hat{E}_{NP}^{tot}$ was found to decrease by only about 1% between $\hat{t} = 0$ and $\hat{t} = 20$, a final time much larger than that required for establishment of the dissipative range. This suggests that, as noted earlier and in [I], there is no forward energy cascade for the NP component in the weak turbulence limit. These results illustrate the distinction between nonlinear energy transfer from large to small scales, which creates such scales, as in figure 1a, and a forward cascade which delivers significant energy to the dissipative ones. Note that random variations, resulting from use of a single realisation, are apparent in figure 1 for the lower values of $k$. Taking the spherical average tends to suppress such fluctuations, the more so as $k$ increases because more Fourier components contribute to the average. Discretisation is also evident.
<figure><img src="image_1.png" /></figure>

Figure 1: Log-log plots of the spherically averaged NP spectrum at times from a) $\hat{t}=0$ to $\hat{t}=3$ in steps of 0.15, b) $\hat{t}=3$ to $\hat{t}=19$ in steps of 2. The dashed line represents the power law $k^{-3}$, the finely dashed ones $k^{-4}$. The arrows indicate the direction of increasing time.

Figure 2 shows the NP spectral flux, $\Phi_{NP}(k)$, at different times. As noted earlier, $\Phi_{NP}(0)=\Phi_{NP}(\infty)=0$ follows from the definition of $\Phi_{NP}$ and conservation of total NP energy by nonlinearity. Furthermore, the initial $\Phi_{NP}(k)$ should be zero because the numerical initialisation using random phases (described in appendix C and in keeping with classical DNS) makes the third-order spectral moments zero. As noted in appendix B, and reflecting the usual turbulence closure problem, the third-order moments determine nonlinear energy transfers between different wave vectors. When they are zero, there is no transfer and hence no energy flux. In figure 2a, $\Phi_{NP}(k)$ is small for $\hat{t}=0$, but not exactly zero. This is because a single realisation, rather than an ensemble average, has been used.

As time increases, figure 2a shows the appearance of regions of both positive and negative $\Phi_{NP}$. The former represent transfer from smaller to larger scales, while the latter result in the formation of the inertial and dissipative ranges found earlier. At the later times of figure 2b, the positive peak drops away and essentially disappears. This leaves only small to large scale transfer at large times and again suggests that there is an inverse, rather than a forward, NP energy cascade. That the NP energy cascade is inverse has been recognised for many years in the geophysical community. As noted in [1], (2.25) is equivalent to the three-dimensional quasi-geostrophic equation. The inverse nature of the energy cascade in quasi-geostrophic turbulence was proposed in the seminal work of Charney (1971) and confirmed by later studies (see e.g. Vallgren and Lindborg (2010)). Also note that Charney (1971) suggested the $k^{-3}$ spectrum found above.
<figure><img src="image_2.png" /><figcaption>Figure 2: Plots of the NP energy spectral flux at times from a) $\hat{t}=0$ to $\hat{t}=1$ in steps of 0.1, b) $\hat{t}=1$ to $\hat{t}=12$ in steps of 1. The horizontal axis is logarithmic, while the vertical one is linear. The arrows indicate increasing time.</figcaption></figure>

Turning attention to the wave component, the scaled results given below are only dependent on $\hat{\beta}$. For each of the $\hat{\beta}$ considered here, a calculation was performed for $-\hat{\beta}$. There are differences between the two, but they are sufficiently small that there is no visible difference for the results given here. The closeness of results for $\hat{\beta}$ and $-\hat{\beta}$ suggests that there may be a precise theoretical connection between them, but, if so, we have been unable to find one. In any case, we only give results for $\hat{\beta} \ge 0$, because those for $\hat{\beta} < 0$ are essentially identical.

Figure 3 shows the time evolution of the scaled overall wave energy, $\hat{E}_w^{tot} / \alpha$, for different values of $\hat{\beta}$. Following a phase in which the energy is nearly constant, there is a transition to
decay. This behaviour is quite different from the near constancy of the NP energy and suggests a forward energy cascade for the wave component.

<figure><img src="image_3.png" /><figcaption>Figure 3: Log-log plots of the scaled wave energy as a function of time for $\hat{\beta}=0, 5, 10, 20$ and $50$. The arrow indicates increasing $\hat{\beta}$.</figcaption></figure>

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Log-log plots of the scaled spherically averaged wave spectrum, $\hat{E}_w / \alpha$, for a) $\hat{\beta}=0$, b) $\hat{\beta}=5$, c) $\hat{\beta}=10$, d) $\hat{\beta}=20$ and e) $\hat{\beta}=50$ at times from $\hat{t}=0$ to $\hat{t}=3$ in steps of 0.15 in the left-hand-side figures and from $\hat{t}=3$ to $\hat{t}=19$ in steps of 2 in the right-hand-side ones. The dashed lines represent power laws. Arrows representing the power laws and the direction of time evolution are given.</figcaption></figure>

Figure 4 shows the scaled spherically averaged wave spectrum, $\hat{E}_w / \alpha$, for different $\hat{\beta}$ at the same times as figure 1. As for the NP component, there is an initial transfer of energy
from large to small scales, followed by the establishment of a dissipative range and slower subsequent evolution. There is very little variation with $\hat{\beta}$ during the first phase, whereas the dependence on $\hat{\beta}$ is more pronounced for the second one, shown on the right-hand side of the figures. During the second phase, the wave spectrum decays with time for all $k$, whereas there are ranges of both increase and decrease of the NP spectrum in figure 1 for the same range of $\hat{t}$. This reflects near energy conservation of the NP component and significant dissipation of the wave one.

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

a

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

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

b

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

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

c

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /></figure>

Figure 5: Plots of the scaled wave energy spectral flux, $\Phi_w(k)/\alpha$, for a) $\hat{\beta}=0$, b) $\hat{\beta}=10$, c) $\hat{\beta}=20$ and d) $\hat{\beta}=50$ at times from $\hat{t}=0$ to $\hat{t}=1$ in steps of 0.1 in the left-hand-side figures and from $\hat{t}=1$ to $\hat{t}=12$ in steps of 1 in the right-hand-side ones. The horizontal axis is logarithmic, while the vertical one is linear. The arrows represent the direction of increasing time.

In figures 4a, 4b and 4c, at large times, $\hat{E}_w$ becomes roughly independent of $k$ in the inertial range, which distributes the wave energy across all scales larger than the dissipative ones. This is different from the NP case and also the classical picture of turbulence in the absence of waves, for which, despite the transfer towards smaller ones, the large scales remain the dominant contributors to the energy. Such a strong transfer towards dissipative scales obviously favours dissipation. In figure 4e, the large-time behaviour is more traditional with the majority of the energy in the large scales, while figure 4d represents a case which is somewhat intermediate. Both figures 4d and 4e show an approximate power law in the inertial range. Such power laws arise at sufficiently large $\hat{\beta}$. Once power laws appear, the exponent is increasingly negative as $\hat{\beta}$ increases. We note that figure 2 of [II] indicates an exponent of about -4 for $\hat{e}_w$ when $\beta=0.7$. This corresponds to a spherically averaged exponent of -2, which may be compared with the value -1.5 of figure 4e, which is for the largest value of $\hat{\beta}$ considered here. Given an increasingly negative spectral exponent for increasing $\hat{\beta}$, this is in accord with the expectation, noted earlier, that the present results approach those of [II] at large $|\hat{\beta}|$.

Figure 5 shows the scaled wave energy flux, $\Phi_w(k)/\alpha$. As for the NP case, $\Phi_w(0)=0$, $\Phi_w(\infty)=0$ and $\Phi_w$ is initially zero. It will be observed that the different values of $\hat{\beta}$ have similar results. Compared to the NP case, a striking feature is that the wave flux is always positive, i.e. there is energy transfer from large to small scales. Furthermore, at large times, the wave energy flux extends through the inertial range before being absorbed in the dissipative range. These results again suggest a forward energy cascade of the wave component.
<figure><img src="image_13.png" /><figcaption>Figure 6: Contour plots of $\hat{e}_{NP}$ for a) $\hat{t}=0$ and b) $\hat{t}=20$. There are ten contour values, which are logarithmically spaced from $10^{-4}$ to 1.</figcaption></figure>

Finally, we give some results concerning anisotropy. Figure 6 shows contour plots of $\hat{e}_{NP}$ in the $k_{\perp}$, $k_{\|}$ plane, where $k_{\perp} = k \sin \theta_k$ and $k_{\|} = k \cos \theta_k$ are wavenumbers perpendicular and parallel to the rotation/vertical axis. Figure 6a is for the initial time and the circular contours reflect the isotropy of the initial conditions. Figure 6b is for the largest time, $\hat{t} = 20$, at which the calculation was performed. The wiggles reflect the fact that these results represent a single realisation, rather than being an ensemble average. Anisotropy favouring the axial direction is apparent. This suggests that NP energy transfer is more effective in that direction than in the transverse direction. These results represent the case $\beta=1$, the value appropriate to the present study. Anisotropy of the NP component for other values of $\beta$ was discussed in section 4.1 of [I]. It was concluded that the NP inertial-range energy distribution is close to isotropic for $\beta=1.4$ and favours the axial direction for lower values of $\beta$ (such as the $\beta=1$ considered here), and the transverse direction for higher values.

Concerning the wave component, figure 7 shows contour plots of $\hat{e}_w / \alpha$ for $\hat{t}=20$ and different values of $\hat{\beta}$. As for the NP component, anisotropy has a preference for the axial direction. Thus, the total energy density, which is the sum of the NP and wave components, also favours that direction. It further appears from figures 7b-d that the energy density decreases near the plane $k_{\|} = 0$. This tendency increases with $\hat{\beta}$ and reflects the conclusion of [II] that wave dissipation is greatest near the given plane. That this effect strengthens with increasing $\hat{\beta}$ is to be expected since $|\hat{\beta}| \to \infty$ should approach the case considered in [II].

## 4. Conclusions

In this article we have considered a particular case, namely small $\beta-1$, of weak turbulence (small Rossby or Froude number, the two being equivalent when $\beta-1$ is small) with rotation and stratification. This case is interesting because it does not permit wave-turbulence
analysis, the wave modes being only weakly dispersive. Furthermore, it completes the studies of weak, rotating, stratified turbulence initiated in [I] and carried further in [II].

<figure><img src="image_14.png" /><figcaption>Figure 7: Contour plots of $\hat{e}_w / \alpha$ for $\hat{t} = 20$ and a) $\hat{\beta} = 0$, b) $\hat{\beta} = 10$, c) $\hat{\beta} = 20$ and d) $\hat{\beta} = 50$. There are ten contour values, which are logarithmically spaced from $10^{-4}$ to 1.</figcaption></figure>

The flow is expressed as a combination of wave and NP components using modal analysis, as described in section 2, where the weak-turbulence equations for the modal amplitudes were derived. Owing to vanishing of the NP-wave coupling coefficient $M_{0s',-s'}(\mathbf{k}, \mathbf{p})$ for $\beta = 1$ and all $\mathbf{k}$ and $\mathbf{p}$, the NP component is found, as in [I] and [II], to evolve independently of the wave one and according to an amplitude equation which is precisely that of [I] and [II] in the limit $\beta \to 1$. As in [II], the NP component induces coupling between wave modes, though such coupling is no longer constrained by the resonance condition of the earlier study, a condition which resulted from wave-turbulence theory. In the absence of the NP component,
nonlinearity has no effect on the wave component, which evolves according to linear
dissipation alone.

The parameter $\hat{\beta} = (\beta - 1) / \varepsilon$ appears in the amplitude equations of the wave component, where $\varepsilon$ is a small parameter (Rossby or Froude number) which measures the weakness of turbulence. $\hat{\beta}$ expresses the departure from $\beta = 1$, scaled in a manner appropriate to the present work. When $|\hat{\beta}|$ is sufficiently large, the analysis of [II] should apply, otherwise the present theory is needed. Note that the NP component is independent of $\hat{\beta}$.

The time variable appearing in the amplitude evolution equations is $\hat{t} = \varepsilon t$. This implies a nonlinear evolution time $O(\varepsilon^{-1})$ times the inertial-gravity wave period, which is comparable with the classical large-eddy turnover time for evolution of turbulence without rotation or stratification. This is short compared with the $O(\varepsilon^{-2})$ periods typical of wave-turbulence theory and which characterised the wave evolution time in [I] and [II]. Thus, here the wave time scale is $O(\varepsilon^{-1})$, whereas it was $O(\varepsilon^{-2})$ before. This suggests that, as $\beta - 1$ decreases from the $O(1)$ values of [II], wave-component evolution accelerates until the associated time reaches the $O(\varepsilon^{-1})$ value of the present work when $\beta - 1 = O(\varepsilon)$.

Energy conservation by nonlinearity was shown to hold for both components individually and hence for their sum, which is the total energy of the flow. Thus, nonlinearity can transfer energy between modes, but not create or destroy it overall. As usual, when dissipation is small at the large scales, which is the usual situation in turbulent flow, nonlinear transfer from large to small scales may lead to a significant energy supply to the smallest ones, which are dissipative. We refer to such a situation as a forward energy cascade. There are however examples, such as the NP component considered here, in which nonlinear transfer from small to large scales occurs, but only small energy dissipation results because the energy cascade is towards the large scales (inverse cascade). On the other hand, we found that the dissipation of the wave component was significant, suggesting a forward cascade for that component.

A numerical scheme for solving the mode-amplitude equations was developed. This scheme is similar to classical, periodic DNS. The initial conditions ((2.33) and (2.34)) were chosen such that the energy distribution is initially isotropic and concentrated in the large scales. Not surprisingly, it was found that transfer to smaller scales lead to the development of inertial and dissipative ranges. This was followed by slower evolution, with rather different consequences for the NP and wave modes.

In the NP case, the result was small energy dissipation and a clear power law for the spherically averaged spectrum in the inertial range. Nonlinear energy transfer from large to small scales plays an important role in the establishment of the inertial range, but the NP component appears to have an inverse, rather than a forward, energy cascade. On the other hand, energy transfer to small scales was much stronger in the wave case, leading to
significant dissipation. We conclude that the wave component is subject to a forward energy cascade. At the same time, it should be recalled that this depends on the existence of the NP component, without which coupling of wave modes does not occur. If the NP amplitude is reduced, nonlinear evolution of the wave component takes longer and longer.

In summary, in both [II] and the present work, the NP component acts as a catalyst for nonlinear coupling of wave modes, leading to a forward energy cascade, hence dissipation, of the wave component. Here we use the term catalyst by analogy with its usage in chemistry because the NP component affects the waves, but is not affected in turn. Assuming small dissipation of the large scales and beginning with both components, the wave component is expected to decay, leaving the NP one. Of course, external forcing of the wave component could maintain it.

Numerical results concerning anisotropy were also obtained. The energy distributions of both the NP and wave components show a preference for the axial direction, hence a trend in that direction for the total flow. As $\hat{\beta}$ was increased from zero, a dip in the intensity of the wave component near the plane perpendicular to the gravity/rotation axis appeared, and intensified at higher $\hat{\beta}$. Such behaviour was perhaps to be expected given the conclusion of [II] that wave dissipation is maximal for the given plane.

Finally, it would be interesting to understand why the results for $-\hat{\beta}$ are nearly the same as those for $\hat{\beta}$.

**Acknowledgements.** Thanks to Anne Cadiou for help with the computational aspects of this work and to Claude Cambon for helpful discussions and suggestions concerning the manuscript.

**Declaration of interests.** The author reports no conflict of interest.

## Appendix A: Derivation of (2.25)

(2.20) gives

$$w_{i,s}^* = \int c_s^*(\mathbf{p}) v_i^{(s)*}(\mathbf{p}) \exp(-i\mathbf{p}\cdot\mathbf{x}) d^3\mathbf{p}. \qquad (A.1)$$

Thus,

$$ 
\begin{aligned}
w_{i,s}^*, w_{j,s}^* &= \iint c_{s^*}^*(\mathbf{p}) c_{s^*}^*(\mathbf{q}) v_i^{(s)*}(\mathbf{p}) v_j^{(s)*}(\mathbf{q}) \exp(-i(\mathbf{p}+\mathbf{q})\cdot\mathbf{x}) d^3\mathbf{q}d^3\mathbf{p} \\
&= \iint c_{s^*}^*(\mathbf{p}) c_{s^*}^*(-\mathbf{k}-\mathbf{p}) v_i^{(s)*}(\mathbf{p}) v_j^{(s)*}(-\mathbf{k}-\mathbf{p}) \exp(i\mathbf{k}\cdot\mathbf{x}) d^3\mathbf{k}d^3\mathbf{p}. \quad (A.2) \\
&= \iint c_{s^*}^*(\mathbf{p}) c_{s^*}^*(-\mathbf{k}-\mathbf{p}) v_i^{(s)*}(\mathbf{p}) v_j^{(s)*}(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p} \exp(i\mathbf{k}\cdot\mathbf{x}) d^3\mathbf{k}
\end{aligned}
$$

It follows that
$$ \widetilde{w}_{i,s'}^{*} \widetilde{w}_{j,s}^{*} = \int c_{s'}^{*}(\mathbf{p}) c_{s}^{*}(-\mathbf{k}-\mathbf{p}) v_{i}^{(s)}(\mathbf{p}) v_{j}^{(s)}(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p}. \quad (A.3) $$

Likewise,

$$ \widetilde{\zeta}_{s',s}^{*} \widetilde{w}_{j,s}^{*} = \int c_{s'}^{*}(\mathbf{p}) c_{s}^{*}(-\mathbf{k}-\mathbf{p}) \eta^{(s)}(\mathbf{p}) v_{j}^{(s)}(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p}. \quad (A.4) $$

Using $w_{i,-s} = w_{i,s}^*$ and $\zeta_{-s} = \zeta_s^*$,

$$ \sum_{s',s'=0,\pm 1} \delta_{s-s'=s^*} \left( v_i^{(s)} \widetilde{w}_{i,s'} \widetilde{w}_{j,s^*} + \eta^{(s)} \widetilde{\zeta}_{s',s} \widetilde{w}_{j,s^*} \right) = \\ \sum_{s',s'=0,\pm 1} \delta_{s+s'=s^*} \left( v_i^{(s)} \widetilde{w}_{i,s'}^* \widetilde{w}_{j,s^*} + \eta^{(s)} \widetilde{\zeta}_{s',s}^* \widetilde{w}_{j,s^*} \right). \quad (A.5) $$

Employing (A.3)-(A.5), (2.24) gives

$$ \frac{\partial c_s(\mathbf{k})}{\partial t} + (D_{ss}(\mathbf{k}) + i2^{-3/2}s(\beta-1)\cos 2\theta_k) c_s(\mathbf{k}) = \sum_{s',s'=0,\pm 1} \delta_{s+s'=s^*} \int N_{ss's^*}^*(\mathbf{k},\mathbf{p}) c_{s'}^*(\mathbf{p}) c_{s^*}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p}, \quad (A.6) $$

where

$$ N_{ss's^*}(\mathbf{k}, \mathbf{p}) = ik_j v_j^{(s^*)}(-\mathbf{k}-\mathbf{p}) (v_i^{(s)}(\mathbf{k})v_i^{(s)}(\mathbf{p}) + \eta^{(s)}(\mathbf{k})\eta^{(s)}(\mathbf{p})) \quad (A.7) $$

represents nonlinear coupling between modes. The change of integration variable $\mathbf{q} = -\mathbf{k} - \mathbf{p}$ gives

$$ \int N_{ss's^*}^*(\mathbf{k}, \mathbf{p}) c_{s'}^*(\mathbf{p}) c_{s^*}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p} = \\ \int N_{ss's^*}^*(\mathbf{k}, -\mathbf{k}-\mathbf{q}) c_{s'}^*(-\mathbf{k}-\mathbf{q}) c_{s^*}^*(\mathbf{q}) d^3\mathbf{q} = \\ \int N_{ss's^*}^*(\mathbf{k}, -\mathbf{k}-\mathbf{p}) c_{s'}^*(-\mathbf{k}-\mathbf{p}) c_{s^*}^*(\mathbf{p}) d^3\mathbf{p} \quad (A.8) $$

Thus, permuting the summation indices according to $s' \leftrightarrow s''$,

$$ \sum_{s',s'=0,\pm 1} \delta_{s+s'=s^*} \int N_{ss's^*}^*(\mathbf{k}, \mathbf{p}) c_{s'}^*(\mathbf{p}) c_{s^*}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p} = \\ \sum_{s',s'=0,\pm 1} \delta_{s+s'=s^*} \int N_{ss's^*}^*(\mathbf{k}, -\mathbf{k}-\mathbf{p}) c_{s'}^*(\mathbf{p}) c_{s^*}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p}. \quad (A.9) $$

It follows that
$$ 
\begin{aligned}
& \sum_{s', s''=0, \pm 1} \delta_{s+s'+s''} \int N_{ss's''}^*(\mathbf{k}, \mathbf{p}) c_{s'}^*(\mathbf{p}) c_{s''}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p} = \\
& \qquad \sum_{s', s''=0, \pm 1} \delta_{s+s'+s''} \int M_{ss's''}^*(\mathbf{k}, \mathbf{p}) c_{s'}^*(\mathbf{p}) c_{s''}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p},
\end{aligned}
\quad (A.10) $$

where

$$ M_{ss's''} (\mathbf{k}, \mathbf{p}) = \frac{1}{2} (N_{ss's''} (\mathbf{k}, \mathbf{p}) + N_{ss's''} (\mathbf{k}, -\mathbf{k}-\mathbf{p})). \quad (A.11) $$

Using (A.10) in (A.6),

$$ 
\begin{aligned}
& \frac{\partial c_s(\mathbf{k})}{\partial t} + (D_{ss}(\mathbf{k}) + i2^{-3/2}s(\beta-1)\cos 2\theta_k) c_s(\mathbf{k}) = \\
& \qquad \sum_{s', s''=0, \pm 1} \delta_{s+s'+s''} \int M_{ss's''}^*(\mathbf{k}, \mathbf{p}) c_{s'}^*(\mathbf{p}) c_{s''}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p}.
\end{aligned}
\quad (A.12) $$

(A.12) is an evolution equation for the mode amplitudes which is equivalent to (2.20), (2.21) and (2.24).

Adopting Cartesian coordinates with $x_3$ directed along the rotation vector and given that $v_i^{(s)}$ and $\eta^{(s)}$ take their values for $\beta=1$, (I.2.12)-(I.2.15), $N = \Omega = 2^{-1/2}$ and $\mathbf{e}=(0,0,1)$ yield

$$ \begin{pmatrix} v_1^{(0)} \\ v_2^{(0)} \\ v_3^{(0)} \end{pmatrix} = \frac{i}{k} \begin{pmatrix} k_2 \\ -k_1 \\ 0 \end{pmatrix}, \qquad \eta^{(0)} = \frac{ik_3}{k} \quad (A.13) $$

and

$$ \begin{pmatrix} v_1^{(s)} \\ v_2^{(s)} \\ v_3^{(s)} \end{pmatrix} = \frac{1}{2^{1/2}kk_{\perp}} \begin{pmatrix} (isk_1 - k_2)k_3 \\ (isk_2 + k_1)k_3 \\ -isk_{\perp}^2 \end{pmatrix}, \qquad \eta^{(s)} = \frac{k_{\perp}}{2^{1/2}k} \quad (A.14) $$

when $s = \pm 1$, where $k = |\mathbf{k}|$ and $k_{\perp} = (k_1^2 + k_2^2)^{1/2} = k \sin \theta_k$. The coefficients, $M_{ss's''}$, in (A.12) follow from (A.7), (A.11), (A.13) and (A.14).

Taking $s=0$, (A.12) yields

$$ 
\begin{aligned}
& \frac{\partial c_0}{\partial t} + D_{NP}c_0 = \int M_{000}^*(\mathbf{k}, \mathbf{p}) c_0^*(\mathbf{p}) c_0^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p} + \\
& \qquad \sum_{s'=\pm 1} M_{0s's''}^*(\mathbf{k}, \mathbf{p}) c_{s'}^*(\mathbf{p}) c_{-s'}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p}
\end{aligned}
\quad (A.15) $$
Using (A.13) and (A.14),

$$v_i^{(0)}(\mathbf{k})v_i^{(s')}(\mathbf{p}) + \eta^{(0)}(\mathbf{k})\eta^{(s')}(\mathbf{p}) = k^{-1}Z_{s'}(\mathbf{k}, \mathbf{p}), \quad (A.16)$$

$$k_j v_j^{(-s')}(\mathbf{-k} - \mathbf{p}) = s' Z_{-s'}(\mathbf{k}, \mathbf{-k} - \mathbf{p}), \quad (A.17)$$

where $s' = \pm 1$,

$$Z_{s'}(\mathbf{k}, \mathbf{p}) = \frac{i}{2^{1/2} pp_{\perp}} (k_2 p_3 (is'p_1 - p_2) - k_1 p_3 (is'p_2 + p_1) + k_3 p_{\perp}^2), \quad (A.18)$$

$p = |\mathbf{p}|$ and $p_{\perp} = (p_1^2 + p_2^2)^{1/2}$. (A.7), (A.16) and (A.17) give

$$N_{0s',-s'}(\mathbf{k}, \mathbf{p}) = is'k^{-1}Z_{s'}(\mathbf{k}, \mathbf{p})Z_{-s'}(\mathbf{k}, \mathbf{-k} - \mathbf{p}), \quad (A.19)$$

hence $M_{0s',-s'}(\mathbf{k}, \mathbf{p}) = 0$ from (A.11). Thus, the second term on the right-hand side of (A.15) is zero, while (A.11) implies

$$\begin{aligned} \int M_{000}^*(\mathbf{k}, \mathbf{p}) c_0^*(\mathbf{p}) c_0^*(-\mathbf{k} - \mathbf{p}) d^3\mathbf{p} = & \\ \frac{1}{2} \left( \int N_{000}^*(\mathbf{k}, \mathbf{p}) c_0^*(\mathbf{p}) c_0^*(-\mathbf{k} - \mathbf{p}) d^3\mathbf{p} + \right. & \\ \left. \int N_{000}^*(\mathbf{k}, \mathbf{-k} - \mathbf{p}) c_0^*(\mathbf{p}) c_0^*(-\mathbf{k} - \mathbf{p}) d^3\mathbf{p} \right) & \end{aligned} \quad (A.20)$$

Using (A.8), (A.20) becomes

$$\begin{aligned} \int M_{000}^*(\mathbf{k}, \mathbf{p}) c_0^*(\mathbf{p}) c_0^*(-\mathbf{k} - \mathbf{p}) d^3\mathbf{p} = & \\ \int N_{000}^*(\mathbf{k}, \mathbf{p}) c_0^*(\mathbf{p}) c_0^*(-\mathbf{k} - \mathbf{p}) d^3\mathbf{p}, & \end{aligned} \quad (A.21)$$

hence

$$\begin{aligned} \int M_{000}^*(\mathbf{k}, \mathbf{p}) c_0^*(\mathbf{p}) c_0^*(-\mathbf{k} - \mathbf{p}) d^3\mathbf{p} = & \\ -ik_j \left( v_i^{(0)*} \overline{w_{i,0} w_{j,0}} + \eta^{(0)*} \overline{\zeta_0 w_{j,0}} \right) & \end{aligned} \quad (A.22)$$

according to (A.3), (A.4) and (A.7). Thus, (A.15) gives (2.25).

## Appendix B: Energy conservation by nonlinearity

The scaled version of (A.6) is
$$ \frac{\partial \hat{c}_s(\mathbf{k})}{\partial \hat{t}} + (\hat{D}_{ss}(\mathbf{k}) + i2^{-3/2} s \hat{\beta} \cos 2\theta_k) \hat{c}_s(\mathbf{k}) = \\ \sum_{s', s''=0, \pm 1} \delta_{s+s'+s''} \int N_{ss's''}^*(\mathbf{k}, \mathbf{p}) \hat{c}_{s'}^*(\mathbf{p}) \hat{c}_{s''}^*(-\mathbf{k}-\mathbf{p}) d^3\mathbf{p} \qquad (B.1) $$

It follows from (2.10) and (2.22) that

$$ \overline{\hat{c}_s^*(\mathbf{k}) \hat{c}_s(\mathbf{k}')} = \hat{A}_{ss}(\mathbf{k}) \delta(\mathbf{k} - \mathbf{k}'). \qquad (B.2) $$

Given statistical homogeneity, the third-order spectral moments have the form

$$ \overline{\hat{c}_s(\mathbf{k}) \hat{c}_{s'}(\mathbf{p}) \hat{c}_{s''}(\mathbf{k}'')} = \Theta_{ss's''}(\mathbf{k}, \mathbf{p}) \delta(\mathbf{k} + \mathbf{p} + \mathbf{k}''). \qquad (B.3) $$

Multiplying the complex conjugate of (B.1) by $\hat{c}_s(\mathbf{k}')$ and taking the ensemble average,

$$ \overline{\frac{\partial \hat{c}_s^*(\mathbf{k})}{\partial \hat{t}} \hat{c}_s(\mathbf{k}')} = \delta(\mathbf{k} - \mathbf{k}') (\tau_s(\mathbf{k}) - (\hat{D}_{ss}(\mathbf{k}) - i2^{-3/2} s \hat{\beta} \cos 2\theta_k) \hat{A}_{ss}(\mathbf{k})), \qquad (B.4) $$

where

$$ \tau_s(\mathbf{k}) = \sum_{s', s''=0, \pm 1} \delta_{s+s'+s''} \int N_{ss's''}(\mathbf{k}, \mathbf{p}) \Theta_{ss's''}(\mathbf{k}, \mathbf{p}) d^3\mathbf{p} \qquad (B.5) $$

and we have used (B.2), (B.3) and the fact that $\hat{D}_{ss}$ is real. Given the Dirac function in (B.4), we have also replaced $\Theta_{ss's''}(\mathbf{k}', \mathbf{p})$ by $\Theta_{ss's''}(\mathbf{k}, \mathbf{p})$ in (B.5). Taking the conjugate of (B.4) and permuting $\mathbf{k} \leftrightarrow \mathbf{k}'$,

$$ \overline{\hat{c}_s^*(\mathbf{k}) \frac{\partial \hat{c}_s(\mathbf{k}')}{\partial \hat{t}}} = \delta(\mathbf{k} - \mathbf{k}') (\tau_s^*(\mathbf{k}') - (\hat{D}_{ss}(\mathbf{k}') + i2^{-3/2} s \hat{\beta} \cos 2\theta_k) \hat{A}_{ss}(\mathbf{k}')) , \qquad (B.6) \\ = \delta(\mathbf{k} - \mathbf{k}') (\tau_s^*(\mathbf{k}) - (\hat{D}_{ss}(\mathbf{k}) + i2^{-3/2} s \hat{\beta} \cos 2\theta_k) \hat{A}_{ss}(\mathbf{k})) $$

where we have used the fact that $\hat{D}_{ss}$ and $\hat{A}_{ss}$ are real. Taking the sum of (B.4) and (B.6), the time derivative of (B.2) implies

$$ \frac{\partial \hat{A}_{ss}}{\partial \hat{t}} = 2(\Re\{\tau_s(\mathbf{k})\} - \hat{D}_{ss}(\mathbf{k})\hat{A}_{ss}(\mathbf{k})). \qquad (B.7) $$

This result shows that the second-order spectral moments $\hat{A}_{ss}$ evolve due to nonlinearity and dissipation. The nonlinear term involves the third-order moments via (B.5), which reflects the usual closure problem of turbulence.

The energy density in spectral space is given by (2.13), hence
$$ \frac{d}{dt} \int \hat{e}(\mathbf{k}) d^3 \mathbf{k} = \Re \left\{ \sum_{s=0, \pm 1} \int \tau_s(\mathbf{k}) d^3 \mathbf{k} \right\} - \sum_{s=0, \pm 1} \int \hat{D}_{ss}(\mathbf{k}) \hat{A}_{ss}(\mathbf{k}) d^3 \mathbf{k} \quad (B.8) $$

for the evolution of the total energy. Using (B.5),

$$ \sum_{s=0, \pm 1} \int \tau_s(\mathbf{k}) d^3 \mathbf{k} = \sum_{s, s', s''=0, \pm 1} \delta_{s+s'+s''} \iint N_{ss's'}(\mathbf{k}, \mathbf{p}) \Theta_{ss's'}(\mathbf{k}, \mathbf{p}) d^3 \mathbf{p} d^3 \mathbf{k}. \quad (B.9) $$

The identity $k_j v_j^{(s)}(\mathbf{k}) = 0$ follows from the definition, (I.2.12) and (I.2.13), of $v_j^{(s)}(\mathbf{k})$. Thus, $(k_j + p_j) v_j^{(s*)}(-\mathbf{k} - \mathbf{p}) = 0$, hence $N_{s'ss'}(\mathbf{p}, \mathbf{k}) = -N_{ss's'}(\mathbf{k}, \mathbf{p})$ from (A.7). On the other hand, $\Theta_{s'ss'}(\mathbf{p}, \mathbf{k}) = \Theta_{ss's'}(\mathbf{k}, \mathbf{p})$ according to (B.3). It follows that the integrand in (B.9) changes sign under the permutation $s \leftrightarrow s'$, $\mathbf{k} \leftrightarrow \mathbf{p}$. Thus,

$$ \sum_{s=0, \pm 1} \int \tau_s(\mathbf{k}) d^3 \mathbf{k} = 0, \quad (B.10) $$

hence nonlinearity conserves the total energy, which evolves according to

$$ \frac{d}{dt} \int \hat{e}(\mathbf{k}) d^3 \mathbf{k} = - \sum_{s=0, \pm 1} \int \hat{D}_{ss}(\mathbf{k}) \hat{A}_{ss}(\mathbf{k}) d^3 \mathbf{k}. \quad (B.11) $$

Next consider the NP component of the energy. As was shown in appendix A, the second term on the right-hand side of (A.15) is zero, while the first term can be expressed using (A.21). Thus,

$$ \frac{\partial \hat{c}_0}{\partial t} + \hat{D}_{NP} \hat{c}_0 = \int N_{000}^*(\mathbf{k}, \mathbf{p}) \hat{c}_0^*(\mathbf{p}) \hat{c}_0^*(-\mathbf{k} - \mathbf{p}) d^3 \mathbf{p}. \quad (B.12) $$

Multiplying the complex conjugate of (B.12) by $\hat{c}_0(\mathbf{k}')$ and taking the ensemble average,

$$ \overline{\frac{\partial \hat{c}_0^*(\mathbf{k})}{\partial t} \hat{c}_0(\mathbf{k}')} = \delta(\mathbf{k} - \mathbf{k}') (\tau(\mathbf{k}) - \hat{D}_{NP}(\mathbf{k}) \hat{A}_{00}(\mathbf{k})), \quad (B.13) $$

where

$$ \tau(\mathbf{k}) = \int N_{000}(\mathbf{k}, \mathbf{p}) \Theta_{000}(\mathbf{k}, \mathbf{p}) d^3 \mathbf{p} \quad (B.14) $$

and we have used (B.2), (B.3) and the fact that $\hat{D}_{NP}$ is real. Given the Dirac function in (B.13), we have also replaced $\Theta_{000}(\mathbf{k}', \mathbf{p})$ by $\Theta_{000}(\mathbf{k}, \mathbf{p})$ in (B.14). Taking the conjugate of (B.13) and permuting $\mathbf{k} \leftrightarrow \mathbf{k}'$,
$$ \begin{aligned} \overline{\hat{c}_0^*(\mathbf{k}) \frac{\partial \hat{c}_0(\mathbf{k}')}{\partial \hat{t}}} &= \delta(\mathbf{k}-\mathbf{k}') (\tau^*(\mathbf{k}') - \hat{D}_{NP}(\mathbf{k}') \hat{A}_{00}(\mathbf{k}')) , \\ &= \delta(\mathbf{k}-\mathbf{k}') (\tau^*(\mathbf{k}) - \hat{D}_{NP}(\mathbf{k}) \hat{A}_{00}(\mathbf{k})) \end{aligned} \qquad (B.15) $$

where we have used the fact that $\hat{D}_{NP}$ and $\hat{A}_{00}$ are real. Taking the sum of (B.13) and (B.15), (B.2) implies

$$ \frac{\partial \hat{A}_{00}}{\partial \hat{t}} = 2(\Re\{\tau(\mathbf{k})\} - \hat{D}_{NP}(\mathbf{k})\hat{A}_{00}(\mathbf{k})). \qquad (B.16) $$

Given the second of equations (2.14),

$$ \frac{d}{d\hat{t}} \int \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k} = \Re\left\{\int \tau(\mathbf{k}) d^3\mathbf{k}\right\} - 2 \int \hat{D}_{NP}(\mathbf{k}) \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k} \qquad (B.17) $$

for the evolution of the NP energy. Using (B.14),

$$ \int \tau(\mathbf{k}) d^3\mathbf{k} = \iint N_{000}(\mathbf{k}, \mathbf{p}) \Theta_{000}(\mathbf{k}, \mathbf{p}) d^3\mathbf{p} d^3\mathbf{k} . \qquad (B.18) $$

It was shown above that $N_{s's's'}(\mathbf{p}, \mathbf{k}) = -N_{s's's'}(\mathbf{k}, \mathbf{p})$ and $\Theta_{s's's'}(\mathbf{p}, \mathbf{k}) = \Theta_{s's's'}(\mathbf{k}, \mathbf{p})$, hence the integrand in (B.18) changes sign under $\mathbf{k} \leftrightarrow \mathbf{p}$. Thus, (B.18) is zero and (B.17) gives

$$ \frac{d}{d\hat{t}} \int \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k} = -2 \int \hat{D}_{NP}(\mathbf{k}) \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k}, \qquad (B.19) $$

i.e. (2.31).

Since $\hat{D}_{00} = \hat{D}_{NP}$ and $\hat{D}_{11} = \hat{D}_{-1,-1} = \hat{D}_{W}$, (2.14) gives

$$ \begin{aligned} \sum_{s=0,\pm 1} \int \hat{D}_{ss}(\mathbf{k}) \hat{A}_{ss}(\mathbf{k}) d^3\mathbf{k} = & \\ 2 \int \hat{D}_{NP}(\mathbf{k}) \hat{e}_{NP}(\mathbf{k}) d^3\mathbf{k} + 2 \int \hat{D}_{W}(\mathbf{k}) \hat{e}_{W}(\mathbf{k}) d^3\mathbf{k} & . \end{aligned} \qquad (B.20) $$

Using (B.11), (B.19), (B.20) and $\hat{e} = \hat{e}_{NP} + \hat{e}_{W}$ yields (2.32).

## Appendix C: The numerical scheme

In what follows, we use Cartesian coordinates, $x_1, x_2, x_3$, where $x_3$ is directed along the rotation vector and the flow is approximated as periodic with respect to $x_1, x_2$ and $x_3$. Thus, (2.27) and (2.28) are replaced by the Fourier series

$$ \hat{w}_{i,s} = \sum_{\mathbf{k}} \hat{c}_s(\mathbf{k}) v_i^{(s)}(\mathbf{k}) \exp(i\mathbf{k}\cdot\mathbf{x}), \qquad (C.1) $$
$$ \hat{\zeta}_s = \sum_k \hat{c}_s(\mathbf{k}) \eta^{(s)}(\mathbf{k}) \exp(i\mathbf{k} \cdot \mathbf{x}), \qquad (C.2) $$

in which the sums are over **k** having components $k_1 = 2\pi l_1 / L_1$, $k_2 = 2\pi l_2 / L_2$, $k_3 = 2\pi l_3 / L_3$, where the $l_i$ take integer values and $L_i$ are the spatial periods. As usual in DNS, the spatial periods, $L_i$, need to be sufficiently large that periodicity is unimportant. Another way of putting this is that the spacing, $2\pi / L_i$, in spectral space should be sufficiently small. The numerical equivalents of (2.29) and (2.30) are

$$ \frac{d\hat{c}_0}{d\hat{t}} + \hat{D}_{NW}\hat{c}_0 = -ik_j \left( v_i^{(0)*} \overline{\hat{w}_{i,0} \hat{w}_{j,0}} + \eta^{(0)*} \hat{\zeta}_0 \overline{\hat{w}_{j,0}} \right), \qquad (C.3) $$

$$ \frac{d\hat{c}_s}{d\hat{t}} + (\hat{D}_W + i2^{-3/2} \hat{\beta} \cos 2\theta_k) \hat{c}_s = -ik_j \left( v_i^{(s)*} \overline{\hat{w}_{i,0} \hat{w}_{j,s}} + \hat{w}_{i,s} \overline{\hat{w}_{j,0}} + \eta^{(s)*} \overline{\hat{\zeta}_0 \hat{w}_{j,s}} + \hat{\zeta}_s \overline{\hat{w}_{j,0}} \right) \quad s = \pm 1, \qquad (C.4) $$

with suitable interpretations of $\overline{\hat{w}_{i,s} \hat{w}_{j,s}}$ and $\overline{\hat{\zeta}_s \hat{w}_{j,s}}$ in terms of discrete Fourier transforms.

The quantities $v_i^{(s)}$ and $\eta^{(s)}$ appearing in (C.1)-(C.4) are given by (A.13), (A.14). There are however problems when $k_\perp = 0$, namely zero-by-zero divisions. We first consider the case $k \neq 0$, for which there is no such problem for $v_i^{(0)}$ and $\eta^{(s)}$. A difficulty arises for $s = \pm 1$ because (A.14) implies $v_i^{(s)}(\mathbf{k})$ has different values in the limit $k_\perp \to 0$, depending on the direction from which the limit is approached. This ambiguity for $k_\perp = 0$ is resolved by choosing a particular limiting direction. Taking the limit $k_\perp \searrow 0$, with $k_2 = 0$ and $k_3 > 0$ constant, (A.14) gives

$$ \begin{pmatrix} v_1^{(s)} \\ v_2^{(s)} \\ v_3^{(s)} \end{pmatrix} = 2^{-1/2} \begin{pmatrix} i s \\ 1 \\ 0 \end{pmatrix}, \quad \eta^{(s)} = 0 \quad s = \pm 1. \qquad (C.5) $$

To respect the identities $v_i^{(s)}(-\mathbf{k}) = v_i^{(-s)*}(\mathbf{k})$ and $\eta^{(s)}(-\mathbf{k}) = \eta^{(-s)*}(\mathbf{k})$, (C.5) is extended to $k_3 < 0$. This provides the values of $v_i^{(s)}(\mathbf{k})$ and $\eta^{(s)}(\mathbf{k})$ used in the numerical implementation for $k_\perp = 0$, $k_3 \neq 0$ and $s = \pm 1$. Given zero-by-zero division in both (A.13) and (A.14), a problem remains when $k=0$. This is resolved by taking $\hat{c}_s(0)=0$. Thus, the $\mathbf{k}=0$ terms in (C.1) and (C.2) are zero and (C.3) and (C.4) are not used for zero $\mathbf{k}$. Note that the right-hand sides of (C.3) and (C.4) are zero when $\mathbf{k}=0$, hence the choice $\hat{c}_s(0)=0$ is consistent with time evolution according to those equations.
In addition to supposing periodicity, the sums in (C.1) and (C.2) are truncated to $|l_i| \le K_i$, where $K_i$ are positive integers, which must be sufficiently large that the smallest dynamically significant length scales of the flow are resolved. Evaluating (C.1) and (C.2) at the discrete points $x_1 = p_1 L_1 / N_1$, $x_2 = p_2 L_2 / N_2$, $x_3 = p_3 L_3 / N_3$, where $N_i > 2K_i$ and $0 \le p_i < N_i$ are integers, the result can be expressed as

$$ \hat{w}_{i,s} = \sum_{q_1=0}^{N_1-1} \sum_{q_2=0}^{N_2-1} \sum_{q_3=0}^{N_3-1} \xi_{i,s} (q_1, q_2, q_3) \exp \left( 2\pi i \left( \frac{q_1 p_1}{N_1} + \frac{q_2 p_2}{N_2} + \frac{q_3 p_3}{N_3} \right) \right), \quad (C.6) $$

$$ \hat{\xi}_s = \sum_{q_1=0}^{N_1-1} \sum_{q_2=0}^{N_2-1} \sum_{q_3=0}^{N_3-1} \xi_s (q_1, q_2, q_3) \exp \left( 2\pi i \left( \frac{q_1 p_1}{N_1} + \frac{q_2 p_2}{N_2} + \frac{q_3 p_3}{N_3} \right) \right), \quad (C.7) $$

in which $\xi_{i,s} = \xi_s = 0$ unless $q_i \le K_i$ or $q_i \ge N_i - K_i$. Otherwise,

$$ \xi_{i,s}(q_1, q_2, q_3) = \hat{c}_s(\mathbf{k}) v_i^{(s)}(\mathbf{k}), \quad (C.8) $$

$$ \xi_s(q_1, q_2, q_3) = \hat{c}_s(\mathbf{k}) \eta^{(s)}(\mathbf{k}), \quad (C.9) $$

where $\mathbf{k} = 2\pi(l_1/L_1, l_2/L_2, l_3/L_3)$, $l_i = q_i$ when $q_i \le K_i$ and $l_i = q_i - N_i$ for $q_i \ge N_i - K_i$. (C.6) and (C.7) are discrete Fourier transforms which are evaluated using FFT.

Having determined $\hat{w}_{i,s}$ and $\hat{\xi}_s$ in physical space, the products $\hat{w}_{i,s'}\hat{w}_{j,s''}$ and $\hat{\xi}_{s'}\hat{w}_{j,s''}$ can be calculated. Performing the discrete Fourier transforms

$$ \chi_{js's''}(q_1, q_2, q_3) = \frac{1}{N_1 N_2 N_3} \sum_{p_1=0}^{N_1-1} \sum_{p_2=0}^{N_2-1} \sum_{p_3=0}^{N_3-1} \hat{w}_{i,s'} \hat{w}_{j,s''} (p_1, p_2, p_3) \exp \left( -2\pi i \left( \frac{q_1 p_1}{N_1} + \frac{q_2 p_2}{N_2} + \frac{q_3 p_3}{N_3} \right) \right), \quad (C.10) $$

$$ \lambda_{js's''}(q_1, q_2, q_3) = \frac{1}{N_1 N_2 N_3} \sum_{p_1=0}^{N_1-1} \sum_{p_2=0}^{N_2-1} \sum_{p_3=0}^{N_3-1} \hat{\xi}_{s'} \hat{w}_{j,s''} (p_1, p_2, p_3) \exp \left( -2\pi i \left( \frac{q_1 p_1}{N_1} + \frac{q_2 p_2}{N_2} + \frac{q_3 p_3}{N_3} \right) \right), \quad (C.11) $$

the nonlinear terms in (C.3) and (C.4) follow from

$$ \widetilde{\hat{w}}_{i,s} \widetilde{\hat{w}}_{j,s''}(\mathbf{k}) = \chi_{js's''}(q_1, q_2, q_3), \quad \widetilde{\hat{\xi}}_{s'} \widetilde{\hat{w}}_{j,s''}(\mathbf{k}) = \lambda_{js's''}(q_1, q_2, q_3), \quad (C.12) $$

where $q_i = l_i$ for $0 \le l_i \le K_i$ and $q_i = N_i + l_i$ for $-K_i \le l_i < 0$. As for (C.6) and (C.7), FFT is used to evaluate (C.10) and (C.11).

Writing the numerical mode-amplitude equations, (C.3) and (C.4), as
$$ \frac{d\hat{c}_s}{d\hat{t}} + \hat{D}_{ss}\hat{c}_s = \gamma_s, \qquad (C.13) $$

a time-stepping scheme is introduced. Time is discretized to the values $\hat{t}_n = n\Delta$, where $\Delta$ is the time step, and (C.13) gives

$$ \hat{c}_s(\hat{t}_{n+1}) = \hat{c}_s(\hat{t}_n)\exp(-\hat{D}_{ss}\Delta) + \int_{\hat{t}_n}^{\hat{t}_{n+1}} \gamma_s(\hat{t}')\exp(\hat{D}_{ss}(\hat{t}' - \hat{t}_{n+1}))d\hat{t}'. \qquad (C.14) $$

Given $\hat{c}_s(\hat{t}_n)$, $\gamma_s(\hat{t}_n)$ is determined using FFT as described above. As a first approximation, $\gamma_s(\hat{t}')$ is taken as constant and equal to $\gamma_s(\hat{t}_n)$, hence

$$ \hat{c}_s^\dagger = \hat{c}_s(\hat{t}_n)\exp(-\hat{D}_{ss}\Delta) + \gamma_s(\hat{t}_n)\frac{1 - \exp(-\hat{D}_{ss}\Delta)}{\hat{D}_{ss}} \qquad (C.15) $$

as an approximate value of $\hat{c}_s(\hat{t}_{n+1})$. To increase the precision to second order, let $\gamma_s^\dagger$ be the value obtained using $(\hat{c}_s(\hat{t}_n) + \hat{c}_s^\dagger)/2$, which approximates the value of $\hat{c}_s(\hat{t})$ at time $\hat{t}_n + \Delta/2$. Thus,

$$ \hat{c}_s(\hat{t}_{n+1}) = \hat{c}_s(\hat{t}_n)\exp(-\hat{D}_{ss}\Delta) + \gamma_s^\dagger \frac{1 - \exp(-\hat{D}_{ss}\Delta)}{\hat{D}_{ss}} \qquad (C.16) $$

completes the time step. The advantage of the above scheme, compared with more traditional ones for integrating differential equations, is that it only assumes that the time scale for significant changes in $\hat{c}_s(\hat{t})$ is long compared with $\Delta$. As $k$ increases, so does $\hat{D}_{ss}(\mathbf{k})$, and, when $\hat{D}_{ss}(\mathbf{k})\Delta$ becomes of $O(1)$, dissipative effects are significant within a time step. At this point, the precision of traditional schemes is degraded, whereas it is maintained for the above scheme provided $\hat{c}_s(\hat{t})$ does not change too much over a single step.

According to (C.1) and (C.2) and the scaled versions of (2.7) and (2.19),

$$ \hat{u}_i = \sum_{s=0,\pm 1} \exp(-ist/2^{1/2}) \sum_{\mathbf{k}} \hat{c}_s(\mathbf{k}) v_i^{(s)}(\mathbf{k}) \exp(i\mathbf{k}\cdot\mathbf{x}), \qquad (C.17) $$

$$ \hat{\eta} = \sum_{s=0,\pm 1} \exp(-ist/2^{1/2}) \sum_{\mathbf{k}} \hat{c}_s(\mathbf{k}) \eta^{(s)}(\mathbf{k}) \exp(i\mathbf{k}\cdot\mathbf{x}), \qquad (C.18) $$

thus
$$ \overline{\hat{u}_i \hat{u}_i} + \overline{\hat{\eta}^2} = \sum_{s,s'=0,\pm 1} \exp(i(s-s')t/2^{1/2}) \overline{\hat{c}_s^*(\mathbf{k}) \hat{c}_{s'}(\mathbf{k}')} \left( \nu_i^{(s)}(\mathbf{k}) \nu_{i'}^{(s')}(\mathbf{k}') + \eta^{(s)*}(\mathbf{k}) \eta^{(s')}(\mathbf{k}') \right) \exp(i(\mathbf{k}'-\mathbf{k})\cdot\mathbf{x}) \quad (C.19) $$

Given homogeneity, $\overline{\hat{c}_s^*(\mathbf{k}) \hat{c}_{s'}(\mathbf{k}')} = 0$ when $\mathbf{k}' \neq \mathbf{k}$, hence

$$ \overline{\hat{u}_i \hat{u}_i + \hat{\eta}^2} = \sum_{s,s'=0,\pm 1} \exp(i(s-s')t/2^{1/2}) \overline{\hat{c}_s^*(\mathbf{k}) \hat{c}_{s'}(\mathbf{k})} \left( \nu_i^{(s)}(\mathbf{k}) \nu_{i'}^{(s')}(\mathbf{k}) + \eta^{(s)*}(\mathbf{k}) \eta^{(s')}(\mathbf{k}) \right) \quad (C.20) $$

(I.2.13) implies that the vectors $\mathbf{e}^{(1)}$ and $\mathbf{e}^{(2)}$ are orthonormal, hence (I.2.12) yields

$$ \nu_i^{(s)*}(\mathbf{k}) \nu_i^{(s')}(\mathbf{k}) = \sum_{l=1}^{2} u_s^{(l)*} u_{s'}^{(l)}. \quad (C.21) $$

Using this result in (C.20), (I.2.16) gives

$$ \frac{1}{2} (\overline{\hat{u}_i \hat{u}_i} + \overline{\hat{\eta}^2}) = \frac{1}{2} \sum_{\mathbf{k}} \sum_{s=0,\pm 1} \overline{|\hat{c}_s(\mathbf{k})|^2}. \quad (C.22) $$

Prior to numerical approximation, the scaled version of (2.12) implies

$$ \frac{1}{2} (\overline{\hat{u}_i \hat{u}_i} + \overline{\hat{\eta}^2}) = \int \hat{e}(\mathbf{k}) d^3 \mathbf{k}. \quad (C.23) $$

(C.22) is a numerical approximation of (C.23), the integral being represented by a sum. Since the density of discrete $\mathbf{k}$ in spectral space is $L_1L_2L_3/8\pi^3$,

$$ \hat{e}(\mathbf{k}) = \frac{L_1 L_2 L_3}{16\pi^3} \sum_{s=0,\pm 1} \overline{|\hat{c}_s(\mathbf{k})|^2} \quad (C.24) $$

approximates the spectral energy density. Comparing the scaled version of (2.13) and (C.24),

$$ \hat{A}_{ss} = \frac{L_1 L_2 L_3}{8\pi^3} \overline{|\hat{c}_s(\mathbf{k})|^2} \quad (C.25) $$

provides a numerical approximation of $\hat{A}_{ss} = A_{ss} / \epsilon^2$.

The initial $\overline{\hat{c}_s(\mathbf{k})}$ must satisfy various constraints. Firstly, $\overline{\hat{c}_s(-\mathbf{k})} = \overline{\hat{c}_{-s}^*(\mathbf{k})}$ because the velocity and buoyancy variable of the flow are real. Secondly, they should be of zero mean and, as noted above, homogeneity requires $\overline{\hat{c}_s^*(\mathbf{k}) \hat{c}_{s'}(\mathbf{k}')} = 0$ when $\mathbf{k}' \neq \mathbf{k}$. Finally,
$$ \overline{|\hat{c}_s(\mathbf{k})|^2} = \frac{8\pi^3}{L_1 L_2 L_3} \hat{A}_{ss} \qquad (C.26) $$

from (C.25). Supposing the initial $\hat{A}_{ss}(\mathbf{k})$ are given and are such that $\hat{A}_{ss}(-\mathbf{k}) = \hat{A}_{-s,-s}(\mathbf{k})$, these conditions are satisfied by taking

$$ \hat{c}_s(\mathbf{k}) = \left( \frac{8\pi^3}{L_1 L_2 L_3} \hat{A}_{ss}(\mathbf{k}) \right)^{1/2} \exp(i\phi_s(\mathbf{k})), \qquad (C.27) $$

where the phases, $\phi_s(\mathbf{k})$, are random variables, uniformly distributed between $-\pi$ and $\pi$ and independent apart from the constraint $\phi_s(-\mathbf{k}) = -\phi_s(\mathbf{k})$. (C.27) also implies statistical homogeneity of the initial flow defined by (C.1), (C.2), another important requirement.

Having described the methods used to initialise and integrate (C.3) and (C.4), we now want to extract information from the results. The problem with using (C.25) to determine $\hat{A}_{ss}(\mathbf{k})$ is that it contains an ensemble average. To calculate such an average numerically would require very many runs, which, given the computational cost, is not attempted. We instead extract information from a single run. In so doing, we assume that, prior to discretization, the flow is statistically axisymmetric. In particular, $\hat{A}_{ss}(\mathbf{k})$ is axisymmetric and thus has the form $\hat{A}_{ss} = \hat{A}_{ss}(k_{\perp}, k_{\|})$, where $k_{\perp} = k \sin\theta_{\perp}$ and $k_{\|} = k \cos\theta_{\|}$ are wavenumbers perpendicular to and parallel to the rotation vector. The aim here is to obtain results which are representative of $\hat{A}_{ss}(k_{\perp}, k_{\|})$ without using ensemble averaging.

For simplicity's sake, from here on we assume $L_1 = L_2 = L_3 = L$ and $K_1 = K_2 = K_3 = K$. Let

$$ k_m = \frac{\pi}{L} (2m - 1) \qquad (C.28) $$

and define axisymmetric regions $D(m_{\perp}, m_{\|})$ in spectral space by $k_{m_{\perp}} \le k_{\perp} < k_{m_{\perp}+1}$ and $k_{m_{\|}} \le k_{\|} < k_{m_{\|}+1}$, where $0 \le m_{\perp} \le K$ and $-K \le m_{\|} \le K$. Let

$$ B_s^{m_{\perp}, m_{\|}} = \frac{L^3}{8\pi^3 N_{m_{\perp}, m_{\|}}} \sum_{\mathbf{k} \in D(m_{\perp}, m_{\|})} |\hat{c}_s(\mathbf{k})|^2, \qquad (C.29) $$

where the sum is over discrete wave vectors inside $D(m_{\perp}, m_{\|})$ and $N_{m_{\perp}, m_{\|}}$ is the number of such wave vectors (which is non-zero given the construction of $D(m_{\perp}, m_{\|})$). Taking the average of (C.29) and using (C.25),
$$ \overline{B_s^{m_{\perp}, m_{\parallel}}} = \frac{1}{N_{m_{\perp}, m_{\parallel}}} \sum_{\mathbf{k} \in D(m_{\perp}, m_{\parallel})} \hat{A}_{ss}(\mathbf{k}). \qquad (C.30) $$

If the region is small in the $k_{\perp} - k_{\parallel}$ plane (large enough $L$), variations of $\hat{A}_{ss}(\mathbf{k})$ across $D(m_{\perp}, m_{\parallel})$ can be neglected and $\overline{B_s^{m_{\perp}, m_{\parallel}}}$ approximates $\hat{A}_{ss}(\mathbf{k})$ throughout the given region, in particular at the point

$$ k_{\perp} = \frac{2\pi}{L} m_{\perp} \qquad 0 \le m_{\perp} \le K, \qquad (C.31) $$

$$ k_{\parallel} = \frac{2\pi}{L} m_{\parallel} \qquad -K \le m_{\parallel} \le K. \qquad (C.32) $$

The problem is that, as noted above, it is not feasible to determine the ensemble average using numerical calculations. Thus, we take $B_s^{m_{\perp}, m_{\parallel}}$, rather than its average, as representative of $\hat{A}_{ss}(k_{\perp}, k_{\parallel})$ at (C.31), (C.32).

The spherically averaged energy spectrum for $k$ given by

$$ k = \frac{2\pi}{L} m \qquad 0 \le m \le K \qquad (C.33) $$

is determined using

$$ \hat{E}(k) = \frac{L}{4\pi} \sum_{s=0, \pm 1} \sum_{\mathbf{k} \in S(m)} |\hat{c}_s(\mathbf{k})|^2, \qquad (C.34) $$

where $S(m)$ denotes the region $k_m \le k < k_{m+1}$ in spectral space. The wave and NP parts of this spectrum are obtained by corresponding restrictions of the sum over $s$. The total energy follows from (C.22) (without averaging) and wave/NP decomposition results from restricting the sum over $s$.

Finally, the spectral fluxes are obtained as follows. The unaveraged wave energy inside the sphere of radius $k$ is

$$ \frac{1}{2} \sum_{|\mathbf{k}'|<k} \sum_{s=\pm 1} |\hat{c}_s(\mathbf{k}')|^2, \qquad (C.35) $$

whose time derivative is

$$ \sum_{|\mathbf{k}'|<k} \sum_{s=\pm 1} \Re(\hat{c}_s^*(\mathbf{k}') \gamma_s(\mathbf{k}')) - \sum_{|\mathbf{k}'|<k} \sum_{s=\pm 1} D_W(\mathbf{k}') |\hat{c}_s(\mathbf{k}')|^2 \qquad (C.36) $$
according to (C.13). Comparison with (3.4) gives

$$ \Phi_W(k) = - \sum_{|\mathbf{k}'|<k} \sum_{s=\pm 1} \Re(\hat{c}_s^*(\mathbf{k}') \gamma_s(\mathbf{k}')). \qquad (C.37) $$

Similarly,

$$ \Phi_{NP}(k) = - \sum_{|\mathbf{k}'|<k} \Re(\hat{c}_0^*(\mathbf{k}') \gamma_0(\mathbf{k}')) \qquad (C.38) $$

for the NP radial flux. (C.37) and (C.38) are applied for the wavenumbers (C.33).

The numerical parameters used to obtain the results described in section 3 were $N_i = 512$, $K_i = 255$, $L_i = 60$, $\Delta = 5 \times 10^{-3}$ and $d = 1.5 \times 10^{-7}$, where, as noted in the main text, the hyperviscous damping factor is $\hat{D}_{ss}(\mathbf{k}) = dk^6$.

## References

* Bartello, P. 1995 Geostrophic adjustment and inverse cascades in rotating stratified turbulence. J. Atmos. Sci., 52, 4410–4428.

* Benney, D.J., Saffman, P.G. 1966 Nonlinear interactions of random waves in a dispersive medium. Proc. Roy. Soc. A 289, 301–320.

* Benney, D.J., Newell, A.C. 1969 Random wave closures. Stud. Appl. Math. 48, 29–53.

* Bourouiba, L., Bartello, P. 2007 The intermediate Rossby number range and two-dimensional–three-dimensional transfers in rotating decaying homogeneous turbulence. J. Fluid Mech. 587, 139–161.

* Cambon, C., Jacquin, L. 1989 Spectral approach to non-isotropic turbulence subjected to rotation. J. Fluid Mech. 202, 295–317.

* Cambon, C., Mansour, N.N., Godeferd, F.S. 1997 Energy transfer in rotating turbulence. J. Fluid Mech. 337, 303–332.

* Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, T.A. 1988 Spectral Methods in Fluid Dynamics. Springer-Verlag.

* Charney, J.G. 1971 Geostrophic Turbulence. J. Atmos. Sci., 28, 1087–1095.

* Coleman, G.N., Ferziger, J.H., Spalart, P.R. 1992 Direct simulation of the stably stratified turbulent Ekman layer. J. Fluid Mech. 244, 677–712. (Corrigendum: J. Fluid Mech. 252, 721.)

* Galtier, S. 2022 Physics of wave turbulence. Cambridge University Press.
Godeferd, F.S., Staquet, C. 2003 Statistical modelling and direct numerical simulations of decaying stably stratified turbulence. Part 2. Large-scale and small-scale anisotropy. J. Fluid Mech. 486, 115-159.

Godeferd, F.S., Cambon, C. 1994 Detailed investigation of energy transfers in homogeneous stratified turbulence. Phys. Fluids 6, 2084-2100.

Haugen, N.E.L, Brandenburg, A. 2004 Inertial range scaling in numerical turbulence with hyperviscosity. Phys. Rev. E 70, 026405.

Liechtenstein, L., Godeferd, F.S., Cambon, C. 2005 Nonlinear formation of structures in rotating stratified turbulence. J. Turbul. 6, 1-18.

Newell, A.C., Rumpf, B. 2011 Wave turbulence. Annu. Rev. Fluid Mech. 43, 59-78.

Orszag, S.A., Patterson, G.S. 1972 Numerical simulation of three-dimensional homogeneous isotropic turbulence. Phys. Rev. Lett. 28, 76-79.

Pedlosky, J. 1987. Geophysical Fluid Dynamics, Springer-Verlag, New York.

Sagaut, P., Cambon, C. 2018 Homogeneous Turbulence Dynamics (second edition). Springer.

Scott, J.F., Cambon, C. 2024 Evolution of weak, homogeneous turbulence with rotation and stratification. J. Fluid Mech. 979, A17.

Scott, J.F. 2025 Evolution of weak, homogeneous turbulence subject to rotation and stratification: Comparable wave and non-propagating components. Phys. Rev. E 111, 035101.

Smith, L. M., Waleffe, F. 1999 Transfer of energy to two-dimensional large scales in forced, rotating three-dimensional turbulence. Phys. Fluids 11, 1608-1622.

Smith, L.M., Waleffe, F. 2002 Generation of slow large scales in forced rotating stratified turbulence. J. Fluid Mech. 451, 145-168.

Thomas, J., Daniel, D. 2020. Turbulent exchanges between near-inertial waves and balanced flows. J. Fluid Mech., 902, A7.

Thomas, J., Daniel, D. 2021 Forward flux and enhanced dissipation of geostrophic balanced energy. J. Fluid Mech. 911, A60.

Thomas, J. 2023 Turbulent wave-balance exchanges in the ocean, Proc. R. Soc. A 479, 20220565.

Vallgren, A., Lindborg, E. 2010 Charney isotropy and equipartition in quasi-geostrophic turbulence. J. Fluid Mech., 656, 448-457.

Waite, M. L., Bartello, P. 2004 Stratified turbulence dominated by vortical motion, J. Fluid Mech., 517, 281-308.
Waite, M. L., Bartello, P. 2006a Stratified turbulence generated by internal gravity waves, J. Fluid Mech., 546, 313–339.

Waite, M. L., Bartello, P. 2006b, The transition from geostrophic to stratified turbulence, J. Fluid Mech., 568, 89-108.

Zakharov, V.E., Lvov, V.S., Falkovich, G.E. 1992 Kolmogorov Spectra of Turbulence I – Wave Turbulence. Springer-Verlag.