Direct Statistical Simulation of the Busse Annulus

We consider direct statistical simulation (DSS) of a paradigm system of convection interacting with mean flows. In the Busse Annulus model zonal jets are generated through the interaction of convectively driven turbulence and rotation; non-trivial dynamics including the emergence of multiple jets and bursting `predator-prey' type dynamics can be found. We formulate the DSS by expanding around the mean flow in terms of equal-time cumulants and arrive at a closed set of equations of motion for the cumulants. Here, we present results using an expansion terminated at the second cumulant (CE2); it is fundamentally a quasi-linear theory. We focus on particular cases including bursting and bistable multiple jets and demonstrate that CE2 can reproduce the results of direct numerical simulation if particular attention is given to symmetry considerations.


Introduction
Turbulence interacts with, and leads to the generation of, mean flows in a wide variety of natural systems. Understanding the nature of these interactions, particularly for systems far from equilibrium, remains a key problem for the description of large-scale flows on planets and stars. The fundamental problem in studying such systems via direct numerical simulation (DNS) is the expense of the calculations, owing to the vast range of spatial and temporal scales that need to be described. One way to make progress is to study the statistics of the flow rather than detailed flow variables themselves ).
However, one should bear in mind that these systems are often inhomogeneous and strongly anisotropic, and such a statistical description should respect this, despite the simplifications that are afforded with assumptions of homogeneity and isotropy.
Rotation is a hallmark of geophysical and astrophysical fluid dynamics. Through the breaking of reflection symmetry, rotating systems generate mean flows such as zonal jets; examples of which include the bands on Jupiter, the jet stream on Earth and the differential rotation in stars (Galperin & Read 2019). In these systems, turbulence that is driven at moderate scales interacts with rotation to produce non-trivial correlations that lead to Reynolds stresses that drive differential rotation and jets. Here, the turbulence is often described as having an anti-frictional character as it moves the system away from solid body rotation; this provides a stern test for theories that seek to describe the statistical properties of the turbulence. We also note that it is common astrophysically for the source of the turbulence to be thermal convection, and that the convection itself is sensitive to the presence or absence of large-scale flows; this leads to non-trivial feedbacks between the mean flows and the turbulence.
Here, we explore the simplest system that describes the interaction of convection with mean flows -the Busse Annulus. The Busse Annulus models rotating convection in an annulus with slanted ends leading to a topographic β effect; this effect leads to vortex stretching and the generation of systematic large-scale flows (see e.g. Busse 1976;Brummell & Hart 1993;Rotvig & Jones 2006). We use this system to test the effectiveness of the statistical description termed CE2. CE2 is a quasilinear approximation of direct statistical simulation (DSS), and is a cumulant expansion truncated at second order. CE2 has been shown to be effective for simple systems where tightly coupled correlations control the dynamics and driving is either stochastic or arises through the instability of a shear flow. However the case of thermal convection, where buoyancy provides the driving in the vorticity equation and there is nonlinearity in both the vorticity and temperature equations, has been much less studied.
Thus, the Busse Annulus system presents an important challenge for CE2, not least because it is known to host multiple solutions at modest Rayleigh numbers (Brummell & Hart 1993). How does a quasilinear statistical theory fare when faced with multiple potential solution basins? In particular, we investigate the sensitivity of CE2 to initial conditions for a case with multiple solution branches.

Model Equations, Formulation, Parity and Numerical Methods
The Busse Annulus we consider is a locally Cartesian model of rotating (at angular velocity Ω = Ωe z ), incompressible Boussinesq fluid with viscosity ν and thermal diffusivity κ. Here, gravity is uniform and in the y-direction. The system is non-dimensionalised, with a fiducial length being the width of the annulus in the y-direction (d), a typical timescale being the viscous timescale d 2 /ν and a typical velocity scale ν/d; a typical temperature is chosen to be ∆T , the temperature difference between the inner and outer walls.
Following Brummell & Hart (1993); Rotvig & Jones (2006), the temperature T = T BS +θ(x, y) is decomposed into a basic state profile T BS satisfying ∇ 2 T BS = 0 and a perturbationθ. After some manipulation, the equation for the vertical (z) component of the vorticity (ζ) is given by where ψ is the streamfunction, (u, v) = (− ∂ ∂y ψ, ∂ ∂x ψ), and ζ is related to ψ by is the Jacobian. The equation for the temperature perturbation θ is given by The system is governed by four dimensionless parameters, β, C, Pr , and Ra (see e.g. Tobias et al. 2018, for the definition of these parameters). Here β measures the degree of vortex stretching engendered by the sloping endwalls and C measures the degree of friction, whilst the more familiar Rayleigh number and thermal Prandtl number have their usual physical interpretations. A zonal CE2 formulation (sometimes termed zCE2) is implemented using a zonal average, to split all dynamical variables into mean and fluctuating components, We then derive evolution equations for the first cumulants c ζ (y) = ζ and c θ (y) = θ and the second cumulants c θθ (ξ, y 1 , y 2 ) = θ (x 1 , y 1 )θ (x 2 , y 2 ) (2.6) with ξ = x 1 − x 2 . Similar definitions are used for the second cumulants c ζζ , c θζ , and c ζθ . Cumulants involving ψ and ζ are related by differential operators (Tobias & Marston 2013). Because the equations require both ψ and ζ cumulants, we solve for c ψ , c θ , c ψψ , c θψ , c ψθ , and c θθ , though we write equations in terms of ζ cumulants. The dynamical equations for the first cumulants are given by  and (2.11) The other second cumulant terms (e.g. c ζθ ) do not need their own evolution equations, as they can be calculated from symmetry considerations. We solve both the dynamical PDE (DNS) system and the cumulant system subject to impenetrable, stress-free, zero temperature perturbation boundary conditions in the y dimension; both systems are periodic in x. This implies that θ, ζ, and ψ all have odd parity in y (and are periodic in x). As the action of the zonal average preserves the parity, we discretise the first cumulants using a sine series in y. Likewise, the second cumulants may be discretised using sine series in y 1 and y 2 and a Fourier basis in x.
We use Dedalus (Burns et al. 2020) to solve both the direct equations (2.1 -2.3) and the CE2 model equations, (2.7)-(2.11). DNS use the SBDF2 timestepper with linear terms implicit, while the CE2 equations use RK222 scheme with all terms evaluated explicitly. We have found fully explicit timestepping to be more stable for CE2. The second cumulants must be positive definite in order to ensure the corresponding probability density function is realisable. At higher Ra, we have found it necessary to remove negative eigenvalues from the second cumulant every 10 timesteps. While this is required for CE2.5 , it is also necessary in these CE2 simulations.

Results
Given that the Busse Annulus system is known to exhibit strong hysteresis, it is key to examine the role of initial conditions when solving the cumulant system. Minimally, we begin each CE2 simulation using one of two extreme initial conditions: maximum ignorance, in which we assume a completely uncorrelated c θθ and all other cumulants (both first and second) zero, and maximum knowledge, in which we initialise all fields using cumulants calculated from the statistics of a DNS solution. This allows the evaluation of CE2's ability to both find solutions and continue solutions already known. In addition, we begin some solutions with biased initial conditions that are designed to ensure solutions are not trapped in symmetry subspaces. In order to facilitate comparison with previous work, we adapt the same run naming scheme as Tobias et al. (2018).
3.1. Run A: 2 Jet/3 Jet solutions The first case we consider (case A) exhibits relatively simple dynamics in DNS. Started from random initial conditions, the driven convection interacts with the βeffect to drive large-scale zonal flows (jets). These takes the form of one prograde   , c θ (bottom) from run A started in DNS and continued with CE2 at t = 2. Attached to the right y axis of the Hovmöller diagrams: first cumulants of cu (top) and c θ (bottom) as a function of y for DNS averaged over 1 t 2 and CE2 averaged from 2 t 3. Right: total kinetic energy. and one retrograde jet as shown in the solutions for the mean zonal velocity u (y, t) and mean temperature θ (y, t) in the Hovmöller plots of Figure 1. Started from maximum ignorance initial conditions, however, CE2 evolves to a 3-jet solution, which is symmetric about the midplane (y = 0.5). These results are similar to the strictly quasilinear run in Tobias et al. (2018), which also produces a 3-jet solution; an unsurprising result given that CE2 is fundamentally a quasilinear theory. Interestingly, generalized quasilinear (GQL) solutions do however yield the correct number of jets Tobias et al. (2018). At these parameters it has been demonstrated in DNS that there is hysteresis between 2-jet and 3-jet solutions, but that 3-jet solutions are unstable to perturbations that break shiftreflect symmetry (Brummell & Hart 1993). To test this, we ran a DNS with initial conditions that respect this symmetry and, owing to our highly accurate spectral methods, the solution remained in the symmetry class and the 3-jet solution remained stable in DNS.
Better comparison with DNS for these parameter values is achieved when the CE2 solution is biased by initialising the first cumulant of the x-velocity, c u to c u (t = 0, y 0 ) = A 0 (λ cos(2π/L y y 0 ) + (1 − λ) cos(π/L y y 0 )) . (3.1) The first term is odd about the center of the domain, while the second is even. We find that CE2 is capable of finding and sustaining a 2-jet solution, with amplitude comparable with the DNS results, if the symmetry of the initial condition is sufficiently biased. This strongly suggests that the transition from 3 to 2 jets in DNS is the product of a subcritical, non-linear transition mediated by eddy-eddy → eddy interactions excluded both from CE2 and our previous quasilinear results. Regardless, this solution remains a fixed point of the CE2 system, suggesting that the selection of a given multi-jet solution is distinct from its maintenence. This hypothesis is supported by our "maximal knowledge" calculation where the CE2 solution is initialised with the cumulants calculated from the saturated state of the DNS run (see Figure 1). CE2 is certainly capable of maintaining the form of these fully nonlinear solutions (even maintaining the slight asymmetry in the temperature distribution) even though the eddy-eddy non-linearity (EENL) term is missing from the system. The degree of success of the quasilinear CE2 models can be investigated further by comparing the second cumulants; it is possible for first cumulants to agree well whilst second cumulants diverge. Figure 2 shows the evolution of the second cumulant c θθ (ξ, y 1 , y 2 = 0.5) under CE2 for the solution started from the statistics of the saturated state of a DNS calculation. The figure shows that this covariance of the DNS is fairly localised in space; we hypothesise that EENL interactions are important in maintaining the locality of these correlations. This hypothesis is validated by the fact that the correlations delocalise in space as the quasilinear CE2 calculation progresses. Eventually the second cumulant is very delocalised and has a wave-like form. The long-range teleconnections are typical of the quasilinear interaction between waves (in this case thermally driven Rossby waves) and a mean flow.

Multiple Jet Solutions
Parameters may also be selected that yield multiple jet solutions; case R, for which the DNS is shown in Figure 3(a), is such a case. Here, after an initial transient that takes the form of a six jet solution, a jet merging leads to the solution settling down into a five jet solution. DNS exhibits a solution for the second cumulant that is very localised in space, which is suggested by the broad power spectrum of the  Figure 4: Power spectra of ζ for run R as a function of kx for DNS, unbiased CE2, and biased CE2. Note that DNS has energy distributed across many k, demonstrating the importance of the EENL. However, both CE2 runs show power only at kx0 and a band from 20 kx 50 despite the fact that the biased solution gets the correct cu and the unbiased one does not. Taken together, this suggests that EENL is not crucial for maintaining the mean flow and may only lead to additional dissipation. Figure 4. For this case, the maximally ignorant quasilinear CE2 model yields an incorrect seven jet solution. Remarkably, if the CE2 calculation is biased initially to yield a three jet solution, the solution evolves via a seven jet solution to form a five jet solution of the correct amplitude, obviously driven by the corresponding evolution of the second cumulants; the basin of attraction for CE2 is clearly very sensitive to the initial condition here. The power spectra of the CE2 calculations shown in Figure 4 demonstrate that power is concentrated in a finite number of k x values indicating a delocalisation of the second cumulant; the second cumulant is again a superposition of wave-like solutions. The precise wavenumbers that receive energy for CE2 are sensitive to the structure of the first cumulant and therefore the initial conditions. Moreover, the "zigzag" structure of the power spectra in CE2 shows that most of the power is trapped within a symmetry subspace and is not efficiently scattered into all modes. Nonetheless we stress that the interactions are sufficient to drive the correct mean flows. As for Case A, the CE2 solution also tracks the DNS solution well if a "maximal knowledge" initial condition is used (not shown).

Bursting Dynamics
The Busse Annulus model is capable of producing extremely complicated nonlinear spatio-temporal dynamics. Perhaps the most nonlinear behaviour exhibited by the DNS model is shown in the top panel Figure 5(a), which show timeseries of the total, zonal, and non-zonal kinetic energies of the flow and Hovmöller plots of the zonally averaged zonal flows for DNS. Here the solution takes the form of a relaxation oscillation; a cycle proceeds as follows. Convection is driven by strong temperature gradients and interacts with the rotation to generate zonal flows (in this case a two-jet solution). As more energy is pumped into the zonal flow, the shear acts so as to switch off the convection (acting as a barrier to transport). This in turn removes the energy source for the zonal flow, which decays on a longer timescale. Once the zonal flow is weak enough convection sets in again and the process is repeated. This type of behaviour has been described in terms of predator-prey dynamics with the convective turbulence taking the role of the prey and the zonal flows acting as a predator. Figure 5(b) shows one of the "maximum knowledge" solutions for CE2. Here the solution is started from the DNS solution when it has reached a peak in the zonal flow energy. Remarkably, CE2 is able to continue relaxation oscillations from this state; we believe that this is the first time a quasilinear model has reproduced such complicated behaviour; previously they have been able to generate either the state at the peak or the trough of such a relaxation oscillation -but not the transition between them (Plummer et al. 2019). If CE2 is started from the trough in the zonal flow energy, similar results are obtained (not shown). The sensitivity to initial conditions for CE2 is, however, demonstrated in Figure 5(c), which shows the evolution from small amplitude "maximal ignorance" initial conditions. This solution clearly does not replicate the correct number of jetsthough it does show some indication of bursting behaviour, as the strength of each zonal jet waxes and wanes in response to the driving.

Rank of Second Cumulants
One important difference between QL and CE2 simulations is that the former explicitly provides only a single realization of the dynamics, while the latter does not; it provides a description of the low-order statistics of the system and hence can be thought of as averaging over a number of simulations of the system. Recently, Nivarti et al. (2022) have shown that an important difference between the two can be quantified as the difference in the ranks of the second cumulantsthat is, the number of non-zero eigenvalues they possess. They have demonstrated that the statistical description CE2 is open to instabilities that can lead to higher rank solutions for the second cumulant than are accessible for QL. Moreover this can lead to different behaviour for CE2 from that found for a single run of a QL system.
Although we do not investigate these instabililities in detail here, it is interesting to note the rank of the modes of the solutions for CE2. Figure 6 shows the rank of the modes of the second cumulant as a function of time, for run A and various initial conditions. For a QL system the rank of the solution should be three (note the counting of the rank here is different from that in Nivarti et al. (2022)). For a maximum knowledge initial condition for Run A, when CE2 gives the correct answer, the solution accesses a higher rank solution than QL; note QL obtains the incorrect answer. The solutions clearly undergoes a rank instability to access the correct answer. However, for a maximum ignorance initial condition both QL and CE2 give the incorrect answer. The CE2 solution eventually returns to a low rank solution and reproduces the QL (rather than the correct) dynamics. Finally, for the biased initial condition case (with the first cumulant correct but the second cumulant incorrect) the second cumulant remains high rank and one gets the correct behaviour for the first cumulant. We stress again here that in all cases the second cumulant is necessarily delocalised for these quasilinear theories and hence incorrect.

Discussion
In this paper we have examined the effectiveness of the statistical quasilinear theory zCE2 in reproducing the statistics of turbulent solutions of a model of the interaction of rotating convection and zonal flows. This Busse Annulus model exhibits complicated spatio-temporal dynamics including the formation of large-scale zonal jets, multiple zonal jets and even "predator-prey" relaxation oscillations between states with strong zonal flows and strong convection.
We show that zCE2 is capable of reproducing even very complicated mean dynamics (such as the driving and decay of the mean flows and modification of the mean temperature profile) though, for this system, it may have multiple stable attractors for any given parameter set. If the system is initiated with maximal ignorance then it is possible to fall into the basin of attraction of an attractor that is not preferred by DNS. It is possible to bias the symmetry of the initial condition in order to achieve solutions that mirror the symmetry of the DNS solutions. Furthermore one may use an maximal knowledge initial condition where the zCE2 is initiated using the statistics from the saturated state of a DNS run. In this case, zCE2 appears to be extremely successful -even reproducing the extremely nonlinear relaxation oscillations of the highly driven system. zCE2 relies on the solution of equations for quantities averaged in the zonal direction, i.e. the mean flows and temperature and the two-point correlation functions. It is a quasilinear theory in that the evolution equation for the twopoint correlation function is linear in the two-point correlation function, though there are terms that involve the product of the means and two-point correlation functions. The results presented here indicate that even the predator-prey dynamics is controlled by the interaction of the mean flows and temperature with the fluctuations, with the eddy/eddy → eddy nonlinearity (EENL) being of secondary importance.
It is interesting to speculate further on the role of the missing physics encoded in the EENL. This term in the equations can be modelled by deriving evolution equations for the higher order cumulants, and truncating the hierarchy at third order. This extension to DSS is known as CE3 -and sometimes CE2.5 when an approximate form of the CE3 equations are utilised . However this is a computationally expensive procedure and solution strategies suffer severely from the "curse of dimensionality". An interesting question then is what is the best strategy for modelling these higher-order interactions. Of course, if the linear operator provided by the mean flows and other fields is sufficiently non-normal (e.g. if for example the shear is sufficiently strong), then the precise form of this modelling term is not critical (since its role will be to act as a driving term for the non-normal modes of the linear operator). In this case, one may as well model these terms via delta correlated white noise. However it is possible to optimise the form of the noise provided to the linear operator so as to best match the true dynamics. Similar optimisation procedures have been utilised for the form of the (coloured) noise in transition problems (Zare et al. 2017). One strategy that appeals is to use data-driven methods to "learn" the optimal form of the term that should be added to CE2 to best reproduce the dynamics. It is to be hoped (though by no means certain) that this optimal form should be relatively insensitive to the form of the large scales since this dependence should be encoded in the interaction between the first and second cumulants of CE2.