2. Dynamics

On a general two-dimensional surface, the Navier–Stokes equations are

(2.1)\begin{equation} \frac{\partial q}{\partial t}+\frac{1}{\sqrt{g}}\frac{\partial (\psi, q)}{\partial (\xi^1, \xi^2)}=D \end{equation}

and

(2.2)\begin{equation} \nabla^2_{LB}\psi \equiv \frac{1}{\sqrt{g}}\frac{\partial }{\partial \xi^i} \left( \sqrt{g} g^{ij} \frac{\partial \psi}{\partial \xi^j} \right) = q ,\end{equation}

where $q$ is the vorticity, $\psi$ is the stream function, $t$ is time and $(\xi ^1,\xi ^2)$ are coordinates covering the surface. Here, $g_{ij}$ is the metric tensor, $g^{ij}$ is its inverse, $g$ is its determinant, and repeated indices are summed from 1 to 2. Additionally, $D$ is the angular-momentum-conserving viscosity recommended by Gilbert, Riedinger & Thuburn (Reference Gilbert, Riedinger and Thuburn2014, hereafter GRT14). Furthermore, $\nabla ^2_{LB}$ is the Laplace–Beltrami operator; it reduces to the ordinary Laplacian on the plane. We also have

(2.3)\begin{equation} \dot{\xi}^i={-}\frac{1}{\sqrt{g}}\epsilon^{ij} \frac{\partial \psi}{\partial \xi^j}, \end{equation}

where the overdot denotes the time derivative following a fluid particle and $\epsilon ^{ij}$ is the permutation tensor. When $D\equiv 0$, the dynamics of (2.1) and (2.2) conserves the energy

(2.4)\begin{equation} E=\frac{1}{2} \iint {\rm d} \xi^1 \, {\rm d} \xi^2 \sqrt{g} \frac{\partial \psi}{\partial \xi^i} g^{ij} \frac{\partial \psi}{\partial \xi^j} ,\end{equation}

and the enstrophy

(2.5)\begin{equation} Z=\frac{1}{2} \iint {\rm d} \xi^1 \, {\rm d} \xi^2 \sqrt{g} q^2 . \end{equation}

On the unit sphere with spherical coordinates $(\xi ^1,\xi ^2)=(\lambda,\theta )$, where $\lambda$ is longitude and $\theta$ latitude, we have $g_{ij}=diag[\cos ^2 \theta,1]$, $g^{ij}=diag[\cos ^{-2}\theta,1]$ and $g=\cos ^2 \theta$. The dynamics of (2.1) and (2.2) take the familiar form

(2.6)\begin{equation} \frac{\partial q}{\partial t}+\frac{1}{\cos \theta}\frac{\partial (\psi, q)}{\partial (\lambda, \theta)}=D \end{equation}

and

(2.7)\begin{equation} \frac{1}{\cos^2 \theta}\frac{\partial^2 \psi}{\partial \lambda^2} + \frac{1}{\cos \theta} \frac{\partial}{\partial \theta} \left( \cos \theta \frac{\partial \psi}{\partial \theta} \right) = q . \end{equation}

On the unit sphere with the stereographic coordinates $(\xi ^1,\xi ^2)=(\xi,\eta )$ used by SP23, $g_{ij}=diag[h^2,h^2]$, $g^{ij}=diag[h^{-2},h^{-2}]$ and $g=h^2$, where

(2.8)\begin{equation} h(\xi,\eta)=\frac{2}{1+\xi^2+\eta^2}.\end{equation}

The dynamics take the simpler form

(2.9)\begin{equation} \frac{\partial q}{\partial t}+\frac{1}{h^2}\frac{\partial (\psi, q)}{\partial (\xi, \eta)}=D \end{equation}

and

(2.10)\begin{equation} \frac{1}{h^2}\left( \frac{\partial^2 \psi}{\partial \xi^2} + \frac{\partial^2 \psi}{\partial \eta^2} \right) = q. \end{equation}

Now let the surface be the ellipsoid (1.1) where $(x,y,z)$ are Cartesian embedding coordinates and $a, b, c$ are constants. We assign coordinates to the ellipsoid as follows. First, we project the point $(x,y,z)$ on the ellipsoid to a point $(X,Y,Z)$ on the unit sphere by requiring that the three points $(x,y,z)$, $(X,Y,Z)$ and $(0,0,0)$ be colinear. Refer to figure 1. This means that

(2.11)\begin{equation} (x,y,z) = R(X,Y,Z) (X,Y,Z)\end{equation}

for some scalar $R(X,Y,Z)$. Since $(x,y,z)$ lies on the ellipsoid,

(2.12)\begin{equation} \frac{(RX)^2}{a^2}+\frac{(RY)^2}{b^2}+\frac{(RZ)^2}{c^2}=1\end{equation}

and hence

(2.13)\begin{equation} R(X,Y,Z)= \left( \frac{X^2}{a^2}+\frac{Y^2}{b^2}+\frac{Z^2}{c^2} \right)^{{-}1/2}. \end{equation}

Next, as in SP23, we cover the unit sphere with two overlapping systems of stereographic coordinates. The first is defined as follows. Let $(\xi,\eta,0)$ be the point at which the line between $(X,Y,Z)$ and $(0,0,1)$ intersects the plane $z=0$. Refer again to figure 1. Then

(2.14a–c)\begin{equation} X=\frac{2\xi}{1+\xi^2+\eta^2}, \quad Y=\frac{2\eta}{1+\xi^2+\eta^2},\quad Z=\frac{\xi^2+\eta^2-1}{1+\xi^2+\eta^2}, \end{equation}

and

(2.15a,b)\begin{equation} \xi=\frac{X}{1-Z}, \quad \eta=\frac{Y}{1-Z}. \end{equation}

For the second system of stereographic coordinates, let $(\hat {\xi },\hat {\eta },0)$ be the point at which the line between $(X,Y,Z)$ and $(0,0,-1)$ intersects $z=0$. We have

(2.16a–c)\begin{equation} X=\frac{2\hat{\xi}}{1+\hat{\xi}^2+\hat{\eta}^2},\quad Y=\frac{2\hat{\eta}}{1+\hat{\xi}^2+\hat{\eta}^2},\quad Z=\frac{1-\hat{\xi}^2-\hat{\eta}^2}{1+\hat{\xi}^2+\hat{\eta}^2},\end{equation}

and

(2.17a,b)\begin{equation} \hat{\xi}=\frac{X}{1+Z}, \quad \hat{\eta}=\frac{Y}{1+Z}. \end{equation}

Figure 1. At the left, the point $(x,y,z)$ on the ellipsoid (1.1) is first projected, along the line containing $(x,y,z)$ and $(0,0,0)$, to the point $(X,Y,Z)$ on the unit sphere, and then, along the line containing $(X,Y,Z)$ and $(0,0,1)$, to the point $(\xi,\eta,0)$ on the $z=0$ plane. This defines the stereographic coordinates $(\xi,\eta )$, which cover the ellipsoid except for the point $(0,0,c)$. At the right, $(x,y,z)$ projects to $(X,Y,Z)$ on the unit sphere, and then, along the line containing $(X,Y,Z)$ and $(0,0,-1)$, to the point $(\hat {\xi },\hat {\eta },0)$ on the $z=0$ plane. This defines the stereographic coordinates $(\hat {\xi },\hat {\eta })$, which cover the ellipsoid except for the point $(0,0,-c)$.

Substituting (2.14a–c) and (2.16a–c) into (2.13) and (2.11), we obtain

(2.18)\begin{equation} (x,y,z) = S(\xi,\eta) (2\xi,2\eta,\xi^2+\eta^2-1) = S(\hat{\xi},\hat{\eta}) (2\hat{\xi},2\hat{\eta},1-\hat{\xi}^2-\hat{\eta}^2),\end{equation}

where

(2.19)\begin{equation} S(\xi,\eta)=\left( \frac{4\xi^2}{a^2}+\frac{4\eta^2}{b^2}+\frac{(\xi^2+\eta^2-1)^2}{c^2} \right)^{{-}1/2}. \end{equation}

In the case of the unit sphere, $(x,y,z)=(X,Y,Z)$, and (2.18) reduces to the coordinate assignment used by SP23.

The mapping (2.18) between $(x,y,z)$ on the ellipsoid and the covering coordinates $(\xi,\eta )$ or $(\hat {\xi },\hat {\eta })$ is one-to-one for all non-vanishing $a,b,c$. It could easily be generalized to an arbitrary, non-ellipsoidal surface provided only that the surface did not contain folds. Both systems of coordinates, $(\xi,\eta )$ and $(\hat {\xi },\hat {\eta })$, cover the ellipsoid except for their respective points at infinity. For the $(\xi,\eta )$ coordinates, the point at infinity corresponds to the point $z=c$, and the area within the unit circle $\xi ^2+\eta ^2<1$ corresponds to the hemi-ellipsoid $z<0$. For the $\hat {\xi },\hat {\eta }$ coordinates, the point at infinity corresponds to the point $z=-c$, and the area within the unit circle $\hat {\xi }^2+\hat {\eta }^2<1$ corresponds to the hemi-ellipsoid $z>0$. Every point on the ellipsoid corresponds to a value of $(\xi,\eta )$ and a value of $(\hat {\xi },\hat {\eta })$, which generally differ. Only at $z=0$ do we have $(\xi,\eta )=(\hat {\xi },\hat {\eta })$ on the ellipsoid. The general relation between the two coordinate systems is

(2.20a,b)\begin{equation} (\xi,\eta)=\frac{(\hat{\xi},\hat{\eta})}{\hat{\xi}^2+\hat{\eta}^2},\quad (\hat{\xi},\hat{\eta})=\frac{(\xi,\eta)}{\xi^2+\eta^2}. \end{equation}

Thus, points within the unit circle on the $\xi,\eta$-plane correspond to points outside the unit circle on the $\hat {\xi },\hat {\eta }$-plane, and vice versa. For analytical developments, it is best to restrict one's attention to a single coordinate system, but numerical computation is a different matter: there is great advantage to using the $\xi,\eta$ coordinates in the domain $\xi ^2+\eta ^2<1$ to cover the ‘southern’ hemi-ellipsoid, and using the $\hat {\xi },\hat {\eta }$ coordinates in the domain $\hat {\xi }^2+\hat {\eta }^2<1$ to cover the ‘northern’ hemi-ellipsoid.

For infinitesimal displacements on the ellipsoid,

(2.21)\begin{equation} {\rm d}s^2={{\rm d}\kern0.7pt x}^2+{{\rm d} y}^2+{\rm d} z^2=g_{\xi \xi}\, {\rm d}\xi^2+g_{\eta \eta}\, {\rm d}\eta^2 +2 g_{\xi \eta}\, {\rm d}\xi \, {\rm d}\eta. \end{equation}

From (2.18), we obtain expressions for the metric components $g_{\xi \xi }$, $g_{\eta \eta }$, $g_{\xi \eta }$ and the determinant $g=g_{\xi \xi }g_{\eta \eta }-{g_{\xi \eta }}^2$ as functions of $\xi,\eta$ or $\hat {\xi },\hat {\eta }$. For example,

(2.22)\begin{equation} g_{\xi \xi}=\left(\frac{\partial x}{\partial \xi}\right)^2 + \left(\frac{\partial y}{\partial \xi}\right)^2 + \left(\frac{\partial z}{\partial \xi}\right)^2 \end{equation}

and

(2.23)\begin{equation} g_{\xi \eta}=\frac{\partial x}{\partial \xi}\frac{\partial x}{\partial \eta}+ \frac{\partial y}{\partial \xi}\frac{\partial y}{\partial \eta}+ \frac{\partial z}{\partial \xi}\frac{\partial z}{\partial \eta}. \end{equation}

Note that the metric coefficients have the same dependence on $\xi,\eta$ as on $\hat {\xi },\hat {\eta }$. For the sphere and for the spheroids, $g_{\xi \eta }= 0$, and hence the coordinates $(\xi,\eta )$ and $(\hat {\xi },\hat {\eta })$ are orthogonal. On the triaxial ellipsoid, $g_{\xi \eta }\neq 0$, and therefore the coordinates are non-orthogonal. The inverse components are $g^{\xi \xi }=g_{\eta \eta }/g$, $g^{\eta \eta }=g_{\xi \xi }/g$ and $g^{\xi \eta }=-g_{\xi \eta }/g$. If $g_{\xi \eta }=0$, then $g^{\xi \xi }=1/g_{\xi \xi }$ and $g^{\eta \eta }=1/g_{\eta \eta }$.

As in SP23, our strategy will be to solve

(2.24)\begin{equation} \sqrt{g} \frac{\partial q}{\partial t}-\frac{\partial (\psi, q)}{\partial (\xi, \eta)}=\sqrt{g}D, \nabla^2_{LB}\psi=q ,\end{equation}

on $\xi ^2+\eta ^2<1$; to solve

(2.25)\begin{equation} \sqrt{g} \frac{\partial q}{\partial t}+\frac{\partial (\psi, q)}{\partial (\hat{\xi},\hat{\eta})}=\sqrt{g}D, \nabla^2_{LB}\psi = q, \end{equation}

on $\hat {\xi }^2+\hat {\eta }^2<1$; and to match the two solutions together at the ‘equator’, $\xi ^2+\eta ^2=\hat {\xi }^2+\hat {\eta }^2=1$. The matching conditions are that $\psi$ and $q$ be continuous at the equator. The sign flip between (2.24) and (2.25) reflects the fact that, as viewed from the exterior of the ellipsoid, $\xi,\eta$ are left-handed coordinates, and $\hat {\xi },\hat {\eta }$ are right-handed coordinates. Our method is designed to be such that the solved values of $\psi$ and $q$ are the same as would be obtained by solving the entire problem in spherical coordinates. However, our method completely avoids the pole problem. For further details, see SP23. For more background on stereographic coordinates, see Needham (Reference Needham1997, Reference Needham2021).

The energy and enstrophy take the forms

(2.26)\begin{align} E&=\frac{1}{2} \iint {\rm d} \xi \, {\rm d} \eta \sqrt{g} \left( \frac{\partial \psi}{\partial \xi} g^{\xi\xi} \frac{\partial \psi}{\partial \xi}+ \frac{\partial \psi}{\partial \eta} g^{\eta\eta} \frac{\partial \psi}{\partial \eta}+ 2 \frac{\partial \psi}{\partial \xi} g^{\xi\eta} \frac{\partial \psi}{\partial \eta} \right) \nonumber\\ &\quad +\frac{1}{2} \iint {\rm d} \hat{\xi} \, {\rm d} \hat{\eta} \sqrt{g} \left( \frac{\partial \psi}{\partial \hat{\xi}} g^{\hat{\xi}\hat{\xi}} \frac{\partial \psi}{\partial \hat{\xi}}+ \frac{\partial \psi}{\partial \hat{\eta}} g^{\hat{\eta}\hat{\eta}} \frac{\partial \psi}{\partial \hat{\eta}}+ 2 \frac{\partial \psi}{\partial \hat{\xi}} g^{\hat{\xi}\hat{\eta}} \frac{\partial \psi}{\partial \hat{\eta}} \right) \end{align}

and

(2.27)\begin{equation} Z=\frac{1}{2} \iint {\rm d} \xi \, {\rm d}\eta \sqrt{g} q^2 + \frac{1}{2} \iint {\rm d}\hat{\xi} \, {\rm d}\hat{\eta} \sqrt{g} q^2 ,\end{equation}

where the integrals are over the interiors of the two unit circles. We define the angular momentum,

(2.28)\begin{equation} \boldsymbol{M}= \iint {\rm d} \xi \, {\rm d} \eta \sqrt{g} ( \boldsymbol{x}(\xi,\eta) \times \dot{\boldsymbol{x}}(\xi,\eta) ) + \iint {\rm d} \hat{\xi} \, {\rm d} \hat{\eta} \sqrt{g} ( \boldsymbol{x}(\hat{\xi},\hat{\eta}) \times \dot{\boldsymbol{x}}(\hat{\xi},\hat{\eta}) ),\end{equation}

whether or not it is conserved. To express the time derivatives $\dot {\boldsymbol {x}}(\xi,\eta )$ and $\dot {\boldsymbol {x}}(\hat {\xi },\hat {\eta })$ in terms of the stream function, we take the time derivative of (2.18), and set

(2.29a,b)\begin{equation} \dot{\xi}=\frac{1}{\sqrt{g}} \frac{\partial \psi}{\partial \eta},\quad \dot{\eta}={-}\frac{1}{\sqrt{g}} \frac{\partial \psi}{\partial \xi}\end{equation}

and

(2.30a,b)\begin{equation} \dot{\hat{\xi}}={-}\frac{1}{\sqrt{g}} \frac{\partial \psi}{\partial \hat{\eta}},\quad \dot{\hat{\eta}}=\frac{1}{\sqrt{g}} \frac{\partial \psi}{\partial \hat{\xi}}.\end{equation}

The sign flip between (2.29a,b) and (2.30a,b) reflects the difference in handedness between $\xi,\eta$ and $\hat {\xi },\hat {\eta }$. Then using (2.28), (2.29a,b) and (2.30a,b), we find that $\boldsymbol {M}=(M_1,M_2,M_3)$ where

(2.31)\begin{gather} M_1=2\iint {\rm d} \xi \, {\rm d} \eta \psi \frac{\partial (y,z)}{\partial (\xi,\eta)} +2\iint {\rm d}\hat{\xi} \, {\rm d}\hat{\eta} \psi \frac{\partial (y,z)}{\partial (\hat{\xi},\hat{\eta})} , \end{gather}

(2.32)\begin{gather}M_2=2\iint {\rm d} \xi \, {\rm d} \eta \psi \frac{\partial (z,x)}{\partial (\xi,\eta)} +2\iint {\rm d}\hat{\xi} \, {\rm d}\hat{\eta} \psi \frac{\partial (z,x)}{\partial (\hat{\xi},\hat{\eta})}, \end{gather}

(2.33)\begin{gather}M_3=2\iint {\rm d} \xi \, {\rm d}\eta \psi \frac{\partial (x,y)}{\partial (\xi,\eta)} - 2\iint {\rm d}\hat{\xi} \, {\rm d}\hat{\eta} \psi \frac{\partial (x,y)}{\partial (\hat{\xi},\hat{\eta})}. \end{gather}

Once again, when $a=b=c$ all three components $M_1, M_2, M_3$ are conserved; when $a=b\neq c$ only the $z$-component, $M_3$, is conserved; and when no two of $a,b,c$ are equal, no component of $\boldsymbol {M}$ is conserved. These conservation laws depend upon the physics ((2.1) and (2.2)) as well as the geometry of the spheroid, as reflected in its metric components.

The viscosity

(2.34)\begin{equation} D=\nu g^{ij}\nabla_i \nabla_j q + 2\nu \nabla_i ( g^{ij} K \nabla_j \psi ),\end{equation}

recommended by GRT14 (equivalent to their (3.61)) conserves angular momentum when the ellipsoid possesses the required rotational symmetry. Here, $\nu$ is the viscosity coefficient, $\nabla _i$ is the covariant derivative and $K(\xi,\eta )$ is the Gaussian curvature. Using tensor identities, we obtain $g^{ij}\nabla _i \nabla _j = \nabla ^2_{LB}$, and hence (2.34) is equivalent to

(2.35)\begin{equation} D=\nu \nabla^2_{LB}q + 2 \nu K q + 2 \nu \frac{\partial K}{\partial \xi^i} g^{ij} \frac{\partial \psi}{\partial \xi^j}. \end{equation}

On the unit sphere, $K\equiv 1$ and (2.35) reduces to the viscosity

(2.36)\begin{equation} D=\nu ( \nabla_{LB}^2 q + 2 q ),\end{equation}

proposed by Willams (Reference Willams1972) and Becker (Reference Becker2001). (SP23 erroneously omitted the last term in (2.36).) As noted by GRT14, in the case of the sphere, it is easy to see that (2.36) conserves angular momentum while dissipating energy and enstrophy: if $q$ is expanded in the spherical harmonics $Y_n^m$, then the coefficients $q_n^m(t)$ obey

(2.37)\begin{equation} \frac{\partial q_n^m}{\partial t}+\cdots = \nu ({-}n(n+1) + 2 )q_n^m, \end{equation}

where the ellipses denote the advective terms. The three conserved components of angular momentum correspond to the three amplitudes $q_1^0, q_1^{\pm 1}$ of degree $n=1$. The viscosity (2.37) vanishes on these three modes, but all higher modes experience dissipation. On the ellipsoid, the Gaussian curvature $K$ generally varies with location. As a function of the embedding coordinates, it is

(2.38)\begin{equation} K=\frac{ 1 } {a^2b^2c^2 \left( \dfrac{x^2}{a^4}+ \dfrac{y^2}{b^4}+ \dfrac{z^2}{c^4} \right)^2 }.\end{equation}

Using (2.18), we obtain the useful formula

(2.39)\begin{equation} K(\xi,\eta)= \frac{ \left( \dfrac{4\xi^2}{a^2}+ \dfrac{4\eta^2}{b^2}+ \dfrac{(\xi^2+\eta^2-1)^2}{c^2} \right)^2 } {a^2b^2c^2 \left( \dfrac{4\xi^2}{a^4}+ \dfrac{4\eta^2}{b^4}+ \dfrac{(\xi^2+\eta^2-1)^2}{c^4} \right)^2 },\end{equation}

from which the derivatives of $K(\xi,\eta )$ can be computed analytically for use in (2.35).

We emphasize that, in our idealized dynamics, the ellipsoid serves only to shape the flow; there is no force between the ellipsoid and the fluid except for the pressure torque that arises from non-sphericity. The fluid ‘rotates’ if it has a conserved component of angular momentum, and the ‘rotation rate’ is constant because angular momentum is conserved. This is in contrast with the Earth-atmosphere system in which the angular momentum of the atmosphere is not by itself conserved, but is leading-order constant because of the torque exerted by the Earth's surface, and because the solid planet has an enormous moment of inertia. In our case, the ‘rotation rate’ is perfectly constant; in the Earth's case, it is nearly constant, as demonstrated by the very small changes in the measured length of day.

Our decision to view all of our results in the inertial reference frame is sensible because of this same lack of frictional torque between the ellipsoid and the fluid. Suppose, on the sphere, that the fluid had non-vanishing angular momentum about the $z$-axis corresponding to solid-body rotation at angular velocity $\varOmega$. If this situation were viewed from a frame rotating with the solid-body flow, then the other two angular momentum components would obey

(2.40a,b)\begin{equation} \frac{{\rm d} M_1}{{\rm d} t}-\varOmega M_2=0,\quad \frac{{\rm d} M_2}{{\rm d} t}+\varOmega M_1=0\end{equation}

(Egger, Weichmann & Hoinka Reference Egger, Weichmann and Hoinka2007, p. 3). In the inertial frame, $M_1$ and $M_2$ are independent constants. (Note that (1.2) vanishes when $n=1$.) Our affection for the inertial frame would certainly disappear if we were to consider an ellipsoid that rotated about an axis that was not an axis of symmetry. This is actually the situation on Earth, if one takes Earth's full surface topography into account.

The presence of angular momentum, whether conserved or not, affects the speed at which energy is transferred between modes. Wave dispersion caused by non-zero $\varOmega$ interferes with the transfer of energy between modes, slowing the inverse cascade in the way first described by Rhines (Reference Rhines1975). Partly for this reason, all of the experiments described in the next section were begun from states of vanishing angular momentum.

It is worth emphasizing that three conservation laws – energy, enstrophy and angular momentum – are fundamentally different constraints. Energy conservation is associated with time independence of the metric components. Enstrophy conservation arises from the fluid-particle relabelling symmetry. Angular momentum conservation arises from the invariance of the metric components with respect to rotation. Specifically, if coordinates $(r,\phi )$ can be found such that a rotation corresponds to a translation in $\phi$ with no change in $r$, and if in these coordinates, the metric components depend only on $r$ and not on $\phi$, then the angular momentum with respect to the axis of rotation is conserved. Angular momentum is an unfamiliar ingredient of two-dimensional turbulence theory, because most planar geometries, including the infinitely periodic box, lack rotational symmetry.

The ellipsoidal harmonics $Y_p(\xi ^1,\xi ^2)$ are defined by

(2.41)\begin{equation} \nabla^2_{LB} Y_p \equiv \frac{1}{\sqrt{g}}\frac{\partial }{\partial \xi^i} \left( \sqrt{g} g^{ij} \frac{\partial Y_p}{\partial \xi^j} \right) ={-}\lambda_p Y_p,\end{equation}

where $p$ is a discrete index that denotes the harmonic, and $\lambda _p$ is the corresponding eigenvalue. Both $Y_p$ and $\lambda _p$ depend parametrically on $a,b,c$. When $a=b=c=1$, the $Y_p$ are the familiar spherical harmonics $Y_n^m$ with eigenvalue $\lambda _{n,m}=n(n+1)$. On the ellipsoid, as on the sphere, the $Y_p$ comprise a complete set of orthogonal functions. Upon normalization, we have

(2.42)\begin{equation} \iint {\rm d} \xi^1\, {\rm d} \xi^2 \sqrt{g} Y_p Y_k=\delta_{pk}.\end{equation}

For a thorough treatment of ellipsoidal harmonics, see Dassios (Reference Dassios2012). The usual development is in terms of Lamé functions, but for ellipsoids close to the sphere, an expansion in terms of the spherical harmonics can be useful. Eswarathasan & Kolokolnikov (Reference Eswarathasan and Kolokolnikov2022) give expressions for the eigenvalues of ellipsoids that are close to the unit sphere. For the spheroid with $a=b=1+\epsilon \alpha$ and $c=1+\epsilon \beta$, they obtain

(2.43)\begin{equation} \lambda=n(n+1)+\epsilon \left[ (\alpha-\beta)\frac{2n(n+1)(2n^2-2~m^2+2n-1)}{(2n+3)(2n-1)} -2\alpha n(n+1) \right] +O(\epsilon^2),\end{equation}

where $n$ and $m$ are integers, and $-n< m< n$. If $\epsilon =0$, there are $2n+1$ distinct harmonics with the same eigenvalue $n(n+1)$. For $\epsilon \neq 0$, the eigenvalues are degenerate only in pairs ($\pm m$).

Let the vorticity and stream function be expanded in ellipsoidal harmonics,

(2.44a,b)\begin{equation} q(\xi^1,\xi^2,t)=\sum_p q_p(t) Y_p(\xi^1,\xi^2),\quad \psi(\xi^1,\xi^2,t)=\sum_p \psi_p(t) Y_p(\xi^1,\xi^2), \end{equation}

where $q_p=-\lambda _p \psi _p$. The dynamics becomes

(2.45)\begin{equation} \frac{{\rm d}q_i}{{\rm d} t}=\frac{1}{2}\sum_j \sum_k \frac{ (\lambda_k-\lambda_j ) }{\lambda_j \lambda_k} q_j q_k \iint {\rm d} \xi^1 \, {\rm d} \xi^2\frac{\partial (Y_j,Y_k)}{\partial(\xi^1,\xi^2)}Y_i.\end{equation}

The energy and enstrophy take the forms

(2.46a,b)\begin{equation} E=\frac{1}{2} \sum_p \lambda_p \psi_p^2,\quad Z=\frac{1}{2} \sum_p \lambda_p^2 \psi_p^2,\end{equation}

just as in the spherical case. However, only on the sphere are the angular momenta proportional to single harmonics. To see this, rewrite (2.31)–(2.33) in arbitrary coordinates as

(2.47)\begin{equation} \boldsymbol{M}=\iint {\rm d} \xi^1 \, {\rm d} \xi^2 \sqrt{g} (Q_1,Q_2,Q_3) \psi, \end{equation}

where

(2.48)\begin{equation} (Q_1,Q_2,Q_3) = \frac{2}{\sqrt{g}} \left( \frac{\partial (y,z)}{\partial (\xi^1,\xi^2)}, \frac{\partial (z,x)}{\partial (\xi^1,\xi^2)} , \frac{\partial (x,y)}{\partial (\xi^1,\xi^2)} \right) ,\end{equation}

and the integration is over the entire ellipsoid. To compute the function $Q_3$, we can use any set of covering coordinates $\xi ^1,\xi ^2$, but the most convenient choice is $(\xi ^1,\xi ^2)=(x,y)$. Then $Q_3=2/\sqrt {g}$, where $g$ is the metric determinant corresponding to the coordinates $x,y$. Setting $a=b$ for the sphere or the spheroid, and using

(2.49) \begin{equation} {\rm d}s^2={{\rm d}\kern0.7pt x}^2+{{\rm d} y}^2+\left({\rm d}\sqrt{1-x^2-y^2}\right)^2=g_{xx}{\,{\rm d}\kern0.7pt x}^2+2g_{xy}\, {\rm d}\kern0.7pt x\, {\rm d}y+g_{yy}{\,{\rm d} y}^2\end{equation}

and $g=g_{xx}g_{yy}-g_{xy}g_{xy}$, we obtain

(2.50)\begin{equation} Q_3=2\sqrt{\frac{1-(x^2+y^2)}{1+\epsilon (x^2+y^2)}}, \end{equation}

where $\epsilon =c^2/a^2-1$. In the case of the unit sphere, $\epsilon =0$ and $Q_3=2z=2\sin \theta$, which is proportional (apart from normalization constant) to the spherical harmonic $Y_1^0$. Similarly, $Q_1=2x$ and $Q_2=2y$ on the sphere.

For the spheriod, $\epsilon \ne 0$. Expanding (2.50) in $\epsilon$ and setting $x^2+y^2=a^2(1-{z^2}/{c^2})$, we obtain an expansion for $Q_3$ in odd powers of $z$. Then, using the definitions of the spherical harmonics in terms of $x,y,z$, we can express $Q_3$ as the weighted sum of spherical harmonics of the form $Y_{2n+1}^0$ with leading term proportional to $Y_1^0$. Thus, $Q_3$ for the spheroid is not proportional to a single spherical harmonic. However, neither is it proportional to a single ellipsoidal harmonic, as can be demonstrated by computing $(\nabla ^2_{LB}Q_3)/Q_3$ and observing that this quotient is not independent of location on the spheroid.

3. Numerics

We consider freely decaying turbulence on four distinct ellipsoids: the unit sphere (S), a prolate spheroid (P), an oblate spheroid (O) and a triaxial ellipsoid (T). The values of $a,b,c$ were chosen such that the ellipsoids, in all four cases, have the same surface area, $4{\rm \pi}$, as the unit sphere. Table 1 summarizes the geometry. Table 1 also gives the eigenvalues of the four lowest modes according to the formulae of Eswarathasan & Kolokolnikov (Reference Eswarathasan and Kolokolnikov2022) for the spheroid (given above as (2.43)) and for the triaxial ellipsoid (their (22)–(25)). The three eigenvalues given as $\lambda _1, \lambda _2, \lambda _3$ correspond to harmonic degree $n=1$. Here, $\lambda _4$ is the smallest of the five eigenvalues corresponding to $n=2$.

Table 1. The geometry of the ellipsoids used for the numerical experiments.

In the case of the sphere, all three of the eigenvalues associated with $n=1$ have the same value, $1(1+1)=2$. The corresponding eigenmodes are the three spherical harmonics, $Y_1^0, Y_1^{\pm 1}$.

The two spheroids, P and O, have a pair of degenerate eigenvalues, $\lambda _1$ and $\lambda _2$. The corresponding two eigenfunctions differ by a 90-degree rotation about the axis of symmetry (the $z$-axis) of the spheroids. These two modes, which have the same latitudinal structure but vary as $\cos \lambda$ and $\sin \lambda$ with longitude, are the lowest modes – with the smallest eigenvalues – on the spheroids, and are hence the modes into which, according to theory, all energy will eventually flow. The mode corresponding to $\lambda _3$ is inhibited from accepting energy by the conservation of angular momentum, but it is not entirely forbidden, because, as explained in § 2, this mode is only the leading order contributor to the conserved component $M_3$. It is interesting that the eigenvalues of $P$ are all higher than the corresponding eigenvalues of $S$, while those of $O$ are lower than those of $S$.

In the case of the general triaxial ellipsoid, no two eigenvalues are equal. No component of angular momentum is conserved, and energy can flow into any mode.

Our numerical method is a generalization of the method described in detail in SP23. The generalization only requires the insertion of metric components in appropriate places. We use a Strang splitting method that alternates between the Jacobian dynamics, i.e. (2.1) with $D=0$, and viscous dynamics, i.e. (2.1) with the Jacobian term omitted and $D$ taken as (2.35). The Jacobian dynamics conserves energy (2.4) and enstrophy (2.5) except for truncation error in the time step; it conserves angular momentum (as geometry permits) except for truncation error in the time step and in the spatial differencing. The viscous dynamics correctly dissipates energy and enstrophy, and it conserves angular momentum except for both temporal and spatial truncation error. Our spatial resolution is equivalent to a harmonic cutoff at degree $n=240$. It corresponds to a grid spacing of approximately 94 km on Earth's surface.

The four experiments begin from nearly identical initial conditions that are assigned as follows. Let $q_S(X,Y,Z)$ be a sum of spherical harmonics in the range $n=5,6,7$ with randomly assigned amplitudes. On all four ellipsoids, we set $q(x,y,z,0)=q_S(X,Y,Z)$ where $(x,y,z)$ is related to $(X,Y,Z)$ by (2.11) and (2.13). This initial condition is then re-scaled in such a way to make the root mean square (r.m.s.) velocity $u_{rms}=1$. This means that a time interval of unity is required for a fluid particle to move a distance equal to the ‘radius’ of the ellipsoid. Figure 2 shows the initial stream function as viewed from $x=\infty$ in the four experiments. In all four experiments, we take $\nu =q_{rms}\varDelta ^2$ where $q_{rms}$ is the r.m.s. vorticity and $\varDelta$ is the grid-spacing on the stereographic plane. As time increases, $u_{rms}$, $q_{rms}$ and $\nu$ all slowly decrease, but viscosity is insignificant at late times because by then, nearly all of the energy is in the lowest modes. By $t=10$, $u_{rms}$ has decreased to approximately $0.90$ in all four experiments, but the energy decay thereafter is quite small; at $t=200$, $u_{rms}$ still ranges between $0.71$ and $0.82$ in the four experiments. Figure 3 shows the vorticity at $t=2.5$, by which time the flow is fully turbulent. The turbulence gradually subsides as the energy reaches the lower modes and Bose–Einstein condensation begins to occur.

Figure 2. The stream function at $t=0$ on the (a) sphere, (b) prolate spheroid, (c) oblate spheroid and (d) triaxial ellipsoid. The initial conditions consist of modes of degrees $n=5,6,7$ with randomly assigned amplitudes. Darker contours correspond to larger values.

Figure 3. The vorticity at $t=2.5$ on the (a) sphere, (b) prolate spheroid, (c) oblate spheroid and (d) triaxial ellipsoid. By this time, the flow is fully turbulent. Darker contours correspond to larger values.

Tables 2 and 3 summarize the results. Table 2 displays the angular momentum components $\mu _i=M_i/4{\rm \pi}$ per unit area, where $M_i$ is computed from (2.31)–(2.33). Division by area makes the ‘scale size’ of $\mu _i$ order one. Thus, $\mu _i \ll 1$ is insignificant. All of our experiments begin with $\mu _i\approx 0$. The asterisks in table 2 denote angular momentum components that are conserved and should therefore remain zero. Differences from that are caused by accumulated numerical truncation error, which is seen to be quite small. The non-conserved angular momentum components generally increase in P, O and T, but the increase is not monotonic, and sign changes even occur. The sign of the non-conserved components of angular momentum is an accident of the initial conditions.

Table 2. Evolution of $\mu _i=M_i/4{\rm \pi}$, the angular momenta per unit area, in the numerical experiments. The asterisks denote components that are conserved by the dynamics.

Table 3. The fraction, $E_n$, of energy in the modes of degree $n$. The energy is initially concentrated in degrees $n=5,6,7$. As time increases, it moves towards the lowest $n$ that is not forbidden by angular momentum conservation.

Table 3 shows the fraction, $E_n$, of total energy in the modes ‘associated’ with degree $n$. By this, we mean the following: the vorticity and stream function at $(x,y,z)$ on the ellipsoid are first projected onto the unit sphere at $(X,Y,Z)$ and then analysed by projection onto spherical harmonics, in a reverse of the method used to assign initial conditions. The better method would be a spectral analysis in terms of ellipsoidal harmonics themselves. These however are very difficult to calculate. Therefore, the data in table 3 must be regarded as an approximation to the true energy spectrum which is good only to the extent that the spherical harmonics resemble the corresponding ellipsoidal harmonics. The eigenvalues recorded in table 1 suggest that this resemblance is adequate: at our values of $a,b,c$, there is no overlap between the eigenvalues $\lambda _1$, $\lambda _2$, $\lambda _3$ corresponding to $n=1$ and the lowest eigenvalue, $\lambda _4$, corresponding to $n=2$.

Table 3 shows that angular momentum conservation on the sphere prevents energy from entering the three lowest modes. By $t=120$ – the time required for a fluid particle to circumnavigate the sphere approximately 20 times – fully 92 % of the remaining energy has condensed into $n=2$, but none has crossed the angular-momentum barrier into $n=1$. For the prolate spheroid, 64 % of the energy has crossed into $n=1$ and most of the remainder resides in $n=2$, waiting to cross. The corresponding figures for the oblate spheroid are 40 % and 48 %. The condensation of energy into $n=1$ on the triaxial ellipsoid proceeds more slowly than on either of the two spheroids. At $t=120$, only 19 % of its energy has reached the lowest mode.

Figures 4(a,b) and 4(c,d) show two views of the stream function on the sphere and on the prolate spheroid, respectively, at $t=120$. The two views are by an infinitely distant observer at (a,c) $y=z=0$ and $x=\infty$, and at (b,d) $y=z=0$ and $x=-\infty$. Thus, the two views offer a complete picture. At this late time, the vorticity strongly resembles the stream function. The sphere exhibits four vortices, which could be anticipated from its dominantly $n=2$ composition (table 3). There is of course no reason why any of these vortices should favour a particular direction on the sphere. At $t=120$, the prolate spheroid exhibits two vortices in keeping with its dominantly $n=1$ composition. As expected, the centres of these two vortices lie in the $xy$-plane. They are analogous to the $Y_1^{\pm 1}$ modes, which, unlike the $Y_1^0$ mode, are not prevented by angular momentum conservation from accumulating turbulent energy. These modes are the lowest modes in the system (table 1) and are therefore the site of Bose–Einstein condensation. The corresponding views of the oblate spheroid and the triaxial ellipsoid (not shown) are less easy to characterize, because the energy of the oblate spheroid at $t=120$ is still about equally divided between $n=1$ and $n=2$ (table 3), and the energy of the triaxial ellipsoid is still dominantly in $n=2$.

Figure 4. Two views of the stream function on the (a,b) sphere and on the (c,d) prolate spheroid at $t=120$. Darker contours correspond to larger values.

Figure 5 shows the same two views of the stream function on the (a,b) oblate spheroid and on (c,d) the triaxial ellipsoid at the final time, $t=200$. At $t=200$, 56 % of the energy has reached $n=1$ on the oblate spheroid, and the stream function (figure 5a,b) shows vortices centred on $z=0$, as predicted by theory. The triaxial ellipsoid, with its complete lack of symmetry, shows no such tendency (figure 5c,d).

Figure 5. Two views of the stream function on the (a,b) oblate spheroid and on the (c,d) triaxial ellipsoid at the final time $t=200.$ Darker contours correspond to larger values.

It is interesting that the prolate spheroid, for which the adjustment to equilibrium is relatively rapid has lowest eigenvalue $2.480$ (table 1) that is displaced upward from the eigenvalues of the sphere, whereas the lowest eigenvalues of the oblate spheroid and triaxial ellipsoid (respectively $1.680$ and $1.838$) are downwardly displaced. One is tempted to say that the adjustment takes longer if it has farther to go. In all cases, except perhaps that of the sphere, this adjustment is extremely slow, but this agrees with previous results; see especially Dritschel et al. (Reference Dritschel, Qi and Marston2015).