## Abstract

The linear instability of divergent perturbations that evolve on a cos^{2} mean steady zonal jet embedded in a zonal channel on the *β* plane and on a rotating sphere is studied for zonally propagating wavelike perturbations of the shallow-water equations. The complex phase speeds result from the imposition of the no-flow boundary conditions at the channel walls on the numerical solutions of the linear differential equations for the wave latitude-dependent amplitude. In addition, the same numerical method is applied to the traditional problem of linear instability of nondivergent perturbations on the *β* plane where results reaffirm the classical, analytically derived, features. For these nondivergent perturbations, the present study shows that the growth rate increases monotonically with the jet maximal speed and that the classical result of a local maximum at some finite westward-directed speed results from scaling the growth rates on the jet’s speed. In contrast to nondivergent perturbations, divergent perturbations on the *β* plane have no short-wave cutoff, and so the nondivergent solution does not provide an estimate for the divergent solution, even when the ocean is 1000 km deep (i.e., when the speed of gravity waves exceeds 10 Mach). For realistic values of the ocean depth, the growth rates of divergent perturbations are smaller than those of nondivergent perturbations, but with the increase in the ocean depth they become larger than those of nondivergent perturbations. For both perturbations, a slight asymmetry exists between eastward- and westward-flowing jets. The growth rates of divergent perturbations on a sphere are similar to those on the *β* plane for the same values of the model parameters, but the asymmetry between eastward and westward jets is more conspicuous on a sphere. The value of *g*′*H*′ (*g*′ is the reduced gravity; *H*′ is the equivalent mean layer thickness), which is filtered out in nondivergent theory, determines for divergent perturbations the relative magnitude of zonal velocity, meridional velocity, and height but has little effect on the growth rates.

## 1. I. Introduction

The classical finding made by Kuo (1949, 1973) regarding the linear instability that results from the shear of a mean zonal current embedded in a channel on the *β* plane is interpreted in terms of the mean current’s absolute vorticity. The linear perturbations that develop on the mean zonal current are assumed to be nondivergent, which enables an analysis of the instability in terms of the perturbation streamfunction. Subsequent theories that were developed in the decades that followed Kuo’s pioneering work have employed both a similar setup and considerations, namely, that the perturbations were nondivergent and the instabilities were interpreted in terms of the conservation of the mean current absolute vorticity. An example for such an extension of Kuo’s paradigm is the theory of the instability of a uniform shear flow on the *f* plane given in Gill (1982) where the zonal channel is replaced by relative vorticity discontinuities but the *f*-plane perturbations are nondivergent (a streamfunction exists).

Our aim in this work is to reexamine Kuo’s theory of nondivergent flows on the *β* plane and to extend it numerically to divergent flows both on the *β* plane and on a rotating sphere in the framework of the shallow-water equations (SWEs). A brief summary of the assumptions and results of the relevant elements of Kuo’s classic theory is now in order.

### a. The *β*-plane barotropic instability theory for nondivergent perturbations on a *cos*^{2} jet

^{2}

In his pioneering original work Kuo (1949) solved the nonlinear eigenvalue equation for the zonal phase speed *C* of nondivergent wavelike perturbations on a mean zonal current *U*(*y*) on the *β* plane. The nonlinear eigenvalue equation for *C* is posed in terms of the linear differential equation for the meridional structure of the perturbation streamfunction *ψ*,

where *k* is the wave zonal wavenumber and subscript “*y*” denotes differentiation.

Equation (1.1) expresses the conservation of absolute (relative plus planetary) vorticity, *Q* = ∇^{2}*ψ* + *f*, of the nondivergent perturbations and such a *ψ* always exists in nondivergent flows. By carefully analyzing the structure of the regular general solution of Eq. (1.1) emanating from the singular point *y _{c}* [defined by

*U*(

*y*) −

_{c}*C*= 0] and imposing the wall boundary conditions where the meridional velocity component,

*υ*=

*ikψ*, vanishes, Kuo (1949) deduced bounds on the growth rate

*kC*and some general properties of the structure of

_{i}*ψ*.

A detailed calculation of the phase speeds and growth rates was carried out for the cos^{2} mean flow |P[i.e., *U*(*y*) = *U*_{max}cos^{2}[*y*/(2*d*)], −*π* < *y*/*d* < *π*; *U*_{max} < 0 for westward flows|P] in the sequel paper, Kuo (1973), whose main findings (Fig. 6 of Kuo 1973) are reproduced here in Fig. 1. Figure 1a shows the real part of the phase speed *C _{r}* (scaled on the maximal speed of the mean current,

*U*

_{max}) as a function of the wavelength

*L*divided by

*d*(where the channel width is 2

*πd*). The coalescence of the two real modes at

*L*/

*d*> 9 and at

*L*/

*d*< 2 to a single mode in the range 2 <

*L*/

*d*< 9 is a clear sign of the existence of a complex root of the eigenvalue equation in this range. Figure 1b shows contours of the growth rate

*kC*on a plane marked by the abscissa

_{i}*b*=

*βd*

^{2}/

*U*

_{max}and the ordinate

*a*=

*kd*= 2

*πd*/

*L*. The main features of the contours are as follows.

A short-wave cutoff is observed at

*kd*= 3^{1/2}/2 for both positive and negative*b*values, consistent with the vorticity considerations raised in Kuo (1949) and Pedlosky (1979).A strong asymmetry between eastward-directed and westward-directed mean currents is evident in these results. For westward-directed mean flows no longwave cutoff exists and the

*kC*= 0 contour is located right on the_{i}*kd*= 0 axis. In contrast, for eastward-directed mean flows the longwave segment of the*kC*= 0 contour is fairly curved (parabolic)._{i}A cutoff at |

*b*| = 0.5 exists, which is uniform for all values of*kd*.The last, and most counterintuitive, point that should be made with regard to the Kuo instability contours is their apparent dependence on the value of the maximal speed

*U*_{max}for westward-directed mean flows (*U*_{max}< 0). The maximal growth rate contour in Fig. 1,*kC*= 0.25, is centered near (_{i}*kd*= 0.45,*b*= −0.15), and the growth rate*decreases*when the value of*b*is increased from −0.15 (for most values of*kd*). Since*b*=*βd*^{2}/*U*_{max}, one faces the strange result that (for fixed*βd*^{2}and for negative*U*_{max}) as the absolute value of*U*_{max }*increases*(and with it the shear of the mean current) the growth rate*decreases*. This result seems to contradict the notion that barotropic instability is generated by the meridional shear of the mean flow since an increase in the shear reduces the growth rate. This odd result does not occur for eastward-directed mean currents for which no closed contours are present, so any increase in*U*_{max}(i.e., any decrease in*b*) brings about an increase of the growth rate,*kC*._{i}

Our aim in the present work is to examine which of these findings prevails for divergent perturbations both in planar and spherical geometries. Accordingly, we formulate our problem as the linear instability analysis of a meridionally varying zonal current embedded in a zonal channel on a rotating sphere. The same numerical search for the instabilities of divergent perturbations will also be done on the *β* plane and the resulting growth rates will be compared to one another and to the Kuo instabilities of nondivergent perturbations on the *β* plane. The paper is organized as follows: The model equations on the *β* plane for both nondivergent and divergent flows are developed, and the resulting instabilities are calculated in section 2. In section 3 we develop the equations on the sphere, and the numerical method used for solving the complex eigenvalue problem is described in section 4. The results of the instability calculations for the cos^{2} mean current are described in section 5. We end with concluding remarks in section 6.

## 2. The eigenvalue problems on the *β* plane and their solutions

Our goal in this section is to solve the barotropic instability problem on the *β* plane as a prelude to its solution on the sphere. We will do so for both nondivergent (where our numerical results can be compared to Kuo’s analysis) and divergent perturbations that will be compared to the instability of divergent perturbations on a sphere in section 5.

In the presence of rotation the vectorial form of the shallow-water equations, which describe the dynamics of a hydrostatic shallow layer of fluid, is

where **V** is the horizontal (i.e., tangential to the earth’s surface) velocity vector, *f* = 2Ωsin*ϕ* is the Coriolis frequency (Ω is the earth’s rotation frequency and *ϕ* is the latitude), **k** is a unit vector pointing positively upward, *g* is the gravitational acceleration, *h* is the height (thickness) of the fluid layer, and time derivatives on the left-hand side (lhs) are Eulerian derivatives; that is,

On the *β* plane, system (2.1) can be linearized about a steady mean state where the zonal velocity *U _{s}*(

*y*) is prescribed as a function of latitude. In the divergent case the mean height

*h*(

_{s}*y*), corresponding to this mean zonal velocity, is given by the geostrophic balance, but for nondivergent perturbations the height field is irrelevant to the perturbation dynamics and is undetermined. Following the Kuo cos

^{2}jet structure we prescribe the mean steady velocity by

where 2*π**d* is the width of the mean current, which is also the width of the channel that bounds the meridional extent of the mean jet and its perturbations (*U*_{max} < 0 for westward jets).

On the *β* plane the Coriolis parameter *f* (*y*) in system (2.1) is given by *f* (*y*) = *f _{c}* +

*βy*, where

*f*= 2Ωsin(

_{c}*ϕ*) (

_{c}*ϕ*is the latitude of the channel center) and

_{c}*β*= 2Ωcos(

*ϕ*)/

_{c}*R*(

*R*earth’s radius). Since on the

*β*plane

*R*appears in the

*β*term system (2.1) is nondimensionalized using R as the length scale and (2Ω)

^{−1}as the time scale, which implies that 2Ω

*R*is the velocity scale. The fluid height

*h*is scaled by the typical (e.g., mean, maximal/minimal) height of the ocean or atmosphere,

*H*. The vectorial form (2.1) is transformed to Cartesian coordinates where

*x*and

*y*point in the local east and north directions, respectively, and

*u*and

*υ*are the corresponding velocity components in these directions. When the (

*u*,

*υ*,

*h*) variables in (2.1) are linearized about a mean steady state

*U*(

_{s}*y*),

*h*(

_{s}*y*) by letting (

*u*′,

*υ*′,

*η*) = (

*u*−

*U*,

_{s}*υ*,

*h*−

*h*) one obtains the linear nondimensional set (upon dropping the primes from

_{s}*u*′ and

*υ*′)

where *f** = sin(*ϕ _{c}*) + cos(

*ϕ*)

_{c}*y*is the nondimensional Coriolis frequency and

*α*=

*gH*/(2Ω

*R*)

^{2}is the nondimensional speed of gravity waves that plays the role of nondimensional gravity in (2.3) (1/

*α*is also known as the Lamb number or the rotational Froude number). As mentioned earlier, the form of

*h*(

_{s}*y*) is different for the nondivergent (Kuo 1973) and divergent perturbations and these two cases are examined separately next. For nondivergent perturbations,

*h*(

_{s}*y*) is set equal to 1.0 in the continuity equation, which implies that

*α*

^{−1}has to equal 0.0 in the geostrophic relation,

that determines *h _{s}*(

*y*) in the divergent case [recall that

*U*(

_{s}*y*) ≠ 0].

### a. Nondivergent perturbations—Revisiting Kuo’s analysis

The standard way of overcoming the subtlety associated with the irregularity of the rhs of the momentum equations in the nondivergent limit is to cross-differentiate these equations so as to eliminate the pressure gradient terms altogether. The nondivergence condition *u _{x}* +

*υ*= 0 implies the existence of a streamfunction

_{y}*ψ*related to the velocity via

**V**=

**k**×

**∇**

*ψ*that satisfies the equation

where the dependence of *ψ*(*x*, *y*, *t*) on *x* and *t* is assumed to be that of a zonally propagating wave −*e ^{ik}*

^{(x−Ct)}. Equation (2.4) is no other than Eq. (1.1) nondimensionalized on the present study’s scales defined earlier. The boundary conditions of vanishing normal flow

*υ*∼

*ψ*at the channel walls are

*ψ*(

*y*= ±

*δϕ*) = 0. For easy comparison between the nondimensional results of the present study (

*y*is scaled on

*R*,

*U*

_{max}and

*C*on 2Ω

*R*,

*β*on 2Ω/

*R*, and

*k*on 2

*π*/

*R*) and those of Kuo (1973) shown in Fig. 1b, one has to identify the Kuo parameters

*a*and

*b*as

where *δϕ* is half-channel width in nondimensional units (in radians).

The details of the numerical method used for calculating the growth rate contours *kC _{i}*(

*b*,

*a*) are given section 4. To enable a comparison with Kuo’s results (Fig. 1b) our results, shown in Fig. 2, were carried out for

*ϕ*= 30° and

_{c}*δϕ*= 0.156 (which correspond to Kuo’s values of

*β*and

*d*). In Fig. 2a these contours are shown on the same (

*b*,

*kd*) plane and in the −0.5 <

*b*< 0.5 range as in Kuo’s work. The growth rate contours,

*kC*(

_{i}*U*

_{max},

*k*), shown in Fig. 2b, conform to the style adopted in this work where

*k*and

*U*

_{max}are the primary model parameters used as the axes in plotting the growth rate contours for fixed values of the other parameters,

*ϕ*and

_{c}*δϕ*.

Our results, shown in Figs. 2a, confirm Kuo’s analytical considerations: There is a uniform (i.e., for all values of −1/2 ≤ *b* ≤ 1/2) cutoff at *a* = *kd* = 3^{1/2}/2 = 0.866 [i.e., *k* = *πa*/(*δϕ*) = 17.44] and there are no instabilities for either *b* > 1/2 or for *b* ≤ −1/2. In addition, the asymmetry between *b* > 0 and *b* < 0 is evident in our results, where the *kC _{i}* = 0.01 contour occupies a much larger region at negative

*b*than for positive

*b*, as in Kuo’s results. However, the main difference between the results of Fig. 2 and those of Kuo relates to the physically unexplainable local maximum at

*b*= −0.15 found by Kuo. In our calculations the maximum is located at

*b*= 0, which implies that, when |

*U*

_{max}| is increased, so do the instability contours, in agreement with the notion that the shear of the mean flow provides the energy source for the growth of the perturbations. This point is further discussed in section 6.

### b. Instabilities of divergent perturbations on the *β* plane

When a zonally propagating wave form, *e ^{ik}*

^{(x−Ct)}is assumed for the (

*x*,

*t*) dependence of

*u*,

*υ*, and

*η*in Eq. (2.3), the

*x*-momentum equation yields

*u*as a linear combination of

*υ*and

*η*, which can be employed to eliminate

*u*from the other two equations. Defining

*V*=

*iυ*/

*k*one then gets the following equations for the meridional variation of the perturbation’ amplitudes:

The results of the search for the instabilities of the divergent case [based on numerical integration of system (2.6); see more details in section 4] are shown in Fig. 3 for three values of *α* that span four orders of magnitude. The instabilities clearly increase with *α*, but even for *α* = 10 (*H* = 1000 km) they bear little resemblance to those of the nondivergent case (Fig. 2b). The most obvious difference is the disappearance of the short-wave cutoff that typifies Fig. 2b. At low *α* values, more than one local maximum in the (*U*_{max}, *k*) plane is encountered.

## 3. The eigenvalue problem on a sphere

In (geographical) spherical coordinates the longitude *λ* points positively eastward and the latitude *ϕ* points positively northward. The corresponding velocity components along these directions are *u* and *υ*, respectively, and the unit vectors in these directions are **i** and **j**, respectively. On a sphere of radius *R* one has

so system (2.1) is written in scalar form and in the (*λ*, *ϕ*) spherical coordinates as

System (3.1) is nondimensionalized using the same scales as in section 2: *R* is the length scale, 1/(2Ω) is the time scale, and 2Ω*R* is the velocity scale. The fluid height *h* is scaled on the typical (e.g., mean, maximal/minimal) height of the ocean or atmosphere, *H*. In nondimensional form system (3.1) is given by

where, as was the case on the *β* plane, *α* = *gH*/(2Ω*R*)^{2}.

As in Kuo’s theory, since a mean zonal flow is assumed, the linear stability analysis of system (3.2) begins by first linearizing it about such mean current, *U _{s}*(

*ϕ*), and its corresponding (via geostrophy) height field. However, the quadratic nonlinearity of Eq. (3.2a) is best handled by replacing the zonal velocity

*u*in system (3.2) by the angular momentum component in the direction parallel to the axis of rotation. In dimensional units this angular momentum component is given by

*M*=

*R*cos

*ϕ*(2Ω

*R*cos

*ϕ*+

*u*) and system (3.2) is transformed to

where *M* = cos*ϕ*(1/2cos*ϕ* + *u*) is the nondimensional angular momentum, so *u* in the advective part of the total time derivative on the lhs of Eqs. (3.3a), (3.3b) is given by *u* = *M*/cos*ϕ* − cos*ϕ*/2.

It is evident that, while the rhs of Eq. (3.2a) is nonlinear and includes both rotation and geometrical terms, its counterpart in Eq. (3.3a) is linear and has the same form as in planar geometry. The nonlinearities on the rhs of Eqs. (3.2b), (3.2c) and of the time derivatives in both equations are identical to those in Eqs. (3.3b), (3.3c).

The steady mean state is defined by a latitude-dependent height field *h _{s}*(

*ϕ*) and vanishing meridional velocity. Following Eq. (3.3b) the angular momentum in the mean state is

from which the mean zonal current *U _{s}*(

*ϕ*) is calculated by

Note that the negative root in Eq. (3.5) corresponds to unphysically large westward (negative) velocities, while the positive root yields the (nondimensional, geostrophic) zonal velocity

in the limit of small *α* (or small height gradient).

System (3.3) is now linearized about the mean state by considering only first-order terms in *m* = *M* − *M _{s}*(

*ϕ*),

*υ*,

*η*=

*h*−

*h*(

_{s}*ϕ*). The linearized variables are then assumed to represent zonally propagating waves with wavenumber

*k*and phase speed

*C*so that in all variables the latitude-dependent amplitude is multiplied by

*e*

^{ik}^{(λ−Ct)}. In Eq. (3.3a) the

*ϕ*derivative occurs only in the term

*υ*

*dM*/

_{s}*d*

*ϕ*, so its linearized version provides

*m*as a linear combination of

*υ*and

*η*. This linear relation is used to eliminate

*m*in Eqs. (3.3b) and (3.3c), so the latitude-dependent amplitude of

*V*=

*iυ*/

*k*and

*η*are solutions of the second-order ordinary differential set

(Note that up to this point the spherical geometry has not been compromised whatsoever.)

The boundary conditions, required in order to close the differential set (3.6), are that *V* (i.e., *υ*) should vanish at the two latitudes where zonal walls bound the channel (again, this is the same setup as in Kuo’s original theory). We thus set the center of the zonal channel at latitude *ϕ _{c}* and its walls at

*ϕ*−

_{c}*δϕ*and

*ϕ*+

_{c}*δϕ*, where 2

*δϕ*is the channel width. The dispersion relation,

*C*(

*k*), for given

*α*and

*h*(

_{s}*ϕ*) [and the corresponding

*U*(

_{s}*ϕ*),

*M*(

_{s}*ϕ*)] is obtained when the no-flow boundary conditions at the channel walls are imposed on numerical solutions of the linear differential set (3.6). Finally, we note that although system (3.6) and its accompanying,

*V*= 0, boundary conditions are linear, the eigenvalue problem for

*C*(

*k*) is nonlinear since

*C*appears nonlinearly there.

## 4. Numerical method for solving the eigenvalue problem and the mean flow

For fixed values of *α*, *k*, *ϕ _{c}*,

*δϕ*, and

*U*

_{max}[which determines

*M*(

_{s}*ϕ*),

*U*(

_{s}*ϕ*) explicitly and

*h*(

_{s}*ϕ*) by integrating the inverse of Eq. (3.5)], system (3.6) is integrated from one wall (e.g., −

*δϕ*), starting with the initial conditions

*V*(

*ϕ*−

_{c}*δϕ*) = 0;

*η*(

*ϕ*−

_{c}*δϕ*) = 1 (the latter condition is an arbitrary normalization of the linear system) to the other wall (e.g.,

*ϕ*+

_{c}*δϕ*). The value of

*C*is then varied by employing a two-dimensional nonlinear zero finder for

*C*and

_{r}*C*, until the

_{i}*V*= 0 boundary condition at the second wall is satisfied. System (3.6) is integrated with a fifth-order Runge–Kutta (RK) method with 10

^{−8}tolerance to ensure sufficient accuracy of the calculated eigenvalues and eigenfunctions. Infinitely many eigenvalues and associated eigenfunctions exist, including real ones, and in this study we looked only for the eigenvalue with the largest growth, that is,

*C*with the largest imaginary part. This method was successfully employed in other instability studies (e.g., Paldor and Nof 1990; Paldor and Ghil 1997; Haza et al. 2004).

An additional numerical issue, specific to the spherical problem, is the determination of the mean state. Equations (3.4) and (3.5) provide straightforward expressions for calculating *M _{s}*(

*ϕ*) and

*U*(

_{s}*ϕ*) by differentiating a given

*h*(

_{s}*ϕ*) distribution. However, we strive to extend Kuo’s original theory to a sphere and therefore wish to employ the same mean current

*U*(

_{s}*ϕ*) used in his

*β*-plane theory. Since the mean current in Kuo’s theory was a cos

^{2}jet structure across the channel, we also chose the mean current to be

Given this choice of the mean current *U _{s}*(

*ϕ*), the calculation of

*M*(

_{s}*ϕ*) follows directly from the definition

*M*(

_{s}*ϕ*) =

*R*cos

*ϕ*[2Ω

*R*cos

*ϕ*+

*U*(

_{s}*ϕ*)], but the corresponding

*h*(

_{s}*ϕ*) can only be calculated by integrating numerically the inverse of Eq. (3.5); that is,

The *h _{s}(ϕ_{c}* −

*δϕ*) = 1 assumption made in Eq. (4.2) does not restrict the generality of our results since it merely changes the value of the height scale

*H*. The numerical integration in Eq. (4.2) to determine the mean height distribution preceded the numerical integration of system (3.6) and was done with the same numerical method (fifth-order RK), but the integration tolerance was 10

^{−5}since higher tolerance did not yield more accurate eigenvalues. For eastward-flowing jets (i.e., for

*U*

_{max}> 0) the mean height field,

*h*(

_{s}*ϕ*), is increasing with

*ϕ*while for westward-flowing mean jets (i.e., for

*U*

_{max}< 0) the height

*h*(

_{s}*ϕ*) is a decreasing function. To avoid circumstances in which

*h*(

_{s}*ϕ*) vanishes at some point inside the channel (for small

*α*) we chose to reverse the direction of integration of Eq. (4.2) (i.e., integrate it from the north wall at

*ϕ*+

_{c}*δϕ*to the south wall

*ϕ*−

_{c}*δϕ*) for westward directed mean currents, which merely changes the scale of the mean height

*h*(

_{s}*ϕ*+

_{c}*δϕ*) = 1 instead of

*h*(

_{s}*ϕ*−

_{c}*δϕ*) = 1.

In the present problem contours of the growth rate, *kC _{i}*, contain troughs in the (U

_{max},

*k*) plane located between local maxima, where

*kC*= 0 (see Figs. 3 and 4). Thus, a continuation method in which the numerical search for the eigenvalue (

_{i}*C*,

_{r}*C*) at any (

_{i}*U*

_{max},

*k*) point begins with the (

*C*,

_{r}*C*) values found at the previous (

_{i}*U*

_{max},

*k*) point is bound to fail near these zero contours. In such cases, when the continuation method failed to yield a nonzero growth rate at some (

*U*

_{max},

*k*) point, we initiated a numerically expensive search for

*C*,

_{r}*C*starting from nine different initial guesses of

_{i}*C*and

_{r}*C*; in each of these initial guesses

_{i}*C*assumed the values of 0.75

_{r}*U*

_{max}, 0.1

*U*

_{max}, and 0.01

*U*

_{max}, while

*C*assumed the values of 0.005, 0.05, and 0.5. The growth rates displayed in Figs. 3 and 4 are the maximal

_{i}*kC*from all searches at that (

_{i}*U*

_{max},

*k*).

## 5. Instabilities of divergent perturbations on a sphere

Our nondimensional, divergent, spherical, model of a cos^{2} zonal jet has five free parameters: *α* [=*gH*/(2Ω*R*)^{2}], *δϕ* (half the channel width), *ϕ _{c}* (the channel’s central latitude),

*k*(the zonal wavenumber), and

*U*

_{max}(the amplitude of the jet’s mean current). Although

*g*= 9.8 m s

^{−2}and

*H*= 4 or 10 km (in the ocean or the atmosphere, respectively), in equivalent-barotropic cases

*g*and

*H*in

*α*should be replaced by the reduced gravity,

*g*′ =

*g*Δ

*ρ*/

*ρ*, and reduced height,

*H*′ =

*H*

_{1}

*H*

_{2}/ (

*H*

_{1}+

*H*

_{2}) so that

*α*is a free parameter in our problem, whose value is changed in each particular application. We will first present our numerical results on contours of instability growth rates

*kC*(

_{i}*U*

_{max},

*k*) and then show the results regarding the associated eigenfunctions. Although only integer values of

*k*are physically acceptable on a sphere towing to the 2

*π*periodicity in

*λ*, the mathematical eigenvalue problem can be solved numerically for noninteger values of

*k*. Thus, as in Fig. 1, our contours are plotted as continuous curves, which are easier to grasp than a set of discrete points.

Figure 4 shows typical growth rate contours for *α* = 0.1 (i.e., *gH* of 10^{5} m^{2} s^{−2}), *ϕ _{c}* = 0.75 rad (∼43°) for −1 <

*U*

_{max}< 1 (in jumps of 0.005) and for 0 <

*k*< 14 (calculated in jumps of 0.1). For the three values of

*δϕ*= 0.1, 0.2, and 0.3 rad in Figs. 4a–c, respectively, the channel widths 2

*δϕ*are about 1270 km (∼11.5°), 2530 km (∼22.9°) and 3800 km (∼34.4°). It is obvious that in contrast to nondivergent perturbations on the

*β*plane (Figs. 1b, 2), divergent perturbations on the sphere have no short-wave [at least up to

*k*= 14, corresponding to a dimensional zonal wavelength [2

*π*

*R*cos(

*ϕ*)/

_{c}*k*∼ 2000 km] cutoff for either

*U*

_{max}> 0 or

*U*

_{max}< 0. A second important difference between the contours of Figs. 1b and those of Fig. 4 is that in the latter (as well as in our results of Fig. 2) there is no finite value of

*U*

_{max}(negative or positive) below which the growth rates vanish, and the contours terminate only on the

*U*

_{max}= 0 ordinate. A third difference is in the maxima of the growth rate contours in all panels of Fig. 4, which are local and occur mostly at large |

*U*

_{max}| values in contrast to the single, global, maximum on the

*β*plane at

*U*

_{max}< 0, which was calculated for |

*U*

_{max}| = 0.01.

On the other hand, the asymmetry between positive and negative values of *U*_{max} that is so clearly evident in Kuo’s nondivergent results (Fig. 1b), where a longwave cutoff exists only for *U*_{max} > 0 and the maximal contour occurs at *U*_{max} < 0, has a much more complex fingerprints in Fig. 4. At some values of *δϕ*, *k*, and |*U*_{max}| eastward directed jets are more unstable than westward directed ones, while at other values of these parameters the opposite occurs. The only general feature that can be drawn from the plots of the *kC _{i}* contours in Fig. 4 is that in the long-wave limit (small

*k*) westward directed jets are more unstable than eastward directed ones. The asymmetries between eastward directed jets and westward directed ones are manifested in the value of the maximal growth rate contour and in the area encompassed by the same

*kC*contours.

_{i}A comparison between the growth rates in Fig. 4 and those for other values of *α*, *ϕ _{c}* shows that the growth rate contours,

*kC*, are slightly affected by even large changes in the value of

_{i}*α*The maximal contours move closer to the

*U*

_{max}= 0 ordinate and for the same values of

*U*

_{max}and

*k*the growth rates at

*α*= 10

^{−4},

*ϕ*= 1.0 are about half those at

_{c}*α*= 0.1,

*ϕ*= 0.75, but they occur at higher values of

_{c}*k*. A possible explanation for the fact that a fairly large change in

*α*(three orders of magnitude) has a small effect on the instability can be sought by realizing that the instabilities are affected by the ratio between the radius of deformation and the channel width; that is, their dependence on

*α*is via

*α*

^{1/2}/[sin(

*ϕ*)

_{c}*δ*

*ϕ*]. A change in the channel central latitude to

*ϕ*= 0.5 rad produced very similar contours to those shown in Fig. 4 (results not shown) but with slightly smaller growth rates and with more spread out contours (for

_{c}*U*

_{max}> 0 the maximal growth rate is slightly lower and occurs at somewhat higher wavenumber).

In contrast to the growth rates that are weakly dependent on *α*, the ratio between the maximal values of *V*(*ϕ*) and that of *η*(*ϕ*) scales like *α* for our normalization of *η* = 1 at one of the channel walls. A sample calculation of the *V*(*ϕ*), *η*(*ϕ*) eigenfunctions for *α* = 10^{−4}; *δϕ* = 0.2; *ϕ _{c}* = 30°;

*U*

_{max}= +0.2;

*k*= 5;

*C*= (4.431, 1.734

*i*) × 10

^{−2}is shown in Fig. 5, which demonstrates that our numerical solution yields associated

*V*(

*ϕ*) and

*η*(

*ϕ*) eigenfunctions that are continuous inside the channel and that V = 0 at the two walls. The amplitude of

*η*(∼100) is four orders of magnitude (i.e., ∼

*α*

^{−1}) larger than that of

*V*(∼0.02). In comparison, for the same values of

*δϕ*= 0.2,

*ϕ*= 30°,

_{c}*U*

_{max}= +0.2,

*k*= 5 but for

*α*= 0.1 and

*C*= (6.151, 4.386

*i*) × 10

^{−2}, the results shown in Fig. 6 demonstrate that in this high

*α*case the amplitude of

*η*(∼3.3) is about one order of magnitude (again ∼

*α*

^{−1}) larger than that of the

*V*(∼0.5). The way

*α*appears in Eq. (3.6b) is consistent with our finding max|

*V*(

*ϕ*)| ∼

*α*max|

*η*(

*ϕ*)|.

## 6. Discussion and concluding remarks

The conspicuous asymmetry between the growth rates on eastward and westward mean flows that dominates the Kuo results (Fig. 1b) exists, but is much weaker in our results (Fig. 2). In particular, the local maximum of the growth rate at *b* = −0.15 in Kuo’s calculations (Fig. 1b) occurs at *b* = 0 in our results (Fig. 2b), in agreement with the intuitive notion that unstable modes grow due to the shear of the mean current (i.e., with the decrease in |*b*| ∼ 1/|*U*_{max}|). The maximum growth rate at finite *U*_{max} (for fixed *β* and *d*) in Fig. 1b results only from the Kuo scaling of *kC _{i}* on

*d*/

*U*

_{max}, which implies that

*kC*can increase with the increase in

_{i}*U*

_{max}when moving to the right of

*b*= −0.15, while

*kC*/(

_{i}*d*/

*U*

_{max}) decreases. This point of misleading scaling by Kuo is evident in Fig. 7, where contours of

*kC*(

_{i}*δϕ*/

*π*)/

*U*

_{max}are plotted using the same

*kC*values as in Fig. 2. When the

_{i}*kC*contours of Fig. 2a are compared with the

_{i}*kC*(

_{i}*δϕ*/

*π*)/

*U*

_{max}contours of Fig. 7, it becomes clear that the apparent maximum at

*b*< 0 in the latter is due only to the scaling of

*kC*in Kuo’s work on U

_{i}_{max}(

*δϕ*/

*π*is the present study’s counterpart of the channel width,

*d*=

*d*′/

*π*, in Kuo). Note that the only difference between Kuo’s analytical results and our numerical results is that we find a maximal contour of 0.23 while Kuo’s maximal contour is 0.25.

Our most surprising finding of the instability of divergent perturbations (SWE) on the *β* plane is that for short waves the classical nondivergent model that yields Kuo’s analytical results does not provide a qualitative approximation to the divergent case even for *α* = 10, which corresponds to an ocean depth of 1000 km. The heuristic rationale behind the rigid-lid approximation, which justifies the nondivergent assumption, is that in a deep ocean the motion of the free surface does not affect the dynamics. Although this argument is correct as a mathematical statement, the relevant depth where divergent solutions are accurately approximated by nondivergent ones is too large for the earth (or for the SWE to apply). Moreover, as is evident from Fig. 3, the growth rates increase with *α* and so any additional increase in the value of *α* will further increase the difference between the maximal growth rates in the divergent and nondivergent cases (cf. Fig. 2b with the top panel of Fig. 3).

As a final note, a comparison between the results of section 2b and section 5 shows that the instabilities on a sphere are similar to those on the *β* plane (cf. Figs. 3 and 4), but the asymmetry between eastward-flowing and westward-flowing jets is more pronounced on a sphere. An increase in the channel width lowers the growth rates, as is evident upon comparing the three panels of Fig. 4.

## Acknowledgments

This work was supported by Grant 579/05 of the Israel Science Foundation to HU. We are grateful to two anonymous reviewers for their comments. Numerous discussions held with A. Sigalov of HU helped us to distill some of the subtle points of the paper.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Nathan Paldor, Department of Atmospheric Sciences, Institute of Earth Sciences, Hebrew University of Jerusalem, Jerusalem 91904, Israel. Email: nathan.paldor@huji.ac.il