Large-scale drift and Rossby wave turbulence

We study drift/Rossby wave turbulence described by the large-scale limit of the Charney–Hasegawa–Mima equation. We define the zonal and meridional regions as Z ≔ { k : ∣ k y ∣ > 3 k x } and M ≔ { k : ∣ k y ∣ < 3 k x } respectively, where k = ( k x , k y ) is in a plane perpendicular to the magnetic field such that kx is along the isopycnals and ky is along the plasma density gradient. We prove that the only types of resonant triads allowed are M ↔ M + Z and Z ↔ Z + Z . Therefore, if the spectrum of weak large-scale drift/Rossby turbulence is initially in Z it will remain in Z indefinitely. We present a generalised Fjørtoft’s argument to find transfer directions for the quadratic invariants in the two-dimensional k -space. Using direct numerical simulations, we test and confirm our theoretical predictions for weak large-scale drift/Rossby turbulence, and establish qualitative differences with cases when turbulence is strong. We demonstrate that the qualitative features of the large-scale limit survive when the typical turbulent scale is only moderately greater than the Larmor/Rossby radius.


Abstract
We study drift/Rossby wave turbulence described by the large-scale limit of the Charney-Hasegawa-Mima equation. We define the zonal and meridional regions as x y is in a plane perpendicular to the magnetic field such that k x is along the isopycnals and k y is along the plasma density gradient. We prove that the only types of resonant triads allowed are « + M M Z and « + Z Z Z. Therefore, if the spectrum of weak large-scale drift/Rossby turbulence is initially in Z it will remain in Z indefinitely. We present a generalised Fjørtoft's argument to find transfer directions for the quadratic invariants in the two-dimensional k-space. Using direct numerical simulations, we test and confirm our theoretical predictions for weak large-scale drift/Rossby turbulence, and establish qualitative differences with cases when turbulence is strong. We demonstrate that the qualitative features of the large-scale limit survive when the typical turbulent scale is only moderately greater than the Larmor/Rossby radius.

Introduction
Drift waves in plasmas and Rossby waves in the ocean and planetary atmospheres, though unrelated at the first sight, have common features in their dynamics and statistics, which at a basic level can be described by the same model-the Charney-Hasegawa-Mima (CHM) equation (2.1). Much of the research using this model has concentrated on its small-scale limit, r  ¥, where ρ is the ion Larmor radius in plasmas or the Rossby radius of deformation in oceans and atmospheres. This limit is applicable to sub-ion-Larmor motions in plasmas and to mesoscale flows in the Earth's atmosphere, where ρ is of the order of a thousand of kilometres. However, the large-scale limit is more relevant to many plasma regimes since the ion Larmor radius is very small, of the order of just a few millimetres, whereas the tokamak diameter is much larger-several metres. Similarly, in a middlelatitude ocean ρ is tens of kilometres, and since Rossby waves can be hundreds of kilometres in length, the largescale limit is more appropriate. On giant planets like Jupiter and Saturn, ρ is not so much larger than on Earth, but the scale of flows is typically much greater.
In the present paper, we will study large-scale drift/Rossby wave turbulence (WT) within the large-scale limit of the CHM model, r  0. First of all, we will discuss a remarkable property of resonant triad interactions related to a new invariant (semi-action) that has recently been found in this limit by Saito and Ishioka (2013). Namely, defining the zonal region of wave vectors > ≔ { | | } Z k k k : 3 y x and the meridional region 3 y x , we prove that the only allowed resonant triad interactions are « + M M Z and « + Z Z Z. This property has profound effects on the nonlinear evolution. In particular, it leads to the claim that the spectrum of weak large-scale drift/Rossby turbulence which is initially in Z will remain in Z indefinitely. We will use this property and the semi-action invariant for revising Fjørtoft's argument of Balk et al (1991) aimed at predicting directions of the turbulent cascades in the two-dimensional (2D) scale space. We test and confirm our theoretical predictions numerically and study their sensitivity to the strength of nonlinearity using an approach previously introduced in the context of the small-scale CHM system in Nazarenko and Quinn (2009). Finally, we study robustness of our results for systems where the large-scale limit is not well satisfied, and find that most qualitative features survive even for flows with typical scales exceeding ρ only by a factor of two.

The CHM model
As already mentioned Rossby waves in geophysical fluids and electron-drift waves in plasmas are frequently discussed together and can both be described by the same partial differential equation known as the Charney equation in the geophysical context (Charney 1948) and the Haseqawa-Mima equation in the plasma context (Hasegawa and Mima 1978)-hence frequently referred to as the CHM equation: , is in a plane perpendicular to the magnetic field such that x is along the isopycnals and y is along the plasma density gradient, = F f gH 2 is the inverse square of the Larmor/Rossby radius, β is a constant proportional to the gradient of the plasma density or Coriolis parameter. Equation (2.1) is similar to the 2D Euler equation and can be written in the form of an advection equation for the potential vorticity. As we shall later see, like the 2D Euler equation, the CHM equation conserves two quadratic invariants-the energy and potential enstrophy. However, the CHM equation differs from the Euler equation by the presence of the linear term b x This extra term means that the CHM equation can support wave motions, unlike the Euler equation, and also has a weakly nonlinear limit since the linear term can be large compared with the nonlinear one. However, the CHM equation is not only used to study weakly nonlinear waves but also strongly nonlinear structures such as solitons and vortices (Petviashvili andPokhotelov 1992, Horton andIchikawa 1996).
The CHM equation has many related models. Often it is argued that the modified CHM equation is more relevant for drift waves in plasma. The modified CHM equation treats purely zonal modes ( = k 0 y ) differently by taking into account absence of an adiabatic response of the electric field. For the derivation of the modified CHM model readers should refer to Dorland et al (1990). However, since in our paper we ignore purely zonal flow and consider imperfect zonal flow with > ¹ k k k , 0 y x y both models, the CHM and the modified CHM, are equivalent with respect to the solutions that we consider. Note that pure zonal flows with = k 0 y will never appear in a system if they are not present initially in the CHM model. For other basic models of the CHM family, including models such as the Hasegawa-Wakatani equations, see Connaughton et al (2015).
For a discussion of the generation and importance of zonal flows, with more of an emphasis on plasma applications, see Smolyakov et al (2000), Diamond et al (2005) and Gürcan and Diamond (2015).
Let the system be a periodic box, with length L in both directions. The Fourier transform of the stream function is: , e , 2 . 2 k k x i with Fourier coefficients: , .
x y x y Equation (2.1) becomes: and zero otherwise. The wave frequency is given by the following dispersion relation: where for short å º å , x x x y y y k 12 1 2 1 2 1 1 2 2 2 2 2 This symmetric form of the interaction coefficient is valid only on the resonant manifold, i.e. such that both the wavenumber triads k k k , , 1 2 in the nonlinear sum satisfy the following resonance conditions: The wave vectors can either be discrete or a continuous limit could be taken. For weakly nonlinear waves in an unbounded domain,  ¥ L , k becomes a continuous vector. Therefore, any k may be a member of infinitely many resonant triads. This is known as the kinetic WT regime (L'vov and Nazarenko 2010, Nazarenko 2011). In bounded domains, where the wave amplitudes are very small the discrete k-space structure remains important. This is the discrete WT regime (L'vov and Nazarenko 2010, Nazarenko 2011). Resonant conditions 2.12 define the dominant nonlinear interactions in both kinetic and discrete WT regimes.

Conservation laws in weak WT
The nonlinear interactions in both kinetic and discrete WT regimes are weak, so in both cases we deal with the so-called weak WT. The discrete WT regime is described by CHM waveaction in which non-resonant terms are discarded from the nonlinear sum.
The main equation governing the kinetic WT regime, the kinetic equation, can be derived from (2.10) by assuming random phases and taking a large-box limit followed by the limit of weak nonlinearity. It is written below in symmetric form as: is the waveaction spectrum (related to the energy spectrum via n n n n n n 2 3 . 3 k k k k k k 12 12 2 12 1 2 1 2 1 2 and R R ,  We can further neglect = k 0 x as purely zonal flows are not weakly nonlinear-the condition required for us to study triad interactions.
In terms of the waveaction density, n , k the energy and enstrophy become: where w k is now the density of the energy and k x the density of the enstrophy. Note that in this case Ω coincides with the x-component of the momentum invariant-it is strictly positive since > k 0 x . The y-component of the momentum is also conserved, but it is not strictly positive since its density k y may change sign.
Generally, one can write for a conserved quantity Λ with density l k : If the spectral density of quantity (3.5) satisfies condition: on the resonant manifold (2.12) then Λ is conserved, i.e. L = const. Schulman 1980, 1988). For the energy and momentum, l k is w k and k respectively, and the resonant condition (3.6) is obviously satisfied, , due to the respective δ-functions in the kinetic equation (3.1), which proves conservation of these quantities. The same condition (3.6) was shown to be necessary and sufficient for the existence of quadratic invariants (3.5) in the case of discrete WT in Harper et al (2013).
For a generic wave system no other invariant besides the energy and momentum have been found to exist in the kinetic WT regime. However, it was discovered in Balk et al (1990), Nazarenko (1990) and Balk et al (1991) that for a system of Rossby waves, one extra conserved quantity exists, the zonostrophy. It was first found for three special cases: large-scale turbulence (rk ). After this it was generalised to all of k-space in Balk (1991) where the zonostrophy invariant was found to be is the Rossby radius of deformation.

The large-scale limit of the CHM equation
Let us turn our attention to the large-scale limit r  k 0 or, for simplicity, r  0. The large-scale dispersion relation can be found by Taylor expanding the general dispersion relation (2.5): In a moving frame of reference, and assuming for simplicity that br = 1, 4 we can replace the dispersion relation (4.1) with a simpler expression: We can do this as in the large-scale limit the problem is scale-invariant-hence these situations with different br 4 can be obtained from each other by rescaling and no generality is restricted by assuming this. The waveaction variable for the large-scale limit is: Substituting this the into (2.4) gives us (2.10) but now with the interaction coefficient: This can be rearranged as: x x x y y x x y y x x 1 2 1 2 2 2 2 2 2 2 1 1 2 1 1 2 Assuming that the frequency resonance condition is satisfied, i.e.
x x x 1 1 2 2 2 2 2 we get: x x x y y y k 12 1 2 1 2 1 1 2 2 2 2 2 Again, we emphasise that this symmetric form of V k 12 is only valid on the resonant manifold, i.e. for weak WT (kinetic or discrete).

Quadratic invariants
In the large-scale limit, the zonostrophy density becomes : The densities of invariants E and Ω are w k (=k k x 2 in this case) and k x respectively. It has recently been discovered in Saito and Ishioka (2013) that expression (4.6) arises in r ( ) O 1 in the Taylor expansion of (3.8) for r  0, whereas another previously unknown invariant arises in r ( ) O . 0 This additional invariant is defined as: Both invariants ϒ and Φ are conserved independently for w = k k x k 2 since this expression for the frequency is valid up to order r 2 corrections in the original expression (2.5). The new expression has since been named semiaction in Connaughton et al (2015) because its density coincides with the one of waveaction in the sector where it is not zero. It is worth mentioning that even though it is tempting to discuss the physical meaning of the zonostrophy and semi-action, it is still unclear and would be premature. This is an interesting problem for the future. Figure 1 shows the distribution of j k in 2D wavenumber space. The dividing line is x which is the zonal region with zonal (Z) modes and j k is equal to one outside this region ( where the modes are meridional (M) modes. The dividing line is very important in the proposition that follows because of the structure of the manifold for the large-scale limit-for example the zonostrophy invariant is singular there. As can be seen from our proof of the proposition in appendix A, the dividing line plays an important role in showing certain triads are prohibited.

Prohibited triads
In order for the semi-action (4.8) to be conserved it must satisfy the resonance condition (3.6). So we have = + 1 1 0for triads of type « + M M Z and = + 0 0 0for « + Z Z Z. In other words, an excitation that is in a M-mode can be transferred to a M-and a Z-mode but not to two M-modes ¹ + ( ) 1 1 1 or two Z-modes ¹ + ( ) 1 0 0 otherwise (3.6) would not be satisfied. On the other hand, an excitation that is in a Z-mode can only move to other Z-modes and not to M-modes. This behaviour follows from the conservation of semi-action, but we would like to prove the following proposition directly.
x be the set of zonal modes of the system with frequency w = k k x k 2 . Then the following triad processes are prohibited: Depending on whether we are considering the kinetic or discrete regime, either one or several of the following triad processes may be realised, « + M M Z and/or « + Z Z Z. This proposition is proven in appendix A. Note that part (4) of the proposition trivially follows from the wavenumber condition alone-it is included here solely for completeness.

Cascade directions
It was shown in Balk et al (1991) for drift/Rossby WT in the case of zonally dominated turbulence, | | | | k k y x  , and in the case of large-scale turbulence, rk 1  , that the presence of an additional quadratic invariant (zonostrophy) allows us to predict the directions of fluxes in the k-space of the three quadratic invariants. This is similar to the famous Fjørtoft's argument for 2D hydrodynamic turbulence where the presence of enstrophy is shown to forbid energy to flow to small scales, and the presence of energy-to forbid enstrophy flow to large scales (Fjørtoft 1953). In Nazarenko and Quinn (2009) this argument was extended to small-scale drift/Rossby turbulence, rk 1  . Cascade boundaries were found for the energy, enstrophy and zonostrophy: it was predicted that each of the invariants was forced by the other two to cascade into its own anisotropic sector of k-space. A numerical study of the CHM equation was carried out which confirmed this prediction.
The most straightforward application of Fjørtoft's argument can be done to strictly positive quadratic invariants. This boils down to saying that an invariant is not allowed to dissipate in parts of the k-space where its spectral density is much less than the spectral density of at least one other invariant-otherwise the latter invariant would have to be dissipated much faster than it is produced at the forcing scales, which is impossible. This splits the entire k-space into non-intersecting sectors to which the considered invariants are allowed to cascade. Under some circumstances Fjørtoft's argument can be extended to invariants whose density may change sign, provided that such invariants remain sign-definite in their own cascade sector of the k-space.
Below, we will return to Fjørtoft's argument for the large-scale drift/Rossby turbulence, and revise the results of the Balk et al (1991) paper by taking into account yet another quadratic invariant-the semi-action Φ. We will modify Fjørtoft's argument (usually formulated for a stationary forced and dissipated system) adopting it to evolving turbulence in a non-dissipative system. This is because our subsequent numerical simulations will be precisely in such an evolving non-dissipative set-up. The picture is qualitatively different depending on which sector the initial spectrum is in-zonal or meridional. Therefore, we will consider these two cases separately.

Initial spectrum in the zonal sector
First, let us put the initial spectrum (at t=0) in the zonal sector near a wavenumber = Î ). Since the semi-action density is zero in the Z-modes, no turbulence is allowed to leave this sector, as this would mean that initially zero semi-action Φ would become finite, contradicting its conservation. (Recall that there is no resonance triads of type + « Z Z M.) But the zonostophy ϒ in the large-scale limit (4.6) is positive in the Z-modes.
Therefore, ϒ can be used in Fjørtoft's argument along with the other two positive invariants, E and Ω (Φ will not be involved in this argument as it is zero in this case). Namely, let us suppose ad absurdum that at > t 0 a significant proportion of a particular invariant (of the order of its total value) has moved into the vicinity of a scale k where the density of this invariant is much less than the density of at least one of the other two invariants. But then the amount of such an invariant with the dominant density would have at the scales around k an amount which is greatly exceeding its total initial value. This would contradict conservation of this invariant and, is therefore, impossible.
Thus, the entire k-space could be divided into three sectors to which respective invariants are allowed to flow. The boundaries of these sectors are 'soft' in a sense that the invariants are allowed to cross into each other's sector but not too deeply. These boundaries can be found by equating the ratios of densities for each pair of invariants to their initial value, i.e.: 3 . The resulting cascade picture is summarised in figure 2. It can been seen that the energy cascades to large k, the enstrophy to small k, both becoming progressively more zonal, and the zonostrophy flows towards the Z/M boundary = k k 3 y x (without crossing it).

Initial spectrum in the meridional sector
Let us now put the initial spectrum (at t=0) in the meridional sector M near a wavenumber ). It is clear that invariant Φ must remain in M because its density is zero outside of this sector. However, E and Ω are free to flow from M to Z. In this case zononstophy ϒ is not signdefinite and, therefore, cannot restrict the fluxes of the other quadratic invariants. On the other hand, semiaction Φ is now non-zero and positive and, therefore, must be used in Fjørtoft's argument along with the other two positive invariants, E and Ω.
In this case, the cascade boundaries obtained by the pairwise equating of the invariant densities are: The W E boundary separates the E-and the Ω-cascades: as before it means that E cannot be transferred to k k 0  and Ω-to k k 0  . Further, since the F E and the WF boundaries are only in M, we have So, because the boundaries are 'soft', we can takek k x , and find that all three boundaries, W E , F E and WF, approximately coincide in M. This means that Φ cannot flow to large k, and, assuming that it should move far from the initial scale k 0 , it must move to modes with small k (while remaining in M). On the other hand, Ω cannot be transferred to neither large nor small wave vectors in M. The only remaining choice for Ω is to flow to small wave vectors in Z. The latter choice does not contradict conservation of Φ because its density is zero in Z. A summary for the allowed cascade directions for this case is given in figure 3.

Numerics
In order to test numerically our theoretical predictions, a pseudo-spectral code using a third-order Runge-Kutta time integration algorithm, originally used for the small-scale case in Nazarenko and Quinn (2009), has been used after adapting it to the large-scale limit-(3.1) with frequency w = k k x k 2 and interaction coefficient (4.4). Similarly to the original set-up, we use initial condition: -a Gaussian spot centred at k 0 with width k * and its mirror image with respect to the k x -axis. Phases f k are are chosen to be random using a random number generator and independent and A is a constant controlling the level of nonlinearity. For simplicity we take the size of the periodic box as p = L 2 . We will consider cases when the initial spectrum is in the zonal sector Z and in the meridional sector M (the Gaussian is suitably truncated to ensure that the initial condition is fully concentrated in one sector only), and cases with both weak and strong nonlinearity. Our theoretical set up is such that interactions are weak. However, more typically nonlinearity is not always weak in all of the k-space and may become strong in some isolated parts of it. Therefore numerics are important to check if our prediction holds when we do not have a purely weakly nonlinear system.
We follow the motion of the invariants in the k-space by tracking the paths followed by their centroids defined as a mean k weighted on the density of the respective invariant, i.e.: Here we took into account (4.3) according to which y = |ˆ| n k x k k 2 (remember that p = L 2 ). The centroids mark the position k around which most of the respective invariant is concentrated in the k-space at time t.

Initial spectrum in the zonal sector
Let us put the initial spectrum in the zonal sector and start with weak nonlinearity. We take the following parameters: 20, 65 , 8, 10, 5 10 0 5 and a resolution of 512 2 . Figure 4 shows the evolution of the total energy, enstrophy and zonostrophy. Here, time t is normalised to the period of the mode k 0 , i.e. p w = ( ) T k 2 0 . It can be seen that the energy, enstrophy and zonostrophy are conserved to within 1.7%, 0.3% and 18% respectively. Conservation of the energy and enstrophy is a good test of the numerical method because these are exact invariants of the underlying equations in the periodic box for any level of nonlinearity. Slightly poorer conservation of the energy is expected since this quantity has a higher contribution from the large k-modes which evolve faster and, therefore, are more sensitive to the error due to a finite time step. On the other hand, zonostrophy is an approximate invariant whose conservation depends on the nonlinearity level.
For the resolution we initially used 256 2 and found that the results are very similar suggesting that there is not much sensitivity to resolution, except for better conservation of the energy and the potential enstrophy at 512 . Although it may seem that this resolution is modest for theses days, to simulate weakly nonlinear systems it is not as they take a long time, evolving very slowly. Figure 5 shows the ratio of the semi-action to the total action in the system. This ratio is chosen as an indicator of the semi-action conservation, obviously, considering the relative change of this invariant based on its initial value would be meaningless as the latter is zero. From this plot we see that the prediction that weak turbulence initially contained if the zonal sector will remain there indefinitely holds with very high accuracy. Indeed, less than 0.6% of action escapes into the meridional sector. Figure 6 shows paths of the centroids of the energy, enstrophy and zonostrophy normalised to the initial position k 0 . One can see that the invariants move into the sectors predicted in section 4.3.1.
The ψ-spectrum in 2D k-space is shown in figure 7. As previously seen in figure 5, the spectrum remains in the zonal sector. The spectrum evolves towards the origin with zonal flow forming at two distinct placesstronger one in a low k region and a weaker one with high kʼs. The latter forms near the zonal scale~( ) k k 0, y 0 in agreement with a theoretical prediction of Nazarenko (1991). In our next run the amplitude was increased to =´-A 1 10 3 so as to make the system strongly nonlinear. Figure 8 shows evolution of the total energy and enstrophy. Like before, we see that the enstrophy is conserved considerably better than the energy. On the other hand, zonostrophy is not conserved at all and, therefore, not shown (its variations exceed 300% of its initial value). Figure 9 shows the ratio of the semi-action to the total action in the system. We can see that the semi-action is not conserved, and turbulence is no longer contained in the zonal sector.
Centroid paths of the energy and enstrophy are shown in figure 10. As the zonostropy is not conserved, it is not expected to restrict the energy and enstrophy fluxes, and, therefore, we are not showing its centroid. On the other hand, the energy and enstrophy still restrict transfer directions of each other in accordance with the predictions of Fjørtoft's argument (which reduces to the standard argument as in 2D turbulence in this case).
Drastic differences with the weakly nonlinear system are particularly evident on the 2D spectra shown in figure 11. Now that nonlinearity is strong and three-wave resonances are no longer dominant, proposition 1 no Figure 4. Evolution of the energy, enstrophy and zonostrophy when the initial spectrum is in the zonal sector and nonlinearity is weak. Figure 5. Plot showing the ratio of action n k in the meridional sector (semi-action) to that in the whole of k-space when the initial spectrum is in the zonal sector and nonlinearity is weak. Figure 6. Centroid paths of the energy, enstrophy and zonostrophy when the initial spectrum is in the zonal sector and nonlinearity is weak. Figure 7. Three frames of the ψ-spectrum in 2D k-space when the initial spectrum is in the zonal sector and nonlinearity is weak. Figure 8. Evolution of the energy and enstrophy when the initial spectrum is in the zonal sector and nonlinearity is strong. longer confines turbulence to the zonal sector only. It can be seen that the spectrum moves into the meridional sector, initially exciting near-meridional modes with wavenumbers close to the sum of the dominant wavenumbers in the initial Gaussian and its image. Subsequent evolution is characterised by isotropisation of the spectra.

Initial spectrum in the meridional sector
Let us now put the initial spectrum in the meridional sector. Namely, we take the initial spectrum as in (5.1) with = Î ( ) M k 40, 20 0 and with a width of = * k 8. In our next run the initial amplitude was chosen as =´-A 1.875 10 4 so as to make the run weakly nonlinear. The resolution was 256 2 . Figure 12 shows the evolution of the energy, enstrophy, zonostrophy and semi-action. The energy, enstrophy and semi-action are well conserved to within 3%, 1% and 7% respectively. The zonostrophy is also conserved initially, but its conservation breaks down rather abruptly in the second half the runtime. This probably happens when the turbulent spectrum reaches the zonal/meridional boundary at which the zonostrophy density (4.6) becomes singular, which has a detrimental effect for the conservation conditions. Figure 13 shows the cascade directions of energy E, enstrophy Ω and semi-action Φ plotted in terms of the centroid paths. It can be seen that each invariant cascades in the direction predicted in section 4.3.2. However, Ω has not reached its designated low-k corner of the zonal sector. Interestingly, to get to this sector the centroid of Ω has to cross the region of low-k meridional scales which is forbidden in a sense that no significant amount of Ω can be present in it at any time. But this means that the transfer occurs nonlocally-passing directly from the initial to the destination scales. At intermediate times this would correspond to two spots in the 2D k-space with Figure 9. Plot showing the ratio of action n k in the meridional sector to that in the whole of k-space when the initial spectrum is in the zonal sector and nonlinearity is strong. Figure 10. Centroid paths of the energy and enstrophy when the initial spectrum is in the zonal sector and nonlinearity is strong. Figure 11. Three frames of the ψ-spectrum in 2D k-space when the initial spectrum is in the zonal sector and nonlinearity is strong. Figure 12. Evolution of the energy, enstrophy, zonostrophy and semi-action when the initial spectrum is in the meridional sector and nonlinearity is weak. Figure 13. Centroid paths of the energy, enstrophy and semi-action when the initial spectrum is in the meridional sector and nonlinearity is weak. significant concentration of Ω, both of which would contribute to the centroid of Ω placing it somewhere in between, i.e. in meridional scales. Figure 14 shows three frames of the ψ-spectrum in 2D k-space at the start, during and at the end of the run. As described in the previous paragraph, we see signs of nonlocal excitation of turbulence bypassing intermediate regions of the k-space. The scales excited are actually on the Z/M boundary, which is allowed by Fjørtoft's argument (recall that the cascade boundaries are 'soft'). Evolution of initially meridional spectrum toward the Z/M boundary was previously reported by Saito and Ishioka (2013).
To see what happens when nonlinearity is increased, we perform a run with initial amplitude =´-A 1 10 3 , and with the rest of the parameters as before. Figures 15-17 show the evolution of the invariants, cascade directions and the ψ-spectrum respectively. The zonostrophy line has been removed from the evolution plot as it quickly looses conservation. The energy, enstrophy and semi-action are conserved to within 3.5%, 1% and 17% respectively. A word of caution is due about the 17% conservation of the semi-action, which might not seem bad at the first sight. Imagine a situation when the final spectrum is fully isotropic: the waveaction loss into the zonal sector would only be 30% in this case.
Secondly, one has to be cautious interpreting the cascade directions in figure 16 which seem to be quite similar to the weakly nonlinear case. In fact, the spectrum evolution seen in figure 17 is drastically different than before: the spectrum evolution is more gradual/local and with a strong tendency to isotropisation, which is expected for strongly nonlinear systems. It is this isotropisation that explains the motion of the centroids in figure 16 (along with the usual tendency for E and Ω to go to the large and small kʼs respectively). Figure 14. Three frames of the ψ-spectrum in 2D k-space when the initial spectrum is in the meridional sector and nonlinearity is weak. Figure 15. Evolution of the energy, enstrophy and semi-action when the initial spectrum is in the meridional sector and nonlinearity is strong.

Finite-ρ effects
So far we have studied the large-scale limit (r  0) for which we proved proposition 1 about the forbidden triads. However, we would also like to see what happens when ρ is small but finite. It turns out that some + « Z Z M triads which are forbidden for r  0 start to appear. They contain wave vectors close to the Z/M boundary for which the deviation of k k y x from 3 shrinks as r  0. In other words, as ρ increases, the angle containing forbidden triads around the Z/M boundary also increases. See below. Fixing q q = ( ) k cos , sin in the resonant conditions (2.12), in figure 18 we plot k k 3 y x for resonant + « Z Z M triads as a function of r 2 for three different values of θ. We see that the closer θ is to the boundary, the bigger the angle containing + « Z Z M triads. We have further found that the triads + « M M Mand + « M Z Z(and obviously + « M M Z) forbidden in r  0 limit, continue to be absent for finite but small ρ.
To check if our numerical results obtained for the large-scale limit are robust when ρ is increased, we performed simulations of the full CHM equation (2.1) with = F 100 000. For our typical modes this corresponds to rk in the range from 0.1 to 0.4, which is not so small for formal validity of the r  k 0 limit. We Figure 16. Centroid paths of the energy, enstrophy and semi-action when the initial spectrum is in the meridional sector and nonlinearity is strong. Figure 17. Three frames of the ψ-spectrum in 2D k-space when the initial spectrum is in the meridional sector and nonlinearity is strong.
use the full CHM equation (2.4) in Fourier space with frequency (2.5) and interaction coefficient (2.6). Again, we considered a set of four cases when the initial spectrum was in the zonal and meridional sectors, and also when the nonlinearity was both weak and strong. The results of this study are reported in appendix B: they reproduce all of the qualitative features of the results obtained for the large-scale limit. This demonstrates that, at least qualitatively, the mechanisms discovered for the r  0 limit remain at work for large-scale turbulence with r  k 0.4.

Conclusion
For the large-scale limit of the CHM equation, a new quadratic invariant has recently been discovered in Saito and Ishioka (2013) for weakly nonlinear systems; it is called the semi-action Φ in the present paper-see (4.7) and (4.8). Prompted by the fact of its conservation, we have proposed and proven that the following resonant triads are prohibited, , and « + Z M M, where Z and M are the zonal and meridional sets of modes defined in proposition 1. This proposition has a drastic consequence for the weakly nonlinear dynamics of the large-scale CHM systems: spectrum initially fully concentrated in the zonal sector Z cannot ever leave this sector (no M-modes can be excited).
Another additional quadratic invariant of the large-scale CHM equation, known since 1990 (Balk et al 1990, Nazarenko 1990-is zonostropy ϒ defined in (3.7) and (4.6). In these papers it was used in a generalised Fjørtoft's argument resulting in a triple cascade picture of anisotropic drift/Rossby turbulence. In the present paper, we have used the newly discovered semi-action invariant to revise Fjørtoft's argument. The latter is now presented in a modified version for evolving non-dissipative systems assuming that the quadratic invariants will eventually move to scales which are greatly separated from the initial scale k 0 . If the initial spectrum is in Z then it remains in Z, but the invariants are transferred among the zonal scales. Namely, the energy cascades to large zonal kʼs, the enstrophy-to small zonal kʼs, and the zonostrophy-towards the the Z/M boundary, = k k 3 y x . If the initial spectrum is in M then energy cascades towards large kʼs, enstrophyto small zonal kʼs, and semi-action-to small meridional kʼs (by definition, the latter cannot leave M).
We have tested and confirmed our theoretical predictions numerically. In particular, we have confirmed conservation of the semi-action and zonostrophy invariant. We also confirmed that when nonlinearity is weak (three-wave interactions dominate) turbulence which is initially in Z remains in Z with remarkable a 0.6% accuracy. In this case, it is redistributed within the zonal sector with the energy, enstrophy and zonostrophy moving in the 2D k-space as predicted by Fjørtoft's argument. The 2D spectrum develops two distinct components: the main spot that spreads to zonal scales with smaller kʼs and a smaller spot forming with large-k zonal scales near~( ) k k 0, y 0 (formation of the latter was predicted in Nazarenko (1991)). When nonlinearity is weak and turbulence is initially in M, the spectrum tends to move towards and concentrate near the Z/M boundary, = k k 3 . y x Such a behaviour was previously reported in Saito and Ishioka (2013). The cascade directions appear to be qualitative predictions of Fjørtoft's argument. Notably, conservation of the zonostrophy, which holds initially, is suddenly broken down when the spectrum reaches the Z/M boundary on which the density of this invariant is singular.
When nonlinearity is increased, as expected, for both zonal and meridional initial conditions, conservation of the semi-action and (especially) zonostropy deteriorates. In particular, initially zonal turbulence is no longer contained in the zonal sector. Turbulence evolution shows tendency to isotropisation, as expected for strongly nonlinear CHM, because the linear term is the only source of anisotropy in this model. In general, the system evolves similar to the classical 2D turbulence, but with reversal of roles of the energy and the enstrophy-they cascade to the small and large scales respectively in our case.
We have also studied sensitivity of our results to increasing values of ρ, so that formally the CHM system is not in the large-scale limit. We showed that even for small values of ρ, some resonant triads of type « + M Z Z (but not the other types forbidden in the large-scale limit) exist. They appear with wavenumbers close to the Z/M boundary, and the sector in which they exist grows as ρ increases. Thus, for finite but small ρ, 'leakage' between the Z and M sectors occurs primarily near the Z/M boundary only. Using direct numerical simulations of the CHM equation, we have found that most of the qualitative and even quantitative features of the large-scale limit survive for our numerical set-ups up to values r  k 0.4.

Acknowledgments
Katie Louise Harper gratefully acknowledges EPSRC DTA funding of her PhD studies.
Appendix A. Proof of proposition 1 2 2 Let us begin by writing out the frequency resonance condition, w w w = + 1 2 , in which we substitute the wavenumber resonance condition, = + k k k 1 2 : Taking > q 0 we can see that > and =´-A 7.5 10 9 . Figures B1-B3 show the conservation and cascade directions of the energy, enstrophy and zonostrophy and the ψ-spectrum respectively.
For the strongly nonlinear run we increased the amplitude to =´-A 1 10 7 . The numerical plots are presented in figures B4-B6.
For the weakly nonlinear run with a meridional initial condition the resolution was b = = = ( ) * k k 256 , 40, 20 , and =´-A 1.875 10 8 . Figures B7-B9 show the conservation of the energy, enstrophy, zonostrophy and semi-action, the cascade directions of the energy, enstrophy and semiaction and the 2D ψ-spectrum respectively.
For the strongly nonlinear run figures B10-B12, we increased the amplitude to =´-A 1 10 6 and kept the remaining parameters the same.
We see a remarkable agreement of the above plots obtained for both weak and strong initial conditions in both zonal and meridional sectors with the respective plots reported in the main text of the large-scale limit simulations. Figure A2. Graph summarising the case when < q p 3 . Figure B1. Evolution of the energy, enstrophy and zonostrophy when the initial spectrum is in the zonal sector, nonlinearity is weak and ρ is small but finite. Figure B2. Centroid paths of the energy, enstrophy and zonostrophy when the initial spectrum is in the zonal sector, nonlinearity is weak and ρ is small but finite. Figure B3. Three frames of the ψ-spectrum in 2D k-space when the initial spectrum is in the zonal sector, nonlinearity is weak and ρ is small but finite.

20
New J. Phys. 18 (2016) 085008 K L Harper and S V Nazarenko Figure B4. Evolution of the energy and enstrophy when the initial spectrum is in the zonal sector, nonlinearity is strong and ρ is small but finite. Figure B5. Centroid paths of the energy and enstrophy when the initial spectrum is in the zonal sector, nonlinearity is strong and ρ is small but finite. Figure B6. Three frames of the ψ-spectrum in 2D k-space when the initial spectrum is in the zonal sector, nonlinearity is strong and ρ is small but finite. Figure B7. Evolution of the energy, enstrophy, zonostrophy and semi-action when the initial spectrum is in the meridional sector, nonlinearity is weak and ρ is small but finite. Figure B8. Centroid paths of the energy, enstrophy and semi-action when the initial spectrum is in the meridional sector, nonlinearity is weak and ρ is small but finite. Figure B9. Three frames of the ψ-spectrum in 2D k-space when the initial spectrum is in the meridional sector, nonlinearity is weak and ρ is small but finite. Figure B10. Evolution of the energy, enstrophy and semi-action when the initial spectrum is in the meridional sector, nonlinearity is strong and ρ is small but finite. Figure B11. Centroid paths of the energy, enstrophy and semi-action when the initial spectrum is in the meridional sector, nonlinearity is strong and ρ is small but finite. Figure B12. Three frames of the ψ-spectrum in 2D k-space when the initial spectrum is in the meridional sector, nonlinearity is strong and ρ is small but finite.