1 Introduction
Observations of Jupiter’s atmosphere reveal a highly turbulent cloud deck containing long-lived coherent vortices, such as the Great Red Spot, that are transported around the planet by an alternating pattern of zonal jets (Vasavada & Showman Reference Vasavada and Showman2005). Following the pioneering non-divergent barotropic model of Williams (Reference Williams1978), shallow water theory has been widely used to model the Jovian atmosphere (Williams & Yamagata Reference Williams and Yamagata1984; Marcus Reference Marcus1988; Dowling & Ingersoll Reference Dowling and Ingersoll1989; Ingersoll Reference Ingersoll1990; Cho & Polvani Reference Cho and Polvani1996a ,Reference Cho and Polvani b ; Iacono, Struglia & Ronchi Reference Iacono, Struglia and Ronchi1999a ; Iacono et al. Reference Iacono, Struglia, Ronchi and Nicastro1999b ; Scott & Polvani Reference Scott and Polvani2007; Showman Reference Showman2007; Scott & Polvani Reference Scott and Polvani2008). The visible cloud deck is treated as a homogeneous dynamically active layer separated from a much deeper and relatively quiescent lower layer by a sharp density contrast, which may be linked to latent heat being released at the water vapour condensation level. The equivalent barotropic approximation (Gill Reference Gill1982; Dowling & Ingersoll Reference Dowling and Ingersoll1989) gives a closed set of shallow water equations for the flow in the cloud deck.
Numerical simulations of the shallow water equations on rotating spheres with Jovian parameter values reproduce a mixture of robust vortices and alternating zonal jets. The latter arise naturally from the nonlinear inverse cascade characteristic of two-dimensional turbulence being arrested through coupling to Rossby waves at the Rhines scale (Rhines Reference Rhines1975) though the quasilinear ‘zonostrophic instability’ mechanism may also have a role (Srinivasan & Young Reference Srinivasan and Young2012, Reference Srinivasan and Young2014, and references therein). The equatorial jets are invariably sub-rotating (retrograde) in both freely decaying and forced-dissipative simulations in the Jovian regime (e.g. Cho & Polvani Reference Cho and Polvani1996a ,Reference Cho and Polvani b ; Iacono et al. Reference Iacono, Struglia and Ronchi1999a ,Reference Iacono, Struglia, Ronchi and Nicastro b ; Scott & Polvani Reference Scott and Polvani2007; Showman Reference Showman2007). By contrast, the mean zonal wind profiles in figure 1, as obtained from tracking visible features in Jupiter’s atmosphere, show a broad, super-rotating equatorial jet. The similarity of the two profiles, taken from the Voyager 2 (Limaye Reference Limaye1986) and Cassini (Porco et al. Reference Porco, West, McEwen, Del Genio, Ingersoll, Thomas, Squyres, Dones, Murray and Johnson2003) missions over 20 years apart, demonstrates the remarkable stability of Jupiter’s zonal winds, and the longevity of the super-rotating equatorial jet.
Figure 1. Observed mean zonal winds versus planetographic latitude for Jupiter. The blue dots show Voyager 2 data from Limaye (Reference Limaye1986) and the red line shows Cassini data from Porco et al. (Reference Porco, West, McEwen, Del Genio, Ingersoll, Thomas, Squyres, Dones, Murray and Johnson2003).
Forced-dissipative shallow water simulations typically include a Rayleigh friction term to absorb the slow inverse cascade of energy past the Rhines scale to the largest scales in the system, which tends to create large coherent vortices in place of the alternating zonal jets (e.g. Vasavada & Showman Reference Vasavada and Showman2005). Scott & Polvani (Reference Scott and Polvani2007, Reference Scott and Polvani2008) added an additional radiative relaxation term with time scale
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
to the continuity equation in their shallow water model:
$$\begin{eqnarray}\displaystyle & h_{t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h\boldsymbol{u})=-(h-h_{0})/\unicode[STIX]{x1D70F}_{\mathit{rad}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{u}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+f\hat{\boldsymbol{z}}\times \boldsymbol{u}=-g^{\prime }\unicode[STIX]{x1D735}h+\boldsymbol{F}-\boldsymbol{u}/\unicode[STIX]{x1D70F}_{\mathit{fric}}. & \displaystyle\end{eqnarray}$$
$h$
is the thickness of the active layer, commonly called the ‘height’,
$\boldsymbol{u}$
is the depth-averaged horizontal velocity,
$f=2\unicode[STIX]{x1D6FA}\sin \unicode[STIX]{x1D719}$
is the Coriolis parameter at latitude
$\unicode[STIX]{x1D719}$
on a planet rotating with angular velocity
$\unicode[STIX]{x1D6FA}$
,
$\unicode[STIX]{x1D735}$
is the horizontal gradient operator,
$\hat{\boldsymbol{z}}$
is a unit vector in the local vertical direction,
$\boldsymbol{F}$
is an isotropic random forcing and
$\unicode[STIX]{x1D70F}_{\mathit{fric}}$
is the time scale for Rayleigh friction. The reduced gravity is
$g^{\prime }=g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}$
, where
$g$
is the actual gravitational acceleration,
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$
is the density contrast between the cloud deck and the lower layer and
$\unicode[STIX]{x1D70C}_{0}$
is the density of the lower layer. These parameters are possibly constrained by observations of what appear to be internal gravity waves radiating from the impact points of Shoemaker–Levy comet debris (Dowling Reference Dowling1995; Ingersoll et al.
Reference Ingersoll, Dowling, Gierasch, Orton, Read, Sánchez-Lavega, Showman, Simon-Miller, Vasavada, Bagenal, Dowling and McKinnon2007), although this interpretation is not without its critics (Walterscheid, Brinkman & Schubert Reference Walterscheid, Brinkman and Schubert2000).
The right-hand side of (1.1a
) models the effects of radiation to space with a Newtonian relaxation of
$h$
towards its constant mean value
$h_{0}$
on the time scale
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
, following earlier shallow water models of the upper ocean mixed layer (Philander, Yamagata & Pacanowski Reference Philander, Yamagata and Pacanowski1984; Hirst Reference Hirst1986) and the terrestrial stratosphere (Juckes Reference Juckes1989; Polvani, Waugh & Plumb Reference Polvani, Waugh and Plumb1995; Thuburn & Lagneau Reference Thuburn and Lagneau1999). Scott & Polvani’s (Reference Scott and Polvani2007, Reference Scott and Polvani2008) simulations of this system with dissipation only by thickness relaxation (no
$\unicode[STIX]{x1D70F}_{\mathit{fric}}$
term) produced a sharply localised super-rotating equatorial jet, as do our subsequent simulations described in § 3, and in Warneford & Dellar (Reference Warneford and Dellar2014).
Thermal shallow water theory was introduced by Lavoie (Reference Lavoie1972) to describe atmospheric mixed layers over frozen lakes. It permits horizontal variations of the thermodynamic properties of the fluid within each layer, which are taken to be uniform in the standard shallow water equations. The density contrast
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$
thus becomes a dynamical variable. However, subsequent developments of thermal shallow water theory have been much more focused on modelling the oceanic mixed layer (Schopf & Cane Reference Schopf and Cane1983; Hirst Reference Hirst1986; McCreary & Kundu Reference McCreary and Kundu1988; McCreary, Fukamachi & Kundu Reference McCreary, Fukamachi and Kundu1991; McCreary & Yu Reference McCreary and Yu1992; Ripa Reference Ripa1993, Reference Ripa1995, Reference Ripa1996a
,Reference Ripa
b
; Røed & Shi Reference Røed and Shi1999).
Warneford & Dellar (Reference Warneford and Dellar2014) simulated Jovian atmospheres using the forced-dissipative thermal shallow water equations
$$\begin{eqnarray}\displaystyle & h_{t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6E9}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D6E9}=-(\unicode[STIX]{x1D6E9}h/h_{0}-\unicode[STIX]{x1D6E9}_{0})/\unicode[STIX]{x1D70F}_{\mathit{rad}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{u}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+f\hat{\boldsymbol{z}}\times \boldsymbol{u}=-h^{-1}\unicode[STIX]{x1D735}\left({\textstyle \frac{1}{2}}\unicode[STIX]{x1D6E9}h^{2}\right)+\boldsymbol{F}-\boldsymbol{u}/\unicode[STIX]{x1D70F}_{\mathit{fric}}. & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D6E9}=g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0}$
for the reduced gravity to emphasise that it is now a function of space and time. The first term on the right-hand side of (1.2c
) is then the usual shallow water pressure gradient term, as rewritten for a spatially varying reduced gravity. The right-hand side of (1.2b
) is a Newtonian cooling term that represents radiative relaxation towards a temperature
$\unicode[STIX]{x1D6E9}_{0}h_{0}/h$
with time scale
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
. We show below that this form of coupling between
$\unicode[STIX]{x1D6E9}$
and
$h$
changes the behaviour of equatorially trapped Rossby waves in almost exactly the same way as the radiative relaxation of
$h$
itself in Scott & Polvani’s (Reference Scott and Polvani2007, Reference Scott and Polvani2008) model. In particular, the damping rates calculated in § 4 are very similar function of wavenumber and
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
. The resulting Reynold stresses calculated in § 5 have very similar spatial structures between the two models, and are proportional to the difference between the dimensionless radiative and frictional decay rate constants in both models.
The thermal shallow water equations (1.2
a
–
c
) conserve both mass and momentum in the absence of the
$\unicode[STIX]{x1D70F}_{\mathit{fric}}$
and
$\boldsymbol{F}$
terms, unlike the Scott and Polvani model (1.1
a
,
b
). Simulations of the thermal shallow water equations for Jovian parameter values (reported in § 3 and Warneford & Dellar Reference Warneford and Dellar2014) reproduce a super-rotating equatorial jet, and more substantial mid-latitude jets than simulations of the Scott & Polvani (Reference Scott and Polvani2007, Reference Scott and Polvani2008) model. Henceforth, for brevity we refer to the radiative relaxation terms as relaxation and the Rayleigh friction term as friction. We also refer to (1.1
a
,
b
) as the ‘standard’ shallow water equations to distinguish them from the thermal shallow water equations (1.2
a
–
c
).
Quasi-geostrophic (QG) theory (Charney Reference Charney1949; Obukhov Reference Obukhov1949) offers a simplified description of slow vortical motions that filters out inertia–gravity waves. Although QG theory is usually employed for mid-latitude beta planes, Matsuno (Reference Matsuno1970, Reference Matsuno1971) presented a global QG theory for fully stratified atmospheres on spheres. Schubert, Taft & Silvers (Reference Schubert, Taft and Silvers2009) and Verkley (Reference Verkley2009) revived this theory for the standard shallow water equations on a sphere in the form
$$\begin{eqnarray}q_{t}+[\unicode[STIX]{x1D713},q]=0,\quad q=f+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}-L_{D}^{-2}\unicode[STIX]{x1D713},\end{eqnarray}$$
where the Jacobian
$[\unicode[STIX]{x1D713},q]=\hat{\boldsymbol{z}}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}q)$
. The Rossby deformation scale is
$L_{D}=c/f$
, where
$c$
is the speed of short shallow water gravity waves. Both
$f$
and
$L_{D}$
vary with latitude when QG theory is used in spherical coordinates. The sole evolving variable is the potential vorticity
$q$
, which is advected by the velocity field derived from a streamfunction
$\unicode[STIX]{x1D713}$
. The elliptic equation relating
$\unicode[STIX]{x1D713}$
to
$q$
involves the spatially varying Coriolis parameter
$f$
. It corresponds to Daley’s (Reference Daley1983) simplest form of geostrophic balance, as justified by assuming that
$\unicode[STIX]{x1D713}$
varies on length scales much smaller than the planetary scales on which
$f$
varies. The system (1.3) reduces to the familiar beta-plane QG equations on approximating
$f$
by
$f_{0}+\unicode[STIX]{x1D6FD}y$
in the first term, and by
$f_{0}$
in the second term. Warneford (Reference Warneford2014) derived a thermal form of this global QG theory, building on the beta-plane thermal QG equations in Ripa (Reference Ripa1996a
) and Warneford & Dellar (Reference Warneford and Dellar2013). Numerical simulations produced super-rotating equatorial jets when the dimensionless radiative relaxation time
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
is sufficiently short, as did simulations of the corresponding global QG form of Scott & Polvani’s (Reference Scott and Polvani2007, Reference Scott and Polvani2008) shallow water model with a relaxing thickness. By contrast, both QG models produced sub-rotating equatorial jets when
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
is sufficiently long. The mechanism responsible for setting the equatorial jet direction must therefore lie within QG theory. This focuses attention on momentum fluxes due to Rossby waves, rather than by inertia–gravity waves, or the equatorially trapped Kelvin and Yanai waves. None of these waves exist in QG theory.
2 Momentum transport by Rossby waves
Dickinson (Reference Dickinson1969) showed that the meridional flux of zonal momentum due to Rossby waves produces a net eastward acceleration in regions where Rossby waves are generated, and a net westward acceleration in regions where Rossby waves are dissipated. This argument was further developed by Green (Reference Green1970) and Thompson (Reference Thompson1971). A detailed version based on ray theory may be found in Held (Reference Held2000) and Vallis (Reference Vallis2006). The QG equations on a mid-latitude beta plane lead to the Rossby wave dispersion relation
$\unicode[STIX]{x1D714}=-\unicode[STIX]{x1D6FD}k/(k^{2}+\ell ^{2}+L_{D}^{2})$
for plane-wave disturbances with streamfunction
$\unicode[STIX]{x1D713}=\hat{\unicode[STIX]{x1D713}}\exp (\text{i}(kx+\ell y-\unicode[STIX]{x1D714}t))$
implies a meridional group velocity
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}\ell }=\frac{2\unicode[STIX]{x1D6FD}k\ell }{(k^{2}+\ell ^{2}+L_{D}^{2})^{2}},\end{eqnarray}$$
while the off-diagonal component of the Reynolds stress (see § 4 for notation) is
$$\begin{eqnarray}\langle u^{\prime }v^{\prime }\rangle =-{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D713}}^{2}k\ell .\end{eqnarray}$$
The radiation condition requires the meridional group velocity to point away from wave sources, so the sign of
$k\ell$
is such as to create a flux of westward (negative) zonal momentum away from sources of Rossby waves towards regions where Rossby waves are dissipated. For example, Schneider & Liu (Reference Schneider and Liu2009) and Liu & Schneider (Reference Liu and Schneider2010) identified Rossby wave sources due to divergent motions associated with the breakdown of geostrophy near the equator in their three-dimensional simulations of Jovian atmospheres. Propagation of these Rossby waves to higher latitudes thus creates a net eastward acceleration at the equator according to the above theory, which is consistent with the formation of super-rotating equatorial jets in their simulations.
However, the focus of our present paper will be on the equatorially trapped Rossby waves that play a key role in near-equatorial dynamics (e.g. Matsuno Reference Matsuno1966; Gill Reference Gill1982; McCreary Reference McCreary1985; Majda & Klein Reference Majda and Klein2003; Khouider, Majda & Stechmann Reference Khouider, Majda and Stechmann2013). These waves are confined by a Gaussian envelope whose extent scales with the equatorial deformation scale
$L_{\mathit{eq}}=\sqrt{c/(2\unicode[STIX]{x1D6FD})}$
. This corresponds to latitudes within about
$30^{\circ }$
of the equator for Jovian parameters (see § 4). We therefore need to compute the latitude dependence of the momentum fluxes inside the scale of this envelope, for which a ray theory approach is inadequate.
The generalised Lagrangian mean (GLM) theory is a powerful approach for addressing this and other problems in wave–mean flow interaction (Andrews & McIntyre Reference Andrews and McIntyre1976a
; Bühler Reference Bühler2000, Reference Bühler2014). It defines the Lagrangian mean of an arbitrary field
$\unicode[STIX]{x1D719}(\boldsymbol{x},t)$
by
$$\begin{eqnarray}\unicode[STIX]{x1D719}^{L}(\boldsymbol{x},t)=\overline{\unicode[STIX]{x1D719}(\boldsymbol{x}+\unicode[STIX]{x1D743}(\boldsymbol{x},t),t)},\end{eqnarray}$$
where
$\overline{(\cdots )}$
is some linear averaging operator and
$\unicode[STIX]{x1D743}(\boldsymbol{x},t)$
is a Lagrangian displacement field such that
$\boldsymbol{x}+\unicode[STIX]{x1D743}(\boldsymbol{x},t)$
is the actual position of the particle whose mean position is
$\boldsymbol{x}$
at time
$t$
. GLM theory establishes circumstances for the validity of a ‘pseudomomentum rule’ (McIntyre Reference McIntyre1981; Bühler Reference Bühler2000, Reference Bühler2014) under which the waves behave for some purposes as if they carried a momentum equal to their pseudomomentum and the background fluid were absent. The pseudomomentum is a disturbance quantity related to the pseudoenergy defined in § 5.
A linearised version of GLM theory for small
$\unicode[STIX]{x1D743}$
was presented in Andrews & McIntyre (Reference Andrews and McIntyre1976a
, Reference Andrews and McIntyre1978). However, the Lagrangian picture of a fluid evolving through displacements of indestructible particles fits awkwardly with the mass source or sink term on the right-hand side of (1.1a
) that models radiative relaxation. In particular, lack of mass conservation requires modification of Bühler’s (Reference Bühler2000) shallow water pseudomomentum rule obtained from GLM theory. Nevertheless, McIntyre (2014, personal communication) has established, firstly, that is possible to reproduce our results below using the linearised GLM theory slightly adapted to apply to (1.1
a
,
b
) and, secondly, that an eastward zonal acceleration at the equator is obtained for two interesting cases, (i) equatorially trapped Rossby waves freely decaying by radiative relaxation – as we compute independently below – and (ii) the same waves, but maintained at constant amplitude by a balance between forcing and dissipation, with the dissipation again by radiative relaxation but with the forcing mimicked by negative Rayleigh friction. The linearised GLM theory shows why negative Rayleigh friction, if spatially uniform, has an effect qualitatively similar to that of free temporal decay – as illustrated by the curves in figure 4(b) of Andrews & McIntyre (Reference Andrews and McIntyre1976a
) with suitable sign changes.
Following the approach of McIntyre & Norton (Reference McIntyre and Norton1990) for stratified fluids, Bühler (Reference Bühler2000) formulated a general approach for shallow water equations of the form
$$\begin{eqnarray}\displaystyle & h_{t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{u}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+f\hat{\boldsymbol{z}}\times \boldsymbol{u}=-g^{\prime }\unicode[STIX]{x1D735}h+\boldsymbol{F}, & \displaystyle\end{eqnarray}$$
$\boldsymbol{F}$
, not necessarily a forcing. The shallow water potential vorticity
$q=(f+\unicode[STIX]{x1D701})/h$
, where
$\unicode[STIX]{x1D701}=\hat{\boldsymbol{z}}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\times \boldsymbol{u})$
is the relative vorticity, evolves according to 
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(hq)+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(hq\boldsymbol{u}-\boldsymbol{F}^{\bot })=0,\end{eqnarray}$$
with an extra flux
$\boldsymbol{F}^{\bot }=\hat{\boldsymbol{z}}\times \boldsymbol{F}$
. Bühler’s (Reference Bühler2000) pseudomomentum rule holds for systems of this form that conserve linear momentum, i.e. those for which
$h\boldsymbol{F}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}$
is the divergence of a stress tensor
$\unicode[STIX]{x1D748}$
.
However, neither of the shallow water models from § 1 fits this form. The shallow water system with relaxation of the thickness does not conserve momentum. In the absence of forcing and friction equations (1.1a ) and (1.1b ) together imply
$$\begin{eqnarray}(h\boldsymbol{u})_{t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(h\boldsymbol{u}\boldsymbol{u}+{\textstyle \frac{1}{2}}g^{\prime }h^{2}\unicode[STIX]{x1D644}\right)+f\hat{\boldsymbol{z}}\times (h\boldsymbol{u})=-(h-h_{0})\boldsymbol{u}/\unicode[STIX]{x1D70F}_{\mathit{rad}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D644}$
is the
$2\times 2$
identity matrix. Moreover, the potential vorticity for the system (1.1
a
,
b
) obeys
$$\begin{eqnarray}(hq)_{t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(hq\boldsymbol{u}+\boldsymbol{u}^{\bot }/\unicode[STIX]{x1D70F}_{\mathit{fric}})=0,\end{eqnarray}$$
where
$\boldsymbol{u}^{\bot }=\hat{\boldsymbol{z}}\times \boldsymbol{u}$
, while the potential vorticity in the thermal shallow water equations (1.2
a
–
c
) obeys
$$\begin{eqnarray}(hq)_{t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(hq\boldsymbol{u}-{\textstyle \frac{1}{2}}h\unicode[STIX]{x1D735}^{\bot }\unicode[STIX]{x1D703}+\boldsymbol{u}^{\bot }/\unicode[STIX]{x1D70F}_{\mathit{fric}}\right)=0.\end{eqnarray}$$
Radiative relaxation thus does not contribute to either potential vorticity equation. However, equation (2.8) contains an extra term involving the perpendicular temperature gradient
$\unicode[STIX]{x1D735}^{\bot }\unicode[STIX]{x1D703}=\hat{\boldsymbol{z}}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
that arises from the curl of the
$(1/2)h\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E9}$
term on the right-hand side of (1.2c
). The thermal shallow water equations thus do not materially conserve potential vorticity, even in the absence of forcing and dissipation. Instead, the total potential vorticity inside a closed isotherm is conserved. This is the thermal shallow water analogue of the potential vorticity impermeability property of isentropic surfaces in three-dimensional stratified fluids (Haynes & McIntyre Reference Haynes and McIntyre1990).
Since the Scott & Polvani (Reference Scott and Polvani2007, Reference Scott and Polvani2008) radiatively relaxed shallow water model lacks both mass and momentum conservation, and hence does not fit naturally into the above general frameworks, in §§ 4 and 5 below we derive from first principles equations for the evolution of the thickness-weighted zonal mean velocity components
$\langle u\rangle ^{\ast }$
and
$\langle v\rangle ^{\ast }$
(see § 5 for notation) caused by the decay of an equatorially trapped Rossby wave by friction and/or radiative relaxation.
3 Numerical experiments
We now show some results from numerical simulations of three different models: the standard shallow water equations with dissipation only through friction (i.e. (1.1
a
,
b
) with no
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
term), the Scott & Polvani (Reference Scott and Polvani2007, Reference Scott and Polvani2008) model with relaxation of the thickness, and our thermal shallow water model. Table 1 gives typical Jovian values for the planetary radius
$a$
, angular velocity
$\unicode[STIX]{x1D6FA}$
, gravitational acceleration
$g$
and internal wave speed
$\sqrt{g^{\prime }h_{0}}$
. The last is deduced from impacts of Shoemaker–Levy comet debris (Dowling Reference Dowling1995; Ingersoll et al.
Reference Ingersoll, Dowling, Gierasch, Orton, Read, Sánchez-Lavega, Showman, Simon-Miller, Vasavada, Bagenal, Dowling and McKinnon2007), and implies that there are 232 polar deformation scales
$L_{D}=\sqrt{g^{\prime }h_{0}}/(2\unicode[STIX]{x1D6FA})$
around the circumference. We used a radiative relaxation time scale
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
of 43 Jovian days, as calculated for a temperature of 120 K at the 25 millibar pressure level, and all simulations employed a friction with a time scale of 1000 Jovian days. Liu & Schneider (Reference Liu and Schneider2010) used the 25 millibar pressure level when comparing their three-dimensional simulations with observations of Neptune. Calculating
$\unicode[STIX]{x1D70F}_{\mathit{rad}}$
at the 25 millibar pressure level for both planets leads to a super-rotating equatorial jet for Jupiter, and a sub-rotating equatorial jet for Neptune, in our thermal shallow water model (Warneford & Dellar Reference Warneford and Dellar2014). The theory in § 5 below suggests that the direction of the equatorial jet is insensitive to these parameters provided the radiative relaxation time scale remains much shorter than the frictional time scale.
Table 1. Parameter values for Jupiter (from Ingersoll Reference Ingersoll1990; Beebe Reference Beebe1994; Warneford & Dellar Reference Warneford and Dellar2014).
Simulations of geostrophic turbulence in the Jovian regime require thousands of rotation periods to reach statistically steady states. The superior memory bandwidth and peak floating point performance of graphical processing units (GPUs) thus becomes attractive, but existing spherical harmonic transform algorithms for GPUs only achieve performance parity with conventional microprocessors (Hupca et al.
Reference Hupca, Falcou, Grigori and Stompor2012). We therefore employed the doubly periodic Cartesian domain sketched in figure 2, which we refer to as the square planet domain. The horizontal axis denotes longitude, while the vertical axis denotes latitude. Starting from the north pole at the top of the domain and moving down, we reach the equator a quarter of the way down, and the south pole half-way down. We then continue to reach the equator again, before returning to the north pole at the bottom of the domain. We imagine following a meridian (constant longitude line) from the north pole to the equator to the south pole, and then following a second meridian with longitude offset by
$180^{\circ }$
back across the equator to the north pole. The Coriolis force is still fully varying in latitude, and is defined in the simulations by
$f(y)=2\unicode[STIX]{x1D6FA}\cos (2\unicode[STIX]{x03C0}y/y_{max})$
, where
$y_{max}$
is the length of the side of the square planet domain, i.e. the circumference of Jupiter (see table 1). The Coriolis force and all its derivatives are thus doubly periodic. This allowed us to exploit a high performance GPU fast Fourier transform library. We discretised the domain using
$1024\times 1024$
Fourier collocation points, and used the Hou & Li (Reference Hou and Li2007) spectral filter to control the buildup of enstrophy at the highest wavenumbers. We integrated the resulting large system of ordinary differential equations using the standard fourth-order Runge–Kutta scheme, with a time step determined dynamically from the Courant–Friedrichs–Lewy stability condition.
Figure 2. Schematic of the doubly periodic Cartesian domain used for the numerical simulations. Each half of the domain represents a full planet, so for ease of interpretation we show simulation results for just the top planet.
Our initial conditions were
$h=h_{0}$
and
$u=v=0$
for all runs, and
$\unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D6E9}_{0}$
for the thermal shallow water simulation. Following Scott & Polvani (Reference Scott and Polvani2007, Reference Scott and Polvani2008) we applied a divergence-free isotropic random forcing
$\boldsymbol{F}$
that is localised to a narrow annulus of wavenumbers
$|\boldsymbol{k}|\in [40,44]$
in Fourier space, counting the longest axis-aligned sinusoidal mode in the domain as wavenumber 1. We forced each mode inside this annulus with amplitude
$\unicode[STIX]{x1D716}_{f}$
using random phases that were
$\unicode[STIX]{x1D6FF}$
-correlated (white) in time. This forcing is widely used in numerical studies of zonal jet formation, and models the energy injected by three-dimensional convection at horizontal length scales comparable to the deformation radius. Our simulations slowly adjusted the amplitude of the forcing
$\unicode[STIX]{x1D716}_{f}$
to reach a prescribed value of the total kinetic energy in the eventual statistical steady state. Further details of the numerical models, parameters, and simulation outputs may be found in Warneford (Reference Warneford2014) and Warneford & Dellar (Reference Warneford and Dellar2014).
Figure 3. Absolute vorticity
$\unicode[STIX]{x1D714}_{a}=v_{x}-u_{y}+f(y)$
in units of
$10^{-4}~\text{s}^{-1}$
for (a) the standard shallow water equations with dissipation by friction alone, (b) the standard shallow water equations with friction and radiative relaxation of the thickness and (c) the thermal shallow water equations with friction and radiative relaxation of the temperature.
Figure 4. Zonal velocity
$u$
in units of
$\text{m}~\text{s}^{-1}$
for (a) the standard shallow water equations with dissipation by friction alone, (b) the standard shallow water equations with friction and radiative relaxation of the thickness and (c) the thermal shallow water equations with friction and radiative relaxation of the temperature.
Figure 5. Mean zonal velocity
$\langle u\rangle$
against latitude for simulations of the standard shallow water equations with dissipation by friction alone, the standard shallow water equations with friction and radiative relaxation of the thickness, and the thermal shallow water equations with friction and radiative relaxation of the temperature.
Figure 3 shows the instantaneous absolute vorticity
$\unicode[STIX]{x1D714}_{a}=v_{x}-u_{y}+f(y)$
after
$2\times 10^{4}$
Jovian days for the top planet in figure 2 for three different simulations: the standard shallow water equations with friction alone, the standard shallow water equations with friction and radiative relaxation of the thickness and the thermal shallow water equations with friction and radiative relaxation of the temperature. Figure 4 shows the corresponding instantaneous zonal velocity and figure 5 shows the instantaneous zonally averaged zonal velocity,
$\langle u\rangle =y_{max}^{-1}\int u(x,y,t)\,\text{d}x$
. All three simulations produced a mixture of vortices, turbulence and multiple zonal jets, with amplitudes that decrease at higher latitudes. The two simulations with radiative relaxation produced a strong super-rotating equatorial jet in line with observations of Jupiter. However, our simulation of the standard shallow water model, with dissipation only by friction, produced a sub-rotating equatorial jet. The simulation with radiative relaxation of the thickness produced very weak jets away from the equator, whereas the other two simulations produced stronger mid-latitude jets in better agreement with observations.
4 Rossby waves on an equatorial beta plane
We now investigate the properties of equatorially trapped waves on an equatorial beta plane in the different shallow water models, as they decay freely due to Rayleigh friction and/or Newtonian radiative relaxation. The equatorial beta plane uses local Cartesian coordinates with
$x$
eastward and
$y$
northward. The underlying spherical geometry appears only through the linearisation of
$f=2\unicode[STIX]{x1D6FA}\sin \unicode[STIX]{x1D719}$
for small latitudes
$|\unicode[STIX]{x1D719}|\ll 1$
into
$f(y)=\unicode[STIX]{x1D6FD}y$
with
$\unicode[STIX]{x1D6FD}=2\unicode[STIX]{x1D6FA}/a$
, where
$a$
is the planetary radius.
4.1 Shallow water equations with thickness relaxation
The linearised form of the unforced (
$\boldsymbol{F}=\mathbf{0}$
) shallow water equations (1.1) with radiative relaxation of the thickness on an equatorial beta plane may be written as (Matsuno Reference Matsuno1966; Gill Reference Gill1982)
$$\begin{eqnarray}\displaystyle & h_{t}^{\prime }+u_{x}^{\prime }+v_{y}^{\prime }=-\unicode[STIX]{x1D705}h^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & u_{t}^{\prime }-{\textstyle \frac{1}{2}}yv^{\prime }+h_{x}^{\prime }=-\unicode[STIX]{x1D6FE}u^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & v_{t}^{\prime }+{\textstyle \frac{1}{2}}yu^{\prime }+h_{y}^{\prime }=-\unicode[STIX]{x1D6FE}v^{\prime }. & \displaystyle\end{eqnarray}$$
$h_{0}$
. We have non-dimensionalised using the internal wave speed
$c=\sqrt{g^{\prime }h_{0}}$
as the velocity scale, the equatorial deformation radius
$L_{eq}=\sqrt{c/2\unicode[STIX]{x1D6FD}}$
as the horizontal length scale, and
$h_{0}$
as the thickness scale. Using the Jovian parameters in table 1, a latitude of
$33^{\circ }$
corresponds to
$y=5$
in the subsequent plots. The dimensionless radiative relaxation and friction rates are
$\unicode[STIX]{x1D705}=T/\unicode[STIX]{x1D70F}_{\mathit{rad}}$
and
$\unicode[STIX]{x1D6FE}=T/\unicode[STIX]{x1D70F}_{\mathit{fric}}$
based on the advective time scale
$T=L_{eq}/c$
. Their values were
$\unicode[STIX]{x1D705}=7.88\times 10^{-3}$
and
$\unicode[STIX]{x1D6FE}=T/\unicode[STIX]{x1D70F}_{\mathit{fric}}=3.42\times 10^{-4}$
in our simulation with radiative relaxation of the thickness.
Following standard practice (e.g. Matsuno Reference Matsuno1966; Gill Reference Gill1982) we seek waves that are harmonic in longitude and time of the form
$h^{\prime }(x,y,t)=\text{Re}\{{\hat{h}}(y)\exp (\text{i}(kx-\unicode[STIX]{x1D714}t))\}$
etc. With no dissipation (
$\unicode[STIX]{x1D705}=\unicode[STIX]{x1D6FE}=0$
) the three equations (4.1) may be combined into a single ordinary differential equation (ODE) for the meridional velocity
$\hat{v}(y)$
,
$$\begin{eqnarray}\frac{\text{d}^{2}\hat{v}}{\text{d}y^{2}}=(Ay^{2}-B)\hat{v},\quad A=\frac{1}{4},B=\unicode[STIX]{x1D714}^{2}-k^{2}-\frac{k}{2\unicode[STIX]{x1D714}},\end{eqnarray}$$
the same equation that governs a quantum harmonic oscillator. The solutions that decay as
$y\rightarrow \pm \infty$
may be written using the Hermite polynomials
$H_{n}(\unicode[STIX]{x1D709})$
as
$$\begin{eqnarray}\hat{v}=H_{n}(\unicode[STIX]{x1D709})\exp (-\unicode[STIX]{x1D709}^{2}/2),\quad \unicode[STIX]{x1D709}=yA^{1/4}.\end{eqnarray}$$
The dispersion relation
$A^{-1/2}B=2n+1$
for
$n=-1,0,1,2,\ldots$
gives a cubic equation for
$\unicode[STIX]{x1D714}$
,
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}-k^{2}-\frac{k}{2\unicode[STIX]{x1D714}}=n+\frac{1}{2}.\end{eqnarray}$$
The Rossby wave branch of solutions is characterised by
$n\geqslant 1$
and
$0<-\unicode[STIX]{x1D714}/k\ll 1$
. The other two roots give inertia–gravity waves with
$|\unicode[STIX]{x1D714}|\geqslant 1$
(the dimensionless inertial frequency). The cases
$n=0$
and
$n=-1$
give the Yanai and Kelvin waves respectively. Figure 6 shows these different branches of the dispersion relation, all of which represent trapped waves localised within a few equatorial deformation scales
$L_{\mathit{eq}}$
of the equator. The inertia–gravity, Kelvin and Yanai wave branches all exceed the inertial frequency
$|\unicode[STIX]{x1D714}|=1$
. The inertia–gravity and Kelvin waves thus do not exist in QG theory, while the frequency of the QG form of the Yanai wave approaches 0 as
$k\rightarrow 0$
like a Rossby wave. The mechanism responsible for setting the direction of equatorial jet in the shallow water models is expected to lie within QG theory, so we only consider the Rossby waves in the remainder of this paper.
The friction and radiative relaxation terms in (4.1
a
–
c
) change the
$A$
and
$B$
coefficients to (Yamagata & Philander Reference Yamagata and Philander1985)
$$\begin{eqnarray}A=\frac{1}{4}\left(\frac{\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE}}\right),\quad B=\unicode[STIX]{x1D714}^{2}-k^{2}-\frac{k}{2\left(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE}\right)}+\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D6FE})-\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FE},\end{eqnarray}$$
so the dispersion relation becomes
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}-k^{2}-\frac{k}{2(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})}+\text{i}\unicode[STIX]{x1D714}(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D6FE})-\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FE}=\left(n+\frac{1}{2}\right)\left(\frac{\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE}}\right)^{1/2}\quad \text{for}~n=1,2,\ldots .\end{eqnarray}$$
Figure 6. Dimensionless dispersion relation (4.4) for equatorially trapped waves in the standard shallow water equations on the equatorial beta plane with no dissipation.
Figure 7. The (a) real and (b) imaginary parts of the dimensionless equatorial Rossby wave frequency
$\unicode[STIX]{x1D714}$
for the standard shallow water equations with no dissipation (green solid lines), dissipation by relaxation of the thickness (red dash-dot lines) and dissipation by friction (blue dotted lines). The three lines appear coincident in (a).
Figure 7 shows the real and imaginary parts of the complex frequency
$\unicode[STIX]{x1D714}$
for the first few equatorial Rossby waves (
$n=1,2,3$
) for different values of
$\unicode[STIX]{x1D6FE}$
and
$\unicode[STIX]{x1D705}$
. The dispersion relation (4.4) is invariant under the transformation
$(k,\unicode[STIX]{x1D714})\mapsto (-k,-\unicode[STIX]{x1D714})$
, so we choose to take
$k>0$
. The Rossby wave frequency
$\text{Re}(\unicode[STIX]{x1D714})$
is now negative, and virtually unchanged by either weak friction (
$\unicode[STIX]{x1D705}=0.001$
,
$\unicode[STIX]{x1D6FE}=0$
) or weak radiative relaxation (
$\unicode[STIX]{x1D705}=0$
,
$\unicode[STIX]{x1D6FE}=0.001$
). The decay rate
$-\text{Im}(\unicode[STIX]{x1D714})$
due to cooling (
$\unicode[STIX]{x1D705}=0.001$
,
$\unicode[STIX]{x1D6FE}=0$
) is largest at
$k=0$
and tends to zero as
$k\rightarrow \infty$
. By contrast, the decay due to friction (
$\unicode[STIX]{x1D6FE}=0.001$
,
$\unicode[STIX]{x1D705}=0$
) is monotonically increasing with
$k$
and
$-\text{Im}(\unicode[STIX]{x1D714})$
tends to
$-\unicode[STIX]{x1D6FE}$
as
$k\rightarrow \infty$
. Figure 7(b) also suggests a symmetry about the line
$\text{Im}(\unicode[STIX]{x1D714})=-\unicode[STIX]{x1D716}/2$
between small radiative relaxation (
$\unicode[STIX]{x1D705}=\unicode[STIX]{x1D716}$
and
$\unicode[STIX]{x1D6FE}=0$
) and small friction (
$\unicode[STIX]{x1D705}=0$
and
$\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D716}$
). We expand the roots of (4.6) as
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{0}+\unicode[STIX]{x1D716}\,\unicode[STIX]{x1D714}_{1,\mathit{rad}}+\unicode[STIX]{x1D716}^{2}\,\unicode[STIX]{x1D714}_{2,\mathit{rad}}+\cdots$
for radiative relaxation, and
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{0}+\unicode[STIX]{x1D716}\,\unicode[STIX]{x1D714}_{1,\mathit{fric}}+\unicode[STIX]{x1D716}^{2}\,\unicode[STIX]{x1D714}_{2,\mathit{fric}}+\cdots$
for friction, where
$\unicode[STIX]{x1D714}_{0}$
is the real root of (4.4). The
$O(\unicode[STIX]{x1D716})$
corrections are both purely imaginary:
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D714}_{1,\mathit{rad}}={\textstyle \frac{1}{2}}\text{i}\unicode[STIX]{x1D714}_{0}(-4\unicode[STIX]{x1D714}_{0}^{2}+2n+1)/(k+4\unicode[STIX]{x1D714}_{0}^{3}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D714}_{1,\mathit{fric}}={\textstyle \frac{1}{2}}\text{i}\unicode[STIX]{x1D714}_{0}(-4\unicode[STIX]{x1D714}_{0}^{2}-2n-1-2k/\unicode[STIX]{x1D714}_{0})/(k+4\unicode[STIX]{x1D714}_{0}^{3}). & \displaystyle\end{eqnarray}$$
$\text{Im}(\unicode[STIX]{x1D714})=-\unicode[STIX]{x1D716}/2$
visible in figure 7(b) is a consequence of the relation 
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{1,\mathit{rad}}+\frac{\text{i}}{2}=\frac{\text{i}}{2}\left(\frac{k+(2n+1)\unicode[STIX]{x1D714}_{0}}{k+4\unicode[STIX]{x1D714}_{0}^{3}}\right)=-\!\left(\unicode[STIX]{x1D714}_{1,\mathit{fric}}+\frac{\text{i}}{2}\right)\end{eqnarray}$$
between the two dispersion relations truncated at
$O(\unicode[STIX]{x1D716})$
.
4.2 Thermal shallow water equations
Linearising the thermal shallow water equations (1.2
a
–
c
) about a rest state with uniform thickness
$h_{0}$
and temperature
$\unicode[STIX]{x1D6E9}_{0}$
, and applying the equivalent non-dimensionalisation to § 4.1 based on the velocity scale
$\sqrt{h_{0}\unicode[STIX]{x1D6E9}_{0}}$
gives
$$\begin{eqnarray}\displaystyle & h_{t}^{\prime }+u_{x}^{\prime }+v_{y}^{\prime }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6E9}_{t}^{\prime }=-\unicode[STIX]{x1D705}(h^{\prime }+\unicode[STIX]{x1D6E9}^{\prime }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & u_{t}^{\prime }-{\textstyle \frac{1}{2}}yv^{\prime }+h_{x}^{\prime }+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6E9}_{x}^{\prime }=-\unicode[STIX]{x1D6FE}u^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & v_{t}^{\prime }+{\textstyle \frac{1}{2}}yu^{\prime }+h_{y}^{\prime }+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6E9}_{y}^{\prime }=-\unicode[STIX]{x1D6FE}v^{\prime }. & \displaystyle\end{eqnarray}$$
$\hat{v}(y)$
, with coefficients 
$$\begin{eqnarray}A=\frac{1}{4}\frac{\unicode[STIX]{x1D714}(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})}{(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705}/2)},\quad B=\frac{\unicode[STIX]{x1D714}(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})}{\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705}/2}-k^{2}-\frac{k}{2(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})},\end{eqnarray}$$
so the dispersion relation for equatorially trapped waves is
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D714}(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})}{\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705}/2}-k^{2}-\frac{k}{2(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})}=\left(n+\frac{1}{2}\right)\left(\frac{\unicode[STIX]{x1D714}(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})}{(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705}/2)}\right)^{1/2}.\end{eqnarray}$$
Figure 8 shows the real and imaginary parts of the complex frequency
$\unicode[STIX]{x1D714}$
for the first few equatorial Rossby waves with
$n\in \{1,2,3\}$
, and three different sets of
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D6FE}$
values. As before, the real part of
$\unicode[STIX]{x1D714}$
is virtually unchanged by weak dissipation (
$\unicode[STIX]{x1D705}=0.002$
or
$\unicode[STIX]{x1D6FE}=0.001$
) while
$\text{Im}(\unicode[STIX]{x1D714})$
becomes negative. The plot of
$\text{Im}(\unicode[STIX]{x1D714})$
in figure 8 closely resembles our earlier figure 7, with the same apparent symmetry, if we now take
$\unicode[STIX]{x1D705}=2\unicode[STIX]{x1D6FE}$
. A similar perturbation analysis to before shows that including this factor of 2 indeed makes
$\text{Im}(\unicode[STIX]{x1D714})$
symmetric about
$-\unicode[STIX]{x1D716}/2$
between the cases of small radiative relaxation (
$\unicode[STIX]{x1D705}=2\unicode[STIX]{x1D716}$
and
$\unicode[STIX]{x1D6FE}=0$
) and small friction (
$\unicode[STIX]{x1D705}=0$
and
$\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D716}$
).
Figure 8. The real and imaginary parts of the dimensionless equatorial Rossby wave frequency
$\unicode[STIX]{x1D714}$
for the thermal shallow water equations with no dissipation (green solid lines), dissipation by relaxation of the temperature (red dash-dot lines) and dissipation by friction (blue dotted lines). The three lines appear coincident in (a).
5 Acceleration of the zonal mean zonal flow
We now investigate the acceleration of the zonal mean zonal flow created by these solutions for decaying equatorially trapped Rossby waves.
5.1 Shallow water equations with thickness relaxation
We consider the shallow water model (1.1) with radiative relaxation of the thickness, and no forcing, on an equatorial beta plane, and apply the non-dimensionalisation from § 4.1 to obtain
$$\begin{eqnarray}\displaystyle & h_{t}+(hu)_{x}+(hv)_{y}=-\unicode[STIX]{x1D705}(h-1), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & (hu)_{t}+\left(hu^{2}+{\textstyle \frac{1}{2}}h^{2}\right)_{x}+(huv)_{y}-{\textstyle \frac{1}{2}}yhv=-\unicode[STIX]{x1D6FE}hu-\unicode[STIX]{x1D705}u(h-1), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & (hv)_{t}+(huv)_{x}+\left(hv^{2}+{\textstyle \frac{1}{2}}h^{2}\right)_{y}+{\textstyle \frac{1}{2}}yhu=-\unicode[STIX]{x1D6FE}hv-\unicode[STIX]{x1D705}v(h-1). & \displaystyle\end{eqnarray}$$
$h_{0}$
becomes unity under this non-dimensionalisation.
Following standard practice, we decompose the velocity and other fields as
$\boldsymbol{u}=\langle \boldsymbol{u}\rangle +\boldsymbol{u}^{\prime }$
into a zonal mean
$\langle \boldsymbol{u}\rangle$
and a deviation
$\boldsymbol{u}^{\prime }$
. The zonal mean of the dimensionless mass conservation equation may then be written as
$$\begin{eqnarray}\langle h\rangle _{t}+\langle hv\rangle _{y}=-\unicode[STIX]{x1D705}(\langle h\rangle -1).\end{eqnarray}$$
Similarly, the zonal means of the two momentum equations are
$$\begin{eqnarray}\langle hu\rangle _{t}+\langle huv\rangle _{y}-{\textstyle \frac{1}{2}}y\langle hv\rangle =-\unicode[STIX]{x1D6FE}\langle hu\rangle -\unicode[STIX]{x1D705}(\langle hu\rangle -\langle u\rangle ),\end{eqnarray}$$
and
$$\begin{eqnarray}\langle hv\rangle _{t}+\langle hv^{2}+{\textstyle \frac{1}{2}}h^{2}\rangle _{y}+{\textstyle \frac{1}{2}}y\langle hu\rangle =-\unicode[STIX]{x1D6FE}\langle hv\rangle -\unicode[STIX]{x1D705}(\langle hv\rangle -\langle v\rangle ).\end{eqnarray}$$
We now consider the form of solutions to the system (5.1a
–
c
) for initial conditions corresponding to a Rossby wave of small amplitude
$O(a)$
superimposed on a rest state with unit thickness. For these initial conditions, and under the subsequent evolution, the thickness and velocity components satisfy the scalings
$$\begin{eqnarray}h=1+\underbrace{h^{\prime }}_{O(a)}+\underbrace{(\langle h\rangle -1)}_{O(a^{2})},\quad u=\underbrace{u^{\prime }}_{O(a)}+\underbrace{\langle u\rangle }_{O(a^{2})},\quad v=\underbrace{v^{\prime }}_{O(a)}+\underbrace{\langle v\rangle }_{O(a^{2})}.\end{eqnarray}$$
We define
$\langle \unicode[STIX]{x1D702}\rangle ^{\ast }=\langle h\unicode[STIX]{x1D702}\rangle /\langle h\rangle$
to be the thickness-weighted zonal mean of a quantity
$\unicode[STIX]{x1D702}$
, as in Thuburn & Lagneau (Reference Thuburn and Lagneau1999). Using
$\langle h\rangle =1+O(a^{2})$
, we can simplify the above to obtain evolution equations correct to
$O(a^{2})$
in the form
$$\begin{eqnarray}\displaystyle & \langle h\rangle _{t}+\langle v\rangle _{y}^{\ast }+\unicode[STIX]{x1D705}(\langle h\rangle -1)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \langle u\rangle _{t}^{\ast }-{\textstyle \frac{1}{2}}y\langle v\rangle ^{\ast }+\unicode[STIX]{x1D6FE}\langle u\rangle ^{\ast }=-\langle u^{\prime }v^{\prime }\rangle _{y}-\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle , & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \langle v\rangle _{t}^{\ast }+{\textstyle \frac{1}{2}}y\langle u\rangle ^{\ast }+\unicode[STIX]{x1D6FE}\langle v\rangle ^{\ast }+\langle h\rangle _{y}=-\langle v^{\prime }v^{\prime }\rangle _{y}-\unicode[STIX]{x1D705}\langle h^{\prime }v^{\prime }\rangle -{\textstyle \frac{1}{2}}\langle h^{\prime 2}\rangle _{y}. & \displaystyle\end{eqnarray}$$
$\langle h\rangle$
,
$\langle u\rangle ^{\ast }$
,
$\langle v\rangle ^{\ast }$
. On the right-hand sides we have only means of quadratic products of the fluctuations
$h^{\prime }$
,
$u^{\prime }$
,
$v^{\prime }$
. These include the meridional derivatives of the Reynolds stress components
$\langle u^{\prime }v^{\prime }\rangle$
and
$\langle v^{\prime }v^{\prime }\rangle$
. The contribution
$-\langle h^{\prime 2}\rangle _{y}/2$
from the fluctuating pressure gradient arises from using (5.5a
) to write 
$$\begin{eqnarray}h^{2}=1+2h^{\prime }+2(\langle h\rangle -1)+h^{\prime 2}+O(a^{3}).\end{eqnarray}$$
There are also radiative relaxation terms proportional to the so-called bolus velocity components
$\langle h^{\prime }u^{\prime }\rangle$
and
$\langle h^{\prime }v^{\prime }\rangle$
. These terms would be absent for the momentum-conserving variant of the shallow water equations with relaxation of the thickness.
The effect of the Rossby wave thus appears as an effective force on the right-hand sides of the evolution equations for
$\langle u\rangle ^{\ast }$
and
$\langle v\rangle ^{\ast }$
. Our use of the thickness-weighted means
$\langle u\rangle ^{\ast }$
and
$\langle v\rangle ^{\ast }$
in place of the unweighted means
$\langle u\rangle$
and
$\langle v\rangle$
eliminates a
$\langle h^{\prime }v^{\prime }\rangle$
term that would otherwise appear in (5.6a
). This is similar to the transformed Eulerian mean (TEM) approach for stratified fluids (Andrews & McIntyre Reference Andrews and McIntyre1976a
, Reference Andrews and McIntyre1978).
A further simplification is possible by assuming that geostrophic balance holds in the meridional momentum equation (5.6c ), since the meridional length scale and velocity component are both smaller than the zonal length scale and zonal velocity component. The same scaling argument is used in the semi-geostrophic theory of atmospheric front formation. The time derivative of this meridional geostrophic balance condition gives
$$\begin{eqnarray}{\textstyle \frac{1}{2}}y\langle u\rangle _{t}^{\ast }+\langle h\rangle _{yt}=-{\textstyle \frac{1}{2}}\langle h^{\prime 2}\rangle _{yt},\end{eqnarray}$$
which can be used in place of (5.6c ) to obtain a closed system with (5.6a ) and (5.6b ).
We now evaluate these right-hand sides using the equatorially trapped Rossby wave solutions in § 4.1. These represent small perturbations about a state of rest with unit thickness in dimensionless variables. The meridional velocity perturbation is
$v^{\prime }=\unicode[STIX]{x1D6FC}\,\text{Re}\{\hat{v}(y)\exp (\text{i}(kx-\unicode[STIX]{x1D714}t))\}$
, with
$\hat{v}(y)$
defined by (4.3) and (4.5), and the constant
$\unicode[STIX]{x1D6FC}$
sets the amplitude. The corresponding zonal velocity and thickness perturbations are
$u^{\prime }=\unicode[STIX]{x1D6FC}\,\text{Re}\{\hat{u} (y)\exp (\text{i}(kx-\unicode[STIX]{x1D714}t))\}$
and
$h^{\prime }=\unicode[STIX]{x1D6FC}\,\text{Re}\{{\hat{h}}(y)\exp (\text{i}(kx-\unicode[STIX]{x1D714}t))\}$
, where
$$\begin{eqnarray}\hat{u} (y)=\text{i}\left(\frac{k(\text{d}\hat{v}/\text{d}y)-{\textstyle \frac{1}{2}}y(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})\hat{v}}{k^{2}-(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})}\right),\quad {\hat{h}}(y)=-\text{i}\left(\frac{(\text{i}\unicode[STIX]{x1D714}-\unicode[STIX]{x1D6FE})\hat{u} +{\textstyle \frac{1}{2}}y\hat{v}}{k}\right),\end{eqnarray}$$
for
$k\neq 0$
. The denominator in the expression for
$\hat{u} (y)$
vanishes when
$k^{2}=(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})$
. This gives the dispersion relation for equatorially trapped Kelvin waves, which are distinguished by having zero meridional velocity (
$\hat{v}=0$
), and frequencies
$\unicode[STIX]{x1D714}=\pm k$
in the absence of dissipation.
To allow a fair comparison between the effective forces generated by waves with different
$n$
and
$k$
values, we normalise the disturbance amplitudes by choosing
$\unicode[STIX]{x1D6FC}$
so that
$$\begin{eqnarray}\frac{1}{2}\int \langle u^{\prime 2}\rangle +\langle v^{\prime 2}\rangle +\langle h^{\prime 2}\rangle \,\text{d}y=a^{2}.\end{eqnarray}$$
The left-hand side of (5.10) is the quadratic approximation to the disturbance pseudoenergy, an exact conserved quantity of the unforced, non-dissipative shallow water equations that takes the above form when expanded for small-amplitude disturbances from a rest state of constant thickness (Shepherd Reference Shepherd1993). Using a standard formula for the average of a product of sinusoidal disturbances represented by the real parts of complex exponentials, we can compute the integrand in (5.10), and the terms on the right-hand sides of (5.6b ) and (5.6c ), using
$$\begin{eqnarray}\langle u^{\prime }v^{\prime }\rangle ={\textstyle \frac{1}{2}}\text{Re}(\hat{u} (y)\hat{v}(y)^{\dagger }),\end{eqnarray}$$
where the superscript
$\text{}^{\dagger }$
denotes a complex conjugate, and similarly for the other quadratic combinations of
$h^{\prime }$
,
$u^{\prime }$
,
$v^{\prime }$
.
Figure 9. The effective zonal forcing terms
$-\langle u^{\prime }v^{\prime }\rangle _{y}$
and
$-\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle$
on the right-hand side of (5.6b
), and their sum, for (a) dissipation solely by radiative relaxation (
$\unicode[STIX]{x1D705}=0.001,\unicode[STIX]{x1D6FE}=0$
) and (b) dissipation solely by friction (
$\unicode[STIX]{x1D705}=0,\unicode[STIX]{x1D6FE}=0.001$
). All terms are calculated for
$n=1$
,
$k=1$
, and
$a=1$
for scaling. In (b) the sum equals the first term as
$\unicode[STIX]{x1D705}=0$
. Here and in subsequent figures the legend applies to both plots.
Figure 10. The effective zonal forcing terms
$-\langle u^{\prime }v^{\prime }\rangle _{y}$
and
$-\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle$
on the right-hand side of (5.6b
), and their sum, for (a) dissipation solely by radiative relaxation (
$\unicode[STIX]{x1D705}=0.001,\unicode[STIX]{x1D6FE}=0$
) and (b) dissipation solely by friction (
$\unicode[STIX]{x1D705}=0,\unicode[STIX]{x1D6FE}=0.001$
). All terms are calculated for
$n=2$
,
$k=1$
, and
$a=1$
for scaling. In (b) the sum equals the first term as
$\unicode[STIX]{x1D705}=0$
.
Figures 9 and 10 show the two quantities
$-\langle u^{\prime }v^{\prime }\rangle _{y}$
and
$-\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle$
on the right-hand side of (5.6b
), and their sum, for waves with
$k=1$
and
$n\in \{1,2\}$
for dissipation by either Rayleigh friction (
$\unicode[STIX]{x1D6FE}>0$
,
$\unicode[STIX]{x1D705}=0$
) or radiative relaxation (
$\unicode[STIX]{x1D6FE}=0$
,
$\unicode[STIX]{x1D705}>0$
). The plots show the quantities divided by
$a^{2}$
, or as given by setting
$a=1$
. The Reynolds stress
$\langle u^{\prime }v^{\prime }\rangle$
vanishes for undamped waves, since
$\unicode[STIX]{x1D714}$
,
$A$
, and
$\hat{v}$
are then all real, while
$\hat{u}$
and
${\hat{h}}$
are purely imaginary. The bolus velocity term
$-\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle$
is explicitly proportional to
$\unicode[STIX]{x1D705}$
, and so vanishes for
$\unicode[STIX]{x1D705}=0$
.
A perturbative analysis similar to that in § 4.1 shows that the Reynolds stress
$\langle u^{\prime }v^{\prime }\rangle$
is proportional to
$\unicode[STIX]{x1D705}-\unicode[STIX]{x1D6FE}$
for small
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D6FE}$
. This explains the antisymmetry visible when comparing case (a) with
$\unicode[STIX]{x1D705}>0$
,
$\unicode[STIX]{x1D6FE}=0$
to case (b) with
$\unicode[STIX]{x1D705}=0$
,
$\unicode[STIX]{x1D6FE}>0$
in each of figures 9 and 10. Increasing
$k$
while holding
$n$
,
$\unicode[STIX]{x1D705}$
,
$\unicode[STIX]{x1D6FE}$
fixed increases the magnitude of
$\langle u^{\prime }v^{\prime }\rangle$
without changing its general shape, as does increasing
$\unicode[STIX]{x1D705}$
or
$\unicode[STIX]{x1D6FE}$
while holding
$k$
and
$n$
fixed.
The Reynolds stress convergence is close to zero at the equator when
$n$
is odd, leaving just the contribution from the bolus velocity visible in figure 9(a). However, the Reynolds stress convergence has a local extremum on the equator when
$n$
is even. Its sign for
$n=2$
is such as to produce eastward acceleration at the equator under radiative relaxation (figure 10
a), and an westward acceleration at the equator under Rayleigh friction (figure 10
b). We can infer the sign of the acceleration directly from the right-hand side of (5.6b
) on the equator. Since
$y=0$
, the Coriolis contribution from the mean meridional velocity
$\langle v\rangle ^{\ast }$
does not contribute to the left-hand side of (5.6b
). This leaves the ordinary differential equation
$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\langle u\rangle ^{\ast }|_{y=0}+\unicode[STIX]{x1D6FE}\langle u\rangle ^{\ast }|_{y=0}=-\text{e}^{-\unicode[STIX]{x1D70E}t}(\langle u^{\prime }v^{\prime }\rangle _{y}+\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle )|_{y=0,t=0},\end{eqnarray}$$
where
$\unicode[STIX]{x1D70E}=-2\text{Im}\,\unicode[STIX]{x1D714}$
is twice the decay rate of the Rossby wave calculated in § 4.1. The quadratic products on the right-hand side of (5.6b
) thus decay in proportion to
$\exp (-\unicode[STIX]{x1D70E}t)$
. The solution of (5.12) for
$\unicode[STIX]{x1D6FE}\neq \unicode[STIX]{x1D70E}$
is
$$\begin{eqnarray}\langle u\rangle ^{\ast }|_{y=0}=\left[\frac{\text{e}^{-\unicode[STIX]{x1D6FE}t}-\text{e}^{-\unicode[STIX]{x1D70E}t}}{\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D6FE}}\right](\langle u^{\prime }v^{\prime }\rangle _{y}+\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle )|_{y=0,t=0}.\end{eqnarray}$$
The time-dependent coefficient in square brackets
$[\cdots ]$
is always positive for
$t>0$
, so the direction of the zonal mean flow
$\langle u\rangle ^{\ast }$
at the equator is set by the sign of the combined Reynolds stress and bolus velocity term evaluated at
$y=0$
and
$t=0$
. This is consistent with the directions of the equatorial jets found in numerical simulations for different rates of Rayleigh friction and radiative relaxation (Scott & Polvani Reference Scott and Polvani2007, Reference Scott and Polvani2008; Warneford & Dellar Reference Warneford and Dellar2014).
However, away from the equator one would need to compute the solution to the full response problem coupling
$\langle h\rangle$
,
$\langle u\rangle ^{\ast }$
,
$\langle v\rangle ^{\ast }$
as three functions of
$y$
and
$t$
to find the mean flow produced by the wave sources on the right-hand sides of (5.6b
) and (5.6c
). Figure 11 shows the three terms
$-\langle h^{\prime 2}\rangle _{y}/2$
,
$-\langle v^{\prime }v^{\prime }\rangle _{y}$
and
$-\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle$
on the right-hand side of (5.6c
), and their sum, for
$k=1$
,
$n\in \{1,2\}$
and dissipation solely by radiative relaxation (
$\unicode[STIX]{x1D705}=0.001,\unicode[STIX]{x1D6FE}=0$
). The first two terms are non-zero even with no dissipation, and typically much larger in magnitude than either the bolus terms or the Reynolds stress convergence in the zonal equation.
Omitting the
$\unicode[STIX]{x1D705}$
terms from the right-hand sides of (5.1b
) and (5.1c
) leads to a shallow water model that relaxes the thickness field but conserves momentum. This model is no less justified than the model that supposes that
$\unicode[STIX]{x1D705}$
should not appear in (1.1b
) for evolving the velocity
$\boldsymbol{u}$
. The contributions
$\langle h^{\prime }u^{\prime }\rangle$
and
$\langle h^{\prime }v^{\prime }\rangle$
from the bolus velocity would then be absent in the right-hand sides of the analogues of (5.6b
) and (5.6c
) for this model. Figures 9(a) and 10(a) show that these contributions are in any case small compared with the contribution from the Reynolds stress convergence. This is consistent with the directions of the equatorial jets in the two models, with and without momentum conservation, showing the same behaviour in simulations (Warneford & Dellar Reference Warneford and Dellar2014).
Figure 11. The effective meridional forcing terms
$-\langle h^{\prime 2}\rangle _{y}/2$
,
$-\langle v^{\prime }v^{\prime }\rangle _{y}$
and
$-\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle$
on the right-hand side of (5.6c
), and their sum, for dissipation solely by radiative relaxation (
$\unicode[STIX]{x1D705}=0.001,\unicode[STIX]{x1D6FE}=0$
). All terms are calculated for
$k=1$
,
$a=1$
for scaling, and either (a)
$n=1$
or (b)
$n=2$
. The bolus velocity term
$-\unicode[STIX]{x1D705}\langle h^{\prime }u^{\prime }\rangle$
is much smaller than the other two terms, which are not proportional to
$\unicode[STIX]{x1D705}$
.
5.2 Thermal shallow water equations
We now apply the same approach to the dimensionless form of the unforced thermal shallow water equations (1.2 a –c) on an equatorial beta plane. The mass and zonal momentum equations are
$$\begin{eqnarray}\langle h\rangle _{t}+\langle v\rangle _{y}^{\ast }=0,\end{eqnarray}$$
and
$$\begin{eqnarray}\langle u\rangle _{t}^{\ast }-{\textstyle \frac{1}{2}}y\langle v\rangle ^{\ast }+\unicode[STIX]{x1D6FE}\langle u\rangle ^{\ast }=-\langle u^{\prime }v^{\prime }\rangle _{y}.\end{eqnarray}$$
There are no
$\unicode[STIX]{x1D705}$
terms in these equations, as radiative relaxation conserves mass and momentum in our thermal shallow water model. The temperature
$\unicode[STIX]{x1D6E9}$
also does not appear in these equations, because the zonal pressure gradient disappears under zonal averaging.
To calculate the thermal contributions to the other equations, it is convenient to decompose the temperature as
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}=1+\unicode[STIX]{x1D703}=1+\underbrace{\unicode[STIX]{x1D703}^{\prime }}_{O(a)}+\underbrace{\langle \unicode[STIX]{x1D703}\rangle }_{O(a^{2})}.\end{eqnarray}$$
This decomposition eliminates
$O(1)$
terms proportional to the reference temperature
$\unicode[STIX]{x1D6E9}_{0}$
, which is scaled to unity by the non-dimensionalisation in § 4.2. Although
$\unicode[STIX]{x1D703}^{\prime }=\unicode[STIX]{x1D6E9}^{\prime }$
, we use
$\unicode[STIX]{x1D703}^{\prime }$
below so that all equations are written using the symbol
$\unicode[STIX]{x1D703}$
.
The meridional momentum equation is
$$\begin{eqnarray}(hv)_{t}+(huv)_{x}+\left(hv^{2}+{\textstyle \frac{1}{2}}h^{2}(1+\unicode[STIX]{x1D703})\right)_{y}+{\textstyle \frac{1}{2}}yhu=-\unicode[STIX]{x1D6FE}hv,\end{eqnarray}$$
and the conservation form of the temperature equation is
$$\begin{eqnarray}(h\unicode[STIX]{x1D703})_{t}+(hu\unicode[STIX]{x1D703})_{x}+(hv\unicode[STIX]{x1D703})_{y}=-\unicode[STIX]{x1D705}(h^{2}(1+\unicode[STIX]{x1D703})-h).\end{eqnarray}$$
The unit term from the decomposition
$\unicode[STIX]{x1D6E9}=1+\unicode[STIX]{x1D703}$
disappears from the left-hand side of (5.18) because
$h$
satisfies the mass conservation equation with no source terms. However, the combination
$h^{2}(1+\unicode[STIX]{x1D703})$
appears in both the pressure gradient in (5.17) and the radiative relaxation term in (5.18). We may reduce this expression to a quadratic product by writing
$$\begin{eqnarray}h^{2}\unicode[STIX]{x1D703}=h\unicode[STIX]{x1D703}+h^{\prime }\unicode[STIX]{x1D703}^{\prime }+O(a^{3}).\end{eqnarray}$$
The zonal mean of the temperature equation may thus be written as
$$\begin{eqnarray}\langle \unicode[STIX]{x1D703}\rangle _{t}^{\ast }+\unicode[STIX]{x1D705}(\langle \unicode[STIX]{x1D703}\rangle ^{\ast }+\langle h\rangle -1)=-\langle v^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle _{y}-\unicode[STIX]{x1D705}(\langle h^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle +\langle h^{\prime 2}\rangle ).\end{eqnarray}$$
As before for the velocity components, we use
$\langle \unicode[STIX]{x1D703}\rangle ^{\ast }$
instead of
$\langle \unicode[STIX]{x1D703}\rangle$
to absorb a contribution from
$\langle \unicode[STIX]{x1D703}^{\prime }h^{\prime }\rangle$
in the time derivative. As
$\unicode[STIX]{x1D703}=O(a)$
under our decomposition, we may neglect factors of
$\langle h\rangle =1+O(a^{2})$
that we could not neglect for
$\langle \unicode[STIX]{x1D6E9}\rangle ^{\ast }$
. Similarly, the zonal mean of the meridional momentum equation is
$$\begin{eqnarray}\langle v\rangle _{t}^{\ast }+{\textstyle \frac{1}{2}}y\langle u\rangle ^{\ast }+\unicode[STIX]{x1D6FE}\langle v\rangle ^{\ast }+\langle h\rangle _{y}+{\textstyle \frac{1}{2}}\langle \unicode[STIX]{x1D703}\rangle _{y}^{\ast }=-\langle v^{\prime }v^{\prime }\rangle _{y}-{\textstyle \frac{1}{2}}\langle h^{\prime 2}\rangle _{y}-{\textstyle \frac{1}{2}}\langle h^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle _{y},\end{eqnarray}$$
so we again have four equations for evolving
$\langle h\rangle$
,
$\langle u\rangle ^{\ast }$
,
$\langle u\rangle ^{\ast }$
,
$\langle \unicode[STIX]{x1D703}\rangle ^{\ast }$
with linear left-hand sides. The right-hand sides are again zonal means of quadratic products of the fluctuating quantities
$h^{\prime }$
,
$u^{\prime }$
,
$v^{\prime }$
,
$\unicode[STIX]{x1D703}^{\prime }$
.
The meridional velocity perturbation for an equatorial Rossby wave is again
$v^{\prime }=\unicode[STIX]{x1D6FC}\,\text{Re}\{\hat{v}(y)\exp (\text{i}(kx-\unicode[STIX]{x1D714}t))\}$
, with
$\hat{v}(y)$
defined by (4.3) and (4.10). The corresponding perturbations to the thickness, temperature and zonal velocity take the same functional forms, with
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\hat{h}}(y)=-\text{i}\left(\frac{\text{i}k\hat{u} +\text{d}\hat{v}/\text{d}y}{\unicode[STIX]{x1D714}}\right),\quad \hat{\unicode[STIX]{x1D6E9}}(y)=-\text{i}\left(\frac{\unicode[STIX]{x1D705}{\hat{h}}}{\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705}}\right),\\ \displaystyle \hat{u} (y)=\text{i}\left(\frac{k(1+\text{i}\unicode[STIX]{x1D705}/(2\unicode[STIX]{x1D714}))\,\text{d}\hat{v}/\text{d}y-{\textstyle \frac{1}{2}}y(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})\hat{v}}{k^{2}(1+\text{i}\unicode[STIX]{x1D705}/(2\unicode[STIX]{x1D714}))-(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})}\right),\end{array}\right\}\end{eqnarray}$$
for
$\unicode[STIX]{x1D714}\neq 0$
and
$\unicode[STIX]{x1D714}\neq -\text{i}\unicode[STIX]{x1D705}$
. The vanishing denominator in the expression for
$\hat{u}$
when
$k^{2}=(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D705})/(1+\text{i}\unicode[STIX]{x1D705}/(2\unicode[STIX]{x1D714}))$
again gives the dispersion relation for the equatorially trapped Kelvin waves. We choose the normalisation constant
$\unicode[STIX]{x1D6FC}$
so that the disturbances
$u^{\prime }$
,
$v^{\prime }$
,
$h^{\prime }$
and
$\unicode[STIX]{x1D703}^{\prime }$
satisfy
$$\begin{eqnarray}\frac{1}{2}\int \langle u^{\prime 2}\rangle +\langle v^{\prime 2}\rangle +\left\langle \left(h^{\prime }+\frac{1}{2}\unicode[STIX]{x1D703}^{\prime }\right)^{2}\right\rangle \,\text{d}y=a^{2}.\end{eqnarray}$$
The left-hand side is the dimensionless quadratic approximation for small-amplitude disturbances to the pseudoenergy
$$\begin{eqnarray}\frac{1}{2}\int h(u^{2}+v^{2})+(h\sqrt{\unicode[STIX]{x1D6E9}}-h_{0}\sqrt{\unicode[STIX]{x1D6E9}_{0}})^{2}\,\text{d}x\,\text{d}y\end{eqnarray}$$
of finite-amplitude disturbances from a rest state of constant thickness
$h_{0}$
and temperature
$\unicode[STIX]{x1D6E9}_{0}$
in the thermal shallow water equations (Ripa Reference Ripa1995; Røed Reference Røed1997; Warneford & Dellar Reference Warneford and Dellar2013)
Figure 12. The effective zonal force
$-\langle u^{\prime }v^{\prime }\rangle _{y}$
on the right-hand side of (5.15) for
$k=1$
,
$n=1$
, and
$a=1$
for scaling. Plot (a) shows dissipation solely by radiative relaxation (
$\unicode[STIX]{x1D705}=0.002,\unicode[STIX]{x1D6FE}=0$
) and (b) shows dissipation solely by friction (
$\unicode[STIX]{x1D705}=0,\unicode[STIX]{x1D6FE}=0.001$
). The zero line is shown dotted. The ratio of
$\unicode[STIX]{x1D705}$
to
$\unicode[STIX]{x1D6FE}$
reflects the symmetry of the dispersion relation in § 4.2.
Figure 13. The effective zonal force
$-\langle u^{\prime }v^{\prime }\rangle _{y}$
on the right-hand side of (5.15) for
$k=1$
,
$n=2$
, and
$a=1$
for scaling. Plot (a) shows dissipation solely by radiative relaxation (
$\unicode[STIX]{x1D705}=0.002,\unicode[STIX]{x1D6FE}=0$
) and (b) shows dissipation solely by friction (
$\unicode[STIX]{x1D705}=0,\unicode[STIX]{x1D6FE}=0.001$
). The zero line is shown dotted. The ratio of
$\unicode[STIX]{x1D705}$
to
$\unicode[STIX]{x1D6FE}$
reflects the symmetry of the dispersion relation in § 4.2.
Figure 14. The effective meridional forcing terms
$-\langle v^{\prime }v^{\prime }\rangle _{y}$
,
$-\langle h^{\prime 2}\rangle _{y}/2$
,
$-\langle h^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle _{y}/2$
on the right-hand side of (5.21), and their sum, for dissipation solely by radiative relaxation (
$\unicode[STIX]{x1D705}=0.002,\unicode[STIX]{x1D6FE}=0$
). All terms are calculated for
$k=1$
,
$a=1$
for scaling, and either (a)
$n=1$
or (b)
$n=2$
.
Figure 15. The effective thermal terms
$-\langle v^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle _{y}$
,
$-\unicode[STIX]{x1D705}\langle h^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle$
,
$-\unicode[STIX]{x1D705}\langle h^{\prime 2}\rangle$
on the right-hand side of (5.20), and their sum, for dissipation solely by radiative relaxation (
$\unicode[STIX]{x1D705}=0.002,\unicode[STIX]{x1D6FE}=0$
). All terms are calculated for
$k=1$
,
$a=1$
for scaling, and either (a)
$n=1$
or (b)
$n=2$
.
Figures 12 and 13 show the Reynolds stress convergences on the right-hand side of (5.15). These all closely resemble those for the previous shallow water equations with radiative relaxation of the thickness. Thus the mean flow at the equator (
$y=0$
) is almost identical to that for the previous model, provided
$\unicode[STIX]{x1D705}$
is rescaled by a factor of 2. A perturbative analysis similar to that in § 4.1 shows that the Reynolds stress
$\langle u^{\prime }v^{\prime }\rangle$
for this model is proportional to
$\unicode[STIX]{x1D705}/2-\unicode[STIX]{x1D6FE}$
for small
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D6FE}$
.
Figure 14 shows the forcing terms in the right-hand side of the mean meridional momentum equation (5.21). The dominant terms,
$-\langle v^{\prime }v^{\prime }\rangle _{y}$
and
$-\langle h^{\prime 2}\rangle _{y}/2$
are very similar to those in figure 11 for the shallow water model with relaxation of the thickness. The thermal term
$-\langle h^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle _{y}/2$
is much smaller. Finally, figure 15 shows the forcing terms in the right-hand side of the mean temperature equation (5.20). These all vanish for dissipation solely by friction (
$\unicode[STIX]{x1D705}=0$
), since the temperature perturbation
$\hat{\unicode[STIX]{x1D703}}$
given by (5.22) then vanishes, so they are only shown for
$\unicode[STIX]{x1D705}>0$
and
$\unicode[STIX]{x1D6FE}=0$
. As for the meridional velocity, the contribution from
$\langle h^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle$
is much smaller than the other contributions.
5.3 Discussion
Our simulation parameters imply a dimensionless radiative relaxation rate
$\unicode[STIX]{x1D705}=7.88\times 10^{-3}$
roughly 23 times larger than our dimensionless frictional rate
$\unicode[STIX]{x1D6FE}=3.42\times 10^{-4}$
. In § 5.1 we established that the Reynolds stress
$\langle u^{\prime }v^{\prime }\rangle$
in the mean zonal momentum equation is proportional to
$\unicode[STIX]{x1D705}-\unicode[STIX]{x1D6FE}$
(for thickness relaxation) or
$\unicode[STIX]{x1D705}/2-\unicode[STIX]{x1D6FE}$
(for temperature relaxation) for small relaxation and friction rates. For
$\unicode[STIX]{x1D705}\gg \unicode[STIX]{x1D6FE}$
, our theory thus suggests that
$\langle u\rangle ^{\ast }>0$
on the equator, consistent with the super-rotating equatorial jets seen in our simulations. Conversely, for the purely frictional case (
$\unicode[STIX]{x1D705}=0$
) our theory gives
$\langle u\rangle ^{\ast }<0$
on the equator, again consistent with the sub-rotating equatorial jet in our simulation. These conclusions hold for the two shallow water models that relax either the thickness or the temperature. The Reynolds stresses in the two models are very similar for small relaxation and friction rates, provided one rescales the radiative relaxation rate
$\unicode[STIX]{x1D705}$
in the thermal model by a factor of 2, as in figures 12(a) and 13(a).
We have performed a similar analysis for decaying equatorially trapped Kelvin, Yanai, and inertia–gravity waves in the two shallow water models with radiative relaxation of the thickness or temperature. The effective zonal forces due to these waves are at least an order of magnitude smaller than those due to decaying Rossby waves of the same pseudoenergy. Thus, as expected from our global QG simulations, these other types of wave make no significant contribution to setting the direction of the equatorial jet.
Our thermal shallow water model conserves both mass and momentum in the absence of friction (
$\unicode[STIX]{x1D6FE}=0$
). This offers a simple demonstration that there is no inconsistency between our calculation of positive (eastward) zonal acceleration at the equator due to radiative dissipation of Rossby waves and the arguments reviewed in § 2 that Rossby waves carry westward pseudomomentum. The evolution equation for
$u$
may be rewritten in conservation form as
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(h(u-y^{2}/4))+\unicode[STIX]{x2202}_{x}(hu(u-y^{2}/4)+\unicode[STIX]{x1D6E9}h^{2}/2)+\unicode[STIX]{x2202}_{y}(hv(u-y^{2}/4))=0.\end{eqnarray}$$
The combination
$u-y^{2}/4$
corresponds to the zonal velocity seen in an inertial frame, such that
$\unicode[STIX]{x2202}_{x}v-\unicode[STIX]{x2202}_{y}(u-y^{2}/4)=y/2+\unicode[STIX]{x2202}_{x}v-\unicode[STIX]{x2202}_{y}u$
is the absolute vorticity (Ripa Reference Ripa1993). The combination
$h(u-y^{2}/4)$
, sometimes called absolute zonal momentum, is the beta-plane analogue of azimuthal angular momentum in spherical geometry.
The zonal mean of (5.25) is
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\langle h(u-y^{2}/4)\rangle +\unicode[STIX]{x2202}_{y}\langle hv(u-y^{2}/4)\rangle =0,\end{eqnarray}$$
which becomes
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\left(\langle u\rangle ^{\ast }-{\textstyle \frac{1}{4}}y^{2}\langle h\rangle \right)+\unicode[STIX]{x2202}_{y}\left(\langle u^{\prime }v^{\prime }\rangle -{\textstyle \frac{1}{4}}y^{2}\langle hv\rangle \right)=0,\end{eqnarray}$$
correct to
$O(a^{2})$
. Integrating (5.27) over
$y$
implies that the global absolute zonal momentum
$$\begin{eqnarray}{\mathcal{M}}=\int _{-\infty }^{\infty }\langle u\rangle ^{\ast }-\frac{1}{4}y^{2}(\langle h\rangle -1)\,\text{d}y\end{eqnarray}$$
is constant in time. We have added the time-independent term
$y^{2}/4$
to the integrand defining
${\mathcal{M}}$
so that this integral converges when
$\langle h\rangle \rightarrow 1$
as
$y\rightarrow \pm \infty$
. Evaluating
${\mathcal{M}}$
at
$t=0$
for the equatorially trapped Rossby wave solution from § 4.2 for which
$\langle h\rangle =1$
and
$\langle u\rangle =0$
initially, we obtain
$$\begin{eqnarray}{\mathcal{M}}|_{t=0}=\int _{-\infty }^{\infty }\frac{1}{2}\text{Re}\{{\hat{h}}(y)\hat{u} (y)^{\ast }\}\,\text{d}y<0,\end{eqnarray}$$
where
${\hat{h}}(y)$
and
$\hat{u} (y)$
are defined in (5.22), and
$\hat{v}(y)$
is defined by (4.3) and (4.10). This negative sign agrees with the arguments reviewed in § 2 that Rossby waves carry westward pseudomomentum.
The mass-weighted mean zonal velocity
$\langle u\rangle ^{\ast }$
generated by the decaying Rossby wave may be positive (eastward) at the equator, but the meridional integral of
$\langle u\rangle ^{\ast }-(1/4)y^{2}(\langle h\rangle -1)$
that defines
${\mathcal{M}}$
is still negative (westward). This is all that is required to reconcile conservation of zonal angular momentum with the requirement that Rossby waves should carry westward zonal pseudomomentum, as established concretely by (5.29).
6 Conclusions
Numerical simulations of shallow water equations on a sphere with isotropic random forcing reliably produce a mix of coherent vortices and alternating zonal jets. These simulations typically include Rayleigh friction to absorb the gradual inverse cascade of energy past the Rhines scale, since the zonal jets would otherwise break up into domain-sized coherent vortices after very long times. Such simulations invariably produce sub-rotating equatorial jets when run with Jovian parameters (Vasavada & Showman Reference Vasavada and Showman2005). However, simulations that model radiative effects by relaxing the thickness field produce a super-rotating equatorial jet (Scott & Polvani Reference Scott and Polvani2007, Reference Scott and Polvani2008), as do simulations of the thermal shallow water equations with a radiative relaxation term in the temperature equation (see § 3 and Warneford & Dellar Reference Warneford and Dellar2014).
We have provided a possible explanation for the different directions of the equatorial jets in different shallow water models with dissipation by Rayleigh friction or Newtonian radiative relaxation. We formulated evolution equations for the zonal mean thickness
$\langle h\rangle$
, and the thickness-weighted zonal mean quantities
$\langle u\rangle ^{\ast }$
,
$\langle v\rangle ^{\ast }$
,
$\langle \unicode[STIX]{x1D703}\rangle ^{\ast }$
. These equations contain effective forces due to zonal means of products of fluctuations, such as the Reynolds stress convergence
$\langle u^{\prime }v^{\prime }\rangle _{y}$
in the evolution equation for
$\langle v\rangle ^{\ast }$
. We evaluated these effective forces for an equatorially trapped Rossby waves superimposed on a background rest state on an equatorial beta plane. Dissipation changes the usual phase relations between the perturbations
$u^{\prime }$
,
$v^{\prime }$
,
$h^{\prime }$
, and permits a non-zero Reynolds stress
$\langle u^{\prime }v^{\prime }\rangle$
that transports zonal momentum in the meridional direction. This Reynolds stress is proportional to
$\unicode[STIX]{x1D705}-\unicode[STIX]{x1D6FE}$
for the shallow water model with relaxation of the thickness, and to
$\unicode[STIX]{x1D705}/2-\unicode[STIX]{x1D6FE}$
for the thermal shallow water model with relaxation of the temperature. The Reynolds stress divergence
$\langle u^{\prime }v^{\prime }\rangle _{y}$
is the sole forcing term in the equation for
$\langle u\rangle _{t}^{\ast }$
in the momentum-conserving shallow water models, and the larger forcing term in the Scott & Polvani (Reference Scott and Polvani2007, Reference Scott and Polvani2008) model with radiative relaxation of both thickness and momentum. In all cases, Newtonian radiative relaxation produces a Reynolds stress with the opposite sign to that due to Rayleigh friction. This suffices to determine that the zonal mean flow
$\langle u\rangle ^{\ast }$
on the equator created by the decaying Rossby wave with
$n=2$
is eastward (super-rotating) for dissipation by Newtonian radiative relaxation, and westward (sub-rotating) for dissipation by Rayleigh friction.
Our work relies upon several major simplifying assumptions. We considered trapped equatorial Rossby waves on a background rest state, which allowed us to use the known analytical expressions for these waves in terms of Hermite polynomials. A more complete theory would calculate the Rossby wave spectrum supported by a zonally symmetric background flow. Secondly, we only calculated the effective forces in the evolution equations for the thickness-weighted zonal mean quantities
$\langle u\rangle ^{\ast }$
,
$\langle v\rangle ^{\ast }$
,
$\langle \unicode[STIX]{x1D703}\rangle ^{\ast }$
. This is only sufficient to determine the evolution of
$\langle u\rangle ^{\ast }$
at the equator. Elsewhere, one would need to solve the system of three or four coupled equations to determine the Coriolis contribution to
$\langle u\rangle ^{\ast }$
from
$\langle v\rangle ^{\ast }$
. However, it is worth noting that the latitude dependence of the Reynolds stress divergence
$\langle u^{\prime }v^{\prime }\rangle _{y}$
for the
$n=1$
equatorially trapped Rossby wave, with two local maxima around
$y=\pm 1$
in figures 9(a) and 12(a), is suggestive of the presence of two ‘horns’, peaks of maximal eastward velocity, at latitudes around
$\pm 7^{\circ }$
near the edges of the broad Jovian equatorial jet, and a slightly smaller velocity precisely on the equator (see figure 1).
Finally, we assumed that the waves are excited by sources outside the equatorial region, and neglected the effect of this excitation on the zonal mean flow. Some support for this second assumption is offered by simulations of the quasilinear version of the Scott & Polvani (Reference Scott and Polvani2007, Reference Scott and Polvani2008) model in spherical geometry by Saito & Ishioka (Reference Saito and Ishioka2015) which show that the zonal mean zonal acceleration due to the ensemble of Rossby waves excited by random forcing of the wavenumbers in the annulus described in § 3, and decaying by either radiative relaxation or friction, produce an equatorial jet in the direction predicted by our arguments, and with an amplitude consistent at early times with their fully nonlinear simulations. They also gave an elegant argument to determine the direction of the zonal acceleration from the tilting of the constant phase lines in the separable Rossby wave solutions, as previously demonstrated by Andrews & McIntyre (Reference Andrews and McIntyre1976b ) for stratified fluids, and later adapted to shallow water models by Mofjeld (Reference Mofjeld1981) and Yamagata & Philander (Reference Yamagata and Philander1985). Although our work falls short of a full wave–mean flow interaction theory, we believe it offers at least a step towards a theoretical explanation of the origins of equatorial super- or sub-rotation in shallow water simulations of Jovian atmospheres.
Acknowledgements
We thank M. E. McIntyre for very detailed and stimulating referee reports that greatly improved our manuscript, and R. Salmon and E. R. Johnson for useful conversations. This work was supported by the Engineering and Physical Sciences Research Council through a Doctoral Training Grant award to E.S.W. and an Advanced Research Fellowship to P.J.D. [grant number EP/E054625/1]. The computations employed the University of Oxford’s Advanced Research Computing facilities (Richards Reference Richards2015), and the Emerald GPU-accelerated High Performance Computer made available by the Science and Engineering South Consortium in partnership with the Science and Technology Facilities Council’s Rutherford–Appleton Laboratory.
