1. Introduction
Fluid turbulence, where nonlinear interactions occur over a wide range of spatial and temporal scales, plays an important role in engineering, geophysical, astrophysical and even biological fluid mechanics. Much turbulence research focuses on the idealised case of homogeneous and isotropic turbulence, despite the canonical situation involving important inhomogeneities and anisotropies. For example, in geophysical and astrophysical situations, rotation and stratification may play important roles in selecting preferred directions, whereas in other cases mean flows and boundaries often lead to both inhomogeneities and anisotropies. For these cases, it is important to develop a framework that builds in inhomogeneity and anisotropy from the outset and turns this ‘bug’ into a ‘feature’. Such a framework involves constructing equations for the evolution of the statistics of the turbulence; it is important to bear in mind that the presence of anisotropy and inhomogeneity often leads to non-trivial low-order statistics; for example, sustained mean flows that interact strongly with fluctuations. Methods designed for describing the evolution of such flows will perform badly for the homogeneous isotropic case, where mean flows are absent. For a description of the many methods that have been developed and the underlying philosophy of this approach, see the review by Marston & Tobias (Reference Marston and Tobias2023).
In general these methods rely on developing equations governing the evolution of the low-order statistics for the flow (often termed the cumulants). Such an approach often requires a closure approximation, where the higher-order statistics are either neglected completely or written as functions of the low-order cumulants. However, if the system exhibits quasilinear dynamics, then the system of low-order statistical equations closes exactly (Herring Reference Herring1963) and no further approximations are needed; such a system of statistical equations is known as CE2 (representing a cumulant expansion at second order). Recent years has seen many systematic investigations of the validity of the quasilinear approximation (QL) in representing the full nonlinear dynamics (see e.g. Tobias & Marston Reference Tobias and Marston2013; Child et al. Reference Child, Hollerbach, Marston and Tobias2016; Marston, Qi & Tobias Reference Marston, Qi and Tobias2019); it has been determined that there are certain circumstances where the QL approximation breaks down and better approximations are needed.
Recent research has focused on a generalisation of the quasilinear approximation (termed GQL) that generalises the definition of the mean flow to include long-wavelength modes in the zonal direction in addition to the spatial mean (Marston, Chini & Tobias Reference Marston, Chini and Tobias2016). A scale separation into long-wavelength modes that interact fully nonlinearly with themselves, and short-wavelength modes with wavenumbers greater than a cutoff 
$\varLambda$ that interact only quasilinearly with the long wavelength modes, has certain advantages in that it allows energy to be scattered among the turbulent eddies through interaction with the generalised mean flow. This eddy scattering often leads to a more faithful representation of the nonlinear dynamics than QL, where these interactions are forbidden.
GQL, like QL, is a conservative approximation but one that systematically interpolates between QL (when 
$\varLambda = 0$) and fully nonlinear (NL). As 
$\varLambda \rightarrow \infty$ the GQL system consists solely of fully interacting low modes and reverts to NL, albeit not necessarily monotonically. In particular, if the cutoff 
$\varLambda$ is large (but not infinite) it is possible that high modes will be stable and (incorrectly) not have any energy (Hernández, Yang & Hwang Reference Hernández, Yang and Hwang2022).
In this paper we derive the statistical representation of the GQL approximation (which we term GCE2) and describe the utility of this approach by comparing it with CE2 for two model fluid dynamics problems (one deterministic and one stochastic) describing the interaction of mean flows with turbulence in two dimensions.
Recently it has been shown that a rank instability can occur in CE2 leading to a divergence in statistics from QL (Nivarti et al. Reference Nivarti, Kerswell, Marston and Tobias2022; see also Oishi et al. Reference Oishi, Burns, Marston and Tobias2022). The existence of this instability has apparently been missed in all prior work going back to Herring (Reference Herring1963), and the resolution presented here now allows us to properly compare GCE2 to statistics obtained from GQL and compare both to the statistics found in QL and NL simulations.
2. The GQL approximation and its closure
The GQL approximation has been studied in a variety of fluid contexts by a number of different groups. We refer the reader to our review article (Marston & Tobias Reference Marston and Tobias2023) for an introduction to the approximation and its physical interpretation.
2.1. The GQL approximation
We consider a system of nonlinear dynamical equations for a state vector 
$\boldsymbol {q}(\boldsymbol {x},t)$ written as
\begin{equation} \partial_t \boldsymbol{q} = \mathcal{L}[\boldsymbol{q}] + \mathcal{N}[\boldsymbol{q},\boldsymbol{q}], \end{equation}
with 
$\mathcal {L} [{\cdot }]$ a linear and 
$\mathcal {N}[{\cdot },{\cdot }]$ the nonlinear (in this case, quadratic) operator. In order to apply the GQL approximation, the state vector 
$\boldsymbol {q}$ is expanded using a spectral basis along the zonal direction (more generally, the direction exhibiting statistical homogeneity). GQL then proceeds (Marston et al. Reference Marston, Chini and Tobias2016) by applying a low-pass filter (projection operator) with cutoff 
$\varLambda$ in the zonal direction, leading to a generalised Reynolds decomposition of the state vector
\begin{equation} \boldsymbol{q} = \boldsymbol{q}_\ell + \boldsymbol{q}_h, \end{equation}
where the subscripts (
$\ell$) and (
$h$) denote low and high zonal wavenumber modes, respectively. For the two-dimensional Cartesian models considered here,
\begin{equation} \boldsymbol{q}_\ell=\sum_{k={-}\varLambda}^{\varLambda}{\boldsymbol{q}}_{k}(y) \,{\rm e}^{{\rm i} k x},\quad \boldsymbol{q}_h = \boldsymbol{q} - \boldsymbol{q}_\ell.\end{equation} This decomposition obeys the usual rules of orthogonality and idempotence, and can be simplified to the conventional Reynolds decomposition (into mean and fluctuation) simply by setting 
$\varLambda = 0$. Note, however, that the conventional Reynolds decomposition obeys the equality 
$\overline {\bar {\boldsymbol {q}} \boldsymbol {q}} = \overline {\bar {\boldsymbol {q}} \bar {\boldsymbol {q}}}$. Note also that we use the overline to indicate a mean computed as a spatial (zonal) average, but other averages are also possible including ensemble and time averages.. In comparison, 
$({\boldsymbol {q}_\ell \boldsymbol {q}})_\ell \ne (\boldsymbol {q}_\ell \boldsymbol {q}_\ell )_\ell$ for 
$\varLambda > 0$ under the generalised Reynolds decomposition into low and high modes; this inequality is remedied after applying certain interaction rules as follows.
Applying the decomposition (2.2) to the 
$\boldsymbol {q}$ in (2.1), gives rise to various classes of nonlinear terms that correspond to different triadic interactions involving low and high modes. The possible triad interactions can be represented using triadic diagrams in which low modes (zonal wavenumber 
$|m| \leq \varLambda$) and high modes 
$|m| > \varLambda$ are denoted using low-frequency and high-frequency wave edges (Marston et al. Reference Marston, Chini and Tobias2016; Marston & Tobias Reference Marston and Tobias2023). GQL equations retain nonlinear self-interactions giving rise to low modes, as well as quasilinear interactions between low and high modes. Nonlinear self-interactions giving rise to high modes and interactions between low and high modes giving rise to low modes are both dropped (see the diagrams in § 2 of Tobias, Oishi & Marston Reference Tobias, Oishi and Marston2018). This elimination of triad interactions, as in the case of triad decimation by pairs (Kraichnan Reference Kraichnan1985), conserves quadratic invariants such as energy and enstrophy. The following GQL equations are obtained by applying these interaction rules:
$$\begin{gather} \partial_t \boldsymbol{q}_\ell = \mathcal{L}[\boldsymbol{q}_\ell] + \mathcal{N}_\ell[\boldsymbol{q}_\ell,\boldsymbol{q}_\ell] + \mathcal{N}_\ell[\boldsymbol{q}_h,\boldsymbol{q}_h], \end{gather}$$
$$\begin{gather}\partial_t \boldsymbol{q}_h = \mathcal{L}[\boldsymbol{q}_h] + \mathcal{N}_h[\boldsymbol{q}_h,\boldsymbol{q}_\ell] + \mathcal{N}_h[\boldsymbol{q}_\ell,\boldsymbol{q}_h], \end{gather}$$
for the low and high modes. Here, 
$\mathcal {N}_\ell [{\cdot }, {\cdot }]$ is the projection of the nonlinearity onto the low modes and 
$\mathcal {N}_h = \mathcal {N}[{\cdot }, {\cdot }] - \mathcal {N}_\ell [{\cdot }, {\cdot }]$, which denotes the high-mode spectral projection of the nonlinear operator 
$\mathcal {N}[{\cdot }, {\cdot }]$. Furthermore, when 
$\varLambda = 0$, we have 
$\boldsymbol {q}_\ell = \bar {\boldsymbol {q}}$ and 
$\boldsymbol {q}_h = \boldsymbol {q}^\prime$, and the GQL equations reduce to the well-known system of QL equations (Marston et al. Reference Marston, Qi and Tobias2019). Similarly, setting 
$\varLambda$ as the full spectral resolution in the zonal direction, we obtain 
$\boldsymbol {q}_\ell = \boldsymbol {q}$ and 
$\boldsymbol {q}_h = 0$; the GQL equations are identical to the NL equations (2.1) in this limit. Hence, the GQL equations interpolate systematically (though not monotonically) between QL and NL dynamics by varying the zonal spectral cutoff 
$\varLambda$. Crucially, the GQL equations lack high-mode nonlinearities which makes GQL amenable to statistical closure as (2.5) is formally linear in 
$\boldsymbol {q}_h$. We shall make use of this property now.
2.2. Deriving the generalised cumulant expansion at second order (GCE2)
Statistically closed equations termed CE2 have been derived for QL equations using cumulant expansions and other methods (see e.g. Farrell & Ioannou Reference Farrell and Ioannou2007; Marston, Conover & Schneider Reference Marston, Conover and Schneider2008; Farrell & Ioannou Reference Farrell and Ioannou2013; Constantinou Reference Constantinou2015). A closure for QL equations is achieved at second order: the equations for the mean and fluctuation terms are used as a starting point for deriving the corresponding equations for the first two cumulants, 
$\bar {\boldsymbol {q}}$ and 
$\overline {\boldsymbol {q}^\prime \boldsymbol {q}^\prime }$, respectively. This strategy allows for direct statistical simulation (DSS) of low-order statistical quantities that correspond to QL dynamics. In a similar manner, statistical closure for GQL dynamics is also achievable at second order using generalised cumulant expansions which employ the notion of the mean implicit within the generalised Reynolds decomposition of (2.2). We thus define the first two generalised cumulants as
\begin{equation} \left.\begin{gathered} c_1 \equiv \boldsymbol{q}_\ell, \\ c_2 \equiv \boldsymbol{q}_h \boldsymbol{q}_h, \end{gathered}\right\} \end{equation}
following the spectral filter (projection operator) notation used earlier. Closed-form equations for these generalised cumulants, termed GCE2, are obtained by following a similar approach to that adopted for deriving CE2 equations from QL equations (Marston et al. Reference Marston, Qi and Tobias2019). The generalised first cumulant is identical to low modes (analogous to 
$\bar {\boldsymbol {q}}$ in CE2). As a result, the equation governing the evolution of the generalised first cumulant 
$c_1$ is the same as that for the low modes 
$\boldsymbol {q}_\ell$ as given by (2.4):
\begin{equation} \partial_t c_1 = \mathcal{L}[c_1] + \mathcal{N}_\ell[c_1, c_1] + \mathcal{N}_\ell[c_2]. \end{equation}
The generalised second cumulant is a field bilinear in high modes (akin to the Reynolds stress 
$C = \overline {\boldsymbol {q}^\prime \boldsymbol {q}^\prime }$ in CE2). The equation for 
$c_2 = \boldsymbol {q}_h \boldsymbol {q}_h$ is thus obtained by multiplying the high-mode equations (2.5) with high modes and subsequently projecting down to the low modes as follows:
\begin{equation} \left.\begin{gathered} \partial_t (\boldsymbol{q}_h \boldsymbol{q}_h) = \{\boldsymbol{q}_h \partial_t \boldsymbol{q}_h\} = \{\boldsymbol{q}_h \mathcal{L}[\boldsymbol{q}_h]\} + \{\boldsymbol{q}_h \mathcal{N}_h[\boldsymbol{q}_h,\boldsymbol{q}_\ell]\}, \\ \implies \partial_t c_2 = \mathcal{L}[c_2] + \{\mathcal{N}_h[c_2, c_1]\}, \end{gathered}\right\} \end{equation}
where the curly braces 
$\{{\cdot } {\cdot } \}$ denote symmetrisation, i.e. 
$\{a b\} = a b + b a$. Taken together, (2.8) and (2.4) form a closed set of equations in terms of generalised cumulants. This allows the implementation of DSS of GQL dynamics, i.e. solve directly for statistics that interpolate between QL and NL dynamics. We note here that similar equation sets have been derived in other contexts (see e.g. Bakas & Ioannou Reference Bakas and Ioannou2013a,Reference Bakas and Ioannoub; Constantinou Reference Constantinou2015; Constantinou, Farrell & Ioannou Reference Constantinou, Farrell and Ioannou2016), though the relationship with GQL models and the possibility of the statistical models diverging from the QL-type models have not been investigated previously.
A recent study has demonstrated that although one can mathematically arrive at the CE2 system starting from the QL equations, the two frameworks are not statistically equivalent (Nivarti et al. Reference Nivarti, Kerswell, Marston and Tobias2022). This is because, as opposed to QL equations, the CE2 system can sustain solutions wherein the second cumulant has rank 
${>}1$ whereas the second cumulant obtained by zonally averaging QL dynamics is always of unit rank. We anticipate therefore that the GCE2 and GQL systems may not be fully equivalent either. However, we emphasise that the ability of GCE2 to predict a generalised second cumulant with rank greater than unity ought to be considered as a feature not a bug, one that motivates further study.
This paper has two aims: (1) to test the predictions of GCE2 against those of GQL, thereby demonstrating that, even for a cutoff 
$\varLambda = 1$, GCE2 is a significant improvement over CE2; and (2) to probe observed divergences between GCE2 and GQL, relating them to the presence of rank instabilities as in the case of CE2 and QL (Nivarti et al. Reference Nivarti, Kerswell, Marston and Tobias2022). We employ only low-resolution models here with a small number of modes to highlight clearly any differences between GQL and GCE2. We have verified that the qualitative results continue to hold at higher resolutions as well.
3. Numerical implementation
We conduct simulations of a rotating, incompressible fluid on a doubly periodic 
$\beta$-plane. The time evolution of the relative vorticity 
$\zeta \equiv \hat {z}\boldsymbol {\cdot } (\boldsymbol {\nabla } \times {\boldsymbol u})$ is given by
\begin{equation} \partial_t \zeta = \beta \partial_x \psi - \kappa\zeta + \nu \nabla^2 \zeta + J[\psi,\zeta] + F, \end{equation}
where 
${\boldsymbol u}$ is the velocity, 
$J[\psi,\zeta ] = \partial _x\psi \partial _y \zeta - \partial _x \zeta \partial _y \psi$ and the streamfunction 
$\psi \equiv \nabla ^{-2}\zeta$. Gradients of rotation are included via the 
$\beta$ term, ‘bottom friction’ via the 
$\kappa$ term and viscosity via the 
$\nu$ term. The forcing term 
$F$ models energy injection into the system. We adopt two different models of forcing (or driving) in this study: (1) a deterministic steady two-scale Kolmogorov-type forcing as used by Tobias & Marston (Reference Tobias and Marston2017), see (4.1); and (2) a white-in-time stochastic model of forcing adapted from Constantinou et al. (Reference Constantinou, Farrell and Ioannou2016), see (4.2). The latter is a commonly used model that mimics thermal driving.
The spectral solver ZonalFlow.jl (written in Julia (Bezanson et al. Reference Bezanson, Edelman, Karpinski and Shah2017) and made available online Nivarti, Marston & Tobias Reference Nivarti, Marston and Tobias2021) is used to obtain direct numerical solutions of (3.1). Timestepping algorithms are imported from the well-tested ecosystem of the DifferentialEquations.jl (Rackauckas & Nie Reference Rackauckas and Nie2017) package. We use different timestepping methods to simulate flows with different types of driving models. For the deterministically driven Kolmogorov flow, we use the explicit 
$5/4$ Runge–Kutta method of Dormand-Prince with a fixed timestep of 
$\Delta t = 0.001$. For the stochastically driven flow, we use the SRIW1 method of order 
$1.5$. To be consistent, we use the same timestepping algorithm to solve the dynamical equations (QL/GQL/NL) and the statistical equations (CE2/GCE2) for a given flow. Unless specified otherwise, the initial condition consists of random noise of mean power 
$10^{-4}$. We employ a 
$[0,2{\rm \pi} ]^2$ grid with resolution 
$M \times N = 10 \times 10$ for the Kolmogorov flow system, and a 
$[0,2{\rm \pi} ]\times [0,{\rm \pi} ]$ grid with a resolution of 
$12 \times 20$ for the stochastically driven system, where 
$|m| < M$ is the wavenumber in the zonal (
$x$) direction and 
$|n| < N$ is the wavenumber in the non-zonal (
$y$) direction.
We have conducted a number of validation tests to ensure that the equations are implemented accurately in our code. Solutions for equations neglecting the dissipative terms have been tested to confirm conservation of energy and enstrophy is satisfied in such regimes. In addition, the agreement of mean equilibrium solutions for a given energy input with dissipation has also been confirmed across the various equation systems. In the presence of nonlinear terms, our tests confirm that setting the spectral cutoff 
$\varLambda = 0$ in GCE2 reproduces CE2 results and setting 
$\varLambda = M - 1$ reproduces the NL results. In addition to the forcing models outlined above, the aforementioned tests are also conducted for a deterministic driving model with relaxation to a point jet (Marston et al. Reference Marston, Conover and Schneider2008); here, an additional validation test is available as the GCE2/CE2 as well as GQL/QL results converge to NL for small relaxation times. All these tests have been automated, such that the Github workflows server conducts them each time a code change is pushed to the online repository (Nivarti et al. Reference Nivarti, Marston and Tobias2021).
In the following, we present results comparing GCE2 against GQL for a spectral cutoff 
$\varLambda = 1$. For each forcing model and specified set of parameters, we compare these results against predictions of the NL system as well as of the CE2 and QL equations. Note that, we use overlines to indicate spatial (zonal) averages (not ensemble or time averages) on QL solutions as in the comparisons against CE2. When time averages are used (such as to facilitate comparison of Hövmöller plots), we explicitly state this along with specifying the averaging window in the accompanying text.
4. Results
4.1. Two-scale Kolmogorov forcing
Figure 1 shows the final time (
$t = 1000$ days) resolved vorticity solution for a flow driven by a deterministic two-scale Kolmogorov forcing. We set 
$F$ in (4.1) to be
\begin{equation} F(y) ={-}\cos(y) -8\cos(4y) \end{equation}
and the parameters 
$\nu = 0.02$ and 
$\beta = \kappa = 0$ as in (3.1). This system is known to lead to non-trivial dynamics (Tobias & Marston Reference Tobias and Marston2017). The forcing term is added to the tendency equation for the first cumulant (2.7).
Figure 1. Snapshot of the vorticity field 
${\zeta (x,y)}$ for NL (a), QL vs CE2 (b,c) and GQL vs GCE2 (d,e) for two-scale Kolmogorov flow with resolution 
$M = N = 10$ at 
$t = 1000$ days. All colour ranges are identical.
The NL solution (figure 1a) consists of a strong band of positive vorticity centred at 
$y = {\rm \pi}$ surrounded by negative vortical regions, as was observed by Tobias & Marston (Reference Tobias and Marston2017). The band of positive vorticity appears to be composed of a strong zonal mean and a sinusoidal zonal wave with wavenumber 
$m = 1$, the latter contributing to the varying width of the band as a function of 
$x$. As noted in Tobias & Marston (Reference Tobias and Marston2017), this vorticity solution corresponds to a rightward flow on the upper half 
$y> {\rm \pi}$ of the domain and leftward flow in the lower half.
The QL solution (figure 1b) captures the central positive vorticity band, albeit with an apparently weaker 
$m = 1$ harmonic resulting in a band of more uniform height and thickness. For presentation of results, we adopt zonal (spatial in 
$x$ direction) averaging so that 
$\bar {\boldsymbol {q}}$ denotes the zonal mean of quantity 
$\boldsymbol {q}$. The zonal mean vorticity 
$\bar {\zeta }$ predicted by CE2 (figure 1c) agrees reasonably well with QL dynamics both in magnitude and distribution. The localised negative vorticity regions predicted by QL (figure 1b) in the upper and lower half of the domain are averaged out in the CE2 prediction (figure 1c). Figure 1(d,e) shows the GQL (figure 1d) and the resolved GCE2 (figure 1e) solutions for a cutoff 
$\varLambda = 1$. These are in very close agreement with the NL solution, improving significantly over the QL/CE2 solutions. Figure 2 contains Hövmöller plots of 
$\zeta (y,t)$ showing the evolution of the central band of vorticity with time for the five different equation systems. All solutions evolve from the same initially random noise field of power 
$10^{-4}$. Time-averaging is conducted in the window 
$500 < t < 1000$ days. Up until time-averaging begins, the QL, GQL, and GCE2 systems exhibit noticeable fluctuations in jet location over a relatively small time scale of 
$t \sim 10$ days. Such fluctuations appear to be much less pronounced in the NL and CE2 systems, wherein the vorticity jet appears to stabilise at its location at around 
$t = 250$ days. The time-averaged solution of all five systems are in excellent qualitative agreement, with the possible exception of QL which appears to predict a slightly weaker jet (lighter red colour in figure 2b).
Figure 2. Hövmöller plots showing 
$\zeta (y,t)$ for NL (a), QL vs CE2 (b,c) and GQL vs GCE2 (d,e) for the Kolmogorov flow case. All colour ranges are identical. Time averaging commences at 
$t = 500$ days.
Focusing on subtle quantitative differences, we show time-averaged energy spectra for the various solutions in figure 3. Each panel plots 
$E(m,n) = \hat \zeta _{m,n}^*\hat \zeta _{m,n}/{(m^2 + n^2)}$ with zonal wavenumbers 
$m$ running along the 
$x$-axis and non-zonal wavenumbers 
$n$ along the 
$y$-axis. Evidently, all equation systems predict energy to be primarily concentrated in two pairs of conjugate modes, namely 
$(m = 0, n = \pm 1)$ and 
$(m = \pm 1, n = 0)$ (note that the mode 
$(m = 0, n = 0)$ denotes the mean across the entire domain and has been set to contain no energy in our solutions). In the NL case (figure 3a), the remaining energy is spread over spectrally-local modes within a relatively narrow band of zonal wavenumbers, with remaining zonal wavenumbers containing small but non-zero energy. In QL and CE2 (figure 3b,c), this zonal spreading of energy is limited to 
$|m| \le 2$, with no energy in modes with larger zonal wavenumbers. On close inspection, minute differences between QL and CE2 can be observed. For instance, QL (figure 3b) contains slightly higher energy in the 
$(m = \pm 1, n = 0)$ harmonic pair (dark orange squares) and a slightly different energy distribution across the zonal mean modes 
$(m = 0,n)$ as compared with CE2 (figure 3c).
Figure 3. Time-averaged energy spectra 
$\overline {E(m,n)}$ at 
$t = 1000$ days for NL (a), QL vs CE2 (b,c) and GQL vs GCE2 (d,e) for the Kolmogorov flow case.
The GQL and GCE2 solutions (figure 3d,e), however, are in excellent agreement with each other. Moreover, they both clearly improve upon the QL/CE2 solution (figure 3b,c) when compared with the NL solution (figure 3a). Both GQL and GCE2 exhibit similar spreading of energy around the four most energetic modes as seen in NL, albeit with slightly less energy in the outer bands of modes with 
$|m| = 4,5$. Crucially, the modes 
$|m| > 5$ contain a small amount of energy in a similar manner to NL, and as opposed to QL/CE2 where energy is completely absent in these modes. These results validate GCE2 as a statistical theory for GQL dynamics, with the ability to improve upon predictions of CE2.
We now proceed to make quantitative comparisons between the various solutions to elucidate on the qualitative differences observed above. We show the energy distribution over zonal wavenumbers 
$m$ in the 
$n = 0$ slice of the time-averaged energy spectra in figure 4(a). NL (black line) contains energy in all non-trivial zonal wavenumbers 
$m \in [1,9]$. Moving away from the strongest wavenumber 
$m = 1$, energy tapers-off rapidly across the wavenumbers 
$m \le 4$ and more gently along the wavenumbers 
$m \ge 5$. GQL (orange line) and GCE2 (orange dots) mimic this distribution of energy reasonably well; however, the 
$m = 1$ mode is slightly stronger and the wavenumbers 
$m \le 4$ have slightly lower energy in GQL/GCE2 than NL. Interestingly, minor departures of energy between GCE2 and GQL are apparent, particularly for 
$m = 5,7$. In comparison, QL (blue line) and CE2 (blue dots) contain most energy in two wavenumbers. Comparing these modes carefully, we note that QL contains significantly more energy in the 
$m = 1$ mode than does CE2, showing that QL and CE2 solutions also have disagreements.
Figure 4. (a) One-dimensional slice of the time-averaged energy spectrum 
$\overline {E(m,0)}$ predicted by the different equation systems for the Kolmogorov flow case. (b) Zonal mean vorticity profile 
$\zeta (y;t = t_\infty )$.
In figure 4(b), we plot the zonal mean vorticity as a function of the 
$y$-coordinate. The centreline vorticity (
$y = {\rm \pi}$) of the NL solution is quite well predicted by CE2, followed closely by GQL and GCE2. However, the QL solution departs conspicuously from the remaining predictions at the centreline. Differences between the solutions away from the centreline are relatively small, although, at 
$y \sim 3{\rm \pi} /2$, the NL solution is over-predicted by both solution pairs GQL/GCE2 and QL/CE2. Again, the QL and CE2 solutions exhibit significant departures from each other. Though we have noted here a difference between QL and CE2 solutions, no major departures are immediately evident between the GQL and GCE2 solutions.
The results above demonstrate that GQL/GCE2 improve upon the predictions of QL/CE2 when compared with the NL solution. Barring the minor differences, these results also validate the predictions of DSS (CE2 or GCE2) against the corresponding dynamical systems (QL or GQL, respectively). Now we focus our attention on investigating the origin of divergences between each solution pair QL/CE2 and GQL/GCE2. As mentioned earlier, Nivarti et al. (Reference Nivarti, Kerswell, Marston and Tobias2022) pointed out recently that divergences can and do arise between the QL and CE2 solutions despite the mathematical correspondence between them. Importantly, they linked such divergences to a rank instability available to the dynamics of the CE2 system and prohibited within the QL system. Divergences occur when the rank of a zonal submatrix 
$C^{(m)}$ in the CE2 second cumulant departs from its initial unity value: the corresponding rank in QL must always remain unity, thus constituting an important source of differences. We define the matrix 
$C^{(m)}$ to be the zonal decomposition of the second cumulant: 
$C^{(m)} = \hat {c}_{\zeta \zeta } (m, n_1; m, n_2) = ({1}/{(2{\rm \pi} )^3})\int _{0}^{2{\rm \pi} }\int _{0}^{2{\rm \pi} }\int _{0}^{2{\rm \pi} }c_{\zeta \zeta }(x_1,y_1; x_2, y_2) \exp ({\textrm {i} m (x_2 - x_1) - \textrm {i} n_1 y_1 - \textrm {i}n_2 y_2})\,\textrm {d} x_1 \,\textrm {d} x_2 \,\textrm {d} y_1 \,\textrm {d} y_2$ for 
$m \in [1,M)$ where 
$M$ is the spectral cutoff in the zonal wavenumber. In order to compare GQL against GCE2, the zonal projection of the GCE2 field can be employed.
Figure 5 simultaneously plots the time evolution of two different quantities for 
$0 \le t \le 20$ days. The left vertical axis quantifies the absolute difference in zonal energy 
$|E(m)_{CE2} - E(m)_{QL}|$ between CE2 and QL for a given zonal mode 
$m$ plotted using blue lines. The right vertical axis quantifies the rank 
$\textrm {rank}(C^{(m)})$ of the cumulant submatrix 
$C^{(m)}$. 
$E(m) = \sum _n \hat {\zeta }_{m,n}^*\hat {\zeta }_{m,n}/(m^2 + n^2)$ where 
$\hat {\zeta }_{m,n}$ is the relevant Fourier coefficient of vorticity 
$\zeta (x,y)$, which is calculated in QL and GQL as 
$\zeta = {\zeta }_\ell + \zeta _h$. In CE2 and GCE2, 
$E(m)$ can be calculated as 
$E(m) = \sum _n |\hat {c}_\zeta (n)|^2/n^2$ for 
$m = 0$ and 
$E(m) = \sum _n C^{(m)}(n,n) /(m^2 + n^2)$ for 
$m \ge 1$. Figure 5(a) compares the QL and CE2 solutions considering the zonal mode 
$m = 1$. As both systems evolve from an identical random noise initialisation, the difference 
$|E(1)_{CE2} - E(1)_{QL}|$ (blue line) is zero initially. The zonal energies begin to differ around 
$t = 5$ days showing an increasing magnitude of differences with time. Simultaneously with the emergence of differences, the rank of second cumulant submatrix 
$C^{(1)}$ in CE2 (orange line) departs from unity, increasing nearly monotonically until 
$t = 20$ days. This is a clear indication of the strong link between the onset of rank instability in CE2 and CE2's divergence from QL as also found by Nivarti et al. (Reference Nivarti, Kerswell, Marston and Tobias2022). In a similar manner, figure 5(b) compares the same quantities for GQL and GCE2 with spectral cutoff 
$\varLambda = 1$ for the high mode 
$m = 2$. The zonal energy difference 
$|E(2)_{GCE2} - E(2)_{GQL}|$ (blue line) shows that the respective solutions begin to diverge even before 
$t = 5$ days. Again, this happens simultaneously with the departure of the corresponding submatrix 
$C^{(2)}$ (orange line) rank from unity. In fact, in this case, we also observe that at times when the solutions have very similar energy (say at 
$t \sim 8$ days), the rank also returns to unity. These results confirm that, akin to divergences observed previously for QL and CE2 (Nivarti et al. Reference Nivarti, Kerswell, Marston and Tobias2022), predictions of GQL and GCE2 may also diverge from one another as a result of the onset of rank instabilities. We have reproduced the GCE2 rank instability using independently written code (see Marston & Tobias Reference Marston and Tobias2023) for the case of Kolmogorov forcing on the sphere.
Figure 5. Difference in energy 
$E(m)$ of zonal mode 
$m$ for (a) the QL and CE2 solutions (
$m = 1$) and (b) the GQL and GCE2 solutions (
$m = 2$), with the corresponding rank 
$C^{(m)}$ shown on the right. The first 
$20$ days of the Kolmogorov flow case are shown with half the timestep size as before. In each comparison, the divergence of zonal energies in the dynamical (QL, GQL) and statistical (CE2, GCE2) solutions appears to be strongly associated with the onset of rank instability in the latter.
4.2. Narrow-band stochastic forcing
We now consider a stochastically-driven system on the rotating 
$\beta$-plane. We adopt the formalism detailed in Constantinou et al. (Reference Constantinou, Farrell and Ioannou2016) which was used for the validation of the statistical state dynamics (SSD) model. The forcing term 
$F$ in the vorticity equation becomes
\begin{equation} F(\boldsymbol{x},t) = {\sqrt{\varepsilon} \xi}, \end{equation}
where 
$\varepsilon$ is the energy injection rate and 
${\xi }$ is the white-in-time stochastic noise with zero mean and covariance 
$Q(\boldsymbol {x})$. Following Constantinou (Reference Constantinou2015, H.4), the covariance 
$Q(\boldsymbol {x})$ is specified on the Fourier domain by the forcing spectrum 
$\hat {Q} (\boldsymbol {k}) = c(m)^2 d^2 \exp ({-n^2d^2})$ which is non-zero for zonal wavenumbers 
$m \in [m_f, m_f + \delta m)$ where 
$d = 0.1$ and 
$c(m)$ is such that the total contribution of 
$\hat {Q}(\boldsymbol {k})$ is unity for each zonal wavenumber (see Constantinou Reference Constantinou2015, H.5), and as a result the total energy injection rate of 
$F(\boldsymbol {x},t)$ in (4.2) is 
$\varepsilon$. Here, we have set 
$m_f = 8$ and 
$\delta m = 2$ (two zonal wavenumbers 
$8, 9$ are forced) on a 
$12\times 20$ grid. Because stochastically driven modes have zonal wavenumbers that exceed the cutoff 
$\varLambda$ the covariance of the stochastic forcing is added to the tendency equation for the second cumulant (2.8).
Figure 6 shows Hövmöller plots 
$\zeta (y,t)$ for a case corresponding to the equatorial 
$\beta$-plane (
$\theta = 0^\circ$) with 
$\beta = 10.0$ and 
$\mu = 0.01$. As is customary (Tobias & Marston Reference Tobias and Marston2013), a hyperviscosity coefficient 
$\nu _4$ is applied such that the largest wavenumber dissipates energy at unit rate. The energy injection rate (described above) via the forced zonal wavenumbers is 
$\varepsilon = 0.02$. The dynamical equations (NL, QL and GQL) are solved for 10 000 days, whereas the statistical equations (CE2 and GCE2) need only be solved for a tenth of that time.
Figure 6. Hövmöller plots showing 
$\zeta (y,t)$ for NL (a), QL vs CE2 (b,c) and GQL vs GCE2 (d,e) for a the stochastically forced case with resolution 
$M = 12$, 
$N = 20$. All colour ranges are identical. Jet migration is captured by GCE2 whereas CE2 fails to capture it.
The NL solution (figure 6a) consists of two opposed jets of vorticity that migrate gently downwards beginning at 
$t = 5000$ days. The migration appears to occur at constant speed (fixed slope in 
$t$–
$y$ coordinates) as has been identified previously by Cope (Reference Cope2020). The QL solution (figure 6b) predicts instead a relocation of the positive vorticity jet from the lower half of the domain to its centreline during the period 
$2500 < t < 7500$ days. After this period, the jet stabilises to a steady position. No vertical movement is recorded in the CE2 solution (figure 6c), which merely captures the dual-jet solution with the sharp inter-jet boundary lying along 
$y = 0 \forall t \in [0,1000]$ so the migration is not captured by a QL theory, as also determined by Cope (Reference Cope2020). The GQL and GCE2 solutions are shown in (figure 6d,e). Whereas GQL (figure 6d) appears to predict a mild degree of jet migration, GCE2 (figure 6e) records the downward jet migration as observed in NL (figure 6a). Note that the slope in GCE2 must be magnified 10 times for an apples-to-apples comparison with that in NL due to the former's much shorter run. Here too then, GCE2 improves significantly over CE2 by capturing jet migration that CE2 cannot.
In order to elucidate on the differences in jet behaviour recorded by GQL and GCE2, we conducted a series of GQL simulations with different random seed values used to generate the initial conditions. This resulted in an ensemble containing a wide range of jet behaviour, including jet migration and steady jet propagation. In the interest of brevity, we show in figure 7 the results of a single GQL run (left) chosen from the ensemble that predicts upwards jet movement with a speed of roughly 
$0.41m^\circ$ per day. We compared this with a GCE2 simulation (left) initialised with maximum ignorance, i.e. with a full-rank second cumulant initialised with power 
$10^{-6}$. GCE2 records the identical upwards jet migration with a speed of roughly 
$0.44m^\circ$ per day. Note that the 
$x$-axis ranges in both figures are significantly different, so the slopes do not appear identical visually. It is interesting that GCE2 with a maximum ignorance initial condition predicts upwards jet migration rather than in the opposite direction as observed in GCE2 with a unity rank IC (see figure 6).
Figure 7. (a) GQL instance from of an ensemble of differently seeded random noise initial conditions. (b) GCE2 initialised with a maximum ignorance initial condition run for a shorter period (colour range is identical to left figure). GCE2 with maximum ignorance captures jet migration with a similar speed (indicated in 
$m^\circ = 10^{-3}$ degrees per day) as seen within a large ensemble of GQL runs.
In figure 8, we show time-averaged energy spectra for the different solution methods for the stochastically forced case shown in figure 6. The NL system (a) contains energy in the 
$m = 1$ zonal mode, primarily via the 
$n = \pm 1$ non-zonal wavenumber (which corresponds to the opposed jet configuration). A significant portion of the energy is also present in the 
$m = 1$ and 
$m = 2$ zonal modes, the remaining being scattered through the entire range of zonal modes but a relatively narrow band of non-zonal modes. In comparison with NL, the QL and CE2 solutions (figure 8b,c) consist of localised bands of excited zonal modes. For instance, the 
$m = 1$ and 
$m = 2$ modes are relatively very weak in the QL and CE2 solutions; thus, also, the modes 
$m = 6$ or 
$m \ge 9$ appear to be absent in comparison with NL. This localisation of energy arises because scattering is unavailable in QL (and, therefore, CE2) owing to the absence of the required nonlinearities; a given zonal mode 
$m$ becomes energetic in QL only via a corresponding instability of the mean flow. Once energetic, the zonal mode 
$m$ can transfer its energy back to the mean via self-interactions, but may not transfer energy to another non-mean mode 
$|m| \ne 0$. On closer examination, it is revealed that QL and CE2 exhibit differences. The differences between energy distribution over non-zonal wavenumbers in the 
$m = 2$ mode are immediately evident between QL (figure 8b) and CE2 (figure 8c). In other words, the CE2 solution obtained is not identical to that of a single realisation of QL. As in the Kolmogorov flow case, we hold that these divergences are linked to the rank instability in CE2.
Figure 8. Time-averaged energy spectra 
$\overline {E(m,n)}$ for NL (a), QL vs CE2 (b,c) and GQL vs GCE2 (d,e) for the stochastically forced case. GCE2 improves considerably over CE2, but also diverges from GQL.
The GQL and GCE2 solutions in figure 8(d,e) improve considerably over the QL/CE2 solutions. The missing zonal modes 
$m = 1,2,6,$ etc. are found to be energetic in GQL/GCE2. This is because nonlinear interactions involving the 
$m = 1$ mode that are available within GQL (and thereby GCE2) allow for scattering of energy leading to a broader spread of energy over the range of zonal wavenumbers. It is remarkable though that GQL with a cutoff 
$\varLambda = 1$ appears to be in excellent agreement with NL (figure 8a). We note that GCE2 diverges from GQL to some extent in the same way CE2 was seen to diverge from the QL solution. Broadly speaking, energy appears to be less spread out in GCE2 as compared with GQL, being more concentrated to small non-zonal wavenumbers. For instance, energy in the 
$m = 1$ mode in GCE2 (figure 8e) is limited to a significantly narrower range of non-zonal wavenumbers 
$|n| \lesssim 3$. In addition, the 
$(\pm 1,0)$ modes are weaker in GCE2 by an order of magnitude.
We show the energy distribution over zonal wavenumbers 
$m$ in the 
$n = 0$ slice of the time-averaged energy spectra for the stochastically forced case in figure 9. NL (black line) is shown to contain energy in all non-trivial zonal wavenumbers 
$m \in [1,11]$ supporting the findings of figure 8. However, figure 9 clarifies that the 
$m = 1$ and 
$m = 2$ zonal wavenumbers are associated with the lowest and highest energies, respectively. Moving away from the wavenumber 
$m = 2$, energy tapers off steadily across the remaining wavenumbers 
$m \le 11$. The energy distributions predicted by GQL (orange line) and GCE2 (orange dots) are in good agreement with the NL energy distribution, save for the 
$m = 1$ mode estimated to be slightly weaker by GCE2. At this wavenumber and at 
$m = 3$, GQL exhibits a slight departure from GCE2, aligning more favourably with NL than GCE2. Contrasting with the GQL and GCE2 solutions, QL (blue line) and CE2 (blue dots) do not contain energy in at least 5 of the 11 wavenumbers 
$m \in [1,11]$. CE2 does not contain energy in the zonal wavenumbers 
$m = {1,2,7,10,11}$. In comparison, QL is additionally deprived of energy in the 
$m = 6$ zonal wavenumber also, for which CE2 predicts an energy level comparable to NL/GQL/GCE2. Thus, the wavenumber 
$m = 6$ is a point of significant departure between the QL and CE2 solutions. The QL and CE2 energy distributions also exhibit a conspicuous departure in the 
$m = 5$ wavenumber and a relatively minor departure in the 
$m = 9$ wavenumber. These quantitative comparisons provide affirmation to the following: (1) GQL and GCE2 solutions better estimate the NL energy distribution; and (2) neither pair of dynamical and statistical solutions, GQL/GCE2 or QL/CE2, is devoid of divergence.
Figure 9. One-dimensional slice of the time-averaged energy spectrum 
$\overline {E(m,0)}$ predicted by the different equation systems for the stochastically forced case.
In light of the divergences shown previously for Kolmogorov forcing, it would suffice here to show a departure in the rank of a second cumulant submatrix to explain the divergences seen above for the stochastically driven case. In figure 10, we show the final-time ranks of second cumulant submatrices. For CE2 (orange bars in figure 10a), the zonal modes 
$m = 8,9$ are full rank by virtue of the stochastic driving; however, we note that a number of additional zonal mode submatrices also depart in rank from unity as is held by the corresponding QL solution (blue bars in figure 10a). Since there is no pathway of interactions that can transfer energy between the various non-mean zonal modes in QL (and CE2 by extension), each of the zonal modes with non-unity rank must have acquired its own rank instability. This contrasts with GCE2 (orange bars in figure 10b) where all high zonal modes 
$|m| > 1$ for the spectral cutoff 
$\varLambda = 1$ are full rank. In essence, scattering of energy allowed within the GQL/GCE2 formalism causes full rank of stochastically driven modes to be transferred to all remaining zonal modes. We term this rank scattering. Thus, we observe that communication matters: moving from CE2 to GCE2 changes significantly the channels of communication so that energy scattering and rank scattering can occur.
Figure 10. Comparisons of ranks 
$C^{(m)}$ in the end point solution for QL and CE2 (a) and for GQL and GCE2 with 
$\varLambda = 1$ (b). In CE2, each zonal mode undergoes its own rank instability, whereas in GCE2, the allowed 
$HL\rightarrow H$ interactions cause ‘rank scattering’, and the (full) rank of forced zonal modes (
$m = 8,9$) is adopted by all other high modes.
5. Conclusions
In this paper we have tested a method of DSS that is obtained as a mathematically exact closure for the GQL equations. This method of DSS, which we term GCE2, adopts generalised cumulant expansions and is capable of systematically interpolating between statistics corresponding to QL and NL equations. We have implemented GCE2 in a numerical code for simulations on the 
$\beta$-plane with two driving models, deterministic and stochastic. Our simulations, which employ a minimal spectral cutoff 
$\varLambda = 1$, confirm that GCE2 improves considerably over CE2, the DSS method corresponding to QL dynamics.
We comment here that GQL (and hence its statistical manifestation GCE2) can be derived using the method of multiple scales in space and time (see e.g. Marston et al. Reference Marston, Chini and Tobias2016). The approximation is made formally rigorous by separating variables into those that undergo rapid dynamics on small scales and more slow dynamics of the large-scale variables. Intriguingly, similar methods have been used for homogeneous isotropic turbulence (see e.g. Laval, Dubrulle & Nazarenko Reference Laval, Dubrulle and Nazarenko2001). That paper showed that the degree of intermittency was crucially determined by the form of the interactions that were included in the model and a new model of turbulence was suggested where non-local interactions coupled large and small scales nonlinearly (with the addition of noise) and the local interactions could be modelled by a turbulent viscosity; it would be interesting to test the applicability of generalisations of such a model to anisotropic flows.
We have also shown that statistics of GCE2 may depart from those of GQL due to the rank instability as found recently for QL and CE2 (Nivarti et al. Reference Nivarti, Kerswell, Marston and Tobias2022). CE2 and QL solutions may diverge at identical parameters despite the fact that CE2 is an exact mathematical closure for QL. Such divergences (also observed by Oishi et al. Reference Oishi, Burns, Marston and Tobias2022) are linked to the emergence of an instability, the rank instability, that is available within the CE2 system but unavailable in QL. The origin of the rank instability can be understood analytically with a simple linear model (Nivarti et al. Reference Nivarti, Kerswell, Marston and Tobias2022). We have found that the GCE2 system admits such instabilities too, and therefore its solutions can and do divergence from solutions of the GQL system. This is a feature rather than a bug. CE2 is a statistical description, whereas any single realisation of QL dynamics is a dynamical one. Likewise GCE2 is a powerful method for self-consistently modelling the full dynamics of a range of spectral scales, coupled with the statistics of the QL dynamics of the smaller spectral scales. An investigation of the GCE2 approximation at high resolution should be illuminating.
Funding
We acknowledge support of funding from the European Union Horizon 2020 research and innovation programme (grant agreement number D5S-DLV-786780). J.B.M. and S.M.T. are supported in part by a grant from the Simons Foundation (grant number 662962, GF). J.B.M. acknowledges funding support from the United States Department of Energy (grant number DE-SC0024572).
Declaration of interests
The authors report no conflict of interest.
Open access
