The linear instability of divergent perturbations that evolve on a cos2 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 cos2 jet
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 yc [defined by U(yc) − C = 0] and imposing the wall boundary conditions where the meridional velocity component, υ = ikψ, vanishes, Kuo (1949) deduced bounds on the growth rate kCi and some general properties of the structure of ψ.
A detailed calculation of the phase speeds and growth rates was carried out for the cos2 mean flow |P[i.e., U(y) = Umaxcos2[y/(2d)], −π < y/d < π; Umax < 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 Cr (scaled on the maximal speed of the mean current, Umax) 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 kCi on a plane marked by the abscissa b = βd2/Umax and the ordinate a = kd = 2πd/L. The main features of the contours are as follows.
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 kCi = 0 contour is located right on the kd = 0 axis. In contrast, for eastward-directed mean flows the longwave segment of the kCi = 0 contour is fairly curved (parabolic).
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 Umax for westward-directed mean flows (Umax < 0). The maximal growth rate contour in Fig. 1, kCi = 0.25, is centered near (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 = βd2/Umax, one faces the strange result that (for fixed βd2 and for negative Umax) as the absolute value of Umax 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 Umax (i.e., any decrease in b) brings about an increase of the growth rate, kCi.
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 cos2 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 Us(y) is prescribed as a function of latitude. In the divergent case the mean height hs(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 cos2 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 (Umax < 0 for westward jets).
On the β plane the Coriolis parameter f (y) in system (2.1) is given by f (y) = fc + βy, where fc = 2Ωsin(ϕc) (ϕc is the latitude of the channel center) and β = 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 Us(y), hs(y) by letting (u′, υ′, η) = (u − Us, υ, h − hs) one obtains the linear nondimensional set (upon dropping the primes from 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 hs(y) is different for the nondivergent (Kuo 1973) and divergent perturbations and these two cases are examined separately next. For nondivergent perturbations, hs(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 hs(y) in the divergent case [recall that Us(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 ux + υy = 0 implies the existence of a streamfunction ψ 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 −eik(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, Umax 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 kCi(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 ϕc = 30° and δϕ = 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, kCi(Umax, k), shown in Fig. 2b, conform to the style adopted in this work where k and Umax are the primary model parameters used as the axes in plotting the growth rate contours for fixed values of the other parameters, ϕc and δϕ.
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 = 31/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 kCi = 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 |Umax| 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, eik(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 (Umax, 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, Us(ϕ), 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 = Rcosϕ(2ΩRcosϕ + u) and system (3.2) is transformed to
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 hs(ϕ) and vanishing meridional velocity. Following Eq. (3.3b) the angular momentum in the mean state is
from which the mean zonal current Us(ϕ) 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 − Ms(ϕ), υ, η = h − hs(ϕ). 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 eik(λ−Ct). In Eq. (3.3a) the ϕ derivative occurs only in the term υdMs/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 hs(ϕ) [and the corresponding Us(ϕ), Ms(ϕ)] 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 Umax [which determines Ms(ϕ), Us(ϕ) explicitly and hs(ϕ) 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 Cr and Ci, until the 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 Ms(ϕ) and Us(ϕ) by differentiating a given hs(ϕ) distribution. However, we strive to extend Kuo’s original theory to a sphere and therefore wish to employ the same mean current Us(ϕ) used in his β-plane theory. Since the mean current in Kuo’s theory was a cos2 jet structure across the channel, we also chose the mean current to be
Given this choice of the mean current Us(ϕ), the calculation of Ms(ϕ) follows directly from the definition Ms(ϕ) = Rcosϕ[2ΩRcosϕ + Us(ϕ)], but the corresponding hs(ϕ) can only be calculated by integrating numerically the inverse of Eq. (3.5); that is,
The hs(ϕ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 Umax > 0) the mean height field, hs(ϕ), is increasing with ϕ while for westward-flowing mean jets (i.e., for Umax < 0) the height hs(ϕ) is a decreasing function. To avoid circumstances in which hs(ϕ) 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 hs(ϕc + δϕ) = 1 instead of hs(ϕc − δϕ) = 1.
In the present problem contours of the growth rate, kCi, contain troughs in the (Umax, k) plane located between local maxima, where kCi = 0 (see Figs. 3 and 4). Thus, a continuation method in which the numerical search for the eigenvalue (Cr, Ci) at any (Umax, k) point begins with the (Cr, Ci) values found at the previous (Umax, 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 (Umax, k) point, we initiated a numerically expensive search for Cr, Ci starting from nine different initial guesses of Cr and Ci; in each of these initial guesses Cr assumed the values of 0.75Umax, 0.1Umax, and 0.01Umax, while Ci assumed the values of 0.005, 0.05, and 0.5. The growth rates displayed in Figs. 3 and 4 are the maximal kCi from all searches at that (Umax, k).
5. Instabilities of divergent perturbations on a sphere
Our nondimensional, divergent, spherical, model of a cos2 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 Umax (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′ = H1H2/ (H1 + H2) 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 kCi(Umax, 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 105 m2 s−2), ϕc = 0.75 rad (∼43°) for −1 < Umax < 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πRcos(ϕc)/k ∼ 2000 km] cutoff for either Umax > 0 or Umax < 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 Umax (negative or positive) below which the growth rates vanish, and the contours terminate only on the Umax = 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 |Umax| values in contrast to the single, global, maximum on the β plane at Umax < 0, which was calculated for |Umax| = 0.01.
On the other hand, the asymmetry between positive and negative values of Umax that is so clearly evident in Kuo’s nondivergent results (Fig. 1b), where a longwave cutoff exists only for Umax > 0 and the maximal contour occurs at Umax < 0, has a much more complex fingerprints in Fig. 4. At some values of δϕ, k, and |Umax| 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 kCi 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 kCi contours.
A comparison between the growth rates in Fig. 4 and those for other values of α, ϕc shows that the growth rate contours, kCi, are slightly affected by even large changes in the value of α The maximal contours move closer to the Umax = 0 ordinate and for the same values of Umax and k the growth rates at α = 10−4, ϕc = 1.0 are about half those at α = 0.1, ϕc = 0.75, but they occur at higher values of 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 ϕc = 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 Umax > 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°; Umax = +0.2; k = 5; C = (4.431, 1.734i) × 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, ϕc = 30°, Umax = +0.2, k = 5 but for α = 0.1 and C = (6.151, 4.386i) × 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/|Umax|). The maximum growth rate at finite Umax (for fixed β and d) in Fig. 1b results only from the Kuo scaling of kCi on d/Umax, which implies that kCi can increase with the increase in Umax when moving to the right of b = −0.15, while kCi/(d/Umax) decreases. This point of misleading scaling by Kuo is evident in Fig. 7, where contours of kCi(δϕ/π)/Umax are plotted using the same kCi values as in Fig. 2. When the kCi contours of Fig. 2a are compared with the kCi(δϕ/π)/Umax contours of Fig. 7, it becomes clear that the apparent maximum at b < 0 in the latter is due only to the scaling of kCi in Kuo’s work on Umax (δϕ/π 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.
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.
Corresponding author address: Nathan Paldor, Department of Atmospheric Sciences, Institute of Earth Sciences, Hebrew University of Jerusalem, Jerusalem 91904, Israel. Email: email@example.com