2. Kinetic energy spectrum on a sphere

In turbulence theory, the kinetic energy wavenumber spectrum is often defined as a Fourier transform of the covariance function of velocity measured at two different points and an assumption of isotropy is often made to simplify the description. Boer (Reference Boer1983) has advanced such an approach to the analysis of two-dimensional incompressible turbulence on a sphere, with the Fourier transform replaced by the spherical harmonic transform. We will adopt a different approach and argue that the constraint of isotropy does not lead to any simplifications but only to restrictions. Fortunately, the spherical geometry offers another attractive mode of description where the spatial mean value over the entire sphere is the only statistical measure that is used and no assumption of isotropy is invoked.

We let $\psi$ be a real scalar function defined on a sphere with radius $a$, $(\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z)$ the unit vectors associated with a Cartesian coordinate system with origin at the centre of the sphere, and $(\theta, \phi )$ the corresponding polar coordinates, where the $z$-axis is the polar axis, $\theta$ is the colatitude and $\phi$ is the longitude. We expand $\psi$ in spherical harmonics as

(2.1)$$\begin{gather} \psi(\theta, \phi) = \sum_{l=0}^{\infty} \psi_{l}(\theta, \phi) , \end{gather}$$

(2.2)$$\begin{gather}\psi_{l}(\theta, \phi) = \sum_{m={-}l}^{l} a_{lm} Y_{l}^{m}(\theta, \phi) , \end{gather}$$

where we have separated the summation into two different equations to emphasise that the $l$- and $m$-numbers have different significance. The spherical harmonics are defined as

(2.3)\begin{equation} Y_l^m(\theta, \phi) = \sqrt{(2l+1) \frac{(l-m)!}{(l+m)!}} \, P_l^m(\cos \theta) {\rm e}^{{\rm i}\phi}, \end{equation}

where $P_l^m$ are the associated Legendre functions. They are eigenfunctions of the Laplace operator and orthogonal with respect to integration over the sphere,

(2.4)$$\begin{gather} \nabla^2 Y_l^m ={-}\frac{l(l+1)}{a^2} Y_l^m , \end{gather}$$

(2.5)$$\begin{gather}\frac{1}{4 {\rm \pi}} \int Y_l^m {Y_s^p}^* \,{\text{d}} \varOmega = \delta_{ls} \delta_{mp} , \end{gather}$$

and have the symmetry property

(2.6)\begin{equation} Y_l^{{-}m} = ({-}1)^m {Y_l^{m}}^{*} , \end{equation}

where $\varOmega$ is solid angle and superscript $*$ denotes the complex conjugate. The fact that $\psi$ is real imposes the constraint

(2.7)\begin{equation} a_{l({-}m)} = ({-}1)^m {a_{lm}}^{*} , \end{equation}

on the expansion coefficients for $\psi$, implying that $\psi _l (\theta, \phi )$ also is a real function.

As pointed out by Tang & Orszag (Reference Tang and Orszag1978), the $l$-number is more significant than the $m$-number since it is only $l$ that enters into the eigenvalue of the Laplace operator. The significance of the $l$-number is more clearly revealed by using a result from group representation theory to show that $\psi _l (\theta, \phi )$ is invariant under a coordinate transformation between systems with different polar axes. A new Cartesian coordinate system is defined as $R [(\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z)] = (\boldsymbol {e}_{x^\prime }, \boldsymbol {e}_{y^\prime }, \boldsymbol {e}_{z^\prime })$, where $R(\alpha, \beta, \gamma )$ is a rotation defined by the Euler angles $(\alpha, \beta, \gamma )$. We let $(\theta ^{\prime }, \phi ^{\prime })$ be the polar coordinates corresponding to the primed system, with the $z^\prime$-axis being the new polar axis. For a given $R(\alpha, \beta, \gamma )$ there is a coordinate transformation $(\theta, \phi ) \rightarrow (\theta ^\prime, \phi ^\prime )$, meaning that $(\theta ^\prime, \phi ^\prime )$ and $(\theta, \phi )$ are two pairs of coordinates referring to the same points on the sphere. In the primed system we can expand $\psi$ as

(2.8)$$\begin{gather} \psi(\theta^\prime, \phi^\prime) = \sum_{l=0}^{\infty} \psi^\prime_{l}(\theta^\prime, \phi^\prime), \end{gather}$$

(2.9)$$\begin{gather}\psi^\prime_{l}(\theta^\prime, \phi^\prime) = \sum_{p={-}l}^{l} b_{lp} Y_{l}^{p}(\theta^\prime, \phi^\prime) . \end{gather}$$

A powerful result of group representation theory (see chapter 15 in Wigner Reference Wigner1959) states that

(2.10)$$\begin{gather} Y_{l}^p (\theta^\prime, \phi^\prime) = \sum_{m={-}l}^{l} D_{pm}^{(l)} Y_{l}^{m} (\theta, \phi) , \end{gather}$$

(2.11)$$\begin{gather}Y_{l}^m (\theta, \phi ) = \sum_{p={-}l}^{l} {D_{pm}^{(l)}}^{*} Y_{l}^{p} (\theta^\prime, \phi^\prime) , \end{gather}$$

where $D_{pm}^{(l)}(\alpha, \beta, \gamma )$ is a unitary matrix of order $2l+1$, the so called ‘Wigner D-matrix’, which is an irreducible representation of the $SO(3)$ group. Using (2.10) and orthogonality we find

(2.12)\begin{align} b_{lp} &= \frac{1}{4{\rm \pi}} \int Y_{l}^{p}(\theta^\prime, \phi^\prime)^{*} \psi\, \text{d}\varOmega\nonumber\\ &= \frac{1}{4{\rm \pi}} \int \sum_{m={-}l}^{l} {D_{pm}^{(l)}}^* Y_{l}^{m} (\theta, \phi)^{*} \psi \,\text{d} \varOmega = \sum_{m={-}l}^{l} a_{lm}{D_{pm}^{(l)}}^{*} . \end{align}

Substituting (2.12) into (2.9) and using (2.11) we obtain

(2.13)\begin{equation} \psi_{l}^\prime(\theta^\prime, \phi^\prime) = \sum_{p={-}l}^{l} \sum_{m={-}l}^{l} a_{lm}{D_{pm}^{(l)}}^{*} Y_{l}^{p}(\theta^\prime, \phi^\prime) = \sum_{m={-}l}^{l} a_{lm} Y_{l}^{m} (\theta, \phi) = \psi_l(\theta, \phi) . \end{equation}

Since $(\theta, \phi )$ and $(\theta ^\prime, \phi ^\prime )$ are two pairs of coordinates referring to the same points on the sphere, (2.13) shows that $\psi _l^\prime$ is exactly the same function as $\psi _l$. In other words, for each $l$, $\psi _{l}$ is a coordinate invariant scalar field.

We will now let $\psi$ play the role of a streamfunction. A standard way of doing this would be to define the velocity field as $\boldsymbol {u} = - \boldsymbol {e}_r \times \boldsymbol {\nabla } \psi$, where ${\boldsymbol {e}}_r$ is the radial unit vector. However, since $\boldsymbol {e}_r$ is not an intrinsic vector of the geometry such an approach would lead us into complications when differentiation is carried out on ${\boldsymbol {u}}$. We will adopt another approach, also used by Eriksson & Nordmark (Reference Eriksson and Nordmark2020), which greatly simplifies differential operations, but comes with a price – the introduction of a notation which is rarely used in fluid mechanics. We let $\boldsymbol {\nabla }$ denote the surface gradient operator acting along the tangent of the sphere and $\star$ the so-called ‘Hodge star operator’ (see e.g. Bishop & Goldberg Reference Bishop and Goldberg1980; Schutz Reference Schutz1980) which produces the Hodge dual of the element that it is operating on. In our case the star is simply rotating a tangent vector on the sphere counter-clockwise a right angle about $\boldsymbol {e}_r$. Equipped with these tools we define the velocity field as

(2.14)\begin{equation} \boldsymbol{u} \equiv{\star} \boldsymbol{\nabla} \psi , \end{equation}

ensuring that it is divergence free, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = \boldsymbol {\nabla } \boldsymbol {\cdot } \star \boldsymbol {\nabla } \psi = 0$. With (2.14) we have the opposite sign convention for $\psi$ as compared with the standard given by $\boldsymbol {u} = - \boldsymbol {e}_r \times \boldsymbol {\nabla } \psi$. The operational rules including $\{ \boldsymbol {\nabla }, \star, \boldsymbol {\cdot } \}$ are given in Appendix A. The advantage of the notation becomes clear when the vorticity field is to be defined. Instead of acting on $\boldsymbol {u}$ with the curl operator and then project on the radial unit vector – a procedure that would bring us into trouble (see Krishnamurthy Reference Krishnamurthy2019) – we define the vorticity as the divergence of the rotated velocity, that is

(2.15)\begin{equation} \omega \equiv \boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u} = \boldsymbol{\nabla} \boldsymbol{\cdot} \star{\star} \boldsymbol{\nabla} \psi ={-} \nabla^2 \psi , \end{equation}

where we have used $\star \star = -1$. Just as in the case of $\psi$, the sign convention of $\omega$, as defined in (2.15), will be the opposite compared with the standard. Readers who feel uncomfortable with this can restore the standard sign conventions of $\psi$ and $\omega$ by introducing a non-standard sign convention of the Hodge star by letting it perform a clockwise instead of a counter-clockwise rotation. That (2.15) is an appropriate definition of vorticity on the sphere is clear from the fact that it permits us to calculate the circulation along any closed curve on the sphere as the corresponding surface integral of $\omega$, in accordance with Stokes’ theorem.

The velocity and vorticity fields can be expanded as

(2.16a,b)\begin{equation} \boldsymbol{u} = \sum_{l=1}^{\infty}\boldsymbol{u}_l , \quad \omega = \sum_{l=1}^{\infty} \omega_l , \end{equation}

where

(2.17a,b)\begin{equation} \boldsymbol{u}_l ={\star} \boldsymbol{\nabla} \psi_l , \quad \omega _l = \frac{l(l+1)}{a^2} \psi_l . \end{equation}

For each $l$, $\boldsymbol {u}_l$ is a coordinate invariant vector field and $\omega _l$ is a coordinate invariant scalar field. That a non-divergent velocity field on a sphere can be completely decomposed as in (2.16a,b) follows from the more general result that vector spherical harmonics form a complete set (Freeden & Gutting Reference Freeden and Gutting2008).

The mean kinetic energy and enstrophy over the sphere can now be expressed as

(2.18)$$\begin{gather} \bar{E} = \frac{1}{4 {\rm \pi}} \int \frac{1}{2} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \,{\text{d}} \varOmega = \sum_{l=1}^{\infty} E(l) , \end{gather}$$

(2.19)$$\begin{gather}\bar{{\mathcal{E}}} = \frac{1}{4 {\rm \pi}} \int \frac{1}{2} \omega^2 \,{\text{d}} \varOmega = \sum_{l=1}^{\infty} \frac{l(l+1)}{a^2} E(l) , \end{gather}$$

where

(2.20)\begin{equation} E(l) = \frac{l(l+1)}{2a^2} \sum_{m={-}l}^{l} a_{lm} a_{lm}^* = \frac{l(l+1)}{2a^2} \sum_{m={-}l}^{l} b_{lm} b_{lm}^* , \end{equation}

is the kinetic energy spectrum. From (2.20) it is clear that $E(l)$ is invariant under coordinate transformations between systems with different polar axes. It should be emphasised that the result is derived without making any assumption regarding the structure of the flow field. For example, the flow may consist of a single localised structure on the sphere, such as a vortex or an equatorial jet. Energy spectra calculated using spherical harmonic systems whose polar axes are oriented differently relative to the structure will all be equal.

2.1. Isotropy

Since the invariance of $E(l)$ with respect to rotations of the polar axis is a property that is satisfied for any flow field, isotropy must constrain the energy distribution even further. The statistics of an isotropic field should contain no information on the direction of the polar axis. Strictly speaking, isotropy therefore implies that the energy distribution among modes with different $m$ and the same $l$ also must be invariant under rotations. It is easy to see that there is no such distribution for a finite $l$ and that isotropy, in a strict sense, therefore can be fulfilled only in the limit $l \rightarrow \infty$. This is analogous to a flow field on a plane with periodic boundary conditions, for which isotropy can be only approximately fulfilled for finite wavenumbers. Likewise, for large but finite $l$, say $l > 100$, the distribution of energy among modes with the same $l$ and different $m$ may become approximately invariant under rotations, with increasing accuracy for larger and larger $l$. According to the analysis of the two-point covariance function by Boer (Reference Boer1983) isotropy implies that the variance of the expansion coefficient should be independent of $m$. Expressed in a different language kinetic energy should be equipartitioned among $m$-modes with the same $l$. Intuitively, it is quite clear that this must be the case. In a neighbourhood of a pole, we can make a flat Earth approximation and pass over to plane polar coordinates, as long as curvature effects are small. The $\theta$-coordinate of the spherical description will transform to a radial coordinate $r = |\theta |a$, while $\phi$ retains the role of an azimuth. Looking at the flow field from the pole, isotropy means statistical independence of the azimuthal coordinate for a fixed radial coordinate, which in spectral space translates to independence of $m$ for a fixed $l$. Since the polar axis can be chosen freely, this reasoning can be applied generally.

It is evident that isotropy cannot be even approximately fulfilled for low-order modes. Consider the $l = 1$ mode. With reference to a specific coordinate system, $(\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z)$, we have

(2.21)\begin{equation} \psi_1(\theta, \phi) = a_{1({-}1)} Y_1^{{-}1}(\theta, \phi) + a_{10} Y_1^0 (\theta, \phi) + a_{11} Y_1^1 (\theta, \phi) . \end{equation}

By (2.10) there exists another coordinate system, $(\boldsymbol {e}_{x^\prime }, \boldsymbol {e}_{y^\prime }, \boldsymbol {e}_{z^\prime })$, where

(2.22)\begin{equation} \psi_1^\prime (\theta^\prime, \phi^\prime) = b_{10} Y_1^0 (\theta^\prime, \phi^\prime) = b_{10} \sqrt{3} \cos \theta^\prime , \end{equation}

with $b_{1(-1)} = b_{11} = 0$. The corresponding velocity field can be written as

(2.23)\begin{equation} \boldsymbol{u}_1 = b_{10} \frac{\sqrt{3}}{a} \sin \theta^\prime \boldsymbol{e}_{\phi^\prime} , \end{equation}

which is a solid body rotation around the $z^\prime$-axis. It is quite clear that $\psi _1$ must be put to zero under the constraint of isotropy, and so must $\psi _2$ and $\psi _3$ and so forth up to a certain $l$ for which it may be assumed that isotropy is approximately fulfilled. The constraint of isotropy would thus force us to refrain from studying large scale dynamics and for this reason we deem it as too restrictive.

3. The spectral energy equation

The incompressible Navier–Stokes equations in a rotating frame of reference on a sphere with radius $a$ can be written as

(3.1)\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} ={-} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} - \boldsymbol{\nabla} (p/\rho) - f \star \boldsymbol{u} + \nu \left (\Delta \boldsymbol{u} + \frac{\boldsymbol{u}}{a^2} \right) , \end{equation}

where $\boldsymbol {\nabla }$ is the gradient operator acting along the tangent of the sphere, $f = 2 {\boldsymbol \varOmega } \boldsymbol {\cdot } \boldsymbol {e}_r$ is the Coriolis parameter, $\boldsymbol{\varOmega}$ is the rotation vector and

(3.2)\begin{equation} \Delta {\boldsymbol{u}} ={-} \star \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) + \frac{\boldsymbol{u}}{a^2} , \end{equation}

is the vector Laplace of $\boldsymbol {u}$ (see Appendix B). The centrifugal force is included in the pressure term. Compared with the standard form of the Navier–Stokes, there is an extra metric term, $\nu {\boldsymbol {u}}/a^2$, on the right-hand side of (3.1). In Appendix B we show how this term arises from the spherical geometry by using tensor analysis on a two-dimensional manifold. At the request of one of the reviewers we also show in Appendix C how the viscous term in (3.1) can be derived from the expression of the viscous term of the three-dimensional equation by considering a fluid layer within a thin spherical shell with stress free boundaries.

The equation for $E(l)$ can be derived by taking the scalar product between each term in (3.1) and $\boldsymbol {u}_l$ and averaging over the sphere. In this way, we obtain the time derivative of $E(l)$ from the left-hand side of (3.1),

(3.3)\begin{equation} \frac{1}{4{\rm \pi}} \int \boldsymbol{u}_l \boldsymbol{\cdot} \sum_{n= 1}^{\infty} \frac{\partial \boldsymbol{u}_n}{\partial t} \,{\text{d}} \varOmega = \frac{1}{4{\rm \pi}} \int \boldsymbol{u}_l \boldsymbol{\cdot} \frac{\partial \boldsymbol{u}_l}{\partial t} \,{\text{d}} \varOmega = \frac{\partial E(l)}{\partial t} . \end{equation}

The spectral term emanating from the nonlinear term in (3.1) we name the transfer term

(3.4)\begin{equation} T(l) ={-}\frac{1}{4{\rm \pi}} \int (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}) \boldsymbol{\cdot} \boldsymbol{u}_l \,{\text{d}} \varOmega , \end{equation}

in accordance with standard vocabulary in isotropic turbulence theory. Energy conservation is expressed as

(3.5)\begin{equation} \sum_{l=1}^{\infty} T(l) ={-}\frac{1}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{u} \frac{1}{2} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \right ) \,{\text{d}} \varOmega = 0 . \end{equation}

Enstrophy is conserved, because vorticity is conserved at a fluid particle. Using the rules in Appendix A we have

(3.6)\begin{align} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} &={-} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) \star \boldsymbol{u} + \tfrac{1}{2} \boldsymbol{\nabla} ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u}) ,\end{align}

(3.7)\begin{align} \boldsymbol{\nabla} \boldsymbol{\cdot} \star ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ) &={-} \boldsymbol{\nabla} \boldsymbol{\cdot} \star [ (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) \star \boldsymbol{u} ] \nonumber\\ &= \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) + (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega . \end{align}

Enstrophy conservation imposes a constraint on $T(l)$ which can be derived in the following way:

(3.8)$$\begin{gather} T(l) ={-}\frac{1}{4{\rm \pi}} \int (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}) \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi_l \,{\text{d}} \varOmega = \frac{1}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} [{\star} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}) \psi_l ] \,{\text{d}} \varOmega \nonumber\\ - \frac{1}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} \star ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ) \psi_l \,{\text{d}} \varOmega ={-} \frac{a^2}{l(l+1)} \frac{1}{4{\rm \pi}} \int (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega ) \omega_l \,{\text{d}} \varOmega , \end{gather}$$

where we have used (3.7), (2.17a,b) and the rules in Appendix A. Multiplying both sides by $l(l+1)$ and summing over $l$ we find

(3.9)\begin{equation} \sum_{l=1}^{\infty} l(l+1) T(l) ={-} \frac{a^2}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{u} \frac{1}{2} \omega^2 \right ) \,{\text{d}} \varOmega = 0 , \end{equation}

which is the expression for enstrophy conservation.

The contribution from the pressure term to the spectral energy equation is trivially zero,

(3.10)\begin{equation} -\frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} (p/\rho) \boldsymbol{\cdot} \boldsymbol{u}_l \,{\text{d}} \varOmega ={-}\frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} (p \boldsymbol{u}_l/\rho ) \,{\text{d}} \varOmega = 0 . \end{equation}

As found in previous investigations based on the vorticity equation (e.g. Tang & Orszag Reference Tang and Orszag1978) the Coriolis term does not make any contribution to the spectral energy equation. We confirm this by direct integration in the system whose polar axis coincides with the rotation axis,

(3.11)\begin{align} &-\frac{1}{4{\rm \pi}} \int f \star \boldsymbol{u} \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi_l \,{\text{d}} \varOmega ={-}\frac{1}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} f \psi_l) \,{\text{d}} \varOmega \nonumber\\ &\qquad+\frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} f \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi \psi_l \,{\text{d}} \varOmega ={-}\frac{1}{4 {\rm \pi}} \int 2 | {\boldsymbol \varOmega} | \sin \theta \boldsymbol{e}_{\theta} \boldsymbol{\cdot} \left( -\frac{\partial_{\phi} \psi \boldsymbol{e}_\theta}{a \sin\theta} +\frac{\partial_{\theta} \psi \boldsymbol{e}_\phi }{a} \right ) \psi_l \,{\text{d}} \varOmega \nonumber\\ &\quad= \frac{1}{4{\rm \pi} a} \int 2 | {\boldsymbol \varOmega} | \sum_{n= 1}^\infty \frac{\partial \psi_n}{\partial \phi} \psi_l \,{\text{d}} \varOmega = 0 . \end{align}

That (3.11) is zero can be seen as follows. For $n = l$, use $\psi _l \partial _\phi \psi _l = 0.5 \partial _\phi (\psi _l^2)$ and integrate in $\phi$. For $n \ne l$, use $\partial _{\phi } Y_{n}^{m} = {\rm i}m Y_{n}^{m}$ and orthogonality of the spherical harmonics. That the Coriolis term does not contribute to the spectral energy equation can also be understood from the fact that it may be incorporated into the $l= 1$ mode, by transforming the equation back to a non-rotating frame of reference. The mean angular momentum with respect to the polar axis of a chosen spherical harmonic system can be calculated as

(3.12)\begin{equation} \overline{L_z} = \frac{1}{4{\rm \pi}} \int_{0}^{2{\rm \pi}} \int_{0}^{\rm \pi} \sqrt{3} \sin^3 \theta a_{10} \,\mathrm{d} \theta \,\mathrm{d} \phi = \frac{2}{\sqrt{3}} a_{10} , \end{equation}

where the contributions from all other modes vanish by Rodrigues’ formula. Clearly, in the absence of drag, $\overline {L_z}$ is a conserved quantity and $a_{10}$ is therefore conserved and so is $\boldsymbol {u}_1$, since the mean angular momentum with respect to any axis is conserved.

The calculation of the contribution from the viscous term is straightforward. First,

(3.13)\begin{align} \Delta \boldsymbol{u}_l &={-}\star \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star{\star} \boldsymbol{\nabla} \psi_l ) + \frac{\boldsymbol{u}_l}{a^2} ={\star} \boldsymbol{\nabla} (\nabla^2 \psi_l ) + \frac{\boldsymbol{u}_l}{a^2} \nonumber\\ &={-}\frac{l(l+1)}{a^2} \star \boldsymbol{\nabla} \psi_l + \frac{\boldsymbol{u}_l}{a^2} ={-} \frac{[l(l+1)-1]}{a^2} \boldsymbol{u}_l . \end{align}

Then, using this expression we find that the contribution from the viscous term in (3.1) to the spectral energy equation is

(3.14)\begin{equation} -\frac{2\nu [l(l+1)-2]}{a^2} E(l) . \end{equation}

Previous studies (e.g. Fjørtoft Reference Fjørtoft1953; Tang & Orszag Reference Tang and Orszag1978; Boer Reference Boer1983) gave an expression of the form $-2\nu l(l+1) E(l)/a^2$. That this cannot be completely accurate can be understood from the fact that there can be no viscous dissipation in the $l=1$ mode.

To summarise, we have derived the spectral energy equation on a sphere

(3.15)\begin{equation} \frac{\partial E(l)}{\partial t} = T(l) -\frac{2\nu [l(l+1)-2]}{a^2}E(l) , \end{equation}

where conservation of angular momentum, energy and enstrophy can be expressed as

(3.16)$$\begin{gather} T(1) = 0 , \end{gather}$$

(3.17)$$\begin{gather}\sum_{l=2}^{\infty} T(l) = 0 , \end{gather}$$

(3.18)$$\begin{gather}\sum_{l= 2}^{\infty} l(l+1) T(l) = 0 . \end{gather}$$

The $l = 1$ mode represents a stationary solid body rotation. A linear Ekman drag force, $-\alpha \boldsymbol {u}$, added to the right-hand side of (3.1), modelling the interaction between the fluid and a solid ground, would give rise to a term $-2\alpha E(l)$ on the right-hand side of (3.15). Such a term would act on all modes, including the $l=1$ mode, attenuating the mean angular momentum of the fluid. The (3.15) is similar to the corresponding equation for the energy spectrum of two-dimensional homogeneous isotropic turbulence on a plane. An important difference is that (3.15) has been derived without invoking assumptions of homogeneity and isotropy.

4. Triad interactions

The transfer term can be written as a sum of interactions involving triad of modes

(4.1)\begin{equation} T(l) ={-}\sum_{s =2}^{\infty} \sum_{n=2}^{\infty} \frac{1}{4{\rm \pi}} \int (\boldsymbol{u}_s \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_n ) \boldsymbol{\cdot} \boldsymbol{u}_l \,{\text{d}} \varOmega . \end{equation}

Using a formula derived by Gaunt (Reference Gaunt1929), involving the integral of the product between three associated Legendre polynomials, Silberman (Reference Silberman1954) brilliantly proved that the nonlinear interaction within a triad of modes, $\{l,n,s\}$, is zero, unless

(4.2)$$\begin{gather} l + n + s = {\text{odd integer}} , \end{gather}$$

(4.3)$$\begin{gather}| l -n | < s < l+n , \end{gather}$$

where it should be noted that the triangle inequality (4.3) is strict. Below, it will be shown that the interaction is also zero if any two of $\{ l, n, s \}$ are identical. Fjørtoft (Reference Fjørtoft1953) derived relations involving energy exchanges within a triad under the assumption that energy exchanges take place only in a single triad, that is at initial time for initial conditions consisting of only three modes. It will now be shown that his relations hold generally, for energy exchanges within any triad of modes independently of all other exchanges.

We define the energy transfer from $n$ to $l$ in the interaction involving $s$ as the third mode, as

(4.4)\begin{equation} T_{E}(s, n, l) ={-} \frac{1}{4{\rm \pi}} \int ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{n} ) \boldsymbol{\cdot} {\boldsymbol{u}}_{l } \, \mathrm{d} \varOmega . \end{equation}

The transfer is antisymmetric with respect to interchange of $n$ and $l$, that is

(4.5)\begin{equation} T_{E}(s, n, l) ={-} T_{E}(s, l, n) , \end{equation}

meaning that the gain of $l$ is the loss of $n$. This follows from

(4.6)\begin{align} 0&= \int \boldsymbol{\nabla} \boldsymbol{\cdot} ({\boldsymbol{u}}_{s} \boldsymbol{u}_{n} \boldsymbol{\cdot} {\boldsymbol{u}}_{l} ) \, \mathrm{d} \varOmega \nonumber\\ &= \int ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{n} ) \boldsymbol{\cdot} {\boldsymbol{u}}_{l} \, \mathrm{d} \varOmega + \int ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{l} ) \boldsymbol{\cdot} {\boldsymbol{u}}_{n} \, \mathrm{d} \varOmega , \end{align}

where we have used $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}_{s} = 0$. In a similar way we define the transfer of enstrophy from $n$ to $l$, in the triad interaction involving $s$ as the third mode, as

(4.7)\begin{align} &T_{{\mathcal{E}}}(s,n,l) ={-}\frac{1}{4{\rm \pi}} \int ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega} _{n} ) {\omega}_{l} \, \mathrm{d} \varOmega \nonumber\\ &\quad ={-}\frac{l(l+1)n(n+1)}{a^4} \frac{1}{4{\rm \pi}} \int \star \boldsymbol{\nabla} \psi_s \boldsymbol{\cdot} \boldsymbol{\nabla} \psi _{n} \psi_{l} \, \mathrm{d} \varOmega . \end{align}

It is clear that

(4.8)$$\begin{gather} T_{\mathcal{E}} (s,n,l) ={-} T_{\mathcal{E}}(s,l,n) , \end{gather}$$

(4.9)$$\begin{gather}T_{\mathcal{E}} (s,s,l) = 0 , \end{gather}$$

implying that the transfer is zero if two modes in a triad are identical.

We define the total gain of energy and enstrophy in unit time of mode $l$ at the expense of modes $n$ and $s$ as

(4.10)$$\begin{gather} G(l, n, s) = T_E (s,n,l) + T_E (n,s,l) , \end{gather}$$

(4.11)$$\begin{gather}H(l, n, s) = T_{\mathcal{E}} (s,n,l) + T_{\mathcal{E}} (n,s,l) . \end{gather}$$

In the next step we show that

(4.12)\begin{equation} H(l, n, s) = \frac{l(l+1)}{a^2} G(l,n,s) . \end{equation}

First, a straightforward calculation gives

(4.13)\begin{equation} H(l,n,s) ={-} \frac{1}{4{\rm \pi}} \frac{ l(l+1)}{a^2} \int ( \boldsymbol{u}_s \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega}_n + \boldsymbol{u}_n \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega}_s ) \psi_l\, \mathrm{d} \varOmega , \end{equation}

where we have used (2.17b). Then, a corresponding calculation gives

(4.14)\begin{align} G(l,n,s) &={-}\frac{1}{4 {\rm \pi}} \int ( {\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} _{n} + \boldsymbol{u}_{n} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{s} ) \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi_l \,\mathrm{d} \varOmega \nonumber\\ & = \frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} [{\star} ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} _{n} + \boldsymbol{u}_{n} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{s} ) \psi_l ] -\boldsymbol{\nabla} \boldsymbol{\cdot} \star ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} _{n} + \boldsymbol{u}_{n} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{s} ) \psi_l \,\mathrm{d} \varOmega \nonumber\\ &= \frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} \star [ (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}_s) \star \boldsymbol{u}_n +(\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}_n) \star \boldsymbol{u}_s ] \psi_l \,\mathrm{d} \varOmega \nonumber\\ &={-} \frac{1}{4{\rm \pi}} \int ( \boldsymbol{u}_n \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega}_s + \boldsymbol{u}_s \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega}_n ) \psi_l \,\mathrm{d} \varOmega , \end{align}

where we have used the rules in Appendix A and $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}_n = \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}_s = 0$. Now, (4.12) follows from (4.13) and (4.14). Energy and enstrophy conservation follow from (4.5) and (4.8),

(4.15)$$\begin{gather} G(l,n,s) + G(n,s,l) + G(s,l,n) = 0 , \end{gather}$$

(4.16)$$\begin{gather}H(l, n,s) + H(n,s,l) + H(s,l,n) = 0 . \end{gather}$$

From (4.12) and (4.16) we also have

(4.17)\begin{equation} l(l+1) G(l,n,s) + n(n+1) G(n,s,l) + s(s+1) G(s,l,n) = 0 . \end{equation}

For a triad of modes, $\{l, n, s \}$, such that $l < n < s$, the Fjørtoft (Reference Fjørtoft1953) relations now follow from (4.15) and (4.17),

(4.18)\begin{equation} \left. \begin{aligned} \frac{G(l,n,s)}{G(s,l,n)} & = \frac{s(s+1) - n(n+1)}{n(n+1)-l(l+1)} > 0 , \\ \frac{G(n,s,l)}{G(l,n,s)} & = \frac{l(l+l) - s(s+1)}{s(s+1)-n(n+1)} < 0 , \\ \frac{G(s,l,n)}{G(n,s,l)} & = \frac{n(n-1) - l(l+1)}{l(l+1)-s(s+1)} < 0 . \end{aligned} \right\} \end{equation}

The relations (4.18) are analogous to the relations for energy exchanges within a triad, $\{k,p,q\}$, of Fourier modes, derived by Kraichnan (Reference Kraichnan1967, equation (2.7)). Apart from the different in notation – our $G(l,n,s)$ corresponds to $T(k,p,q)$ of Kraichnan – the difference is in the eigenvalue of the Laplace operator, which is $k^2$ in Fourier space and $l(l+1)$ in spherical harmonics space.

Nowhere in the derivation of (4.18) have we used that the elementary mode for each $l$ has been defined as a sum over the contributions from all $m$, as in (2.16a,b). The relations (4.18) would therefore hold also for a triad of single spherical harmonic modes. Such a triad can be specified as $\{ (l,m), (n,p), (s,q) \}$, corresponding to spherical harmonics $Y_l^m$, $Y_n^p$ and $Y_s^q$. However, the distribution of energy between modes with the same $l$ and different $m$ depends on the chosen coordinate system and for this reason it appears to be of limited interest that (4.18) actually hold even though no summation over $m$ is carried out. A related but more interesting result, is that for a chosen coordinate system there is no net enstrophy or energy exchange within a triad where all three modes have the same $l$. With a slight generalisation of the notation in (4.7) we can define the enstrophy gain of mode $(l,m)$ at the expense of mode $(n,p)$ in the triad interaction involving $(s, q)$ as the third mode as

(4.19)$$\begin{gather} T_{{\mathcal{E}}}[(s,q), (n,p), (l,m) ]={-}\frac{1}{4{\rm \pi}} \int ({\boldsymbol{u}}_{sq} \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega} _{np} ) {\omega}_{lm} \, \mathrm{d} \varOmega \nonumber\\ ={-}\frac{l(l+1)n(n+1)}{a^4} \frac{1}{4{\rm \pi}} \int \star \boldsymbol{\nabla} \psi_{sq} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi _{np} \psi_{lm} \, \mathrm{d} \varOmega . \end{gather}$$

Clearly, we have

(4.20)$$\begin{gather} T_{{\mathcal{E}}}[(s,q), (l,m), (n,p) ] ={-} T_{{\mathcal{E}}}[(s,q), (n,p), (l,m) ] , \end{gather}$$

(4.21)$$\begin{gather}T_{{\mathcal{E}}}[(s,q), (s,p), (l,m) ] ={-} T_{{\mathcal{E}}}[(s,p), (s,q), (l,m) ] . \end{gather}$$

The enstrophy gain of mode $(l,m)$ at the expense of the two other modes is

(4.22)\begin{equation} H[ (l,m), (s,q), (n,p) ] = T_{{\mathcal{E}}}[(s,q), (n,p), (l,m) ]+ T_{{\mathcal{E}}}[(n,p), (s,q), (l,m) ] . \end{equation}

It now follows that

(4.23)\begin{equation} H[(l,m),(l,q), (l,p)] = 0 , \end{equation}

which shows that there is no net enstrophy or energy exchange within a triad where all three modes have the same $l$.

5. Vanishing kinetic energy dissipation and cascade directions

Apart from the simultaneous conservation of energy and enstrophy by the nonlinear term there is a rigorous result which is encapsulating the difference between two- and three-dimensional turbulence, and that is the vanishing of kinetic energy dissipation in the limit $\nu \rightarrow 0$. The mean kinetic energy and enstrophy equations can be written as

(5.1)$$\begin{gather} \frac{\mathrm{d} \bar{E}}{\mathrm{d} t } ={-}\epsilon , \end{gather}$$

(5.2)$$\begin{gather}\frac{\mathrm{d} \bar{\mathcal{E}}}{\mathrm{d} t} ={-}\eta , \end{gather}$$

where

(5.3)$$\begin{gather} \epsilon \equiv \sum_{l = 2}^{\infty} \frac{2\nu ( l(l+1)-2)}{a^2} E(l) , \end{gather}$$

(5.4)$$\begin{gather}\eta \equiv \sum_{l=2}\frac{2\nu ( l(l+1)-2)l(l+1)}{a^4} E(l) , \end{gather}$$

are the mean kinetic dissipation and mean enstrophy dissipation, respectively. We now have

(5.5)\begin{equation} \epsilon < 2 \nu \bar{\mathcal{E}} , \end{equation}

implying that there can be no Richardson energy cascade in two-dimensional turbulence.

In order to say something more on energy and enstrophy transfer in wavenumber space some additional assumption is needed. A deep insight of Fjørtoft (Reference Fjørtoft1953) was that the inequalities (4.18) imply that there are two types of triad interactions. The first type is where the middle wavenumber $n$ acts as a source for each of the other two modes $l$ and $s$, and the second type is where it acts as a sink. He assumed that first type interactions are dominant and argued that as a consequence more energy will flow upscale than downscale. As pointed out by Merilees & Warn (Reference Merilees and Warn1975) his arguments were partly flawed, since they are based on the assumption that more energy is transferred to the smallest wavenumber mode $l$ than to the largest wavenumber mode $s$ in a triad interaction, which is not generally true. Kraichnan (Reference Kraichnan1967) developed the notion of two separated infinitely broad constant flux ranges – a downscale enstrophy flux range at large wavenumbers where the energy flux is zero and an upscale energy flux range at small wavenumbers where the enstrophy flux is zero. In both ranges, first type interactions are dominant, in accordance with the assumption of Fjørtoft (Reference Fjørtoft1953). Most subsequent studies (e.g. Eyink Reference Eyink1996; Tran & Shepherd Reference Tran and Shepherd2002; Gkioulekas & Tung Reference Gkioulekas and Tung2007) have used an assumption of stationarity in the analysis of cascade directions, which can be invoked only if some extra terms are added to (3.1). Usually, the added terms are forcing acting in a narrow wave number band and dissipation acting at small and large wave numbers. Of particular interest is an inequality derived by Gkioulekas & Tung (Reference Gkioulekas and Tung2007) under the assumption of stationarity and forcing in a wavenumber band $[k_1, \, k_2]$. In Fourier space, the ‘Danilov inequality’, named by Gkioulekas & Tung (Reference Gkioulekas and Tung2007) after the researcher who communicated it to them, reads

(5.6)\begin{equation} \varPi_{{\mathcal{E}}}(k) - k^2 \varPi_{E}(k) > 0 , \end{equation}

where $\varPi _{{\mathcal {E}}}(k)$ and $\varPi _{E}(k)$ are the spectral enstrophy and energy fluxes respectively, and $k$ is the magnitude of the wavenumber vector. As shown by Gkioulekas & Tung (Reference Gkioulekas and Tung2007) (5.6) holds for a stationary flow at all wavenumbers outside $[k_1, \, k_2]$. It will now be shown that a corresponding inequality can be derived in spherical harmonic space under the assumption of dominance of first type of interactions for which it holds at all wavenumbers.

In spherical harmonic space, we define spectral energy and enstrophy fluxes as

(5.7)$$\begin{gather} \varPi_{E}(k) \equiv \sum_{l = k}^{\infty} T(l) ={-} \sum_{l=2}^{k-1} T(l) , \end{gather}$$

(5.8)$$\begin{gather}a^2\varPi_{\mathcal{E}}(k) \equiv \sum_{l = k}^{\infty} l(l+1)T(l) ={-} \sum_{l=2}^{k-1} l(l+1)T(l) . \end{gather}$$

There are some basic mathematical facts that deserve to be pointed out. First, it is evident that the fluxes cannot be completely constant in a range unless $T(l)$ is identically zero in the same range. The notion of a constant flux range is therefore an idealisation which should be adopted with caution. Second, the maxima and minima of $\varPi _{E}(k)$ coincide with the maxima and minima of $\varPi _{\mathcal {E}}(k)$. If we permit ourselves to think about $T(l)$ as a continuous function, the maxima and minima of the fluxes are the zeros of $T(l)$.

Following a similar derivation as carried out by Kraichnan (Reference Kraichnan1967) we can use the conservation and symmetry properties of $G(l,n,s)$ to rewrite the fluxes as

(5.9)$$\begin{gather} \varPi_{E} (k) = \sum_{A_1(k,l,n,s)} G(s, n, l)- \sum_{A_2(k,l,n,s)} G(l,n,s) , \end{gather}$$

(5.10)$$\begin{gather}a^2 \varPi_{\mathcal{E}}(k) = \sum_{A_1(k,l,n,s)} s(s+1)G(s, n, l) - \sum_{A_2(k,l,n,s)} l(l+1)G(l,n,s) , \end{gather}$$

where

(5.11)$$\begin{gather} A_1(k,l,n,s) = \{ \{l,n,s\} ; 1 < l < n < k \le s < l+n ; \; l+n+s = {\text{odd}} \} , \end{gather}$$

(5.12)$$\begin{gather}A_2(k,l,n,s) = \{ \{l,n,s\} ; 1 < l < k \le n < s < l+n ; \; l+n+s = {\text{odd}} \} . \end{gather}$$

A detailed derivation of (5.9) is carried out in Appendix D and a derivation of (5.10) can be carried out in the same way. We note that in the summation over $A_1$, $G(s,n,l)$ is the gain of the largest mode $s$ and in the summation over $A_2$, $G(l, n, s)$ is the gain of the smallest wavenumber mode $l$. If the first type of interaction is dominant for all wavenumbers, the sums over $A_1$ in (5.9) and (5.10) give positive contributions to the fluxes while the sums over $A_2$ give negative contributions. For such a flow we will have

(5.13)\begin{equation} a^2 \varPi_{\mathcal{E}}(k) - k(k+1)\varPi_E(k) \ge 0 , \end{equation}

for all $k$. The inequality (5.13) implies that $T(l)$ must be positive both in the limit of small and large $l$. Therefore, $T(l)$ must have at least two maxima, one at small and one at large $l$. Assume that there are only two maxima. Then, there is a middle range where $T(l)$ is negative and has a minimum. If we permit ourselves to think about the fluxes as continuous functions, they will both have a single zero. The enstrophy flux will be positive at the zero of the energy flux and the energy flux will be negative at the zero of the enstrophy flux. With narrow band forcing and a sufficiently large separation of length scales it can also be shown that the amount of energy which is flowing to the small wavenumber region is larger than the amount of energy flowing to the large wavenumber region and that the opposite is true for enstrophy (Eyink Reference Eyink1996; Gkioulekas & Tung Reference Gkioulekas and Tung2007).

If the second type of interaction were dominant the opposite inequality

(5.14)\begin{equation} a^2 \varPi_{\mathcal{E}}(k) - k(k+1)\varPi_E(k) \le 0 , \end{equation}

would hold for all $k$. In this case, $T(l)$ must be negative in the limit of small as well as large wavenumbers and have at least two minima. Assume that there are just two minima. Then there is a middle wavenumber region where $T(l)$ is positive and has a maximum. Energy and enstrophy will flow from large and small wavenumbers into the middle region. Even though it is possible that such a transfer function can be generated for a short period of time from some initial conditions – for example an initial $E(l)$ having a broad minimum in a middle wavenumber range – it must be highly improbable. In particular, it seems extremely unlikely that enstrophy would flow from the region of very large wavenumbers where viscosity is acting into the middle wavenumber region. Two-dimensional turbulence cannot relax towards a stationary or quasi-stationary state as long as the transfer is dominated by the second type of interaction.

Thus, it seems that the Fjørtoft (Reference Fjørtoft1953) assumption of dominance of the first type of interaction may be the key to a deeper understanding of the problem of cascade directions. Whether the assumption can be rooted in some deeper principle is a very interesting question. Kraichnan (Reference Kraichnan1967) suggested that first type interactions represent ‘a statistical plausible spreading of the excitation in wave number: out of the middle wave number into the extremes’. Waleffe (Reference Waleffe1992) developed a stability argument suggesting that energy is generally flowing from the mode in which the gain/loss is largest. In two dimensions this is always the middle mode. In three dimensions, on the other hand, the dominant triad interaction is the one where the smallest wavenumber mode is loosing energy to the two larger ones. The argument is appealing, because the upscale energy transfer in two dimensions and the downscale transfer in three dimensions are derived from the same principle.

6. The double cascade theory expressed in spherical harmonic space

Kraichnan (Reference Kraichnan1967) formulated the equations of motions in an infinitely large quadratic domain with period boundary conditions in which case the discrete Fourier description passes over to an integral description. Formally, the description is valid only in the limit $Lk \rightarrow \infty$, where $L$ is the side length of the domain and $k$ is the magnitude of the wave vector. Under the assumption of isotropy he carried out a similarity analysis of the expressions for the energy fluxes, in which he made ingenious use of the conservation laws for triadic interactions. Based on this analysis he predicted that there are two possible similarity ranges – an upscale energy flux range in which the enstrophy flux is zero and a downscale enstrophy flux range where the energy flux is zero. It should be realised that such ranges are mathematical idealisations which only can be approached in the limit of large separation of length scales. In the constant energy flux range, the Fourier energy spectrum scales as (1.1) and in the constant enstrophy flux range it scales as (1.2), with a possible logarithmic correction.

From the similarity between the spectral energy equation (3.15), the Fjørtoft (Reference Fjørtoft1953) relations (4.18) and the flux expressions (5.9) and (5.10), on the one hand, and the corresponding expression in Fourier space, on the other hand, it is clear that there is a very close correspondence between the two descriptions in the limit of large wavenumbers. A heuristic demonstration of how they are related in this limit can be developed in the following way. Assume that the spherical harmonic spectrum peaks at $l \gg 1$ and the Fourier spectrum peaks at $k \gg 1/L$, and both spectra fall off exponentially at large wavenumbers. Then, for $n \ge 0$ we have

(6.1)\begin{equation} \sum_{l=1}^{\infty}a^{{-}2n} [l(l+1)]^n E(l) \approx \sum_{l=1}^{\infty} (l/a)^{2n} E(l) \approx \int_{0}^{\infty} k^{2n} E_F(k) \,{\text{d}} k , \end{equation}

with increasing accuracy for increasing $n$. The reason why (6.1) must hold is that the two sides are expressions of the same physical quantity. For $n = 0$, (6.1) is the mean energy, for $n =1$ the mean enstrophy and for $n= 2$ the mean palenstrophy. With a broad band spectrum peaking at $l \gg 1$, we can approximate the sum in the middle term of (6.1) with an integral. There is only one straightforward way to do this, and that is

(6.2)\begin{equation} \sum_{l=1}^{\infty} (l/a)^{2n} E(l) \approx \int_{0}^{\infty} k^{2n} E(ak) a\,\text{d} k \approx \int_{0}^{\infty} k^{2n} E_F(k) \,{\text{d}} k . \end{equation}

Since this must hold for all finite $n$, it can be argued that the spherical harmonic spectrum must converge to the Fourier spectrum in the limit of large $l$, where we must have

(6.3)\begin{equation} E(l) = \frac{E_F(l/a)}{a} . \end{equation}

For the same reason it can be argued that $G(l,n,s)$ must converge to the corresponding triadic transfer function (see Kraichnan Reference Kraichnan1967), $T(k,p,q)$, in Fourier space in the limit of large wavenumbers, where we must have

(6.4)\begin{equation} G(l,n,s) = \frac{T(l/a, n/a, s/a)}{a^3} . \end{equation}

The triangle inequality (4.3) of Silberman (Reference Silberman1954), which is a necessary condition for non-zero triad interactions in spherical harmonic space, corresponds exactly to the triangle inequality $| k-p| < q < k+p$, which is the condition in Fourier space. Thus, it is quite clear that the double cascade theory can equally well be formulated in spherical harmonic space as in Fourier space. If we use the expressions (4.18) in (5.9) and (5.10), make the approximation $l(l+1) \approx l^2$ and let the sums pass over to integrals, the similarity analysis of Kraichnan (Reference Kraichnan1967) can be repeated in every detail in spherical harmonic space. It may be objected that the approximation $l(l+1) \approx l^2$, which must be made in order for the similarity assumption to be consistent, indicates that the theory cannot be as accurately expressed in spherical harmonic space as in Fourier space. However, the assumptions of isotropy and similarity are only justified in the limit of large wavenumbers and in this limit the approximation $l(l+1) \approx l^2$ is also justified. The theory predicts that the spherical harmonic spectrum should scale as

(6.5)\begin{equation} E(l) = \frac{C |\varPi_E|^{2/3} (l/a)^{{-}5/3}}{a} , \end{equation}

in the constant energy flux range and

(6.6)\begin{equation} E(l) = \frac{\kappa \varPi_{{\mathcal{E}}}^{2/3} (l/a)^{{-}3}}{a} , \end{equation}

in the constant enstrophy flux range.

6.2. The condensate may be reached in two different ways

As pointed out in the introduction, a number of numerical investigations show that the energy spectrum in the constant energy flux range can be considerably steeper than $k^{-5/3}$, sometimes as steep as $k^{-3}$, as in the simulation of Boure (Reference Boure1994). Common features of these simulations are that viscous damping is weak, both in the enstrophy and the energy cascade range, allowing for built up of non-local interactions, that a multitude of vortices is formed in the upscale cascade range and that the peak of the energy spectrum is moving very slowly towards small wavenumbers. Although the characteristic radius of the coherent vortices is much smaller than the side, $L$, of the computational domain, the state is in several ways similar to the end state in the simulation by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007) of a condensate, as also pointed out by Boffetta & Ecke (Reference Boffetta and Ecke2012). The condensate consists of two counter-rotating vortices of size ${\sim }L$ with energy spectrum of the form $k^{-3}$. Of course, the latter state has a higher degree of regularity than the former. However, as shown by Burgess & Scott (Reference Burgess and Scott2017, Reference Burgess and Scott2018), the former has also a quite high degree of regularity, which reflects itself in the distribution of vortices with respect to radius, which can be predicted using scaling arguments.

These observations indicate that the condensate is not necessarily reached through a local cascade as predicted by Kraichnan (Reference Kraichnan1967) but may also be reached through a non-local cascade with a flow field dominated by coherent vortices and a steep energy spectrum, which slowly builds up in the direction from small to large wavenumbers, in contrast to the steep spectrum in the simulation by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007), which was built up in the opposite way. The high degree of regularity and the steep energy spectrum of the final state suggest that it is dominated by highly non-local interactions, which is also in line with the fact that the number of local triads decreases with decreasing wavenumber. Consider the nonlinear transfer to the lowest-order interacting mode, $l= 2$, in spherical harmonic space, which is also equal to the negative energy flux through the mode $l=3$,

(6.7)\begin{equation} T(2) ={-} \varPi_{E}(3) = \sum_{n=3}^{\infty} G(2,n,n+1) . \end{equation}

If we use the threshold $n/2 \le 5$ for locality, there are eight local triads, while the rest are non-local. It would be surprising if local triad interactions were dominant.

7. Applications

In Earth's atmosphere, there is a relatively localised kinetic energy source at $l \in [10, 20]$ (see Augier & Lindborg Reference Augier and Lindborg2013), where available potential energy (Lorenz Reference Lorenz1955) is transferred to kinetic energy in a process involving baroclinic instability (see Vallis Reference Vallis2006). The subsequent energy transfer in spherical harmonic space is complicated by the fact that the velocity field is not strictly rotational (or non-divergent). Although the rotational field is totally dominant in terms of magnitude, the divergent field is transferring a non-negligible part of the energy in the downscale direction. On the other hand, the rotational field is transferring energy only in the upscale direction, which has been shown in a number of studies (e.g. Boer & Sheperd Reference Boer and Sheperd1983; Augier & Lindborg Reference Augier and Lindborg2013; Brune & Becker Reference Brune and Becker2013; Burgess, Erler & Shepherd Reference Burgess, Erler and Shepherd2013). Energy is transferred upscale at a rate ${\sim }1 \, {\textrm {W}}\,{\textrm {m}}^{-2}$, integrated over the height of the troposphere. This is the mechanism by which the general circulation of the atmosphere is maintained, as also predicted by Lorenz (Reference Lorenz1955), based on an analysis using a mean-field/eddy-field decomposition. There can be little doubt that this transfer is caused by first type interactions as predicted by Fjørtoft (Reference Fjørtoft1953). As shown by analyses of Boer & Sheperd (Reference Boer and Sheperd1983) and Burgess et al. (Reference Burgess, Erler and Shepherd2013), the rotational field is also transferring a substantial amount of enstrophy downscale, which is the other signature of first type interactions.

The fact that the source is localised at such low wavenumbers, $l \sim 10-20$, implies that there is no room for a broad upscale energy cascade range to develop. As for the enstrophy transfer, the approximate $l^{-3}$-spectrum which is observed in a quite narrow range in analyses of general circulation models, may be a sign of something like an enstrophy cascade range, although the strong non-locality of enstrophy transfer makes such an interpretation questionable. In any case, three-dimensional effects are becoming important already at $l \sim 50$, so for Earth's atmosphere the theory of two-dimensional turbulence cannot be used to make much more predictions than determining the energy and enstrophy transfer directions for the large scale rotational field. As for large scale vortical structures, such as cyclones, they do not resemble the very long lived coherent vortices of two-dimensional turbulence, but are highly three-dimensional structures carrying helicity (Lilly Reference Lilly1986).

As for the atmosphere of the giant gas planets, Saturn and Jupiter, there is growing evidence from observations as well as simulations that the theory of two-dimensional turbulence is highly relevant. Evidence of an upscale energy cascade range has been extracted from analyses of Jupiter's tropospheric winds using Cassini near-infrared imaging (Choi & Showman Reference Choi and Showman2011; Galperin et al. Reference Galperin, Young, Sukoriansky, Dikovskaya, Read, Lancester and Armstrong2014; Young & Read Reference Young and Read2017). Since the dynamics is highly anisotropic and also involves Rossby waves, the assumption of isotropy which is often used in studies of two-dimensional turbulence, is not applicable. An extended concept of ‘zonostrophic turbulence’ accounting for Rossby waves and the formation of jets has been developed by Sukoriansky, Galperin & Dikovskaya (Reference Sukoriansky, Galperin and Dikovskaya2002). According to this theory, the zonal mean spectrum, which is the spectrum associated with all spherical harmonic modes with zonal wavenumber zero, should scale as $l^{-5}$, while the spectrum made up by all other modes conforms to the predictions of the Kraichnan (Reference Kraichnan1967) theory. Recently, it has been shown (Cabanes, Spiga & Young Reference Cabanes, Spiga and Young2020) that the dynamics can be successfully simulated and the spectra as well as the spectral transfers show all signs of two-dimensional turbulence, including both an upscale energy cascade and a downscale enstrophy cascade. However, the prediction of a $l^{-5/3}$-spectrum in the upscale cascade range is only partially fulfilled. In the low wavenumber end of the spectrum, there is indeed a range with $l^{-5/3}$-scaling, but for somewhat larger wavenumbers there is a local maximum (see figure 6 in Cabanes et al. Reference Cabanes, Spiga and Young2020) of the same type as seen in simulations where a multitude of vortices are formed in the upscale cascade range. It remains to be investigated if the $l^{-5/3}$-spectrum can be reproduced in different set-ups and to what degree the observed upscale energy cascade is consistent with the theory of Kraichnan (Reference Kraichnan1967) of a local cascade.