## Abstract

Depth-dependent barotropic instability of short mixed Rossby–gravity (MRG) waves is proposed as a mechanism for the formation of equatorial zonal jets. High-resolution primitive equation simulations show that a single MRG wave of very short zonal wavelength and small to moderate amplitude is unstable and leads to the development of a largely barotropic, zonally symmetric flow, featuring a westward jet at the equator and extra-equatorial jets alternating in direction with latitude. At higher but still moderate amplitude, westward flow still prevails at the equator at depths of maximum horizontal velocity amplitude in the initial wave, but the long-term equilibrated state can also feature eastward “superrotating” jets at the equator near the depths of zero horizontal velocity in the initial wave. The formation of the superrotating jets in the simulations is found to be sensitive to the inclusion of the nontraditional Coriolis force in the equations of motion. A linear theory is used to demonstrate the existence of exponentially growing horizontally nondivergent perturbations with a dominant zonally symmetric zonal velocity component. An argument for the sense and positioning of the jets relative to the equator is given in terms of inertial instability and the meridional mixing of planetary vorticity by the small zonal-scale components of the linearly unstable modes. In the long time evolution of the flow, if the amplitude of the westward equatorial jet becomes too great, zonally symmetric inertial instability limits the growth of the jets, and inertial adjustment leads to the homogenization of potential vorticity in latitude and depth around the equator.

## 1. Introduction

Understanding the formation and evolution of zonally symmetric zonal flows is a central problem in the theory of atmospheric and oceanic circulation. Examples include the equatorial jets on Jupiter (Vasavada and Showman 2005) and the zonal jets in the circulations of all three equatorial oceans (Ollitrault et al. 2006; Richards et al. 2006). In the latter case, the jets have a rich latitudinal and vertical structure not easily explained by traditional arguments.

The formation of zonal jets by the destabilization of short low-frequency planetary waves was first studied, if indirectly, by Lorenz (1972) and Gill (1974). Building on the result of Lorenz, Gill examined the stability of barotropic Rossby waves on the midlatitude *β* plane and showed that a short Rossby wave is prone to a longer wave instability that favors the creation of low-frequency zonal flows if the base wave flow involves predominantly meridional motions. In the absence of the *β* effect, his results are related to those of Meshalkin and Sinai (1961) for the Kolmogorov sinusoidal base flow in the two-dimensional Navier–Stokes equations, with instability occurring preferentially for perturbations of longer wavelength. Such a large-scale instability is deemed to be the fundamental mechanism underlying the well-known inverse energy cascade of two-dimensional turbulence (Sivashinsky 1985). This idea was applied by Manfroi and Young (1999), who added friction and nonlinearity to the problem of Gill and showed how the destabilization of the Rossby wave leads via an inverse energy cascade to a steady zonal flow consisting of narrow eastward jets and broad westward jets whose widths and number depend on *β* and the amount of friction.

The method of Gill has been adapted recently to the problem of the formation of equatorial jets through the destabilization of mixed Rossby–gravity (MRG) waves by d’Orgeville et al. (2007) and Hua et al. (2008). They note that the time variability of motions in the western boundary region of the equatorial oceans is such as to excite zonally short, eastward-propagating wave trains that are prone to barotropic instability of the type identified by Gill for short barotropic Rossby waves. The present article is a continuation of those studies. Using high-resolution numerical simulations of the primitive equations on an equatorial *β* plane, we investigate the instability of MRG waves of even shorter wavelength and concentrate on waves of the lowest baroclinic mode. The resulting zonal jets are more barotropic and of narrower meridional scale than in the earlier studies, leading to the emergence of various interesting phenomena not considered before, including the formation of eastward “superrotating” equatorial jets.

A theme of the present work is the role of equatorial inertial instability in the evolution and equilibration of the zonal jets. Dunkerton (1981) and Stevens (1983) showed how meridionally sheared, zonally symmetric zonal flows at the equator are inertially unstable, leading to so-called Taylor vortices in the vertical–meridional plane. The physical origin of the instability is the decrease of angular momentum in the absolute reference frame with perpendicular distance from the axis of rotation (Rayleigh 1917). As we will see, this instability can also come about from the growth of a narrow enough westward jet at the equator, precisely the type of flow generated by the barotropic instability of the zonally short MRG wave. The nonlinear equilibration of inertial instability has been studied numerically by Hua et al. (1997) and Griffiths (2003), who find that the adjustment toward inertial stability leads to the homogenization of potential vorticity (PV) extending over a latitude interval wider than the initial instability. Kloosterziel et al. (2007) show that when inertial instability acts alone, the region of eventually mixed PV can be accurately estimated by assuming global conservation of angular momentum between the unstable and equilibrated states. Potential vorticity homogenization resulting in the formation of “PV staircases” is related to zonal jets by global conservation and stability requirements, as shown by Dunkerton and Scott (2008).

Equatorial superrotation, observed in some of our simulations, is unusual in the sense that it implies the concentration of zonal-mean angular momentum at a point, thus requiring, in the absence of friction, the redistribution of angular momentum by eddies (defined here as departures from zonally symmetric flow) to act nondiffusively (Hide 1969). It is, however, observed in various contexts, notably the circulation in the atmospheres of Jupiter and Venus. It can spontaneously emerge in idealized numerical models of the troposphere (see Kraucunas and Hartmann 2005, and references therein), although it does not occur in the real troposphere. Explanations for equatorial superrotation usually depend on convergence of planetary wave momentum flux at the equator (Holton and Lindzen 1972) or forcing by a prescribed (zero-mean) equatorial heat or momentum source (Suarez and Duffy 1992; Saravanan 1993; Williams 2003).

The superrotation observed in our simulations turns out to be especially sensitive to the inclusion of the Coriolis force due to the poleward component of planetary rotation. The neglect of these terms, known as the “traditional approximation” (Eckart 1960), is a part of the derivation of the hydrostatic primitive equations. Wave solutions to the equatorial *β*-plane equations linearized about a state of rest have been calculated (Fruman 2009) and show that the effect of the nontraditional terms is small like (Ω*/*N**)^{2} in the relative amplitudes and phases of the velocity components, where Ω* is the rate of planetary rotation and *N** is the buoyancy frequency. This is consistent with arguments by Phillips (1968) that the terms are negligible when Ω*/*N** ≪ 1. In the simulations considered here, Ω*/*N** ≈ 0.1, so a priori, one would not expect a large effect. For more general applications, White and Bromley (1995) estimate that the neglected terms should be considered important when 2Ω**H**/*U** ≪ 1, where *H** is the vertical scale and *U** the velocity scale of the motions of interest, citing as significant the acceleration of a fluid parcel by orders of tens of centimeters per second per kilometer of angular-momentum-conserving vertical displacement. Colin de Verdiere and Schopp (1994) argue that for equatorial flows that vary on a horizontal scale of less than *H***a**, where *a** is the planetary radius, the nontraditional Coriolis force is not negligibly small. For barotropic motions in the equatorial ocean, *H***a** ∼ 200 km, which is comparable with the widths of the equatorial jets observed in the simulations described here. A growing literature on the subject of the traditional approximation is summarized in Gerkema et al. (2008).

The paper is organized as follows: In section 2 the numerical model and simulations are presented along with a description of typical results. In section 3, the instability of a low-frequency MRG wave is discussed, beginning with a linear instability analysis explaining the formation of quasi-barotropic zonal jet structures. The predictions of the linear model are compared to the simulations, with the particular physics of the equatorial region being used to explain the differences. Finally, in section 4 the arrest of the growth of the zonal jets and the equilibration of the flow is discussed, with potential vorticity homogenization across the equator being the characteristic feature of the final state and inertial instability (when the amplitude of the equatorial westward jet becomes too great) indicated as the controlling mechanism.

## 2. Numerical simulations

### a. Model description and configuration

The primitive equations code we have used for all the numerical simulations is the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams 2005) which is based on third-order advection schemes whose accuracy proved to be particularly important for our problems. The code is both highly parallelized and vectorized. Nonlinear Smagorinsky lateral diffusion was used. None of the principal results are dependent on the Smagorinsky scheme, as has been demonstrated by experiments without it. Minimal biharmonic lateral diffusion and viscosity (with a damping time scale of 3 days on the grid scale), which act as a very high-pass filter for small spatial scales, were also used. Laplacian diffusion and viscosity with values close to those of molecular diffusion were used in the vertical.

The simulations integrate the initial value problem of a finite-amplitude MRG wave in an equatorial *β*-plane channel allowed to evolve freely in time. Boundary conditions are free-slip at the bottom, north, and south, with a free upper surface and periodicity in the zonal direction. The parameters varied between simulations are the horizontal wavelength and the amplitude of the basic state wave. We focus on the low-frequency, very short zonal wavelength limit (that is, with the zonal wavelength being short compared to 2*π* times the equatorial radius of deformation for the lowest baroclinic mode) and relatively low amplitudes (with the maximum horizontal velocity being small compared to the phase speed of the lowest baroclinic-mode Kelvin wave).

The problem can be considered in nondimensional terms, using

where *u**, *υ**, and *w** are the components of velocity in the *x**, *y**, and *z** directions; is the vertical component of relative vorticity; *t** is time; *β** is the gradient of planetary vorticity at the equator; *ν** ≡ *νπ*/*H** (with *ν* being an integer and *H** the depth of the channel) is the vertical wavenumber of the MRG wave to be considered; and *c** ≡ *N**/*ν** is the phase speed of the Kelvin wave of the same vertical mode (with *N** the constant Brunt–Väisälä frequency).

The crucial nondimensional parameters involved are the meridional velocity amplitude of the wave (the Froude number) *V*_{0} and the nondimensional zonal wavenumber *k* ≡ −2*πL*** _{R}*/

*λ** (where

*λ** is the dimensional zonal wavelength,

*L**

*is the equatorial deformation radius, and the negative sign refers to westward phase propagation of the wave). We consider waves with the lowest available vertical mode*

_{R}*ν*= 1, nondimensional zonal wavenumbers

*k*between −6 and −17, and meridional velocity amplitudes

*V*

_{0}between 0.05 and 0.15. In this wavelength regime, the MRG wave has characteristics similar to a short, barotropic, zonally propagating Rossby wave, the case considered by Gill (1974).

The results apply equally to ocean- and atmosphere-scale motions, provided the configurations respect the nondimensional parameter ranges required for the barotropic instability mechanism to act. The vertical boundary conditions, while not central to the phenomena being studied, can represent either the ocean bottom and surface or a rigid (land) surface and the tropical tropopause. The simulations presented below were performed with an idealized configuration representing the equatorial Atlantic because the model and the subgrid-scale schemes are optimized for such a case. Similar results (not shown) have been obtained for cases with the corresponding atmospheric parameter values. Table 1 summarizes the nondimensional ranges of parameters used and the corresponding values representative of both oceanic and tropospheric MRG waves.

The simulations described in the next section were run in a 30° wide, 5000-m-deep equatorial *β*-plane channel, with horizontal resolution 0.1° × 0.1° and between 100 and 200 vertical levels. By comparison, the radius of deformation of the lowest baroclinic mode is 3.4° and the Kelvin wave phase speed 3.2 m s^{−1}. For the stratification used in the simulations (*N** = 2 × 10^{−3} s^{−1}), the waves have periods between 52 and 143 days.

Simulations were also performed with the model modified to include the nontraditional Coriolis force terms.

### b. Simulation results

We observe that in the low-frequency limit the waves are unstable, and their destabilization leads to the formation of zonal jets, alternating in direction with latitude and largely barotropic in the depth ranges around the maximum horizontal velocity amplitude of the initial wave. At such depths, there are westward jets centered on the equator. In certain cases, we observe the formation of vertically short superrotating (eastward) equatorial jets between the depths of the westward jets and the depths of the zero horizontal velocity of the initial wave.

Figure 1 shows equatorial cross sections of initial meridional velocity and meridional and zonal velocity as the wave destabilizes. The wave in this case had zonal wavelength *λ** = 2.5° and amplitude 0.32 m s^{−1}, corresponding to a nondimensional zonal wavenumber *k* = −8.4 and nondimensional amplitude *V*_{0} = 0.1. Observe the weakening of the wave and the development of a strong vertically sheared equatorial zonal flow (and the shearing of the wave by the zonal flow). This case features superrotating jets just above and below middepth. For atmospheric-scale parameter settings, with an initial wave amplitude of 10 m s^{−1} the westward equatorial jet reaches a maximum amplitude of about 10 m s^{−1} and the superrotating jet reaches up to 5 m s^{−1}. Also shown is the total kinetic energy as a function of time, together with the contributions from the meridional and zonal velocity components as well as from the zonal-mean zonal velocity component. The kinetic energy is initially dominated by the contribution from the meridional velocity and in the long time limit by the contribution from the zonal-mean zonal velocity (the zonal jets).

Figure 2 shows depth versus latitude of zonal-mean zonal velocity in the long time equilibrated state for three experiments with waves of amplitude 0.32 m s^{−1} and zonal wavelengths of 3.3°, 2.5°, and 1.7° (*k* = −6.3, −8.4, and −12.6). Also shown is depth versus latitude of early-time meridional velocity at a longitude where its amplitude is maximum. Note the largely barotropic zonal jets alternating in direction with latitude that form in the latitudes where the initial wave is strongest. For a shorter initial wave, the meridional scale of the jets is shorter, the extra-equatorial jets are more pronounced, and superrotating jets appear at the equator near middepth. The formation of superrotating jets depends also on the amplitude of the initial wave; it does not occur for waves below a threshold amplitude.

Simulations initialized with higher vertical mode MRG waves and low meridional mode equatorial Rossby waves were also performed and similarly reveal the formation of zonal jets. In the higher vertical mode simulations, when superrotating jets form, they occur with the same vertical periodicity as the basic state, indicating that vertical boundary conditions are probably not involved in their formation. Here we concentrate on the limit of large zonal wavenumbers, the susceptibility of the jets to inertial instability, and the implications for potential vorticity homogenization and will not discuss these other variations of the problem any further.

We note that the MRG wave is not an exact time-dependent solution to the nonlinear equations, unlike the Rossby wave to the barotropic vorticity equation, and thus part of the short time behavior observed in the simulations is the adjustment due to the nonlinear self-interaction of the initial wave. However, as will be shown in the following section, there is significant cancellation between the nonlinear terms so that in the parameter range considered here, the MRG wave is a very good approximate solution to the nonlinear equations, and the growth of the zonal jets results from barotropic instability due to the horizontal shear in the horizontal velocity of the initial condition. The nonlinear terms are important only insofar as they project onto the fastest-growing instability modes of the linearized perturbation equation.

This is supported by the results of Hua et al. (2008), where in simulations with smaller |*k*| and larger *V*_{0} than are considered here, such that the total horizontal shear is weaker, the initial wave takes much longer to break down (if it does at all). Whereas in the simulations discussed here the initial wave is largely replaced with zonal jets inside one wave period (as can be seen in the meridional velocity component in Fig. 1 where the period of the wave is 0.2 yr), with *k* = −3 an integration of more than three periods of the initial wave was required for the zonal jets to grow to dominate the flow, and with *k* ≳ −2 the shear instability mechanism does not predominate and the initial wave persists throughout the simulations, not producing jets and adjusting only very slowly through nonlinear triad interactions.

## 3. Barotropic instability of an MRG wave

Although the MRG wave is inherently baroclinic, with a sinusoidal structure in the vertical, the numerical simulations suggest that in the short wavelength, low-frequency limit, its destabilization leads to largely barotropic structures. As such we propose that its destabilization is predominantly due to barotropic instability of the type studied by Gill (1974) for the case of a barotropic Rossby wave. We show below that in the given limit, assuming horizontally nondivergent perturbations, the perturbation vorticity equation reduces to almost the same equation as in the case considered by Gill, with the depth dependence of the basic state entering as a parameter affecting the growth rate (or frequency) of perturbations.

In the earlier paper of Hua et al. (2008), where the configurations of the experiments were not zonally periodic so that different vertical modes of the perturbation were important in different parts of the domain, the approximation of horizontally nondivergent perturbations could not be made (although it was considered as a special case). In addition, Hua et al. considered zonally longer initial waves and so the time variation of the basic state had to be included in the perturbation equation. Thus, the problem considered here has two advantages that simplify the analysis, which can be framed here in terms of a single scalar streamfunction.

Consider the equation for the vertical component of vorticity on an equatorial *β* plane (*β* being unity in the nondimensional representation):

The MRG wave is a solution to (3.1) linearized about a state of rest. It has velocity components and vorticity:

where *V*_{0} is constant (the Froude number for the perturbation equation to be derived below) and where the frequency *ω* and wavenumber *k* satisfy the dispersion relation

with the branches of *ω*(*k*) chosen such that *ω* > 0 for all *k* and *k* < 0 for westward-propagating waves.

Although the MRG wave (3.2) is only a linear solution to the governing equations, the nonlinear terms are small, particularly in the low-frequency limit, and it is a very good approximation to an exact solution to the nonlinear equations. Substituting (3.2) into the nonlinear terms in the vorticity Eq. (3.1) gives

There is significant cancellation between the horizontal advection of relative vorticity terms, much like in the case of the barotropic vorticity equation. Whereas in that case the cancellation is complete and the Rossby wave is an exact nonlinear solution, here the nonlinear terms are nonzero. But when *k* ≪ −1 and *V*_{0} ∼ |*k*|^{−1}, they are small like *V*_{0}|*k*|^{−1} compared to the leading linear terms (*ζ _{t}* +

*υ*). In the parameter range of the simulations presented,

*V*

_{0}≈ 0.1 and

*k*is on the order of 10, so

*V*

_{0}|

*k*|

^{−1}is on the order of 0.01.

### a. Short-wave limit of the perturbation vertical vorticity equation

We consider the vorticity equation [(3.1)], linearized about the MRG wave in the limit of *k* ≪ −1, where the dispersion relation tends to *ω* = −*k*^{−1}.

Anticipating growing modes of a spatial scale on the order of the wavelength of the basic state wave and a growth rate on the order of its maximum horizontal velocity shear, we rescale the independent variables according to

We introduce the perturbation velocities *u*′ and *υ*′ and relative vorticity *ζ̃*′ ≡ *υ*′_{ξ} − *u*′_{η} = *k*^{−1}*ζ*′ and seek solutions with *w*′ ≡ 0 (hence “quasi-barotropic”) and with *u*′ and *υ*′ allowed to vary with depth. Linearizing (3.1) about (3.2) gives

where q.t. refers to terms quadratic in disturbance quantities.

In the *k* ≪ −1 limit, the MRG wave velocity components and vorticity are approximately

where we have used the identity

The *O*(*V*_{0}*k*^{2}) terms in the perturbation vorticity equation thus give

with the “*β* term” affecting the *O*(1) vorticity; advection of the perturbation by the basic state *U* component and the time dependence of the basic state only enter at *O*(*V*_{0}*k*). Here *E* represents the self-interaction of the basic state. When *V*_{0}|*k*|^{−1} ≪ 1, *E* acts much like an initial seed perturbation for the instability, as shown below. Equation (3.10) is similar to Eq. (6.20) in Hua et al. (2008), with the notable exception of the right-hand side (rhs) terms. The rhs terms are not crucial for understanding the high-vertical-mode equatorial deep jets in ocean basin simulations that were the focus of that paper.

Because we have assumed horizontally nondivergent perturbations, we can define a streamfunction *ψ*′ such that *u*′ = −*ψ*′* _{η}* and

*υ*′ =

*ψ*′

_{ξ}. Inserting (3.7) into (3.10) and approximating

*V*and

*Z*by their leading terms from (3.7) gives

Following the approach of Gill (1974) for the barotropic Rossby wave problem, Eq. (3.11) is solved using a Floquet series technique. The solution *ψ*′(*ξ*, *η*, *z*, *τ*) and the rhs *E* are expanded in the series

where in the present case the only nonzero terms in the expansion of *E* are *E*_{−2} and *E*_{2}. By construction, the form (3.12) has the same periodicity as the coefficients on the left, which must be the case because the rhs *E* also has that same periodicity. In the language of Floquet theory, the Floquet exponent that would normally appear in (3.12), here amounting to a wavenumber shift applied to all terms in the series, must be an integer and has been chosen equal to zero without loss of generality.

An approximate solution is obtained by substituting a truncated form of (3.12) into (3.11), where the sum is over −*M* < *n* < *M*, and setting the coefficients of *e ^{inξ}* on the left-hand side (lhs) (where again −

*M*<

*n*<

*M*) equal to their counterparts on the right-hand side. This yields the system of 2

*M*+ 1 coupled PDEs:

where

Because there are no *z* derivatives in (3.14), *z* may be considered as a parameter. The system may be solved by finite-difference discretization in *η* to yield the system of ODEs:

where ** ψ** and

**E**are vectors whose elements represent, respectively, the values of the 2

*M*+ 1

*ψ*′

*and*

_{n}*E*at each of the discrete

_{n}*η*points; 𝗔 and 𝗕 are (2

*M*+ 1)

*N*× (2

_{p}*M*+ 1)

*N*matrices, where

_{p}*N*is the number of points used in the discretization of

_{p}*η*.

Equation (3.16) is solved by projecting onto the eigenvectors of the matrix 𝗕^{−1}𝗔. We denote the eigenvectors of 𝗕^{−1}𝗔 by **w*** _{j}* and the corresponding eigenvalues by

*μ*. Let

_{j}**and 𝗕**

*ψ*^{−1}

**E**be written as linear combinations of the

**w**

*according to*

_{j}where *ψ̂ _{j}* and , the coefficients of the

*j*th eigenvector of 𝗕

^{−1}𝗔 in the expansions of the streamfunction and of

*Ê*

^{(B)}, are functions of

*τ*and

*z*. The general solution to (3.16) may be written

where is the projection of the perturbation onto **w*** _{j}* at time

*τ*= 0.

It is evident from (3.18) that the rhs *E* acts much like an initial value of the perturbation and that the growth of small perturbations to large amplitude depends on the existence of eigenvalues of 𝗕^{−1}𝗔 with nonzero real parts. Also evident is that the growth rate of such perturbations is proportional to cos(*z*), so that the vertical profiles of perturbation fields become deformed with time.

Where cos(*z*) = 0, (3.18) implies the growth of mode *ψ̂ _{j}* linear in time with rate . The part of

*E*that is nonzero where cos(

*z*) = 0 is

*O*(

*V*

_{0}

*k*

^{−3}), and its effect is very small when

*τ*is

*O*(1).

The eigenvalues and eigenvectors of 𝗕^{−1}𝗔 for the case *k* = −8.4 have been calculated for a nine-term truncation of the Floquet series (*M* = 4), with a discretization of 201 points in −50 < *η* < 50. Eigenvalues are either pure real or pure imaginary and come in matching positive and negative pairs with corresponding eigenvectors of the same spatial structure.

Figure 3 shows the total streamfunction for the fastest-growing even mode. Notice that the next-largest terms in the Floquet series after the zonal jets feature meridional motion across the axes of the jets, in particular cross-equatorial motion. Also shown in Fig. 3 are the zonal-mean zonal velocity components of the three fastest-growing even modes (even in terms of their symmetry about the equator). Eigenvectors come in even and odd symmetric pairs with nearly identical growth rates, and modes with higher growth rates have fewer and wider jets.

The latitudes of the minima in the zonal-mean zonal velocity of the most unstable mode are equally spaced with peaks separated by 3.36*π*, or 1.68 times the *zonal* wavelength of the basic state. This is close to the value of 1.55 found by Gill for the wavelength in the transverse direction of the leading term in the fastest-growing disturbance to a barotropic Rossby wave. Each of the next-fastest-growing modes simultaneously exhibits meridional scales both larger and smaller than that of the most unstable mode, as is evident in Fig. 3 from the oscillations with roughly equal amplitude near the equator whose wavelength decreases with decreasing growth rate and the isolated global maxima whose distance from the equator increases with decreasing growth rate. Note that the spacing between adjacent extrema in even the slowest-growing unstable modes is wider than 2*π* (the nondimensional zonal wavelength of the basic state). This transfer of energy to larger spatial scales is to be anticipated because of the need to globally conserve enstrophy in barotropic dynamics, as discussed by Gill (1974). Oscillatory modes, with pure imaginary eigenvalues, have meridional scales smaller than the zonal scale of the basic state.

Where cos(*z*) > 0, the mode with the largest negative eigenvalue should dominate, and where cos(*z*) < 0, the mode with the largest positive eigenvalue should dominate. The symmetry of the eigenvectors is such that if the zonal-mean components of the eigenvectors corresponding to eigenvalues *μ* and −*μ* are identical, the leading zonally varying components are equal in magnitude and shape but opposite in sign. This is consistent with the observation in the simulations that whereas the zonal-mean zonal velocity is symmetric with respect to middepth (because of the bias of the physics toward westward flow at the equator), the zonally varying signal (in particular the meridional velocity) is antisymmetric with respect to middepth. Figure 3 (lower right) shows a schematic of the vertical structure of the zonal velocity at the equator at four equally spaced times assuming an rhs *E* in (3.10) with amplitude 0.001 cos^{2}(*z*). Plotted is the sum of the respective solutions (3.18) for *μ* and −*μ*.

### b. Comparison with simulations

Figure 4 compares the most unstable even linear barotropic instability mode to the zonal-mean zonal velocity near the free upper surface in the simulations (the flow near the rigid bottom is almost identical because of the vertical symmetry of the initial condition and the barotropic nature of the instability). Observe that the jet curvature at the equator is well predicted by the linear theory, and the number of zeros and hence the number of extra-equatorial jets increases with |*k*| as predicted.

The selection of the ultimately dominant mode in the simulations can be partially explained by appealing to the special physics of the equatorial region.

First, although unstable modes come in even and odd pairs with almost identical growth rates, odd modes cannot grow at the equator because they are inertially unstable. The meridional shear in zonal velocity at the equator implicit in odd modes implies the displacement of the maximum in zonal-mean absolute zonal angular momentum from the equator.

Second, the fact that the simulations produce a westward jet at the equator and not an eastward jet, although the eigenvectors of the linear problem are obviously sign reversible, can be explained by the fact that the rectification of the zonally variable component of the velocity field, which features meridional motion across the jet axes (see Fig. 3), will lead to potential vorticity mixing within the jets. This encourages westward zonal mean flow at the equator and is self-reinforcing because the weakening of PV gradients facilitates the meridional movement of fluid parcels, a phenomenon referred to as the PV Phillips effect by Dritschel and McIntyre (2008). An eastward jet would, on the other hand, inhibit meridional motion and inhibit the growth of the instability. Although it is probably a weak effect at early times, given that the *β* term is small compared to the relative vorticity advection in the large |*k*| limit, this might be enough to create the preference for westward flow at the equator.

Referring to Fig. 4, notice that although jet positioning is well predicted at the lower values of |*k*| (Fig. 4, top), at the higher values, the extra-equatorial jets are farther from the equator in the simulations than in the linear solutions. This could be attributable to the tendency of the system to preserve barotropic stability and thus to favor a slower-growing combination of modes with weaker vorticity gradients.

The growth rate is well predicted only at low to moderate |*k*| (where |*k*| is still greater than unity) and low *V*_{0} [see Hua et al. (2008) for more on that regime]. At the higher |*k*| investigated here and at higher *V*_{0}, the growth rate begins to diverge from the prediction of the linear theory. That the jets in the simulations tend to grow more slowly than predicted by the linear theory, especially for larger values of |*k*|, might also have been expected because for higher *|k|,* the narrower jet centered on the equator is more prone to inertial instability if its axis drifts off the equator. The resulting adjustment through small-vertical-scale cells would invalidate the assumption of quasi-barotropic perturbations and draw energy from the basic state that is no longer available to the growing nearly barotropic jets.

The vertical structure of the fields is difficult to compare with the linear theory because the specific structure depends on the initial conditions and the projection of the small terms in **B**^{−1}**E** on the eigenvectors. The difficulty is related to the fact that the time and *z* dependence are not separable. Nevertheless, we note that the vertical structure of the zonal velocity at the equator, for example, develops in a qualitatively similar way to what is shown in Fig. 3 based on *E* ∝ cos^{2}(*z*).

When the vertical viscosity and diffusion are increased in the simulations (by as much as three orders of magnitude), the horizontal structure of the perturbations near the top and bottom is virtually unchanged, supporting the theory that the primary instability is barotropic.

The linear theory for the instability of an MRG wave in the large-*k*, small-*V*_{0} limit, assuming horizontally nondivergent perturbations, thus predicts the growth of signals dominated by zonal jets with a meridional scale on the order of but larger than the zonal wavelength of the initial wave, with vertical shear exponentially growing with time. The requirement of inertial stability selects modes with a jet centered on the equator. The next-largest contribution to the unstable perturbations after the zonal jet features meridional motion with the zonal wavelength of the initial wave and aligned in latitude with the axes of the jets. This zonally variable part of the perturbation does not persist in the simulations, and its rectification is likely responsible for mixing potential vorticity and biasing the response toward westward jets at the equator and eastward extra-equatorial jets. The initial vertical structure of the perturbations is inherited from the self-interaction of the initial wave. The linear theory agrees very well in most respects with the short-time behavior in the simulations, although growth rates and extra-equatorial jet positioning agree less well at higher *k* and *V*_{0}.

## 4. Inertial instability and equilibration of zonal jets

To explain the arrest of the growth of the zonal jets in the simulations, we again appeal to the requirement for inertial stability. At a certain amplitude, the westward equatorial jet will have curvature at the equator greater than unity (curvature greater than *β** in dimensional terms), implying that at that moment, which will be at a different time at each level, the zonal-mean absolute angular momentum,

with the overbar indicating zonal mean, will take its maximum values off the equator, and the flow will violate the inertial stability condition *y**m*_{y} < 0 in a latitude range surrounding the equator. The instability leading to jet growth will then compete with a process of inertial adjustment, leading in the long time limit to a state with homogenized *m* and PV. The negative of the gradient of *m*,

is the nondimensional zonal-mean absolute vorticity. An equivalent statement of the inertial stability condition is that *q* (or, in dimensional form, *q** ≡ *β***y* − _{y*}) will have the same sign as latitude.

Note that we are referring to the “parcel-type” zonally symmetric inertial instability associated with angular momentum shear confined to the region where angular momentum increases with latitude. While instability is implied if the jet curvature at the equator is greater than *β** and thus when there are local maxima in |*q**| near the equator, this is not the “wave-type” barotropic instability associated with a maximum in absolute vorticity. In the present case, although the criteria for the two instabilities are both satisfied, the simulations show evidence that inertial instability is particularly active. This is consistent with Stevens and Ciesielski (1986), who found that among all unstable modes of a near-equatorial zonal jet, the axisymmetric inertial instability modes with smallest vertical scale have the fastest growth rate.

Nonlinear studies of inertial instability (Hua et al. 1997; Griffiths 2003; Kloosterziel et al. 2007) show that equilibration, and therefore the ultimate width of the interval of zero potential vorticity, extends in latitude beyond the interval in which the flow was initially unstable. This behavior is observed in the simulations, as shown in Fig. 5, where *q* from the simulation with *k* = −8.4 and *V*_{0} = 0.11 is plotted at various times. Initially equal to the planetary contribution (because the zonal mean of the initial wave is zero), the jet formation leads eventually to the creation of anomalous *q* simultaneously on either side of the equator. The long time evolution of the flow leads to the homogenization of *q* across a latitude interval roughly twice the width of the interval over which the flow was ever unstable. Figure 6 shows the *q* profile at the moment of strongest inertial instability and the subsequent homogenized profile for different values of *k* and shows that the width of the interval of homogenized *q* in the equilibrated state is approximately independent of *k*. For the range of parameters used in our simulations, it is also insensitive to the amplitude of the initial wave (assuming enough energy is present to induce inertial instability). This can be seen from Fig. 7 where the equilibrated *q* is plotted for *k* = −8.4 and *V*_{0} equal to 0.075, 0.1, and 0.11.

Figure 8 shows both the maximum westward zonal velocity at the equator and *q* integrated over regions where it is anomalous as functions of time, showing the arrest of the growth of the westward jet coinciding in time with the appearance of anomalous PV. Subsequently, episodes of jet weakening coincide with removal of anomalous PV, and the steady-state jet amplitude is achieved once anomalous PV is exhausted and inertial stability restored. In some instances there is a slight time lag between the two curves, perhaps due to the destabilization of the initial wave continuing even after inertial instability has been triggered.

Also consistent with the notion of inertial instability being the mechanism responsible for halting the growth of the equatorial jet is the fact (evident from Fig. 8) that the peak and the equilibrated westward velocities are greater for runs with smaller |*k*|. This seems counterintuitive because the instability is due to the shear in the basic wave, which increases with |*k*|; however, the larger |*k*| is, the narrower the resulting equatorial jet will be, and the lower the amplitude at which its curvature at the equator will exceed *β** and anomalous PV will begin to appear (also evident from Fig. 8).

### a. Inertial instability of a narrow westward equatorial jet

Normal modes of inertial instability due to an oversharpened westward equatorial jet can be calculated by linearizing the zonally symmetric primitive equations about a depth-independent zonal flow with the shape of the zonally symmetric component of the most unstable even linear mode, calculated in section 3, and amplitude such that its curvature at the equator is greater than unity.

Still using the nondimensionalizations of section 3, let such a solution be *U*(*y*), and let perturbations to it be represented by the streamfunction *ϕ*′(*y*, *z*, *t*), such that *υ*′ = *ϕ*′* _{z}* and

*w*′ = −

*ϕ*′

*are the meridional and vertical components of the velocity perturbation, where*

_{y}*ϕ*′ satisfies the partial differential equation

Equation (4.3) has solutions of the form

where *μ _{ii}* and

*ν*are constants and the function

_{ii}*ϕ̂*′(

*y*) satisfies the ordinary differential equation

which is the equation found by Dunkerton (1981) for the case where *U* is linear in *y*.

Equation (4.5) is solved numerically. The most unstable solutions resemble the Taylor vortices familiar from inertial instability due to linear meridional shear in zonal velocity at the equator (Dunkerton 1981), but they are double-peaked in latitude and can be symmetric or antisymmetric about the equator (as opposed to the shifted latitude of symmetry in the linear shear problem). The smallest vertical scale modes have the highest growth rates, with the maximum growth rate tending to a constant value as the vertical mode number tends to infinity. The maximum growth rate in the limit of infinite vertical mode number increases with the curvature of the basic state. Figure 9 shows the streamfunction for the overturning components of the fastest-growing odd meridional mode for vertical mode *ν _{ii}* = 75 (here, “odd” refers to the symmetry of

*ϕ*′ about the equator; even modes are less interesting because they lead to meridional shear in zonal velocity at the equator, which is itself an inertially unstable condition). The chosen basic state (Fig. 9, left) is the most unstable linear barotropic instability mode for the case

*k*= −8.4 (cf. Fig. 3) at nondimensional amplitude 0.09 (0.3 m s

^{−1}in terms of the ocean parameters of the simulations). Notice that the linear instability is strongly confined to the latitude interval of instability (indicated on the graph).

Evidence of inertial instability in the simulations is the occurrence of and removal of anomalous potential vorticity already described (Fig. 8) and the transient formation of zonally symmetric structures in the *υ* and *w* fields (Fig. 10), which appear neither in the basic state nor in the linear barotropic instability modes of the initial wave (nor indeed in modes of barotropic instability of the equatorial jet). Figure 10 compares qualitatively with Fig. 9. The pattern of overturning cells is sensitive to the vertical diffusion in the model, becoming weaker, less persistent, and lacking in smaller-scale features with increasing diffusion. Presumably wavelike barotropic instability plays a more significant role in the equilibration of the flow when inertial instability is suppressed through vertical diffusion.

A remaining question is not answered here: What are the consequences of a vertically variable basic state? In particular, because at a given time the flow will be unstable at some depths and stable at others, will there be mixed propagating/growing modes that could potentially redistribute angular momentum vertically (in the long time-equilibrated state)? On the other hand, because the most unstable solutions are of smallest vertical scale, and because the instability is due to the horizontal gradient of angular momentum, the dynamics of the instability should be well represented locally by the solution linearized about a depth-independent state. Griffiths (2000) found that for the problem of shear inertial instability with depth-dependent shear, the most unstable modes are concentrated near the depths of maximum instability.

### b. Superrotation

In the simulations with higher values of |*k*| and still moderate values of *V*_{0}, the destabilization of the MRG wave leads to the development of superrotating (eastward) equatorial jets near the depths of zero horizontal velocity in the initial wave (see Figs. 1 and 2).

As mentioned in the introduction, atmospheric superrotation has been explained in terms of the convergence of planetary-wave momentum flux at the equator [the mechanism described in, e.g., Holton and Lindzen (1972), with respect to the QBO] or forcing by a prescribed zero-mean equatorial heat or momentum source (Suarez and Duffy 1992; Saravanan 1993; Williams 2003). Here, there are elements of both processes, with inertial instability providing a source of planetary-scale waves at the equator and the basic-state wave and its harmonics potentially serving as a zonally variable but zero zonal-mean forcing.

As is the case for the instability of the initial wave, the formation of the superrotating jet also scales with *k* and *V*_{0}. Figure 11 charts the occurrences of superrotating jets in the simulations. For extremely high vertical diffusion, the superrotating jets do not form.

How might the formation of the superrotating jets be explained? One suggestion is that the mixing leading to the meridional homogenization of potential vorticity extends in the vertical direction, outside of the depth range in which the flow was unstable. If so, we might expect horizontal PV mixing aligned in latitude with unstable modes (with a structure as in Fig. 9) but at depths where the flow is inertially stable. This organization of mixing weakens and widens an inertially unstable westward jet, but acting on the rest state, it would create PV steps on either side of the equator and a mixing barrier at the equator. If global angular momentum were conserved relative to the rest state, the mixing barrier would correspond to a superrotating jet [see Dunkerton and Scott (2008) for a discussion of jet patterns associated with PV staircases].

Whatever the mechanism for the creation of the superrotating jets, it must involve the violation of the assumption of eddies acting diffusively, the central assumption behind Hide’s theorem excluding the possibility of superrotation. This usually comes about through the nonlocal transfer of zonal-mean angular momentum by waves: vertical separation of regions of excitation and absorption of waves implies a vertical flux of angular momentum that depends on the intrinsic zonal phase speed of the waves (Bretherton and Garrett 1969). Here the modes of inertial instability are zonally symmetric (hence zero zonal phase speed), so no momentum flux should ensue. For that matter, in the absence of dissipation, angular momentum is materially conserved by zonally symmetric motions, so angular momentum (and potential vorticity) homogenization should be impossible. Therefore, the momentum transfer required for angular momentum homogenization probably involves the interaction of eddies with zonally symmetric inertial instability, which seems to be nonetheless the driving mechanism behind the adjustment.

### c. The effect of the nontraditional Coriolis force terms

Simulations were also performed with the Coriolis force terms due to the poleward component of planetary rotation included in the model. The nontraditional terms couple the zonal and vertical velocity components, and their inclusion has important consequences for the material conservation of angular momentum by zonally symmetric motions, especially for fluid parcels displaced vertically. The materially conserved quantity in the nontraditional system is

where *γ* ≡ 2Ω*/*N**. It is the spatial distribution of *m*_{γ} that is critical for inertial instability in the nontraditional system (Fruman and Shepherd 2008), and inertial adjustment should act in such a way as to redistribute *m*_{γ} to create a stable flow. While the instabilities in the present study are due to horizontal shear (zonal shear in *υ* in the basic state and meridional shear in *u* in the equatorial jet), as suggested in the previous section, the redistribution of *m*_{γ} and the corresponding zonal-mean potential vorticity

where *ρ* is density, may be in the meridional *and* vertical directions when the inertially unstable flow depends on *z*.

Indeed, when the simulations are run with and without the nontraditional Coriolis force terms included, there is a distinct vertical symmetry breaking in the nontraditional runs, particularly with respect to the superrotating jet. The upper jet is intensified and the lower is weakened (Fig. 12). This suggests that if the jets are created and maintained by vertically propagating disturbances, there is evidence for conservation of total angular momentum for disturbance wave packets (such as exists for zonally symmetric parcel displacements). If that were the case, upward-propagating disturbances would lose some of their relative angular momentum depositing potential because they gain planetary angular momentum, whereas the potential of downward-propagating disturbances would be increased. This would be an interesting question to address in a future study.

## 5. Summary

Using high-resolution numerical simulations, the instability of low-vertical mode, low-frequency, short zonal wavelength mixed Rossby–gravity (MRG) waves and the subsequent equilibration of the flow have been studied. The simulations reveal that westward-propagating MRG waves of sufficiently short wavelength and sufficiently large amplitude are unstable and lead to the growth of zonal jets alternating in direction with latitude, with a westward jet centered on the equator. The equatorial jet becomes inertially unstable at moderate amplitudes due to its very narrow meridional scale, and the ensuing inertial adjustment leads to the homogenization of absolute angular momentum and potential vorticity (PV) about the equator.

A linear model for the initial instability, based on the study of Gill (1974) for a barotropic Rossby wave, has been presented, and the argument is put forward that in the limit being considered, the instability of the intrinsically baroclinic MRG wave is dominated by horizontal shear. The most unstable mode of the vorticity equation for horizontally nondivergent perturbations is dominated by a zonally symmetric zonal velocity component, thus explaining the pattern of zonal jets alternating in direction with latitude observed in the simulations.

Inertial instability and PV mixing at the equator were invoked to explain several features of the simulations not explained by the linear model. First, the jets form such that a jet is centered on the equator (whereas the linear model predicts modes symmetric and antisymmetric about the equator) because zonal velocity shear at the equator is inertially unstable. Second, the fact that the equatorial jet observed in the simulations is always westward, even though the linear solutions are sign-reversible, is explained by noting that the unstable mode consists of a zonally symmetric term (the jets) and an infinity of terms of short zonal scale that dissipate in the long time limit. Because the largest terms of short zonal scale have motion across the jet axes, their dissipation will lead to potential vorticity mixing across the equator, which favors the formation of a westward jet over an eastward jet. Third, in the very short wavelength limit, the jets grow more slowly than predicted by the linear model, perhaps because the resulting narrower equatorial jet is more prone to being deflected off the equator, exciting highly baroclinic inertial instability cells that consume energy that would otherwise feed the barotropic instability of the wave.

The arrest of the growth of the zonal jets and the equilibration of the flow were attributed to the tendency of the westward equatorial jet toward inertial instability due to its narrowness, which implies anomalous angular momentum gradients when it reaches a certain amplitude. The appearance and removal of PV with the opposite sign as latitude, the well-known signature of inertial instability, coincide in time with, respectively, the peak strength and the weakening of the equatorial jet. The most striking feature of the equilibrated state is the homogenization of PV across the equator. The width of the latitude interval over which PV is homogenized seems to be independent of the wavelength and amplitude of the initial MRG wave (in the range of parameter space sampled).

For certain parameter values, the simulations exhibit small-vertical-scale superrotating equatorial jets straddling the depth of zero horizontal velocity in the initial wave. A qualitative interpretation was proposed in terms of the mixing of angular momentum in both latitude and depth to most efficiently stabilize the inertially unstable equatorial jet. A clue to the importance of vertical momentum redistribution is the vertical symmetry breaking due to the nontraditional Coriolis force terms on the structure of the superrotating jets.

Figure 13 summarizes the transitions that occur in the simulations. Wave motion dominated by meridional motion varying strongly with longitude is converted through barotropic instability to largely barotropic zonal jets of short meridional scale alternating in direction with latitude. The jets amplify until they become inertially unstable when the curvature of the westward jet at the equator becomes large enough that the angular momentum increases away from the equator, leading to the widening and weakening of the jet and homogenization of potential vorticity. Finally, in certain cases, superrotating jets form outside the depths of inertial instability, implying potential vorticity mixing on either side of the equator.

## Acknowledgments

Simulations described in this work were performed at l’Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) under Grant 81777 and at the Earth Simulator Center under IFREMER-CNRS-ES MOU. The authors thank Sylvie Le Gentil for her work adapting ROMS for this project and Marc d’Orgeville and Claire Ménesguen for useful discussions during the preparation of the manuscript.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

* Current affiliation: Institut für Atmosphäre und Umwelt, Goethe-Universität Frankfurt, Frankfurt am Main, Germany.

*Corresponding author address:* B. L. Hua, Laboratoire de Physique des Océans, IFREMER Centre de Brest, 29280 Plouzané, France. Email: lien@ifremer.fr