Zeitlin Truncation of a Shallow Water Quasi‐Geostrophic Model for Planetary Flow

In this work, we consider a Shallow‐Water Quasi Geostrophic equation on the sphere, as a model for global large‐scale atmospheric dynamics. This equation, previously studied by Verkley (2009, https://doi.org/10.1175/2008jas2837.1) and Schubert et al. (2009, https://doi.org/10.3894/james.2009.1.2), possesses a rich geometric structure, called Lie‐Poisson, and admits an infinite number of conserved quantities, called Casimirs. In this paper, we develop a Casimir preserving numerical method for long‐time simulations of this equation. The method develops in two steps: first, we construct an N‐dimensional Lie‐Poisson system that converges to the continuous one in the limit N → ∞; second, we integrate in time the finite‐dimensional system using an isospectral time integrator, developed by Modin and Viviani (2020, https://doi.org/10.1017/jfm.2019.944). We demonstrate the efficacy of this computational method by simulating a flow on the entire sphere for different values of the Lamb parameter. We particularly focus on rotation‐induced effects, such as the formation of jets. In agreement with shallow water models of the atmosphere, we observe the formation of robust latitudinal jets and a decrease in the zonal wind amplitude with latitude. Furthermore, spectra of the kinetic energy are computed as a point of reference for future studies.


Introduction
In the context of Geophysical Fluid Dynamics, the rotating stratified Euler equations constitute an accurate model governing the dynamics of oceanic and atmospheric motions.Due to their complexity, simpler models are often used to describe features of geophysical flows.Important examples of these are the Primitive Equations (PE) and shallow water models.In particular, the use of the Rotating Shallow Water (RSW) equations on the sphere has played a pivotal role in advancing our understanding of large-scale oceanic and atmospheric flows.The RSW equations describe a single shallow layer of fluid on a rotating sphere under hydrostatic pressure.The equations can be derived from the PE under the assumption of constant density and columnar motion and the subsequent vertical averaging of pressure gradients.The relation between these models can be found in more detail in Zeitlin [2018] and Luesink [2021].
From the one-layer global shallow-water model of the atmosphere, a simplified equation has been independently derived by Verkley [2009] and Schubert et al. [2009], based on previous works of Charney and Stern [1962] and Kuo [1959].This equation, which we will refer to as the Balanced Shallow Water (BSW) equation, is derived from RSW using a non-divergence assumption for the horizontal velocity field, together with a linear balance equation linking the stream function and the geopotential.This simplifying procedure leads to an equation for the potential vorticity that is identical in form to the rotating Euler equation for the relative vorticity, except for an additional f 2 -dependent term, where f is the Coriolis parameter.To validate the BSW model, Verkley considers the dynamics of linearized Rossby waves and shows that these waves accurately reproduce the westward-propagating waves of the second class of the original shallow-water equations, studied in Longuet-Higgins [1968].Moreover, the associated 2D turbulence is studied in Verkley [2009] and in Schubert et al. [2009], where the anisotropic Rhines barrier is derived.
The geometric structure is a property that the BSW equation and the rotating Euler equations share.As we will show, both equations form a Lie-Poisson system and admit an infinite number of independent conserved quantities, called Casimirs.These Casimirs can be expressed as integrated powers of potential vorticity, of which the integrated squared potential vorticity is called the enstrophy.
In a dynamical system, the presence of conserved quantities has profound physical consequences.If we consider the role of enstrophy in decaying 2D turbulence in the context of fluid dynamics, the well-known physical consequence is the inverse energy cascade.This phenomenon is opposite to the forward energy cascade that is observed in 3D turbulence, where enstrophy is not conserved (for the inviscid dynamics).
The presence of conserved Casimirs for 2D flows similarly has implications on the long-time behaviour of the fluid and for the formation of large-scale coherent structures (Bouchet and Venaille [2012]).These phenomena are sometimes explained through the equilibrium statistical theory.This theory typically makes use of some conserved quantities of the underlying (inviscid) dynamics to make predictions.For example, Kraichnan [1975], Salmon et al. [1976] and Carnevale and Frederiksen [1987], derive equilibrium statistical states for the 2D Euler and the quasi-geostrophic equations using a Fourier-spectral truncation of the equations of motion.This spectral truncation preserves circulation, energy and enstrophy, but not the Casimirs of order higher than two.Therefore, it is of interest to investigate whether higher-order Casimirs are statistically relevant for geophysical equations.This topic has been addressed in Abramov and Majda [2003] for the barotropic two-dimensional flow with topography on a flat geometry (see also Dubinkina and Frank [2010]).In this work the authors use a Casimir preserving numerical method to integrate the equations of motion.In particular, they prove the statistical relevance of the third order Casimir C 3 by showing an increasing discrepancy with the energy-enstrophy statistical theory (which predicts collinearity between the time-averaged streamfunction and potential vorticity) as a function of C 3 .We remark that, in order to study numerically the statistical relevance of high-order Casimirs in geophysical flows, it is necessary to have a numerical method that preserves the Casimirs in the discrete model.
Motivated by these considerations, the present work introduces a Casimir preserving numerical method for the BSW model on the sphere, based on Zeitlin discretization Zeitlin [1991Zeitlin [ , 2004]], and on an isospectral time-integrator developed by Modin and Viviani [2020].This method was recently applied successfully to underpin Kraichnan's double scaling hypothesis for 2D forced turbulence Cifani et al. [2022] and is also adopted here.In this paper, we show that the combination of Zeitlin discretization and isospectral time-integration for the BSW flow on a sphere captures two important physical phenomena: the first is the double cascade in the kinetic energy spectrum as is typical for 2D turbulence; the second is the emergence of zonal jets across the sphere.In particular, the intensity of these zonal jets is not uniform across latitudes.We observe that the strong equatorial jets are attenuated towards the poles, and we will show that the intensity of this attenuation is regulated by a single parameter in this model, called Lamb's parameter.We remark that the attenuation of zonal jets, which is seen in physical systems such as the atmospheres of giant gaseous planets, has been observed in the simulations of RSW equations, but not in simulations of the 2D Navier-Stokes equations (see the work of Scott and Polvani [2007] for the simulations of RSW where the attenuation of zonal jets is regulated by changing the Rossby deformation radius, similar to our case).For this reason, the BSW model may be a valuable tool for studying the formation and attenuation of zonal jets.
The paper is organized as follows.In Section 2 the mathematical model for BSW flow on a sphere is presented.Section 3 introduces the Zeitlin truncation and the Lie-Poisson time-integration method.At the end of this section we simulate an inviscid unforced BSW flow to demonstrate the numerical conservation of energy and the Casimirs.Simulations of forced BSW flow in the presence of dissipation are performed and discussed in Section 4 and concluding remarks are collected in Section 5.

BSW equations on the sphere
The one-layer rotating shallow water equation on the sphere is a model that, under the assumptions of constant density and columnar motion, describes the motion of a thin layer of fluid on the surface of a rotating sphere.The fluid is described by a two-dimensional horizontal (or tangential) velocity field u(ϕ, θ) = (u, v) and a local height h, where ϕ and θ are the latitude and longitude respectively.This system admits a materially conserved quantity called potential vorticity: where ω := r • (∇ × u) is the radial component of the relative vorticity and f = 2Ω sin ϕ is the Coriolis parameter.Moreover, Ω denotes the angular velocity of the planet.
Starting from the one-layer RSW equations on the sphere, Verkley [2009] and Schubert et al. [2009] derived the simplified BSW model, by adding three assumptions as follows: Assumption 1: We denote the local height of the fluid column by where H is a uniform average height and η is a variable surface elevation.The quantity b denotes the topography and in the following will be assumed to be zero (b(x) = 0).The first assumption consists of assuming that the variable surface elevation is much smaller than the average height, i.e., η/H ≪ 1 (3) Assumption 2: Motivated by this first assumption and by the continuity equation for RSW, we assume that the divergence of the horizontal velocity can be approximated as zero.Using the Helmholtz decomposition, we write the horizontal velocity in terms of a stream function ψ such that u = r × ∇ψ (4) ∂θ θ the horizontal gradient on the surface of the sphere with radius R.
Assumption 3: In order to derive a closed dynamical system from the material conservation of the potential vorticity (1), we propose a type of balance that relates the velocity field u to the local height h, or, equivalently to the surface elevation η.Among the linear balance relations, we choose the simplest one, which relates the stream function and the variable height as follows: where g is the gravitational acceleration.This relation is inspired by quasi-geostrophic systems, in which there exists an approximate balance between the Coriolis force and pressure forces.It was first considered by Daley [1983], as a simplification of the balance relation ∇ Using Assumptions 1, 2 and 3 we can derive the BSW system from the RSW point of reference.
Since we assume a small free surface elevation, an expansion up to first order in η/H yields so that the potential vorticity becomes In the regime of rapid rotation, the dominant contribution to the potential vorticity according to Verkley [2009] reads defining the central dynamic variable q to which we turn next.By using the linear balance relation (5) we write the potential vorticity q as: In this approximation, we can formulate the material conservation law for the potential vorticity (1) as where where ∆ := ∇ • ∇ is the Laplace-Beltrami operator on the surface of the sphere.
Equations ( 10) and ( 11) are the central formulation of the BSW model on the sphere.Anticipating the Zeitlin truncation Zeitlin [2004] in the next section, we non-dimensionalize the equations using the radius of the sphere R as the length scale, and Ω −1 as the time scale.Furthermore, we denote the convective term in (10) as The binary operation {•, •} forms a Poisson bracket on the space of smooth functions on the sphere (i.e. the operation is bilinear, skew-symmetric and respects Jacobi and Leibniz identities Arnold and Khesin [1992]).Thus, equations ( 10) and ( 11) can be reformulated as q = {ψ, q} where µ = sin ϕ and γ is Lamb's parameter, which expresses the ratio between the radius of the sphere, and the typical size of quasi-geostrophic vortices Vallis [2019] given by the Rossby deformation length R d : We now turn to the geometric properties of BSW equations in (13).As anticipated, the BSW model possesses a rich geometric structure, which it inherits from the RSW model.Indeed, equations ( 13) constitute an infinite-dimensional Lie-Poisson system Marsden and Weinstein [1983], Arnold and Khesin [1992] on the space of smooth functions on the sphere C ∞ (S 2 ).In other words, the BSW equations have a Hamiltonian structure with an additional Lie algebra property.
The Hamiltonian structure is defined as follows.Firstly, the space C ∞ (S 2 ), i.e. the space of vorticity fields, forms the infinite dimensional phase space for (13).Secondly, the Poisson bracket in ( 12) defines a new infinite dimensional Poisson bracket {•, •} LP that acts on functionals defined on the phase space as follows: {F, G} LP (q) := for smooth functionals F, G : Here dA = − cos ϕdϕdθ is the infinitesimal surface element, and δ(•)/δq denotes the variational derivative.When, as in Verkley [2009], we define the Hamiltonian of BSW as the potential vorticity equation ( 13) can be written in equivalently as In fact, in the case F(q) = q(⃗ x), for ⃗ x ∈ S 2 a given point on the sphere, we recover (13).
The additional Lie algebra property arises from the following observation.Since the phase space C ∞ (S 2 ) is a vector space, the couple (C ∞ (S 2 ), {•, •}) forms an infinite dimensional Lie algebra Olver [1993].This additional property implies the existence of special functionals These functionals are called Casimirs and, due to equation ( 17), are constants of motion.The Casimir functions for the BSW model are an infinite family and are given, for any smooth function ξ ∈ C ∞ (R), by In practice, the functions ξ(q) can be chosen to be monomials, so the corresponding Casimirs read The conservation of Casimirs for BSW can also be verified via direct computation: where the second equality uses equations ( 4) and (10).Notably, the conservation of the Casimirs is independent of the particular relation between q and ψ.In other words, their conservation is independent of the choice of the Hamiltonian, as equation ( 18) suggests.This reflects the underlying Lie-Poisson structure.In the next section, we will derive a special finite-dimensional truncation of equations ( 13) that preserves this Lie-Posson structure and, in particular, that preserves all Casimirs of the discrete model.
We conclude by commenting on the geometric structure of BSW.As mentioned, this structure is inherited from the RSW model, which also exhibits a Lie Poisson structure and the conservation of infinitely many Casimirs (see, for example, Salmon [2004] and Dellar and Salmon [2005] for a geometric formulation of RSW).This fact is not exceptional, as many other approximate dynamical equations exhibit geometric properties analogous to the exact parent cases (see Salmon [1988] and references therein).For this reason, it would be interesting to investigate the possibility of deriving the BSW model utilizing the Hamiltonian formulation of the RSW equations, and this topic will be left for future research.

Numerical method for Zeitlin truncation
In this section, we follow the approach by Zeitlin to construct a numerical scheme with which the geometric structure of the BSW equations can be preserved optimally -this procedure is sometimes referred to as 'quantization', contrasting the familiar 'discretization'.Specifically, we show how to construct and integrate a finite-dimensional analogue of equations ( 13) while preserving all the associated independent Casimirs.The first part of Subsection 3.1 is quite technical, and may be skipped upon first reading, continuing after the bullet points.

Finite truncation of the Poisson bracket
We begin our construction of the Casimir-preserving numerical method by constructing a finitedimensional analogue of ( 13).This will be referred to as the Zeitlin truncation.Fundamental to the derivation, we notice that the couple (C ∞ (S 2 ), {•, •}) forms an infinite dimensional Lie algebra and the existence of Casimirs is a direct result of this Lie-Poisson structure.Therefore, following Zeitlin, we look for a sequence of N-dimensional Lie algebras (g N , [•, •] N ) that constitute an approximation, in a sense specified later, of (C ∞ (S 2 ), {•, •}) in the limit N → ∞.Here, g N , denotes an N-dimensional vector space and [•, •] N refers to the bilinear form of the Lie algebra.The family of finite-dimensional Lie algebras (g N , [•, •] N ) is referred to as L N -approximationBordemann et al. [1991].
In Bordemann et al. [1994], it is shown that the L N -approximation for (C ∞ (S 2 ), {•, •}) can be obtained using the Lie algebra (gl(N ), [•, •]), where gl(N ) is the set of complex N × N matrices and [•, •] is the matrix commutator.Indeed, we can construct a surjective projection Π N : This allows us to construct the finite-dimensional analogue of equations ( 13).Using the projection Π N , we associate with every smooth function on the sphere S 2 a complex N × N matrix in the Lie algebra gl(N ).The action of the projection Π N is linear and explicitly known: to each spherical harmonics Y l,m ∈ C ∞ (S 2 ) we associate a matrix i T l,m ∈ gl(N ).Here i is the imaginary unit and the explicit expression of the matrices T l,m is reported in Appendix A. Thus, using the spherical harmonics as a basis, we associate with any function given by the matrix Notice that for a real-valued function q, the spherical harmonics coefficients satisfy q l,m = (−1) m q l,−m .This translates to the matrix condition Q + Q † = 0, where Q † is the conjugate transpose of Q.For this reason, we can restrict to the subalgebra of skew-Hermitian matrices Q ∈ u(N ) ⊂ gl(N ).In terms of the discrete basis, our system of equations ( 13) takes the following form: where Q is usually referred to as the potential vorticity matrix and P as the stream matrix.
Here, have replaced the Poisson brackets {•, •} in ( 13) with the scaled matrix commutator •] as given by the aforementioned projection of the Poisson bracket Modin and Viviani [2020].Moreover, we have replaced the Laplace-Beltrami operator ∆ in (13) with the quantized Laplace-Beltrami operator ∆ N , that satisfies ∆ N T l,m = −l(l + 1) T l,m , as further detailed in A.

Computation of the stream function
The scheme presented in equations ( 24) involves computation of the stream matrix P .Determining the stream matrix involves quantizing (i.e.projecting using Π N ) the functions µ and µ 2 ψ.For the projection of the former, we use the expansion in terms of spherical harmonics: The discrete representation is therefore given by the following projection: which can be calculated explicitly.The term µ 2 ψ is more involved, as it requires the projection of a product of functions.In B, we discuss a method for approximating the projection of a product of two functions on the sphere, together with its accuracy.In practice, given two real functions f and g and their associated matrices in u(N ) given by F = Π N f and G = Π N g, we use the following approximation: which associates multiplication with the anti-commutator on u(N ).In our specific case, using the expansion: we can define which enables writing: The final step which allows us to solve the system for P is the observation that the matrix harmonics T N l,0 only have non-zero entries on the diagonal (see appendix A).If we denote the diagonal of S as the vector s ∈ C N such that s i = S ii , we can write: where • is the Hadamard, or entry-wise product.
Finally, we arrive at the discrete system by substituting the above results into equations ( 24): where the matrices M and S are defined in equations ( 26) and ( 31) respectively.The equation for P can be solved efficiently using a clever splitting of the problem along diagonals of the matrix P .The implementation is based on the methodology used by Cifani et al. [2023] and is further detailed in C.
The time-evolution of Q in equation ( 32) is isospectral since the spectrum of the matrix is preserved by the matrix commutator Viviani [2020].This implies that traces of powers of the quantized potential vorticity are conserved.In other words, the quantities are conserved by (32).This represents the discrete analogue of the conserved Casimirs Cifani et al.
[2023] in equation ( 19).In this case the conservation law for the Casimirs follows directly from the cyclic property of the trace.
Additionally, the system is Lie-Poisson with the Hamiltonian H N given by: System (32) provides a finite-dimensional analogue of equations ( 13), relating the potential vorticity matrix to the stream matrix.In the next section, we integrate (32) by means of an isospectral time integrator that preserves all the Casimirs in equation (33).

Isospectral time-integration
In order to preserve the discrete Casimirs (33), we use an isospectral Lie-Poisson time-integrator following Modin and Viviani [2020].This retains the basic discrete Lie-Poisson dynamics, and, moreover, provides near-conservation of the Hamiltonian as obtained through backward error analysis Hairer et al. [2006].The time integrator is based on the implicit midpoint rule and can be written as follows: The first equation is solved for Q using fixed point iteration that creates a sequence { Qk }, with Q0 = Q n , which converges to Q as k → ∞.This contraction was found to yield rapid convergence and only a small number of iterations was found necessary to reach a set tolerance.In test calculations we may typically reach an L ∞ -norm of Qk+1 − Qk of O(10 −12 ) after only 3 iterations.
To complete a time-step, the sufficiently converged approximation QK of Q for suitable K is used once again to evaluate the full right-hand side in the second equation in (35).
The computational effort associated with the integration scheme is dominated by the evaluation of the matrix commutator.As shown in C , the stream matrix can be evaluated very efficiently from the vorticity matrix, requiring only O(N 2 ) operations.Evaluating the commutator, however, involves dense matrix multiplications, requiring O(N 3 ) operations.This shows a comparable computational complexity to that of typical pseudo-spectral methods since transformations involving spherical harmonics require at least O(N 2 log 2 (N )) operations, though typical implementations display O(N 3 ) scaling Schaeffer [2013].Additionally, the current scheme allows for an efficient MPI parallelisation as shown in Cifani et al. [2023].Similar to pseudo-spectral methods, the memory requirements for this scheme are O(N 2 ), associated with the storage of several dense complex matrices or vectors of spherical harmonics coefficients.Hence, the complexity of the new method is comparable to popular methods in literature, but in addition, the new method retains the geometric structure of the equations.

Numerical performance
To test the performance of the numerical method, we set up a numerical experiment in which the flow develops freely from an initial condition.In the absence of any forcing or dissipation, the Hamiltonian energy and all Casimirs remain constant at their respective initial values.We solve the system in equations ( 32) on the unit sphere using the integration scheme in equations ( 35).
The angular velocity of the sphere is set to Ω = 250, corresponding to approximately 40 revolutions per second, and the Lamb parameter is set at γ = 10 3 .The rotation axis is placed vertically.We simulate the system using N = 512, which represents a truncation in the expansion of spherical harmonics up to degree l = 511.The system is solved up to t = 500, corresponding to around 20 000 revolutions of the unit sphere.The time step is set at ∆t = 4 • 10 −4 , corresponding to around 60 time steps per revolution of the sphere, or, 'day'.
Figure 1 shows the evolution of the kinetic energy spectrum in time.As an initial condition, all modes between l = 40 and l = 60 are given a random amplitude around the reference value of 50/l(l + 1).As the flow evolves in time, energy spreads to other modes through nonlinear interaction, with energy going to the larger scales as well as smaller scales.After this initial mixing phase, the spectrum reaches a statistically stationary state, as in Modin and Viviani [2022].As a reference, the dashed line shows the -5/3 scaling for energy cascade in two-dimensional flows.
For any discretization of two-dimensional flows, the forward cascade of enstrophy is limited by the smallest length scales represented by the discrete system, leading in many cases to spectral blocking Boyd [2001].In this case, where the spatial resolution is determined by truncation of spherical harmonic modes on the sphere, spectral blocking occurs at the modes of the highest degree.At long simulation times, this does not seem to lead to instabilities of the numerical method and we find a similar trend as was observed for a Zeitlin truncation of the Euler equations by Modin and Viviani [2022].They report that after an initial strong cascade of enstrophy to small scales, the system develops to statistically steady dynamics.Furthermore, they show that the resulting flow can be decomposed into large-scale dynamics over a noisy background, and they study the long-time behaviour of these large-scale structures (see also Modin and Viviani [2020]).This noisy background can be linked to the presence of the aforementioned spectral blocking effect.
Figure 1: Kinetic energy spectrum of freely decaying turbulence on the sphere at several times.For reference, the slope -5/3 is shown as a dashed line.
In Figure 2 the relative drift of the Hamiltonian and the first 8 Casimirs of the discrete system are shown.The left panel shows the evolution of the error Hamiltonian, normalized by its initial value, given as H = (H(t) − H(0))/H(0).As expected, the Hamiltonian is approximately conserved, i.e., it fluctuates close to the exact value with minimal linear drift.On the right panel, the relative error in the first 8 Casimirs are shown as defined in equation ( 33) , where the relative error is given by Although exact Casimir conservation is by the isospectrality of the time-integrator, it is in practice limited by the tolerance set for the iterative solver for equations (35), which is set to 10 −12 for this simulation.The figure shows that the relative error in the Casimirs stays below 10 −10 even up to long times, with negligible drift, with Casimirs with an even index even showing conservation up to 10 −14 .
These results show that simulations of freely decaying turbulence of the BSW equation on a sphere with an energy and Casimir preserving method are indeed feasible and numerically stable.The effect of spectral blocking of the forward enstrophy cascade will be the subject of future research.
In the next section, we will focus on forced turbulence in the presence of viscous dissipation.

Simulation of zonal jet formation
To illustrate the Zeitlin truncated BSW model with Casimir-preserving iso-spectral time-integration, we simulate a particular rotating sphere case.Results obtained using the BSW equation will be compared with findings based on the Navier-Stokes equations on a rotating sphere from Cifani et al. [2022].Therefore, we consider the BSW equation with additional forcing, dissipation and (Rayleigh) friction: q = {ψ, q} + ν(∆q where ν is the dimensionless viscosity, F is the forcing and αq is the dimensionless Rayleigh friction to avoid an accumulation of energy at large scales due to the inverse energy cascade in 2D turbulence.Moreover, 2νq arises from the spherical geometry Lindborg and Nordmark [2022].Following  Cifani et al. [2022], the forcing is time-dependent, derived from a Wiener process uncorrelated to the time scales of the flow.The forcing is localized in a narrow band in spherical harmonic space around degree l f = 100.
The viscous dissipation, damping and forcing are integrated using a second-order Crank-Nicolson scheme, while the isospectral integrator is used for the convective term.This ensures that the solution departs from the level sets of the Casimirs only due to the non-conservative terms in the equation and not by numerical inaccuracies in the discretization of convection.By denoting with Φ iso,h the isospectral map for convection and by Φ CN,h the Crank-Nicolson map for the remaining terms on a time interval h, the time integration is obtained using the second-order Strang splitting: where h is the time step and n the time level.
In the following, we report the flow generated at different values of Lamb's parameter γ.In order to appreciate the qualitative changes in the flow, we chose to display three cases: γ = 0 (which corresponds to the limit R d → ∞ in which we recover the two-dimensional Navier-Stokes equations), γ = 10 3 and γ = 10 4 .
Figures 3 show a snapshot of the zonal component of the velocity field at a time at which the flow is in a statistically stationary state for each choice of parameter γ.The spatial resolution is set to N = 2048, which was found to accurately represent the smallest scales in the flow.
As a reference, we consider the solution obtained at γ = 0 as shown in the left snapshot of Figure 3.The solution displays a clear zonal structure, with the formation of a large number of zonal jets.
In particular, we notice how these jets appear to have an intensity that is quite independent of the latitude.The situation is different when we solve the BSW equation at γ 1 = 10 3 (middle snapshot in Fig. 3).Here, we see that while the number of zonal jets appears unaffected, there is an apparent graded intensity of the jets which ranges from comparably weak jets near the poles to strong jets near the equator.This trend is even stronger when Lamb's parameter is increased to γ = 10 4 (right image in Fig. 3), which corresponds to a smaller Rossby deformation length.
This effect is illustrated quantitatively in Figure 4, where the zonal velocity, averaged along longitudes, is displayed as a function of latitude.Here we clearly see how the width and intensity of jets near the equator are almost unaffected by the change in the Lamb parameter.However, we observe that the intensity of jets close to the poles gets attenuated strongly with increasing γ.This phenomenon of jet attenuation is a physical mechanism which is observed in giant gaseous planets, like Jupiter.The attenuation of zonal jets, while not present in simulations of 2D Navier-Stokes equations (as we see in Figure 4), is captured by numerical simulations of RSW equations.For example, in the work of Scott and Polvani [2007] it is shown how the equatorial confinement is present in the RSW model and it is regulated by the Rossby deformation radius, similar to the current simulations.
In Figure 5 we show the kinetic energy spectrum for the three values of Lamb's parameter, in terms of the degree l of the spherical harmonics.The spectra have been split into the zonal (m = 0) and the non-zonal (m ̸ = 0) energy contributions.First of all, we notice how, in all spectra, the non-zonal energy undergoes a double cascade, with a slope of −5/3 towards large scales and of −3 towards small scales, as expected for 2D turbulence on the sphere Lindborg and Nordmark [2022].This double cascade was also reported recently for isotropic turbulence in Cifani et al. [2022].Also, we notice how the most energetic large-scale modes come from the zonal contribution, with a peak around l = 20 for all three cases.
As we increase the value of γ, we observe two trends in the kinetic energy spectra which relate to the phenomenon of equatorial confinement of zonal jets.First, we notice that, while the degree l of the peak in the energy spectrum of the zonal modes appears unaffected by γ, the energy contained in these modes clearly diminishes as γ grows.Second, in the case where γ = 0, we notice that there is a rapid decay in the energy of the non-zonal modes as we move towards low wavenumbers l.This effect is known as the anisotropic Rhines barrier effectRhines [1975], and it is believed to be responsible for the formation of zonal jets Vallis and Maltrud [1993].This effect takes place at low wavenumbers l with m ̸ = 0, i.e. in the spectral region where the flow is dominated by waves rather than turbulence Huang and Robinson [1998].When γ ̸ = 0, the decay of energy in non-zonal modes becomes less rapid, i.e., they contain more energy than their counterparts at γ = 0.As a consequence, even if at low wavenumbers the spectrum is still dominated by zonal modes, the non-zonal modes at low l become important when γ ̸ = 0.
We conclude by noticing that a similar phenomenon has been studied in the works of Theiss Theiss [2004, 2006].In the first of these works, a generalized quasi-geostrophic equation that allows a latitude-dependent Rossby deformation length is considered.Similarly to what we find in our simulations of the BSW model, they show an equatorial confinement of zonal jets, together with an attenuation of the Rhines barrier, which can be bypassed at certain latitudes.The similarity between these numerical results can be explained by noticing that in the BSW model, the effect of a variable Rossby deformation length is naturally encoded in the factor γµ 2 in (12).

Conclusions
In this work, we developed a Casimir preserving numerical method for the BSW equations ( 13) on the sphere.Simulations of unforced turbulence in the absence of dissipative terms show that this method approximately conserves energy and conserves Casimirs with a relative error of 10 −10 .The BSW model has been simulated in the presence of dissipation and forcing.In this case, we observed that the BSW model is able to capture the phenomenon of jet attenuation, which is shown to be regulated by Lamb's parameter.
The presented method for the BSW equations can be readily extended to include a non-trivial bottom topography, since the topography only enters the model as an additional source of potential vorticity in the diagnostic equation ( 11) for the streamfunction (see Verkley [2009]).One application of the present method involves examining the statistical relevance of high-order Casimirs for the BSW model, and their influence on the coarse-grained large-scale flow, following the work of Abramov and Majda [2003].Furthermore, this modelling approach may be applied to multilayered models such as the multi-layered quasi-geostrophic equations, which will be the subject of future research.
The important property of the twisted product is that where the limit above is taken uniformly, i.e. we have uniform convergence of the sequence of functions on the l.h.s. to the function on the r.h.s Rios and Straume [2014].We remark that this asymptotic (N → +∞) expansion of the twisted product is invalid without the assumption l, l ′ ≪ N .
So we proceed as follows: we approximate the product sin 2 ϕψ by means of (41); then we project through Π N .We obtain sin 2 ϕψ → i 2 (SP + P S) where P = Π N ψ ∈ u(N ), S = Π N sin 2 ϕ ∈ u(N ).Now, we know that the Poisson brackets are approximated by the rescaled commutator N 3/2 √ 16π [•, •] up to an error O(1/N 2 ).But how well is the product sin 2 ϕψ approximated by −i 2 (P S + SP )?To compute it we need to use some asymptotics of the 6j Wigner symbol.Taking the inner product of our function product with a basis function, we use the following identity: l 1 l 2 l 3 0 0 0 On the other hand, we have for our renormalized basis Hoppe and Yau [1998] T r( T l1m1 T l2m2 T l3m3 ) = (−1) 2s (2l 1 + 1) (2l 2 + 1) (2l 3 + 1) Which indicates that the right scaling for the approximation is i 2 N 4π (−1) 2s (SP + P S), with an error of O(1/N ).As a final comment, we notice that in case we want to approximate the Poisson Brackets {Y l1,m1 , Y l2,m2 }, we have an exact expression in terms of commutators when the l 1 , l 2 , l 3 in (43) have even sum l 1 + l 2 + l 3 (See Hoppe and Yau [1998]).When the sum is odd, the first order in the expansion (45) is zero, due to the 3j Wigner symbol property, and for this reason the commutator [•, •] is eventually rescaled by a factor N 3/2 .

C Tridiagonal splitting
The discrete Laplacian ∆ N admits a decomposition in tridiagonal blocks, which enables an efficient implementation of the associated eigenvalue problemCifani et al. [2023].In this appendix, we show that the same splitting applies when solving the for the stream matrix P in equation 11.Since the matrices T lm are eigenvectors of the discrete Laplacian and have entries only in the m th diagonals,

Figure 2 :
Figure 2: Time-evolution of the relative error H in the Hamiltonian as a function of time (left) and the error in the first 8 Casimirs of the discrete system, also normalized by their respective initial values (right).The Casimirs are grouped by the parity of their index.Casimirs with an odd index are clustered at the top of the figure, while Casimirs with an even index appear at the bottom of the figure.

Figure 3 :
Figure 3: Zonal velocity of statistically stationary BSW flow at γ = 0 (left), γ = 10 3 (middle) and γ = 10 4 (right).The colour scale is the same in each simulation and ranges between red and blue, representing prograde and retrograde motion respectively

Figure 5 :
Figure 5: Kinetic energy spectra for the cases γ = 0 (left), γ = 10 3 (middle) and γ = 10 4 (right).The spectra are shown as a function of the spherical harmonics number l, split into the zonal modes (m = 0) in grey, and the non-zonal modes (m ̸ = 0) in black.For reference, the Kraichnan scalings −5/3 and −3 are shown.The dotted vertical line indicates the modes at which the forcing is centred.
s = N −1 2 and {:::} denotes the 6j Wigner symbol.For the latter, there is a known asymptotic formulaFlude[1998]: given a, b, c satisfying the triangle inequality and R ≫ a, b, c we have