2.1. Fourier transforms and modes

Defining the Fourier transforms

(2.4)\begin{gather}{\tilde{u}_i}(\boldsymbol{k}) = \frac{1}{{8{{\rm \pi} ^3}}}\int {{u_i}(\boldsymbol{x})\exp ( - \textrm{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{x}} ,\end{gather}

(2.5)\begin{gather}\tilde{\eta }(\boldsymbol{k}) = \frac{1}{{8{{\rm \pi} ^3}}}\int {\eta (\boldsymbol{x})\exp ( - \textrm{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{x}} ,\end{gather}

with similar definitions of $\tilde{p}$, $\widetilde {{u_i}{u_j}}$ and $\widetilde {\eta {u_j}}$, (2.1)–(2.3) yield

(2.6)\begin{gather}\frac{{\partial {{\tilde{u}}_i}}}{{\partial t}} + \varOmega {\varepsilon _{ijk}}{e_j}{\tilde{u}_k} + N\tilde{\eta }{e_i} + \textrm{i}{k_i}\tilde{p} ={-} {\rm i}{k_j}\widetilde {{u_i}{u_j}} - {D_u}{k^2}{\tilde{u}_i},\end{gather}

(2.7)\begin{gather}{k_i}{\tilde{u}_i} = 0,\end{gather}

(2.8)\begin{gather}\frac{{\partial \tilde{\eta }}}{{\partial t}} - N{e_i}{\tilde{u}_i} ={-} \textrm{i}{k_j}\widetilde {\eta {u_j}} - {D_\eta }{k^2}\tilde{\eta },\end{gather}

where $k = |\boldsymbol{k}|$.

Appendix A examines the consequences of (2.6)–(2.8). In the absence of nonlinearity and visco-diffusion, the right-hand sides of (2.6) and (2.8) are zero. There are then three linearly independent solutions, referred to as modes. These solutions have time dependence $\exp ( - \textrm{i}s\omega (\boldsymbol{k})t)$, where s takes one of the three values $s = 0,\; \pm 1$,

(2.9)\begin{equation}\omega (\boldsymbol{k}) = \frac{{{{({N^2}k_ \bot ^2 + {\varOmega ^2}k_\parallel ^2)}^{1/2}}}}{k} = {({N^2}{\sin ^2}{\theta _{\boldsymbol{k}}} + {\varOmega ^2}{\cos ^2}{\theta _{\boldsymbol{k}}})^{1/2}},\end{equation}

${k_ \bot } = |\boldsymbol{k} \times \boldsymbol{e}|$ is the transverse wavenumber, ${k_\parallel } = \boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{e}$ is the axial wavenumber, and $0 \le {\theta _{\boldsymbol{k}}} \le {\rm \pi}$ is the angle between $\boldsymbol{k}$ and $\boldsymbol{e}$. The modes with $s ={\pm} 1$ represent inertial-gravity waves, for which (2.9) is the dispersion relation, while that with $s = 0$ will be referred to as the non-propagating (NP) mode because, when regarded as a wave, it has zero group velocity. The modal solutions of (2.6)–(2.8) (with zero right-hand sides) are

(2.10)\begin{gather}{\tilde{u}_i} = v_i^{(s)}(\boldsymbol{k})\exp ( - \textrm{i}s\omega (\boldsymbol{k})t),\end{gather}

(2.11)\begin{gather}\tilde{\eta } = {\eta ^{(s)}}(\boldsymbol{k})\exp ( - \textrm{i}s\omega (\boldsymbol{k})t),\end{gather}

where

(2.12)\begin{equation}v_i^{(s)}(\boldsymbol{k}) = \sum\limits_{l = 1}^2 {u_s^{(l)}(\boldsymbol{k})e_i^{(l)}(\boldsymbol{k})} ,\end{equation}

(2.13a,b)\begin{equation}{\boldsymbol{e}^{(1)}} = \frac{{\boldsymbol{k} \times \boldsymbol{e}}}{{{k_ \bot }}},\quad {\boldsymbol{e}^{(2)}} = \frac{{\boldsymbol{k} \times {\boldsymbol{e}^{(1)}}}}{k},\end{equation}

are unit vectors orthogonal to each other and to $\boldsymbol{k}$,

(2.14)\begin{equation}\left( {\begin{array}{@{}c@{}} {u_s^{(1)}}\\ {u_s^{(2)}}\\ {{\eta^{(s)}}} \end{array}} \right) = \frac{1}{{{2^{1/2}}{{({N^2}k_ \bot ^2 + {\varOmega ^2}k_\parallel ^2)}^{1/2}}}}\left( {\begin{array}{@{}c@{}} { - \varOmega {k_\parallel }}\\ {{\rm i}s{{({N^2}k_ \bot^2 + {\varOmega^2}k_\parallel^2)}^{1/2}}}\\ {N{k_ \bot }} \end{array}} \right),\end{equation}

when $s ={\pm} 1$ and

(2.15)\begin{equation}\left( {\begin{array}{@{}c@{}} {u_0^{(1)}}\\ {u_0^{(2)}}\\ {{\eta^{(0)}}} \end{array}} \right) = \textrm{i}\frac{1}{{{{({N^2}k_ \bot ^2 + {\varOmega ^2}k_\parallel ^2)}^{1/2}}}}\left( {\begin{array}{@{}c@{}} {N{k_ \bot }}\\ 0\\ {\varOmega {k_\parallel }} \end{array}} \right).\end{equation}

It follows from (2.12)–(2.15) and the definitions of ${k_ \bot }$ and ${k_\parallel }$ that $v_i^{(s)}( - \boldsymbol{k}) = v_i^{{{( - s)}^\ast }}(\boldsymbol{k})$ and ${\eta ^{(s)}}( - \boldsymbol{k}) = {\eta ^{{{( - s)}^\ast }}}(\boldsymbol{k})$, where * denotes complex conjugation. Note the orthonormality relation

(2.16)\begin{equation}\sum\limits_{l = 1}^2 {u_s^{{{(l)}^\ast }}u_{s^{\prime}}^{(l)}} + {\eta ^{{{(s)}^\ast }}}{\eta ^{(s^{\prime})}} = {\delta _{ss^{\prime}}},\end{equation}

where ${\delta _{ss^{\prime}}}$ is the Kronecker delta. Note also that the $s = 0$ modes, referred to as NP here, are often described in the geophysical literature as ‘vortical modes’ or ‘potential vorticity (PV) modes’.

2.2. Mode-amplitude equation

The modes form a complete set for ${\tilde{u}_i}$, $\tilde{\eta }$ satisfying (2.7). Thus,

(2.17)\begin{gather}{\tilde{u}_i} = \sum\limits_{s = 0, \pm 1} {{b_s}(\boldsymbol{k},t)v_i^{(s)}(\boldsymbol{k})} ,\end{gather}

(2.18)\begin{gather}\tilde{\eta } = \sum\limits_{s = 0, \pm 1} {{b_s}(\boldsymbol{k},t){\eta ^{(s)}}(\boldsymbol{k})} ,\end{gather}

where, as shown in Appendix A, the modal coefficients evolve according to

(2.19)\begin{equation}\frac{{\partial {b_s}}}{{\partial t}} + \textrm{i}s\omega {b_s} ={-} \textrm{i}{k_j}(v_i^{{{(s)}^\ast }}\widetilde {{u_i}{u_j}} + {\eta ^{{{(s)}^\ast }}}\widetilde {\eta {u_j}}) - \sum\limits_{\hat{s} = 0, \pm 1} {{D_{s\hat{s}}}{b_{\hat{s}}}} ,\end{equation}

in which

(2.20)\begin{equation}{D_{s\hat{s}}}(\boldsymbol{k}) = {k^2}\left( {{D_u}\sum\limits_{l = 1}^2 {u_s^{{{(l)}^\ast }}(\boldsymbol{k})u_{\hat{s}}^{(l)}(\boldsymbol{k})} + {D_\eta }{\eta^{{{(s)}^\ast }}}(\boldsymbol{k}){\eta^{(\hat{s})}}(\boldsymbol{k})} \right)\end{equation}

is a Hermitian, positive definite matrix expressing viscosity and diffusion. Real ${u_i}$ and $\eta$ requires that ${\tilde{u}_i}( - \boldsymbol{k}) = \tilde{u}_i^\ast (\boldsymbol{k})$ and $\tilde{\eta }( - \boldsymbol{k}) = {\tilde{\eta }^\ast }(\boldsymbol{k})$, conditions which, according to (2.17), (2.18), $v_i^{(s)}( - \boldsymbol{k}) = v_i^{{{( - s)}^\ast }}(\boldsymbol{k})$ and ${\eta ^{(s)}}( - \boldsymbol{k}) = {\eta ^{{{( - s)}^\ast }}}(\boldsymbol{k})$, are met provided ${b_s}( - \boldsymbol{k}) = b_{ - s}^\ast (\boldsymbol{k})$.

Using (2.17), (2.18), and the inverse transforms of (2.4) and (2.5), the total flow can be expressed as the sum of wave and NP components:

(2.21a,b)\begin{equation}{u_i} = u_i^W + u_i^{NP},\quad \eta = {\eta ^W} + {\eta ^{NP}},\end{equation}

where

(2.22)\begin{gather}u_i^W = \sum\limits_{s ={\pm} 1} {\int {{b_s}(\boldsymbol{k})v_i^{(s )}(\boldsymbol{k})\exp (\textrm{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{k}} } ,\end{gather}

(2.23)\begin{gather}{\eta ^W} = \sum\limits_{s ={\pm} 1} {\int {{b_s}(\boldsymbol{k}){\eta ^{(s)}}(\boldsymbol{k})\exp (\textrm{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{k}} } ,\end{gather}

(2.24)\begin{gather}u_i^{NP} = \int {{b_0}(\boldsymbol{k})v_i^{(0)}(\boldsymbol{k})\exp (\textrm{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{k}} ,\end{gather}

(2.25)\begin{gather}{\eta ^{NP}} = \int {{b_0}(\boldsymbol{k}){\eta ^{(0)}}(\boldsymbol{k})\exp (\textrm{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{k}} .\end{gather}

Given ${b_s}( - \boldsymbol{k}) = b_{ - s}^\ast (\boldsymbol{k})$, $v_i^{(s)}( - \boldsymbol{k}) = v_i^{{{( - s)}^\ast }}(\boldsymbol{k})$ and ${\eta ^{(s)}}( - \boldsymbol{k}) = {\eta ^{{{( - s)}^\ast }}}(\boldsymbol{k})$, both components are real, the $s ={\pm} 1$ contributions to the wave component being complex conjugates.

If nonlinearity, viscosity and diffusion were neglected, the solution of (2.19) would have the expected modal form ${b_s} \propto \exp ( - \textrm{i}s\omega t)$. The resulting oscillations of ${b_{ {\pm} 1}}$ due to waves can be supressed by defining ${a_s} = {b_s}\exp (\textrm{i}s\omega t)$, which evolves according to

(2.26)\begin{equation}\frac{{\partial {a_s}}}{{\partial t}} ={-} \textrm{i}{k_j}(v_i^{{{(s)}^\ast }}\widetilde {{u_i}{u_j}} + {\eta ^{{{(s)}^\ast }}}\widetilde {\eta {u_j}})\exp (\textrm{i}s\omega t) - \sum\limits_{\hat{s} = 0, \pm 1} {{D_{s\hat{s}}}\exp (\textrm{i}(s - \hat{s})\omega t){a_{\hat{s}}}} .\end{equation}

In the absence of nonlinearity and visco-diffusion, the mode amplitudes ${a_s}$ are time-independent, whereas, when we later specialise to weak turbulence and small visco-diffusion, they evolve slowly with time. Here, ${b_s}( - \boldsymbol{k}) = b_{ - s}^\ast (\boldsymbol{k})$ and the definition of ${a_s}$ imply ${a_s}( - \boldsymbol{k}) = a_{ - s}^\ast (\boldsymbol{k})$.

As shown in Appendix A, the nonlinear term in (2.26) can be expressed as

(2.27) \begin{align}& -\! \textrm{i}{k_j}(v_i^{{{(s)}^\ast }}\widetilde {{u_i}{u_j}} + {\eta ^{{{(s)}^\ast }}}\widetilde {\eta {u_j}})\exp (\textrm{i}s\omega t)\nonumber\\ & \quad = \sum\limits_{{s_p},{s_q} = 0, \pm 1} {\int {N_{s{s_p}{s_q}}^\ast (\boldsymbol{k},\boldsymbol{p})a_{{s_p}}^\ast (\kern1pt \boldsymbol{p})a_{{s_q}}^\ast ( - \boldsymbol{k} - \boldsymbol{p})\exp (\textrm{i}{F_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p})t)\,{\textrm{d}^3}\boldsymbol{p}} },\end{align}

where

(2.28)\begin{equation}{N_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = \textrm{i}{k_j}v_j^{({s_q})}( - \boldsymbol{k} - \boldsymbol{p})(v_i^{(s)}(\boldsymbol{k})v_i^{({s_p})}(\kern1pt \boldsymbol{p}) + {\eta ^{(s)}}(\boldsymbol{k}){\eta ^{({s_p})}}(\kern1pt \boldsymbol{p}))\end{equation}

represents nonlinear coupling between modes and

(2.29)\begin{equation}{F_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = s\omega (\boldsymbol{k}) + {s_p}\omega (\kern1pt \boldsymbol{p}) + {s_q}\omega ( - \boldsymbol{k} - \boldsymbol{p}).\end{equation}

A symmetrised version of (2.27), namely

(2.30) \begin{align} & -\! \textrm{i}{k_j}(v_i^{{{(s)}^\ast }}\widetilde {{u_i}{u_j}} + {\eta ^{{{(s)}^\ast }}}\widetilde {\eta {u_j}})\exp (\textrm{i}s\omega t)\nonumber\\ & \quad = \sum\limits_{{s_p},{s_q} = 0, \pm 1} {\int {M_{s{s_p}{s_q}}^\ast (\boldsymbol{k},\boldsymbol{p})a_{{s_p}}^\ast (\kern1pt \boldsymbol{p})a_{{s_q}}^\ast ( - \boldsymbol{k} - \boldsymbol{p})\exp (\textrm{i}{F_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p})t)\,{\textrm{d}^3}\boldsymbol{p}} },\end{align}

where

(2.31)\begin{equation}{M_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = {\textstyle{1 \over 2}}({N_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) + {N_{s{s_q}{s_p}}}(\boldsymbol{k}, - \boldsymbol{k} - \boldsymbol{p}))\end{equation}

is also derived in Appendix A. Using either (2.27) or (2.30), (2.26) provides an evolution equation for ${a_s}$. Both versions are employed in what follows. Note the symmetries ${M_{s{s_p}{s_q}}}(\boldsymbol{k}, - \boldsymbol{k} - \boldsymbol{p}) = {M_{s{s_q}{s_p}}}(\boldsymbol{k},\boldsymbol{p})$ and ${F_{s{s_p}{s_q}}}(\boldsymbol{k}, - \boldsymbol{k} - \boldsymbol{p}) = {F_{s{s_q}{s_p}}}(\boldsymbol{k},\boldsymbol{p})$. Note also that the wave vectors $\boldsymbol{p}$ and $\boldsymbol{q} ={-} \boldsymbol{k} - \boldsymbol{p}$, which appear in (2.27) and (2.30), satisfy the usual condition, $\boldsymbol{k} + \boldsymbol{p} + \boldsymbol{q} = 0$, for formation of a triad with $\boldsymbol{k}$.

2.3. The spectral matrix and energy

Here and henceforth, the random nature of turbulent flow is recognised and ensemble averaging introduced. The flow is assumed statistically homogeneous, i.e. its statistical properties are the same at all spatial locations, in particular, any one-point average is independent of position. Averaging (2.1) and (2.3) eliminates the nonlinear, pressure and visco-diffusive terms, while the average of (2.2) is automatically satisfied. We have in mind that there is no mean flow, i.e. $\overline {{u_i}} = 0$, where the overbar denotes ensemble averaging, hence the average of (2.1) gives $\bar{\eta } = 0$. Here, $\overline {{u_i}} = \bar{\eta } = 0$ leads to $\overline {{a_s}} = \overline {{b_s}} = 0$, while (2.22)–(2.25) imply that the wave and NP components are both of zero mean. Those components also inherit the statistical homogeneity of the total flow.

The ensemble-averaged, non-dimensional energy density in physical space is the sum of kinetic and potential energies:

(2.32)\begin{equation}\underbrace{{{\textstyle{1 \over 2}}\overline {{u_i}{u_i}} }}_{{Kinetic}} + \underbrace{{{\textstyle{1 \over 2}}\overline {{\eta ^2}} }}_{{Potential}}.\end{equation}

Given statistical homogeneity,

(2.33a,b)\begin{equation}\overline {\tilde{u}_i^\ast (\boldsymbol{k}){{\tilde{u}}_i}(\boldsymbol{k^{\prime}})} = {\varPhi _u}(\boldsymbol{k})\delta (\boldsymbol{k} - \boldsymbol{k^{\prime}}),\quad \overline {{{\tilde{\eta }}^\ast }(\boldsymbol{k})\tilde{\eta }(\boldsymbol{k^{\prime}})} = {\varPhi _\eta }(\boldsymbol{k})\delta (\boldsymbol{k} - \boldsymbol{k^{\prime}}),\end{equation}

in which $\delta$ represents the Dirac function. Hence, using the inverse transform of (2.4),

(2.34) \begin{align} \dfrac{1}{2}\overline {{u_i}{u_i}} &= \dfrac{1}{2}\int {\int {\overline {\tilde{u}_i^\ast (\boldsymbol{k}){{\tilde{u}}_i}(\boldsymbol{k^{\prime}})} \exp (\textrm{i}(\boldsymbol{k^{\prime}} - \boldsymbol{k})\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{k^{\prime}}} \,{\textrm{d}^3}\boldsymbol{k}}\nonumber \\ &= \dfrac{1}{2}\int {{\varPhi _u}(\boldsymbol{k})\,{\textrm{d}^3}\boldsymbol{k}} , \end{align}

which indicates that the kinetic-energy density in spectral space is ${e_K}(\boldsymbol{k}) = {\varPhi _u}/2$. Similarly, the potential-energy density in spectral space is ${e_P}(\boldsymbol{k}) = {\varPhi _\eta }/2$.

Again, using statistical homogeneity,

(2.35)\begin{equation}\overline {a_s^\ast (\boldsymbol{k}){a_{s^{\prime}}}(\boldsymbol{k^{\prime}})} = {A_{ss^{\prime}}}(\boldsymbol{k})\delta (\boldsymbol{k} - \boldsymbol{k^{\prime}}),\end{equation}

where ${A_{ss^{\prime}}}$ is the spectral matrix and s, $s^{\prime}$ run over the values $- 1$, $0$ and $+ 1$; ${A_{ss^{\prime}}}$ is Hermitian and positive semi-definite. In particular, the diagonal elements of ${A_{ss^{\prime}}}$ are real and non-negative. Also, given ${a_s}( - \boldsymbol{k}) = a_{ - s}^\ast (\boldsymbol{k})$, ${A_{ss^{\prime}}}( - \boldsymbol{k}) = {A_{ - s^{\prime}, - s}}(\boldsymbol{k})$. Using (2.12), (2.17), (2.18), (2.35), $e_i^{(l)}(\boldsymbol{k})e_i^{(l^{\prime})}(\boldsymbol{k}) = {\delta _{ll^{\prime}}}$ and ${b_s} = {a_s}\exp ( - \textrm{i}s\omega t)$ to evaluate the averages in (2.33),

(2.36)\begin{gather}{e_K} = \frac{1}{2}\sum\limits_{s,s^{\prime} = 0, \pm 1} {{A_{ss^{\prime}}}(u_s^{{{(1)}^\ast }}u_{s^{\prime}}^{(1)} + u_s^{{{(2)}^\ast }}u_{s^{\prime}}^{(2)})} \exp (\textrm{i}(s - s^{\prime})\omega t),\end{gather}

(2.37)\begin{gather}{e_P} = \frac{1}{2}\sum\limits_{s,s^{\prime} = 0, \pm 1} {{A_{ss^{\prime}}}{\eta ^{{{(s)}^\ast }}}{\eta ^{(s^{\prime})}}\exp (\textrm{i}(s - s^{\prime})\omega t)} .\end{gather}

Note that (2.16) implies

(2.38)\begin{equation}e = {e_K} + {e_P} = \frac{1}{2}\sum\limits_{s = 0, \pm 1} {{A_{ss}}} ,\end{equation}

for the total energy density. Thus, ${A_{ss^{\prime}}}$ determines the energy densities in spectral space, but it contains considerably more statistical information than that: any second-order, two-point moment involving ${u_i}$ and $\eta$ can be obtained knowing ${A_{ss^{\prime}}}$.

The diagonal elements of ${A_{ss^{\prime}}}$ can be interpreted as the energy densities of the different modes and will often be referred to as spectra, whereas the off-diagonal ones represent correlations between modes. As usual, it can be shown that the nonlinear term in (2.26) conserves the total energy, which decays due to visco-diffusion according to

(2.39)\begin{equation}\frac{\textrm{d}}{{\textrm{d}t}}\int {e(\boldsymbol{k})\,{\textrm{d}^3}\boldsymbol{k}} ={-} \textrm{Re} \biggl( {\int {\sum\limits_{s,\hat{s} = 0, \pm 1} {D_{s\hat{s}}^\ast (\boldsymbol{k})\exp (\textrm{i}(\hat{s} - s)\omega (\boldsymbol{k})t){A_{\hat{s}s}}(\boldsymbol{k})} \,{\textrm{d}^3}\boldsymbol{k}} } \biggr).\end{equation}

Finally, a spherically averaged spectrum, $E(k)$, can be defined by averaging $e(\boldsymbol{k})$ over the sphere $|\boldsymbol{k}|= k$, then multiplying by $4{\rm \pi} {k^2}$. When integrated over k, $E(k)$ gives the total energy. It represents the energy distribution in spectral space, without regard for anisotropy. Obviously, $E(k)$ contains less information than $e(\boldsymbol{k})$.

The energy density in spectral space, given by (2.38), can be split into wave and NP contributions as $e = {e_W} + {e_{NP}}$, where

(2.40a,b)\begin{equation}{e_W} = \frac{1}{2}\sum\limits_{s ={\pm} 1} {{A_{ss}}} ,\quad {e_{NP}} = \frac{1}{2}{A_{00}}.\end{equation}

The spherically averaged spectrum can also be split into wave and NP contributions as $E = {E_W} + {E_{NP}}$, where ${E_W}$ and ${E_{NP}}$ are obtained from ${e_W}$ and ${e_{NP}}$ in the same way that E follows from e.

3.1. Evolution of the NP component

Applying (2.26) with $s = 0$,

(3.1)\begin{equation}\frac{{\partial {a_0}}}{{\partial t}} ={-} \textrm{i}{k_j}(v_i^{{{(0)}^\ast }}\widetilde {{u_i}{u_j}} + {\eta ^{{{(0)}^\ast }}}\widetilde {\eta {u_j}}) - \sum\limits_{\hat{s} = 0, \pm 1} {{D_{0\hat{s}}}\exp ( - \textrm{i}\hat{s}\omega (\boldsymbol{k})t){a_{\hat{s}}}} .\end{equation}

Using the decomposition into wave and NP components, (2.21),

(3.2)\begin{gather}{u_i}{u_j} = \underbrace{{u_i^{NP}u_j^{NP}}}_{{NP\textrm{-}NP}} + \underbrace{{u_i^Wu_j^W}}_{{wave \textrm{-} wave}} + \underbrace{{u_i^Wu_j^{NP} + u_i^{NP}u_j^W}}_{{wave \textrm{-} NP}},\end{gather}

(3.3)\begin{gather}\eta {u_j} = \underbrace{{{\eta ^{NP}}u_j^{NP}}}_{{NP \textrm{-} NP}} + \underbrace{{{\eta ^W}u_j^W}}_{{wave \textrm{-} wave}} + \underbrace{{{\eta ^W}u_j^{NP} + {\eta ^{NP}}u_j^W}}_{{wave \textrm{-} NP}}.\end{gather}

Thus, there are three contributions to the nonlinear term in (3.1): NP-NP, wave-wave and wave-NP.

The right-hand side of (3.1) represents nonlinearity and visco-diffusion. As discussed above, both are negligible over time spans of $O(1)$, but can have significant cumulative effects over the longer time scales considered here. The wave component has oscillations of periods $O(1)$, whereas the NP one is steady on such time scales. Thus, the wave-NP contributions to (3.2) and (3.3) are oscillatory and their cumulative effect following evolution according to (3.1) remains small and is neglected. Furthermore, the $\hat{s} ={\pm} 1$ contributions to the visco-diffusive term are also oscillatory and hence negligible. Thus, we obtain the approximation

(3.4)\begin{equation}\frac{{\partial {a_0}}}{{\partial t}} + {D_{00}}{a_0} ={-} \textrm{i}{k_j}(v_i^{{{(0)}^\ast }}\widetilde {u_i^{NP}u_j^{NP}} + {\eta ^{{{(0)}^\ast }}}\widetilde {{\eta ^{NP}}u_j^{NP}}) + f,\end{equation}

where the first term on the right-hand side expresses nonlinear interactions between NP modes and

(3.5)\begin{equation}f(\boldsymbol{k},t) ={-} \textrm{i}{k_j}(v_i^{{{(0)}^\ast }}\widetilde {u_i^Wu_j^W} + {\eta ^{{{(0)}^\ast }}}\widetilde {{\eta ^W}u_j^W})\end{equation}

represents forcing of the NP mode by the wave component.

Following the procedure which led to (2.30), but without the NP contributions,

(3.6)\begin{equation}f(\boldsymbol{k},t) = \sum\limits_{{s_p},{s_q} ={\pm} 1} {\int {M_{0{s_p}{s_q}}^\ast (\boldsymbol{k},\boldsymbol{p})a_{{s_p}}^\ast (\kern1pt \boldsymbol{p})a_{{s_q}}^\ast ( - \boldsymbol{k} - \boldsymbol{p})\exp (\textrm{i}{F_{0{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p})t)\,{\textrm{d}^3}\boldsymbol{p}} } .\end{equation}

Using (2.29), ${F_{0{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = {s_p}\omega (\kern1pt \boldsymbol{p}) + {s_q}\omega ( - \boldsymbol{k} - \boldsymbol{p})$. Since $\omega > 0$, the exponential in (3.6) is oscillatory when ${s_p} = {s_q}$ and such terms are therefore neglected. Assuming $|\varOmega - N|t$ is large, if ${s_p} ={-} {s_q}$, the exponential is oscillatory with period $O(1)$ for $\boldsymbol{p}$ away from the surface $\omega (\kern1pt \boldsymbol{p}) = \omega ( - \boldsymbol{k} - \boldsymbol{p})$. Neglecting such $\boldsymbol{p}$, we focus on $\boldsymbol{p}$ close to the surface, where the exponential oscillates slowly, potentially allowing significant cumulative effects at long times. However, although we have been unable to show it analytically, numerical calculations with different values of $\beta \ne 1$, $\boldsymbol{k}$ and $\boldsymbol{p}$ show that ${M_{0,s, - s}}(\boldsymbol{k},\boldsymbol{p})$ is zero to IEEE double precision when $s \ne 0$ and $\omega (\kern1pt \boldsymbol{p}) = \omega ( - \boldsymbol{k} - \boldsymbol{p})$. Assuming this result, which is a priori far from obvious, is exactly true, slow oscillations are suppressed and we suppose the forcing term in (3.4) is negligible.

Recall from (2.40) that the energy density in spectral space can be split into wave and NP contributions and that the NP contribution is

(3.7)\begin{equation}{e_{NP}} = {\textstyle{1 \over 2}}{A_{00}}.\end{equation}

Appendix D (Supplementary material is available at https://doi.org/10.1017/jfm.2023.1046) shows that the first term on the right-hand side of (3.4), representing nonlinear interactions between NP modes, conserves the total energy of the NP component, which, neglecting the forcing term, decays due to visco-diffusive dissipation according to

(3.8)\begin{equation}\frac{\textrm{d}}{{\textrm{d}t}}\int {{e_{NP}}(\boldsymbol{k})\,{\textrm{d}^3}\boldsymbol{k}} ={-} 2\int {{D_{00}}(\boldsymbol{k}){e_{NP}}(\boldsymbol{k})\,{\textrm{d}^3}\boldsymbol{k}} .\end{equation}

Let ${\varepsilon _{NP}}$ be a small parameter measuring the amplitude of the NP component, thus $\hat{u}_i^{NP} = u_i^{NP}/{\varepsilon _{NP}}$ and ${\hat{\eta }^{NP}} = {\eta ^{NP}}/{\varepsilon _{NP}}$ are $O(1)$. Equations (2.24) and (2.25) with ${b_0} = {a_0}$ give

(3.9)\begin{gather}\hat{u}_i^{NP}(\boldsymbol{x}) = \int {{{\hat{a}}_0}(\boldsymbol{k})v_i^{(0)}(\boldsymbol{k})\exp (\textrm{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{k}} ,\end{gather}

(3.10)\begin{gather}{\hat{\eta }^{NP}}(\boldsymbol{x}) = \int {{{\hat{a}}_0}(\boldsymbol{k}){\eta ^{(0)}}(\boldsymbol{k})\exp (\textrm{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,{\textrm{d}^3}\boldsymbol{k}} ,\end{gather}

where ${\hat{a}_0} = {a_0}/{\varepsilon _{NP}}$. Without the forcing term, (3.4) yields

(3.11)\begin{equation}\frac{{\partial {{\hat{a}}_0}}}{{\partial \hat{t}}} + {\hat{D}_{00}}{\hat{a}_0} ={-} \textrm{i}{k_j}(v_i^{{{(0)}^\ast }}\widetilde {\hat{u}_i^{NP}\hat{u}_j^{NP}} + {\eta ^{{{(0)}^\ast }}}\widetilde {{{\hat{\eta }}^{NP}}\hat{u}_j^{NP}}),\end{equation}

where $\hat{t} = {\varepsilon _{NP}}t$ is time, scaled appropriately for evolution of the NP component, and ${\hat{D}_{00}} = {D_{00}}/{\varepsilon _{NP}}$. Thus, the NP component evolves according to (3.9)–(3.11). Assuming visco-diffusion is sufficiently small that it does not kill the turbulence before nonlinearity intervenes, the time scale for NP evolution is $O(\varepsilon _{NP}^{ - 1})$. This time scale is generally distinct from that of $O({\varepsilon ^{ - 2}})$, which, as we will see later, characterises the wave component according to wave-turbulence theory.

Before going further, we should discuss the close relationship between the theory of the NP component described here and quasi-geostrophic (QG) theory (see Pedlosky Reference Pedlosky1987), which is one of the cornerstones in the study of atmospheric and oceanic flows since its development by Charney (Reference Charney1948, Reference Charney1971). NP modes represent geostrophic flows (i.e. the Coriolis and pressure-gradient terms in the horizontal momentum equation balance) and the NP component of the real flow can be regarded as its mathematical projection onto such idealised flows. With this in mind, it can be shown that, for the present problem, (3.9)–(3.11) are equivalent to the three-dimensional QG approximation, thus providing support for the present approach.

Embid & Madja (Reference Embid and Madja1998) give asymptotic analysis of (2.1)–(2.3) for small Froude number and consider two cases. In the first, the Rossby number is small, and they conclude that the flow consists of oscillatory waves and a slowly varying component which evolves according to the quasi-geostrophic approximation, in agreement with the present results. In the second case, the Rossby number is of order one, corresponding to a small value of $\beta$. This makes $N\sim 1$ and $\varOmega \sim \beta$, thus the wave frequency, (2.9), is small for small ${\theta _{\boldsymbol{k}}}$, i.e. waves having wave vectors near the axis ${k_ \bot } = 0$ are slowly varying. Given such modes, one might question the neglect of the wave-NP contribution to (3.1), which was based on oscillatory waves and led to (3.4). For Rossby numbers of order one, Embid and Majda supposed the slow, horizontal component of the flow to be the sum of two parts (see their (3.31)), one of which is the so-called VSHF (vertically sheared horizontal flow), which is independent of ${x_1}$ and ${x_2}$, and can be considered as a combination of wave modes with ${k_ \bot } = 0$. They derived evolution equations for both parts of the flow (see their (3.33) and (3.34)) and found that the VSHF component entered into the equation for the other part. This corresponds to coupling of the wave and NP components, and neglect of the wave-NP contribution to (3.1) does not capture such coupling. However, the assumed form of the flow places significant wave energy precisely on the axis, ${k_ \bot } = 0$, which is not the situation we have in mind. Instead, we envisage a continuous distribution of energy over wave vectors. In that case, we believe that the neglect of the wave-NP contribution to (3.1) is justified at small $\beta$ by the smallness of the region in ${\theta _{\boldsymbol{k}}}$ for which the wave frequency is small. Whether or not the flow will eventually evolve towards a state close to that assumed by Embid and Majda is unclear.

Equation (3.11) is integrated numerically. Equations (3.9) and (3.10) yield $\hat{u}_i^{NP}(\boldsymbol{x})$ and ${\hat{\eta }^{NP}}(\boldsymbol{x})$ as inverse Fourier transforms. Then, $\hat{u}_i^{NP}\hat{u}_j^{NP}$ and ${\hat{\eta }^{NP}}\hat{u}_j^{NP}$ can be calculated in physical space, while forward transformation gives $\widetilde {\hat{u}_i^{NP}\hat{u}_j^{NP}}$ and $\widetilde {{{\hat{\eta }}^{NP}}\hat{u}_j^{NP}}$ for use in (3.11). This resembles classical DNS of homogeneous turbulence and the numerical method employed here is based on that approach and is described in Appendix F (Supplementary material). However, as noted in the introduction, it differs significantly from classical DNS because of the projection onto the NP modes at each time step, a projection which is implicit in (3.11).

3.2. Evolution of the wave component

Terms in the sum of (2.26) with $\hat{s} \ne s$ are oscillatory and hence negligible for the long-time evolution of ${a_s}$. Adopting this approximation and using (2.30), (2.26) becomes

(3.12) \begin{align} & \dfrac{{\partial {a_s}(\boldsymbol{k})}}{{\partial t}} + {D_{ss}}(\boldsymbol{k}){a_s}(\boldsymbol{k})\nonumber\\ & \quad = \sum\limits_{{s_p},{s_q} = 0, \pm 1} {\int {M_{s{s_p}{s_q}}^\ast (\boldsymbol{k},\boldsymbol{p})a_{{s_p}}^\ast (\kern1pt \boldsymbol{p})a_{{s_q}}^\ast ( - \boldsymbol{k} - \boldsymbol{p})\exp (\textrm{i}{F_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p})t)\,{\textrm{d}^3}\boldsymbol{p}} } . \end{align}

Appendix B describes wave-turbulence analysis based on (3.12). It is assumed that $|\varOmega - N|t$ is large. Furthermore, to obtain closed equations for the wave component, it is supposed that the amplitude of the NP component is small compared with that of the wave component. The result is (B52) for evolution of the wave elements (${A_{ss^{\prime}}}$, $s,s^{\prime} \ne 0$) of the spectral matrix. Since the wave-component spectral energy density is given by the first of (2.40), the most interesting application of (B52) is $s^{\prime} = s$, hence (B53), which involves Cauchy principal-value integrals. However, somewhat remarkably, it is found that the sum of such contributions is zero, as is the sum of ${s_q} = 0$ contributions to the first term on the right-hand side of (B53). The final result is

(3.13) \begin{align} & \dfrac{{\partial {{\hat{A}}_{ss}}(\boldsymbol{k})}}{{\partial T}} + 2\hat{D}(\boldsymbol{k}){{\hat{A}}_{ss}}(\boldsymbol{k}) = 2{\rm \pi} \sum\limits_{{s_p},{s_q} ={\pm} 1} {\int_{{S_{s{s_p}{s_q}}}(\boldsymbol{k})} {\dfrac{{{{\hat{A}}_{{s_p}{s_p}}}(\kern1pt \boldsymbol{p})}}{{|{s_p}{\boldsymbol{c}_g}(\kern1pt \boldsymbol{p}) + {s_q}{\boldsymbol{c}_g}(\boldsymbol{k} + \boldsymbol{p})|}}} } \nonumber\\ & \quad\times ({\lambda _{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}){{\hat{A}}_{{s_q}{s_q}}}( - \boldsymbol{k} -\kern1pt\boldsymbol{p}) + \textrm{Re} ({\zeta _{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p})){{\hat{A}}_{ss}}(\boldsymbol{k}))\,{\textrm{d}^2}\boldsymbol{p}, \end{align}

for $s \ne 0$, where $T = {\varepsilon ^2}t$ is time, scaled appropriately for evolution of the wave component, ${\hat{A}_{ss}} = {A_{ss}}/{\varepsilon ^2}$ are the $O(1)$ scaled wave spectra, $\hat{D} = D/{\varepsilon ^2}$,

(3.14)\begin{equation}D(\boldsymbol{k}) = {k^2}\frac{{{D_u}({N^2}k_ \bot ^2 + 2{\varOmega ^2}k_\parallel ^2) + {D_\eta }{N^2}k_ \bot ^2}}{{2({N^2}k_ \bot ^2 + {\varOmega ^2}k_\parallel ^2)}}\end{equation}

is the damping coefficient of wave modes (i.e. ${D_{ss}}(\boldsymbol{k}) = D(\boldsymbol{k})$ for $s \ne 0$),

(3.15)\begin{equation}{\boldsymbol{c}_g}(\boldsymbol{k}) = {\nabla _{\boldsymbol{k}}}\omega (\boldsymbol{k})\end{equation}

gives the group velocity of the wave modes as $s{\boldsymbol{c}_g}$,

(3.16)\begin{gather}{\lambda _{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = 2|{M_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}){|^2},\end{gather}

(3.17)\begin{gather}{\zeta _{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = 4{M_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p})M_{{s_q}s{s_p}}^\ast ( - \boldsymbol{k} - \boldsymbol{p},\boldsymbol{k}),\end{gather}

and ${S_{s{s_p}{s_q}}}(\boldsymbol{k})$ is the surface in $\boldsymbol{p}$-space defined by ${F_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = 0$. As is apparent from (2.29), the surface ${S_{s{s_p}{s_q}}}(\boldsymbol{k})$ is such that $s\omega (\boldsymbol{k}) + {s_p}\omega (\kern1pt \boldsymbol{p}) + {s_q}\omega ( - \boldsymbol{k} - \boldsymbol{p}) = 0$. This condition represents triadic resonances, triadic since it involves three wave vectors, $\boldsymbol{k}$, $\boldsymbol{p}$ and $- \boldsymbol{k} - \boldsymbol{p}$ which sum to zero, and resonant because it says that the sum of modal frequencies is zero. For this reason, ${S_{s{s_p}{s_q}}}(\boldsymbol{k})$ will be referred to as the resonant surface. If there are no such resonances for given $\boldsymbol{k}$, s, ${s_p}$ and ${s_q}$ (i.e. ${F_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = 0$ has no solutions in $\boldsymbol{p}$-space), ${S_{s{s_p}{s_q}}}(\boldsymbol{k})$ is the empty set and the surface integral in (3.13) should be interpreted as zero.

It is perhaps interesting to discuss (3.13) in the context of EDQNM. Eddy damping does not appear (for the case of pure rotation, it was shown in [B] that it vanishes in the limit of weak turbulence). Quasi-normality is a consequence of wave-turbulence analysis and is apparent in (3.13) via the products of spectra on the right-hand side which represent fourth-order spectral moments. The fact that, according to (3.13), the spectrum at time t evolves according to its values at the same instant justifies the final, Markovian, closure hypothesis of EDQNM. An important difference between wave-turbulence theory and EDQNM is that the nonlinear term in (3.13) is a surface integral over resonant triads, whereas it is a volume integral over all triads, including non-resonant ones, according to EDQNM. For the case of pure rotation, it was shown in [B] that the volume integral of EDQNM is dominated by resonant triads in the limit of weak turbulence, hence the wave-turbulence result is approached.

Appendix E (Supplementary material) analyses the energetics of the wave component based on (3.13). It is shown that the right-hand side of (3.13) conserves the total wave-component energy, which evolves according to

(3.18)\begin{equation}\frac{\textrm{d}}{{\textrm{d}t}}\int {{e_W}(\boldsymbol{k})\,{\textrm{d}^3}\boldsymbol{k}} ={-} 2\int {D(\boldsymbol{k}){e_W}(\boldsymbol{k})\,{\textrm{d}^3}\boldsymbol{k}} .\end{equation}

Taking the sum of (3.8) and (3.18),

(3.19)\begin{equation}\frac{\textrm{d}}{{\textrm{d}t}}\int {e(\boldsymbol{k})\,{\textrm{d}^3}\boldsymbol{k}} ={-} 2\int {({D_{00}}(\boldsymbol{k}){e_{NP}}(\boldsymbol{k}) + D(\boldsymbol{k}){e_W}(\boldsymbol{k}))\,{\textrm{d}^3}\boldsymbol{k}} .\end{equation}

This result may be compared with the exact energy equation (2.39). In the weak-turbulence limit considered here, terms in the sum of (2.39) with $\hat{s} \ne s$ are oscillatory and hence negligible. Dropping these terms gives (3.19). Thus, as regards the total energy, the present approximations agree with the weak-turbulence limit of the exact result.

Appendix C concerns the existence of solutions of ${F_{s{s_p}{s_q}}}(\boldsymbol{k},\boldsymbol{p}) = 0$, where s, ${s_p}$ and ${s_q}$ take one of the values ${\pm} 1$. It is shown that the resonant surface does not exist if $s = {s_p} = {s_q}$ or $1/2 \le \beta \le 2$ (a result in agreement with Smith & Waleffe (Reference Smith and Waleffe2002), § 6.1). Thus, when $1/2 \le \beta \le 2$, the right-hand side of (3.13) is zero according to the present theory. Nonlinear interactions between wave modes are then absent at the order to which we are working, suggesting the need to go to higher order. However, that lies beyond the scope of the present work. When $\beta < 1/2$ or $\beta > 2$, the resonant surface exists for ${s_p} = {s_q} \ne s$ provided

(3.20)\begin{equation}\omega (\boldsymbol{k}) > 2\min(\varOmega ,N),\end{equation}

and for ${s_p} ={-} {s_q}$ when

(3.21)\begin{equation}\omega (\boldsymbol{k}) < |\varOmega - N|.\end{equation}

If $\beta < 1/3$ or $\beta > 3$, one or other of (3.20) or (3.21) is satisfied for all $\boldsymbol{k}$, hence resonant surfaces can always be found for such values of $\beta$. However, when $1/3 < \beta < 1/2$ or $2 < \beta < 3$, there are directions of $\boldsymbol{k}$ for which neither (3.20) nor (3.21) hold, hence no resonant surfaces. Such wave vectors are decoupled from all others according to (3.13), which predicts that these modes simply decay under visco-diffusive dissipation. Thus, spectral space can be divided into two regions. In the first (referred to as A), either (3.20) or (3.21) applies and there is nonlinear coupling between wave modes within that region. In the second (B), such coupling is absent. This distinction expresses different dynamics for the two types of modes, but is only significant for $1/3 < \beta < 1/2$ or $2 < \beta < 3$, otherwise all modes are of type A. Note that the extent of the region in which inter-mode coupling occurs shrinks down to zero as either $\beta \nearrow 1/2$ or $\beta \searrow 2$, which is the way in which nonlinearity of the order considered here disappears as those boundaries are approached.

Consider a wave vector for which there is one or more resonant surfaces and suppose visco-diffusive dissipation small enough $(D \le O({\varepsilon ^2}))$ that it does not kill the wave component before nonlinearity intervenes, the time scale for evolution of the wave component implied by (3.13) is $t\sim O({\varepsilon ^{ - 2}})$, the usual one for wave turbulence. If there is no resonant surface, we expect the nonlinear evolution time to be longer. As noted earlier, the time scales for evolution of the wave and NP components are generally distinct.

Appendix G (Supplementary material) describes the numerical method used to solve (3.13) for the statistically axisymmetric case. The procedure is essentially that of Bellet (Reference Bellet2003), who used the wave-turbulence equations for the case of pure rotation. For the applications described in the next section, the initial ${A_{ss}}(\boldsymbol{k})$ are not just axisymmetric, but also symmetric under reflection in the plane ${k_\parallel } = 0$. This corresponds to statistical symmetry of the underlying flow under reflection in any plane perpendicular to the rotation axis, a symmetry which, like axisymmetry, is preserved by time evolution according to the governing equations. As a result, ${A_{ss}}(\boldsymbol{k})$ remains reflection symmetric. Together, ${A_{ss}}( - \boldsymbol{k}) = {A_{ - s, - s}}(\boldsymbol{k})$, axisymmetry and reflection symmetry imply ${A_{ss}}(\boldsymbol{k}) = {A_{ - s, - s}}(\boldsymbol{k})$, hence ${A_{11}}(\boldsymbol{k}) = {A_{ - 1, - 1}}(\boldsymbol{k}) = {e_W}(\boldsymbol{k})$ according to the first of (2.40).

Although, for the sake of simplicity in the analysis, the case $N = 0$ (pure rotation) was earlier excluded, it may be interesting to consider that case because, to our knowledge, it is the only one for which wave-turbulence theory has previously been applied to the present problem. According to (2.10), (2.12), (2.14) and (2.15), if $N = 0$, the velocity field is carried solely by the wave modes, while the scalar field, $\eta$, consists of NP modes alone, i.e. there is a precise correspondence between velocity and wave modes, and $\eta$ and NP ones. As discussed towards the beginning of § 2, when $N = 0$, the velocity field decouples from the scalar $\eta$, hence a corresponding decoupling of the wave modes from the NP ones. Under these circumstances, it is natural to leave the NP component to one side and consider the wave component alone. Furthermore, since the wave component is exactly decoupled from the NP one, the wave-turbulence equations are closed without the need for the assumption ${\varepsilon _{NP}} \ll \varepsilon$ made in the general case.

Focusing on the wave component, as $N \to 0$, (3.13) approaches the wave-turbulence equations used by Galtier (Reference Galtier2003) and [B] for the case of pure rotation. Numerical integration of those equations by [B] indicate that the wave energy density, ${e_W}(\boldsymbol{k})$, develops an infinite, but integrable, singularity at the plane ${k_\parallel } = 0$. This is the point of view of wave-turbulence theory, i.e. $\varepsilon = 0$. However, for small, but non-zero $\varepsilon$, there is energy transfer towards the plane ${k_\parallel } = 0$ (see Waleffe Reference Waleffe1991), leading to large, but bounded, values of ${e_W}$ for small ${k_\parallel }$. Thus, we see how the wave-turbulence result is approached in the limit $\varepsilon \to 0$. That the singularity of ${e_W}(\boldsymbol{k})$ is integrable in that limit is significant because it implies that the contribution of small ${k_\parallel }$ to the total energy is small, rather than dominant, as suggested by some authors (e.g. Hossain Reference Hossain1994).

4.1. DNS results for the NP component

Let us begin with figure 1, which gives log-log plots of the scaled NP spectra, ${\hat{A}_{00}} = {A_{00}}(k,{\theta _{\boldsymbol{k}}})/\varepsilon _{NP}^2$, where ${\theta _{\boldsymbol{k}}}$ is the angle between $\boldsymbol{k}$ and $\boldsymbol{e}$, as a function of k for different values of $\beta$ and ${\theta _{\boldsymbol{k}}}$. Perhaps the most striking feature is the wiggles in the curves. These are due to two traditional limitations of DNS. First, the results are based on a single run of DNS, not the many runs which would be needed to accurately compute the ensemble average required by the definition of ${A_{00}}$, runs which would render the calculations too computationally expensive. Thus, there are random fluctuations about the smooth curves which would presumably result from averaging. Second, DNS discretises the wave vector $\boldsymbol{k}$, leading to the considerable jumping around apparent for smaller values of k. These problems decrease in importance as k increases and we focus on k having larger values.

Figure 1. Log-log plots of the scaled NP spectra for $\hat{t} = 8$ and (a) $\beta = 0$, (b) $\beta = 0.25$, (c) $\beta = 0.5$, (d) $\beta = 0.7$, (e) $\beta = 1.4$, (f) $\beta = 2$, (g) $\beta = 4$ and (h) $\beta = 10$. Each figure shows curves for eight equally spaced values of ${\theta _{\boldsymbol{k}}}$ between ${\rm \pi} /32$ and $15{\rm \pi} /32$. The dashed line corresponds to the power law ${\hat{A}_{00}} \propto {k^{ - 5}}$. If required, N and $\varOmega$ follow from $\beta$ using $N = {({\beta ^2} + 1)^{ - 1/2}}$ and $\varOmega = \beta {({\beta ^2} + 1)^{ - 1/2}}$.

Figures 1(a)–1(h) represent the scaled, NP spectra for $\hat{t} = 8$ (recall that $\hat{t} = {\varepsilon _{NP}}t$ is the scaled time for the NP component) for the different values of $\beta$ given in the figure caption. Each figure shows results for eight equally spaced values of ${\theta _{\boldsymbol{k}}}$ between ${\rm \pi} /32$ and $15{\rm \pi} /32$. The dashed line represents the power law ${k^{ - 5}}$ and has the same location in all figures. For those readers more used to plots of spherically averaged spectra, we recall that the ${k^{ - 5}}$ for the spectral density shown here corresponds to behaviour like ${k^{ - 3}}$ of ${E_{NP}}(k)$. The symmetry ${A_{00}}(k,{\theta _{\boldsymbol{k}}}) = {A_{00}}(k,{\rm \pi} - {\theta _{\boldsymbol{k}}})$ allows the results for values of ${\theta _{\boldsymbol{k}}}$ above ${\rm \pi} /2$ to be deduced from those given.

In most cases, an inertial range with power-law exponent close to $- 5$ is apparent. In particular, this applies for the values of $\beta$ represented by figure 1(b–g) for all ${\theta _{\boldsymbol{k}}}$ and for figure 1(a) (pure stratification), apart from ${\theta _{\boldsymbol{k}}} = {\rm \pi}/32$. Such generic near constancy of the exponent is remarkable. However, for $\beta = 10$ (near to pure rotation), the situation is rather different and harder to interpret, since the power law observed in the other cases is no longer apparent. Indeed, it is hard to discern large-$k$ power laws from these results and the limit of small N merits further study of the NP component. Furthermore, the lack of a ${k^{ - 5}}$ range in figure 1(a) for ${\theta _{\boldsymbol{k}}} = {\rm \pi}/32$ suggests that the case of small $\varOmega$ and small ${\theta _{\boldsymbol{k}}}$ is also worth further investigation.

In addition to the inertial-range exponent, figure 1 gives information on anisotropy of the small scales. For the $\beta$ of figures 1(a)–1(d), the large-$k$ spectra decrease as ${\theta _{\boldsymbol{k}}}$ increases, thus the small-scale NP energy is concentrated towards the poles, ${\theta _{\boldsymbol{k}}} = 0$ and ${\theta _{\boldsymbol{k}}} = {\rm \pi}$. When $\beta = 1.4$ (figure 1e), the energy distribution is almost isotropic, whereas for larger $\beta$, the energy is concentrated towards the equator, ${\theta _{\boldsymbol{k}}} = {\rm \pi}/2$. This is in agreement with the general consensus that dominant stratification favours energy transfer towards the poles, whereas dominant rotation sends it towards the equator. However, it should again be recalled that we are only considering one component of the flow here.

Figure 2 provides another view of anisotropy in which contours of ${\hat{A}_{00}}$ are plotted in the ${k_ \bot }-{k_\parallel }$ plane for the cases $\beta = 0.1$ and $\beta = 10$. We once again see that close to pure stratification, the NP energy density tends towards the axial direction, whereas it concentrates near the equator as the case of pure rotation is approached.

Figure 2. Contour plots of ${\hat{A}_{00}}$ as functions of ${k_ \bot }$ and ${k_\parallel }$ for $\hat{t} = 8$ and (a) $\beta = 0.1$, (b) $\beta = 10$. There are ten contours, whose heights are logarithmically spaced from ${10^{ - 8}}$ to $1$.

As regards time evolution, figure 3 shows a representative example $(\beta = 0.5)$ of the spherically averaged NP energy spectrum ${E_{NP}}(k)$, defined earlier, at $17$ equally spaced times from $\hat{t} = 0$ to $\hat{t} = 4$. It will be seen that, at early times, there is energy transfer from large to small scales, thus forming the inertial and dissipative ranges. At later times, the spectral peak moves slowly towards smaller k. The inertial range has approximate power-law behaviour close to ${k^{ - 3}}$, corresponding to the ${k^{ - 5}}$ of figure 1(c) and agreeing with the spectral power law proposed by Charney (Reference Charney1971).

Figure 3. Log-log plots of the scaled, spherically averaged NP spectra for $\beta = 0.5$ and times $\hat{t} = 0,0.25,0.5, \ldots ,3.75,4$. The arrow indicates the direction in which the large-$k$ spectral curves move with increasing time in the early stages.

The total NP energy decreases with time by an amount which depends on the value of $\beta$. The decrease from $\hat{t} = 0$ to $\hat{t} = 8$ is $0.7\,\%$ for $\beta = 0.7$ and $\beta = 1.4$, $1\,\%$ when $\beta = 0.5$ and $\beta = 2$, $3\,\%$ for $\beta = 0.25$ and $\beta = 4$, $11\,\%$ when $\beta = 0.1$, $13\,\%$ for $\beta = \infty$, $16\,\%$ when $\beta = 10$ and $33\,\%$ for $\beta = 0$. Considering that these small to moderate decreases correspond to a time, $\hat{t} = 8$, much greater than the establishment time of the dissipative range, it appears that there is no energy NP energy cascade. This conjecture is reinforced by log-log plots (not shown here) of the total NP energy as a function of time which give no indication of the power laws which might be expected if there were an energy cascade. The absence of a cascade is in agreement with Charney (Reference Charney1971).

In conclusion, our results for the NP component are consistent with the theoretical predictions of Charney (Reference Charney1971), both in terms of power laws and the lack of a cascade. They also provide detailed information on spectral anisotropy and its variation with $\beta$. It has been suggested (see e.g. Herbert, Pouquet & Marino Reference Herbert, Pouquet and Marino2014) that the absence of a cascade of energy to large k is because the cascade is towards small k (an inverse cascade). This is related to the existence of two QG invariants, namely energy and potential enstrophy, in the absence of visco-diffusion. The theoretical basis of the inverse cascade for three-dimensional QG turbulence is analogous to that of the classical two-dimensional case. Several DNS studies (e.g. Marino et al. Reference Marino, Mininni, Rosenberg and Pouquet2013) have claimed that the inverse cascade is strongest in the range $1/2 < \beta < 2$, though why this should be, given that this condition appears to refer to the wave component, is unclear.