[

# [

[
27 August 2013, in revised form 22 January 2014, 28 June 2014, 31 March 2015, 30 December 2016
###### Abstract

Numerical simulations of the shallow water equations on rotating spheres produce mixtures of robust vortices and alternating zonal jets, as seen in the atmospheres of the gas giant planets. However, simulations that include Rayleigh friction invariably produce a sub-rotating (retrograde) equatorial jet for Jovian parameter regimes, whilst observations of Jupiter show a super-rotating (prograde) equatorial jet that has persisted over several decades. Super-rotating equatorial jets have recently been obtained in shallow water simulations that include a Newtonian relaxation of perturbations to the layer thickness to model radiative cooling to space, and in simulations of the thermal shallow water equations that include a similar relaxation term in their temperature equation. Simulations of global quasigeostrophic forms of these different models produce equatorial jets in the same directions as the parent models, suggesting that mechanism responsible for setting the direction lies within quasigeostrophic theory. We provide such a mechanism by calculating the effective force acting on the thickness-weighted zonal mean flow due to the decay of an equatorially trapped Rossby wave. Decay due to Newtonian cooling creates an eastward zonal mean flow at the equator, consistent with the formation of a super-rotating equatorial jet, while decay due to Rayleigh friction leads to a westward zonal mean flow at the equator, consistent with the formation of a sub-rotating equatorial jet. In both cases the meridionally integrated zonal mean of the absolute zonal momentum is westward, consistent with the standard result that Rossby waves carry westward pseudomomentum, but this does not preclude the zonal mean flow being eastward on and close to the equator.

Super- and sub-rotating equatorial jets: Newtonian cooling versus Rayleigh friction]Super- and sub-rotating equatorial jets in shallow water models of Jovian atmospheres: Newtonian cooling versus Rayleigh friction E. S. Warneford and P. J. Dellar]Emma S. Warneford and Paul J. Dellarthanks: Email address for correspondence: dellar@maths.ox.ac.uk

## 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 2005). Following the pioneering non-divergent barotropic model of Williams (1978), shallow water theory has been widely used to model the Jovian atmosphere (Marcus 1988; Dowling & Ingersoll 1989; Williams & Yamagata 1984; Cho & Polvani 1996a, b; Iacono et al. 1999a, b; Showman 2007; Scott & Polvani 2007, 2008; Ingersoll 1990). The visible cloud deck is treated as a homogeneousm 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 1982; Dowling & Ingersoll 1989) 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 (1975) scale, though the quasilinear “zonostrophic instability” mechanism may also have a role (Srinivasan & Young 2012, 2014, 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 1996a, b; Iacono et al. 1999a, b; Showman 2007; Scott & Polvani 2007). 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 1986) and Cassini (Porco et al. 2003) missions over 20 years apart, demonstrates the remarkable stability of Jupiter’s zonal winds, and the longevity of the super-rotating equatorial jet.

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 2005). Scott & Polvani (2007, 2008) added an additional radiative relaxation term with timescale to the continuity equation in their shallow water model:

 ht+∇⋅(hu) = −(h−h0)/τrad, (1.0a) ut+(u⋅∇)u+f^z×u = −g′∇h+F−u/τfric, (1.0b)

Here is the thickness of the active layer, commonly called “height”, is the depth-averaged horizontal velocity, is the Coriolis parameter at latitude on a planet rotating with angular velocity , is the horizontal gradient operator, is a unit vector in the local vertical direction, is an isotropic random forcing, and is the timescale for Rayleigh friction. The reduced gravity is , where is the actual gravitational acceleration, is the density contrast between the cloud deck and the lower layer, and is a reference density. 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 1995; Ingersoll et al. 2007), although this interpretation is not without its critics (Walterscheid et al. 2000).

The right hand side of (1) models the effects of radiation to space with a Newtonian relaxation of towards its constant mean value on the time scale , following earlier shallow water models of the upper ocean mixed layer (Philander et al. 1984; Hirst 1986) and the terrestrial stratosphere (Juckes 1989; Polvani et al. 1995; Thuburn & Lagneau 1999). Scott & Polvani’s (2007; 2008) simulations of this system with dissipation only by thickness relaxation (no term) produced a sharply localised super-rotating equatorial jet, as do our subsequent simulations described in §3, and in Warneford & Dellar (2014).

Thermal shallow water theory introduced by Lavoie (1972) 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 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 1983; Hirst 1986; McCreary & Yu 1992; McCreary & Kundu 1988; McCreary et al. 1991; Røed & Shi 1999; Ripa 1993, 1995, 1996b, 1996a).

Warneford & Dellar (2014) simulated Jovian atmospheres using the thermal shallow water equations

 ht+∇⋅(hu) =0, (1.0a) Θt+(u⋅∇)Θ =−(Θh/h0−Θ0)/τrad, (1.0b) ut+(u⋅∇)u+f^z×u =−h−1∇(12Θh2)+F−u/τfric. (1.0c)

We introduce the symbol 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) is then the usual shallow water pressure gradient term, as rewritten for a spatially varying reduced gravity. The right hand side of (1) is a Newtonian cooling term that represents radiative relaxation towards a temperature with time scale . We show below that this form of coupling between and changes the behaviour of equatorially trapped Rossby waves in almost exactly the same way as the radiative relaxation of itself in Scott & Polvani’s (2007; 2008) model. In particular, the damping rates calculated in §4 are very similar function of wavenumber and . 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 (1a,b,c) conserve both mass and momentum in the absence of the and terms, unlike the Scott and Polvani model (1a,b). Simulations of the thermal shallow water equations for Jovian parameter values (reported in §3 and Warneford & Dellar 2014) reproduce a super-rotating equatorial jet, and more substantial mid-latitude jets than simulations of the Scott & Polvani (2007, 2008) model. Henceforth, for brevity we refer to the radiative relaxation terms as relaxation and the Rayleigh friction term as friction. We also refer to equations (1a,b) as the “standard” shallow water equations to distinguish them from the thermal shallow water equations (1a,b,c).

Quasi-geostrophic (QG) theory (Charney 1949; Obukhov 1949) 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 (1970, 1971) presented a global QG theory for fully stratified atmospheres on spheres. Verkley (2009) and Schubert et al. (2009) revived this theory for the standard shallow water equations on a sphere in the form

 qt+[ψ,q]=0,q=f+∇2ψ−L−2Dψ, (1.0)

where the Jacobian . The Rossby deformation scale is , where is the speed of long surface gravity waves. Both and vary with latitude when QG theory is used in spherical coordinates. The sole evolving variable is the potential vorticity , which is advected by the velocity field derived from a streamfunction . The elliptic equation relating to involves the spatially varying Coriolis parameter . It corresponds to Daley’s (1983) simplest form of geostrophic balance, as justified by assuming that varies on lengthscales much smaller than the planetary scales on which varies. The system (1) reduces to the familiar beta-plane QG equations on approximating by in the first term, and by in the second term. Warneford (2014) derived a thermal form of this global QG theory, building on the beta-plane thermal QG equations in Ripa (1996a) and Warneford & Dellar (2013). Numerical simulations produced super-rotating equatorial jets when the dimensionless radiative relaxation time is sufficiently short, as did simulations of the corresponding global QG form of Scott & Polvani’s (2007; 2008) shallow water model with a relaxing thickness. By contrast, both QG models produced sub-rotating equatorial jets when 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 (1969) 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 (1970) and Thompson (1971). A detailed version based on ray theory may be found in Held (2000) and Vallis (2006). The QG equations on a mid-latitude beta plane lead to the Rossby wave dispersion relation for plane-wave disturbances with streamfunction implies a meridional group velocity

 ∂ω∂ℓ=2βkℓ(k2+ℓ2+L2D)2, (2.0)

while the off-diagonal component of the Reynolds stress (see §4 for notation) is

 ⟨u′v′⟩=−12^ψ2kℓ. (2.0)

The radiation condition requires the meridional group velocity to point away from wave sources, so the sign of 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 (2009) and Liu & Schneider (2010) 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 1966; Gill 1982; McCreary 1985; Majda & Klein 2003; Khouider et al. 2013). These waves are confined by a Gaussian envelope whose extent scales with the equatorial deformation scale . This corresponds to latitudes within about 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 1976a; Bühler 2000, 2014). It defines the Lagrangian mean of an arbitrary field by

 ϕL(x,t)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕ(x+ξ(x,t),t), (2.0)

where is some linear averaging operator, and is a Lagrangian displacement field such that is the actual position of the particle whose mean position is at time t. GLM theory establishes circumstances for the validity of a “pseudomomentum rule” (McIntyre 1981; Bühler 2000, 2014) 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 linearized version of GLM for small was presented in Andrews & McIntyre (1976a, 1978). 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) that models radiative relaxation. In particular, lack of mass conservation requires modification of Bühler’s (2000) 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 (1a,b) and, secondly, that an eastward zonal acceleration at the equator is obtained for two interesting cases, (1) equatorially trapped Rossby waves freely decaying by radiative relaxation – as we compute independently below – and (2) 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 4b of Andrews & McIntyre (1976a) with suitable sign changes.

Following the approach of McIntyre & Norton (1990) for stratified fluids, Bühler (2000) formulated a general approach for shallow water equations of the form

 ht+∇⋅(hu) = 0, (2.0a) ut+(u⋅∇)u+f^z×u = −g′∇h+F (2.0b)

with a general term , not necessarily a forcing. The shallow water potential vorticity , where is the relative vorticity, evolves according to

 ∂t(hq)+∇⋅(hqu−F⊥)=0, (2.0)

with an extra flux . Bühler’s (2000) pseudomomentum rule holds for systems of this form that conserve linear momentum, those for which is the divergence of a stress tensor

However, neither of 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, (1) and (1) together imply

 (2.0)

where is the identity matrix. Moreoever, the potential vorticity for the system (1a,b) obeys

 (hq)t+∇⋅(hqu+u⊥/τfric)=0, (2.0)

where , while the potential vorticity in the thermal shallow water equations (1a–c) obeys

 (hq)t+∇⋅(hqu−12h∇⊥θ+u⊥/τfric)=0. (2.0)

Radiative relaxation thus does not contribute to either potential vorticity equation. However, (2) contains an extra term involving the perpendicular temperature gradient that arises from the curl of the term on the right hand side of (1). 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 1990).

Since the Scott & Polvani (2007, 2008) 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 and (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 numerical simulations of the standard shallow water equations with dissipation only through friction (i.e. (1a,b) with no term), the Scott & Polvani (2007, 2008) model with relaxation of the thickness, and our thermal shallow water model. Table 1 gives typical Jovian values for the planetary radius , angular velocity , gravitational acceleration , and internal wave speed . The last is deduced from impacts of Shoemaker–Levy comet debris (Dowling 1995; Ingersoll et al. 2007), and implies that there are polar deformation scales around the circumference. Simulations with radiative relaxation used a timescale of 43 Jovian days, as calculated for a temperature of K at the 25 millibar pressure level, and all simulations employed a friction with a time scale of 1000 Jovian days. Liu & Schneider (2010) used the 25 millibar pressure level when comparing their three-dimensional simulations with observations of Neptune. Calculating 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. The theory in §5 below suggests that the direction of the equatorial jet is insensitive to these parameters provived the radiative relaxation timescale remains much shorter than the frictional timescale.

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. 2012). 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 back across the equator to the north pole. The Coriolis force is still fully varying in latitude and is defined in the simulations by , where 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 discretized the domain using Fourier collocation points, and used the Hou & Li (2007) spectral filter to control the build-up 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.

Our initial conditions were and for all runs, and for the thermal shallow water simulation. Following Scott & Polvani (2007, 2008) we applied a divergence-free isotropic random forcing that is localised to a narrow annulus of wavenumbers in Fourier space, counting the longest sinusoidal mode in the domain as wavenumber . We forced each mode inside this annulus with amplitude using random phases that were -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 lengthscales comparable to the deformation radius. Our simulations slowly adjusted the amplitude of the forcing 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 (2014) and Warneford & Dellar (2014).

Figure 3 shows the instantaneous absolute vorticity after Jovian rotation periods 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 plots, and figure 5 shows the instantaneous zonally averaged zonal velocity, . All three simulation runs exhibit a mixture of vortices, turbulence and multiple zonal jets, with amplitudes that decrease at higher latitudes. Both 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 eastward and northward. The underlying spherical geometry appears only through the linearisation of for small latitudes into with , where is the planetary radius.

### 4.1 Shallow water equations with thickness relaxation

The linearised form of the unforced () shallow water equations (1) with radiative relaxation of the thickness on an equatorial beta-plane may be written as (Matsuno 1966; Gill 1982)

 h′t+u′x+v′y = −κh′, (4.0a) u′t−12yv′+h′x = −γu′, (4.0b) v′t+12yu′+h′y = −γv′. (4.0c)

The dashes denote small perturbations about a rest state with uniform depth . We have non-dimensionalised using the internal wave speed as the velocity scale, the equatorial deformation radius as the horizontal length scale, and as the thickness scale. Using the Jovian parameters in table 1, a latitude of corresponds to in the subsequent plots. The dimensionless radiative relaxation and friction rates are and based on the advective time scale . Their values were and in the simulation with radiative relaxation of the thickness.

Following standard practice (e.g.   Matsuno 1966; Gill 1982) we seek waves that are harmonic in longitude and time of the form etc. With no dissipation () the three equations (4.1) may be combined into a single ordinary differential equation (ODE) for the meridional velocity ,

 d2^vdy2=(Ay2−B)^v,A=14,B=ω2−k2−k2ω, (4.0)

the same equation that governs a quantum harmonic oscillator. The solutions that decay as may be written using the Hermite polynomials as

 ^v=Hn(ξ)exp(−ξ2/2),ξ=yA1/4. (4.0)

The dispersion relation for gives a cubic equation for ,

 ω2−k2−k2ω=n+12. (4.0)

The Rossby wave branch of solutions is characterised by and . The other two roots give inertia-gravity waves with (the dimensionless inertial frequency). The cases and 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 of the equator. The inertia-gravity, Yanai, and Kelvin wave branches all exceed the inertial frequency , and so do not exist in QG theory. 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.1a–c) change the and coefficients to (Yamagata & Philander 1985)

 (4.0)

so the dispersion relation becomes

 ω2−k2−k2(ω+iγ)+iω(κ+γ)−κγ=(n+12)(ω+iκω+iγ)1/2 for n=1,2,…. (4.0)

Figure 7 shows the real and imaginary parts of the complex frequency for the first few equatorial Rossby waves () for different values of and . The dispersion relation (4.1) is invariant under the transformation , so we choose to take . The Rossby wave frequency is now negative, and virtually unchanged by either weak friction (, ) or weak radiative relaxation (, ). The decay rate due to cooling (, ) is largest at , and tends to zero as . By contrast, the decay due to friction (, ) is monotonically increasing with , and tends to as . Figure 7b also suggests a symmetry about the line between small radiative relaxation ( and ) and small friction ( and ). We expand the roots of (4.1) as for radiative relaxation, and for friction, where is the real root of (4.1). The corrections are both purely imaginary:

 ω1,rad = 12iω0(−4ω20+2n+1)/(k+4ω30), (4.0a) ω1,fric = 12iω0(−4ω20−2n−1−2k/ω0)/(k+4ω30). (4.0b)

The symmetry about visible in figure 7b is a consequence of the relation

between the two dispersion relations truncated at .

### 4.2 Thermal shallow water equations

Linearising the thermal shallow water equations (1a–c) about a rest state with uniform depth and temperature and applying the equivalent nondimensionalisation based on the velocity scale gives

 h′t+u′x+v′y = 0, (4.0a) Θ′t = −κ(h′+Θ′), (4.0b) u′t−12yv′+h′x+12Θ′x = −γu′, (4.0c) v′t+12yu′+h′y+12Θ′y = −γv′. (4.0d)

The four equations (4.2) may again be combined into the single ODE (4.1) for , with coefficients

 A=14ω(ω+iκ)(ω+iγ)(ω+iκ/2),B=ω(ω+iγ)(ω+iκ)ω+iκ/2−k2−k2(ω+iγ), (4.0)

so the dispersion relation for equatorially trapped waves is

 ω(ω+iγ)(ω+iκ)ω+iκ/2−k2−k2(ω+iγ)=(n+12)(ω(ω+iκ)(ω+iγ)(ω+iκ/2))1/2. (4.0)

Figure 8 shows the real and imaginary parts of the complex frequency for the first few equatorial Rossby waves with , and three different sets of and values. As before, the real part of is virtually unchanged by weak dissipation ( or ) while becomes negative. The plot of in figure 8 closely resembles our earlier figure 7, with the same apparent symmetry, if we now take . A similar perturbation analysis to before shows that including this factor of indeed makes symmetric about between the cases of small radiative relaxation ( and ) and small friction ( and ).

## 5 Acceleration of the zonal mean zonal flow

We now investigate the acceleration of the zonal mean zonal flow caused by these solutions for decaying equatorially trapped Rossby waves.

### 5.1 Shallow water equations with thickness relaxation

We consider the shallow water model (1) with radiative relaxation of the thickness, and no forcing, on an equatorial beta-plane, and apply the nondimensionalisation from §4.1 to obtain

 ht+(hu)x+(hv)y =−κ(h−1), (5.0) (hu)t+(hu2+12h2)x+(huv)y−12yhv =−γhu−κu(h−1), (5.0) (hv)t+(huv)x+(hv2+12h2)y+12yhu =−γhv−κv(h−1). (5.0)

The mean thickness becomes unity under this nondimensionalisation.

Following standard practice, we decompose the velocity and other fields as into a zonal mean and a deviation . The zonal mean of the dimensionless mass conservation equation may then be written as

 ⟨h⟩t+⟨hv⟩y=−κ(⟨h⟩−1). (5.0)

Similarly, the zonal means of the two momentum equations are

 ⟨hu⟩t+⟨huv⟩y−12y⟨hv⟩=−γ⟨hu⟩−κ(⟨hu⟩−⟨u⟩), (5.0)

and

 ⟨hv⟩t+⟨hv2+12h2⟩y+12y⟨hu⟩=−γ⟨hv⟩−κ(⟨hv⟩−⟨v⟩). (5.0)

We now consider the form of solutions to the system (5.1) to (5.1) for initial conditions corresponding to a Rossby wave of small amplitude 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

 h=1+h′O(a)+(⟨h⟩−1)O(a2),u=u′O(a)+⟨u⟩O(a2),v=v′O(a)+⟨v⟩O(a2). (5.0)

We define to be the thickness-weighted zonal mean of a quantity , as in Thuburn & Lagneau (1999). Using , we can simplify the above to obtain evolution equations correct to in the form

 ⟨h⟩t+⟨v⟩∗y+κ(⟨h⟩−1) =0, (5.0) ⟨u⟩∗t−12y⟨v⟩∗+γ⟨u⟩∗ =−⟨u′v′⟩y−κ⟨h′u′⟩, (5.0) ⟨v⟩∗t+12y⟨u⟩∗+γ⟨v⟩∗+⟨h⟩y =−⟨v′v′⟩y−κ⟨h′v′⟩−12⟨h′2⟩y. (5.0)

On the left hand sides we have only the three mean quantities , , . On the right hand sides we have only means of quadratic products of the fluctuations , , . These include the meridional derivatives of the Reynolds stress components and . The contribution from the fluctuating pressure gradient arises from using (5.1) to write

 h2=1+2h′+2(⟨h⟩−1)+h′2+O(a3). (5.0)

There are also radiative relaxation terms proportional to the so-called bolus velocity components and . 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 and . Our use of the thickness-weighted means and in place of the unweighted means and eliminates a term that would otherwise appear in (5.1). This is similar to the transformed Eulerian mean (TEM) approach for stratified fluids (Andrews & McIntyre 1976a, 1978).

A further simplification is possible by assuming that geostrophic balance holds in the meridional momentum equation (5.1), since the meridional lengthscale and velocity component are both smaller than the zonal lengthscale and zonal velocity component. The same scaling argument is used in the semigeostrophic theory of atmospheric front formation. The time derivative of this meridional geostrophic balance condition gives

 12y⟨u⟩∗t+⟨h⟩yt=−12⟨h′2⟩yt, (5.0)

which can be used in place of (5.1) to obtain a closed system with (5.1) and (5.1).

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 depth in dimensionless variables. The meridional velocity perturbation is , with defined by (4.1) and (4.1), and the constant sets the amplitude. The corresponding zonal velocity and thickness perturbations are and , where

 ^u(y)=i⎛⎝k(d^v/dy)−12y(ω+iκ)^vk2−(ω+iγ)(ω+iκ)⎞⎠,^h(y)=−i⎛⎝(iω−γ)^u+12y^vk⎞⎠, (5.0)

for . The denominator in the expression for vanishes when . This gives the dispersion relation for equatorially trapped Kelvin waves, which are distinguished by having zero meridional velocity (), and frequencies in the absence of dissipation.

To allow a fair comparison between the effective forces generated by waves with different and values, we normalise the disturbance amplitudes by choosing so that

 12∫⟨u′2⟩+⟨v′2⟩+⟨h′2⟩dy=a2. (5.0)

The left hand side of (5.1) 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 depth (Shepherd 1993). 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.1), and the terms on the right hand sides of (5.1) and (5.1), using

 ⟨u′v′⟩=12Re(^u(y)^v(y)†) (5.0)

where the superscript denotes a complex conjugate, and similarly for the other quadratic combinations of , , .

Figures 11 and 11 show the two quantities and on the right hand side of (5.1), and their sum, for waves with and for dissipation by either Rayleigh friction (, ) or radiative relaxation (, ). The plots show the quantities divided by , or equivalently correspond to setting . The Reynolds stress vanishes for undamped waves, since , , and are then all real, while and are purely imaginary. The bolus velocity term is explicitly proportional to , and so vanishes for .

A perturbative analysis similar to that in §4.1 shows that the Reynolds stress is proportional to for small and . This explains the antisymmetry visible when comparing panel (a) with , panel (b) with , in each of figures 11 and 11. Increasing while holding , , fixed increases the magnitude of without changing its general shape, as does increasing or while holding and fixed.

The Reynolds stress convergence is close to zero at the equator when is odd, leaving just the contribution from the bolus velocity visible in figure 11(a). However, the Reynolds stress convergence has a local extremum on the equator when is even. Its sign for is such as to produce eastward acceleration at the equator under radiative relaxation (figure 11(a)), and an westward acceleration at the equator under Rayleigh friction (figure 11(b)). We can infer the sign of the acceleration directly from the right hand side of (5.1) on the equator. Since , the Coriolis contribution from the mean meridional velocity does not contribute to the left hand side of (5.1). This leaves the ordinary differential equation

 ddt⟨u⟩∗∣∣y=0+γ⟨u⟩∗∣∣y=0=−e−σt(⟨u′v′⟩y+κ⟨h′u′⟩)∣∣y=0,t=0, (5.0)

where is twice the decay rate of the Rossby wave calculated in §4.1. The quadratic products on the right hand side of (5.1) thus decay in proportion to . The solution of (5.1) for is

 ⟨u⟩∗∣∣y=0=[e−γt−e−σtσ−γ](⟨u′v′⟩y+κ⟨h′u′⟩)∣∣y=0,t=0, (5.0)

The time-dependent coefficient in square brackets is always positive for , so the direction of the zonal mean flow at the equator is set by the sign of the combined Reynolds stress and bolus velocity term evaluated at and . 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 2007, 2008; Warneford & Dellar 2014).

However, away from the equator one would need to compute the solution to the full response problem coupling , , as three functions of and to find the mean flow produced by the wave sources on the right hand sides of (5.1) and (5.1). Figure 11 shows the three terms , and on the right hand side of (5.1), and their sum, for , and dissipation solely by radiative relaxation (). The first two terms are non-zero even with no dissipation, and typically much larger in magntitude than either the bolus terms or the Reynolds stress convergence in the zonal equation.

Omitting the terms from the right hand sides of (5.1) and (5.1) leads to a shallow water model that relaxes the thickeness field but conserves momentum. This model is no less justified than the model that supposes that should not appear in the equation (1) for evolving the velocity . The contributions and from the bolus velocity would then be absent in the right hand sides of the analogues of (5.1) and (5.1) for this model. Figures 11(a) and 11(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 2014).

### 5.2 Thermal shallow water equations

We now apply the same approach to the dimensionless form of the unforced thermal shallow water equations (1) on an equatorial beta plane. The mass and zonal momentum equations are

 ⟨h⟩t+⟨v⟩∗y=0, (5.0)

and

 ⟨u⟩∗t−12y⟨v⟩∗+γ⟨u⟩∗=−⟨u′v′⟩y. (5.0)

There are no terms in these equations, as radiative relaxation conserves mass and momentum in our thermal shallow water model. The temperature 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

 Θ=1+θ=1+θ′O(a)+⟨θ⟩O(a2). (5.0)

This decomposition eliminates terms proportional to the reference temperature , which is scaled to unity by the nondimensionalisation in §4.2. Although , we use below so that all equations are written using .

The meridional momentum equation is

 (hv)t+(huv)x+(hv2+12h2(1+θ))y+12yhu=−γhv, (5.0)

and the conservation form of the temperature equation is

 (hθ)t+(huθ)x+(hvθ)y=−κ(h2(1+θ)−h). (5.0)

The unit term from the decomposition disappears from the left hand side of (5.2) because satisfies the mass conservation equation with no source terms. However, the combination appears in both the pressure gradient in (5.2) and the radiative relaxation term in (5.2). We may reduce this expression to a quadratic product by writing

 h2θ=hθ+h′θ′+O(a3). (5.0)

The zonal mean of the temperature equation may thus be written as

 ⟨θ⟩∗t+κ(⟨θ⟩∗+⟨h⟩−1)=−⟨v′θ′⟩y−κ(⟨h′θ⟩+⟨h′2⟩). (5.0)

As before for the velocity components, we use instead of to absorb a contribution from in the time derivative. As under our decomposition, we may neglect factors of that we could not neglect for . Similarly, the zonal mean of the meridional momentum equation is

 ⟨v⟩∗t+12y⟨u⟩∗+γ⟨v⟩∗+⟨h⟩y+12⟨θ⟩∗y=−⟨v′v′⟩y−12⟨h′2⟩y−12⟨h′θ′⟩y, (5.0)

so we again have four equations for evolving , , , with linear left hand sides. The right hand sides are again zonal means of quadratic products of the fluctuating quantities , , , .

The meridional velocity perturbation for an equatorial Rossby wave is again , with defined by (4.1) and (4.2). The corresponding perturbations to the thickness, temperature, and zonal velocity take the same functional forms, with

 ^h(y) =−i(ik^u+d^v/% dyω),^Θ(y)=−i(κ^hω+iκ), ^u(y) (5.0)

for and . The vanishing denominator in the expression for when again gives the dispersion relation for the equatorially trapped Kelvin waves. We choose the normalisation constant so that the disturbances , , and satisfy

 12∫⟨u′2⟩+⟨v′2⟩+⟨(h′+12θ′)2⟩dy=a2. (5.0)

The left hand side is the dimensionless quadratic approximation for small amplitude disturbances to the pseudoenergy

 12∫h(u2+v2)+(h√Θ−h0√Θ0)2dxdy (5.0)

for finite amplitude disturbances from a rest state of constant depth and temperature in the thermal shallow water equations (Ripa 1995; Røed 1997; Warneford & Dellar 2013)

Figures 13 and 13 show the Reynolds stress convergences on the right hand side of (5.2). These all closely resemble those for the previous shallow water equations with radiative relaxation of the thickness. Thus the mean flow at the equator () is almost identical to that for the previous model, provided is rescaled by a factor of . A perturbative analysis similar to that in §4.1 shows that the Reynolds stress for this model is proportional to for small and .

Figure 14 shows the forcing terms in the right hand side of the mean meridional momentum equation (5.2). The dominant terms, and are very similar to those in figure 11 for the shallow water model with relaxation of the thickness. The thermal term is much smaller. Finally, figure 15 shows the forcing terms in the right hand side of the mean temperature equation (5.2). These all vanish for dissipation solely by friction (), since the temperature perturbation given by (5.2) then vanishes, so they are only shown for and . As for the meridional velocity, the contribution from is much smaller than the other contributions.

### 5.3 Discussion

Our simulation parameters imply a dimensionless radiative relaxation rate roughly 23 times larger than our dimensionless frictional rate . In §5.1 we established that the Reynolds stress in the mean zonal momentum equation is proportional to (for thickness relaxation) or (for temperature relaxation) for small relaxation and friction rates. For , our theory thus suggests that on the equator, consistent with the super-rotating equatorial jet in simulations. Conversely, for the purely frictional case () our theory gives on the equator, again consistent with the sub-rotating equatorial jet in simulations. These conclusions hold for both shallow water models, with either the thickness or the temperature being relaxed. The Reynolds stresses in the two models are very similar for small relaxation and friction rates, provided one rescales the radiative relaxation rate in the thermal model by a factor of 2, as in figures 13(a) and 13(a).

We have performed a similar analysis for decaying equatorially trapped Kelvin, Yanai, and inertia-gravity waves in our 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 (). This offers a simple demonstration that there is no inconsistency between our calculations of positive (eastward) zonal accelerations at the equator due to radiative dissipation of Rossby waves and the arguments of §2 that Rossby waves should carry westward pseudomomentum. The evolution equation for may be rewritten in conservation form as

 ∂t(h(u−y2/4))+∂x((hu(u−y2/4)+Θh2/2)+∂y(hv(u−y2/4))=0. (5.0)

The combination corresponds to the zonal velocity seen in an inertial frame, such that