## Abstract

The baroclinic instability of a meridional current in a north–south channel is investigated in a two-layer model for the case when the current has no horizontal shear. The vertical shear of the current provides a potential vorticity gradient in the zonal direction while the beta effect provides a potential vorticity gradient in the meridional direction. The normal modes of the two-layer baroclinic flow are found both numerically and analytically.

In contrast to the situation when the current is in the zonal direction there seems to be no minimum shear required for instability in spite of the active presence of the planetary vorticity gradient, *β,* although the growth rates of the instability are reduced as the shear is weakened. Also, the horizontal structure of the unstable mode is a strong function of the parameters, and weakly growing modes exhibit a boundary layer structure and are compressed to a narrow region near the western edge of the channel. The unstable modes are connected to the neutral Rossby modes that exist in the channel in the absence of shear and an analysis shows how those modes provide information about the long-wave stability threshold of the flow. The short-wave cutoff, which coincides with the Eady cutoff for zero *β* moves to *higher* along-channel wavenumber as *β* increases, thus expanding the range of instability on the short-wave side of the instability interval in wavenumber. For weak shears (or large values of *β*) the weakly unstable modes are very oscillatory in the zonal direction. Each mode occupies a small interval in meridional wavenumber and a connection to the neutral Rossby mode is inferred. In the special case when the two layers of the model have equal thicknesses, symmetries of the basic equations allow modes with the same growth rate but differing phase speeds and, hence, *linear* vacillating modes result when two such modes are superposed.

The results of the modal analysis are checked against a direct numerical integration of the initial value problem with excellent agreement and provides a point of contact with earlier numerical calculations by other authors. These findings support the hypothesis that baroclinic instability of midocean flows may represent a significant source of eddy energy.

## 1. Introduction

The stability properties of baroclinic zonal flows on the *β* plane have been long studied as a consequence of their relevance to the explanation of the development of wavelike disturbances in atmospheric flow whose wavefree state would be essentially zonal. Necessary conditions for instability, bounds on growth rates, and phase speeds of unstable normal modes are well known [see Pedlosky (1987) for a review]. Additional theorems speaking to the range of wavenumbers, in particular the maximum allowable wavenumber for instability for the modes, can also be developed in some cases (Pedlosky 1964). Particular examples of the normal modes of instability of such flows have formed the conceptual paradigm most researchers rely upon in thinking about the general stability properties of atmospheric and oceanic flows.

While such models are certainly pertinent to those oceanic flows that are largely zonal, the ocean provides examples of broad, baroclinic flows that deviate substantially from zonality. Indeed, in the eastern portions of the major subtropical gyres the flow, constrained by the presence of the eastern boundary of its basin, is more meridional than zonal. The shears of such Sverdrup interior flows are relatively weak compared to the western boundary currents but they contain substantial stores of available potential energy and it is natural to anticipate local instabilities growing on that potential energy stored in the east–west slope of the supporting density surfaces in the meridional shear.

If the flow were horizontally uniform and laterally unbounded, such currents would always be unstable. This is possible since disturbances consisting of plane waves with crests oriented zonally would release the potential energy of the current without suffering the stabilizing effects of *β* because the perturbation velocities would be strictly zonal. Such an idealized limit is, of course, never realized. The finite size of the gyre, the variation of the basic current, and the inevitable nonuniformity of initial conditions will always impose a zonal variation of the structure of the disturbance. It is natural then to ask what the role of the planetary vorticity gradient will be in the development of the instabilities of meridional currents of laterally finite extent.

Spall (1994, 2000) has addressed some of these issues in papers dealing with the variability of the circulation of the eastern North Atlantic. The emphasis in his papers cited above is the nonlinear development of disturbances from initial data, their eventual amplitude equilibration, and their effect on tracer fields such as salinity. The disturbances in his numerical calculation are permitted to develop and leave the domain of the flow through the device of an absorbing layer near the western boundary of the region of calculation. The flow is strongly unstable and secondary instabilities spontaneously develop before the disturbances are able to evolve to anything resembling a normal mode structure.

Kamenkovich and Pedlosky (1996) examined the linear instability of nonzonal jets and showed that the combination of nonzonality and wave radiation could destabilize otherwise stable flows. Indeed, it is the nonzonality that is responsible for the radiation. The study concentrated on jets in unbounded regions but did not confront directly the case of purely meridional flow.

We believe it would be of value to understand the underlying normal mode structure of such flows, their characteristic growth rates, and propagation speeds as well as the parameter domain in which they would be observed and the relation of such modes to both the initial value problem and the subsequent nonlinear evolution and equilibration of the disturbances. The normal mode problem is an important element in understanding the basic physics of the instability of such flows, even when the individual modes only eventuate for very long times. The spectrum of the modes presents us with bounds on growth rates and useful information of preferred structures.

The present study forms a first step in that effort and presents an analysis of the normal modes of an idealized channel flow oriented in the north–south direction. This idealization is meant to substitute for the longitudinal variation a current of finite width would impose on the perturbation without the added complexity of dealing with both horizontal and vertical shear. The flow in this study has only vertical shear and the channel width can be of arbitrary size with respect to the deformation radius. This is, therefore, the Phillips model (1954) turned 90 degrees.

We find that in contrast to the zonal flow case there is no minimum shear required for instability. Similarly, we find that there is instability even for channel widths less than the deformation radius. The modal structure in the cross-channel direction is a strong function of the parameters and the mode's growth rate. We infer from the sum of our results a strong connection between the instability occurring at weak shears and the existence of Rossby normal modes although we have not been able to prove the connection conclusively. Although the growth rates of the cases with weak shears are small, there is no reason to believe that the resulting amplitudes in an equilibrated nonlinear state need be small. Indeed, the analysis of Shepherd (1988) shows how an a priori bound on the amplitude of the equilibrated state can be given in terms of the degree of supercriticality of the flow. Absent a threshold for instability, the supercriticality is always *O*(1) and it is natural to assume that the resulting equilibrated amplitudes will be consequently larger than in the zonal flow case. That will depend on the efficacy of the equilibration process, which itself may be weak for broad meridional flows, as Spall (2000) suggests; therefore, the existence of such weak-shear modes may be oceanographically very significant.

Section two describes the basic model and its formulation. The symmetry properties of the model when the layer thicknesses are equal are described, as well as the symmetry between northward and southward flowing currents. Section 3 gives a preliminary overview of our detailed calculation and provides a starting point for an analytical prediction of both the long-wave and short-wave cutoffs of the model. We also present in section 4 an asymptotic analysis of the instabilities for a very broad channel that illustrates the tendency for western intensification of the modes of instability. Section 5 is a detailed look at the modal calculations in which the phase speeds and growth rates for the modes are presented and contrasted with the zonal flow case. Section 6 briefly describes the model with unequal layer thicknesses in which certain symmetries of the basic model are broken. Section 7 describes the results of initial value calculations that provide a link between the modal theory and numerical calculations of Spall, at least in the early linear stages. We conclude in section 8 with a discussion of our results and their significance.

## 2. The model

We employ throughout the two-layer quasigeostrophic model on the *β* plane (Pedlosky 1987). Since the basic flow is nonzonal, this implies a potential vorticity (PV) source in each layer for which the basic flow is nonzonal. We will restrict attention to cases in which the basic flow is a uniform flow in only the upper of the two layers so that the PV forcing is only in the upper layer. The reader may think of it as a uniform wind stress curl acting on the upper layer alone.

Thus the mean flow will be taken to be

where subscripts refer to the relevant layer: 1 for the upper layer and 2 for the lower layer.

The flow is contained in a channel on the *β* plane. The disturbance streamfunctions in each layer will be denoted by *ϕ*_{n}(*x, **y, **t*), where *x* and *y* are the zonal and meridional coordinates and *t* is the time. The spatial coordinates are scaled with the deformation radius. When each layer has the same thickness, *H,* this scale is

while the time is scaled on the advective time, *L*_{scale}/*V.* If velocities are also scaled with *V* the linearized perturbation equations for *ϕ*_{n}(*x, **y, **t*) are

In the above equations the mean meridional velocity in the upper layer is unity in dimensionless units. The parameters *F*_{n} are

where *H*_{n} are the layer depths. The choice (2.2) for the length scale allows these parameters each to be set equal to unity when, as here, *H*_{1} = *H*_{2} = *H.* The single parameter remaining is the ratio of the planetary vorticity gradient in the *y* direction to the PV gradient in the *x* direction associated with the thickness variation produced by the vertical shear of the meridional flow,

where each term on the right-hand side of (2.5) is dimensional. We note that (2.5) is also the ratio of the meridional PV gradient to the zonal PV gradient.

It is important to note that the major physical difference between the problem of the instability of a meridional shear and the zonal equivalent is that, as can be seen from (2.3), the potential vorticity gradients associated with *β* and the potential vorticity gradient of the basic flow are not colinear. It is therefore impossible for the beta effect to completely determine the sign of the potential vorticity gradient of the ambient flow as seen by the perturbation field. This is the fundamental physical novelty of the meridional instability problem.

We will search for wavelike solutions of the linear perturbation equations in the form:

Instability corresponds to solutions with complex *c* such that the imaginary part of *c* > 0. The *y* wavenumber *l* is real. The *x*-dependent amplitudes *A*_{n} satisfy the ordinary differential equations, obtained by placing (2.6) in (2.3a,b), here given for the case of equal layer depths:

subject to the condition that each *A*_{n} vanish on *x* = 0, *L* the (*L* is nondimensional and scaled with the deformation radius) position of the channel's side walls. We note that the orthogonal character of the two PV gradients has introduced complex coefficients into the perturbation equations.

Try as we might, we have not been able, using the usual integral methods, to deduce a useful necessary condition for instability, that is, one that would give us a maximum value of *β* beyond which the flow would be definitely stable. Nor have we been able to deduce useful bounds on phase speeds and growth rates. The presence of the complex term in (2.7) has thwarted such attempts, and we have come to believe this is a reflection of the absence of a well-defined stability threshold, a result which is further suggested by our calculations that follow.

There are some useful symmetry properties of the equations worth noting. If we write

then it is possible to show directly from (2.7a,b) that, if

is a solution, then from the invariant properties of the equations under complex conjugation,

is also a solution.

This implies that a second solution with the same growth rate and a real phase speed placed, with respect to the first, symmetrically about *c* = 0.5 will also be a solution *if the wave amplitudes of the two layers are exchanged after complex conjugation.* This in turn implies that any spatial slope of the phase lines in the *x*–*y* plane in each layer of one mode,

will be reversed from one mode to another.

We will also find solutions in which the real part of *s* is identically zero so that the perturbation is moving exactly with the mean speed of the current. The two “modes” obtained in that case from the symmetry property are really the same mode but the invariance in that case implies that the phase lines in the two layers will be sloping in opposite directions.

Similarly, under a change in sign of the basic flow it is possible, using the invariant properties of the equations, to show that a solution with positive *V* implies the same solution for the same negative *V* with a change in sign of the real phase speed and an interchange of wave functions between the two layers with the same growth rate.

These precise symmetries are broken when the layer thicknesses are unequal, as will be discussed in section 6.

## 3. Overview of results and long- and short-wave cutoffs

### a. Method of solution

Solutions to the eigenvalue problem (2.3a,b) are rendered technically complicated by the presence of the first derivative in *x* in the equations as well as the complex coefficient of that term. We have proceeded in two ways to find the eigensolutions. In the first, semi-analytical method, we search, for each trial value of *c,* for solutions of the set in the form,

This leads to a fourth-order polynomial equation for *k*_{j}. For a given *c,* the roots can be found numerically or, since the quartic is analytically tractable, the solutions can be found directly from the standard solution for quartics (Abramowitz and Stegun 1972). For each root the ratio *a*_{1j}/*a*_{2j} is determined by either of (2.3a,b). Application of the four boundary conditions at the channel boundaries yields a set of four homogeneous algebraic equations whose determinant must vanish for nontrivial solutions. Only for the complex eigenvalues of the problem will the determinant vanish and this condition, implemented numerically, provides us with the eigenvalues and the eigenfunctions from (3.1). These solutions were checked by simply discretizing the ordinary differential equations (2.3a,b), which turns the problem into a large matrix problem XA = *c*A with *c* as the eigenvalue and X is the matrix resulting from the discretization. We have used the MATLAB function [*c, **A*] = eig (*X*) to find the solutions. This second method becomes unwieldy when the growth rates are small, that is, for large *β* where the eigenfunctions are very oscillatory in *x.* In some cases we have checked both methods by time stepping the time-dependent version of the perturbation equations, assuming an initial disturbance of a particular *y* wavenumber. The asymptotic state of the time-stepping problem yields the most unstable mode at that *y* wavenumber in agreement with the two methods described above (which yield *all* cross-channel modes for a given *y* wavenumber).

### b. A preliminary sample of results

Before proceeding with further analytical results it is useful to examine in a preliminary way the results of the numerical calculations outlined above since it is those results that suggested the analytical approaches to be presented and justify what might seem to be arbitrary assumptions about the solutions.

Figure 1 shows the eigenvalues *c* for the case where *L* = 10, *β* = 0.5 and for equal layer thicknesses. The eigenvalues of the first two modes are shown. Starting at the right-hand side of the figure, at short waves, we see two branches of the first mode with real phase speeds symmetrically centered around the mean flow speed of *c* = 0.5. As the *y* wavenumber is reduced, the two branches coalesce at the phase speed *c* = 0.5. This is found to be a general birthplace for unstable modes. This occurs first at the point M_{1}. At still shorter wavenumbers the phase speed becomes complex. The positive values of imag (*c*) = *c*_{i} are shown with dark crosses, while its negative, less pertinent root is shown in faint gray. We were astonished to find, in contrast to the zonal flow case, that further decrease in wavenumber, while leading to increased growth rate, proceeded with the same value of real (*c*) = *c*_{r} = 0.5. In the interval between points M_{1} and M_{2} there is a single unstable mode, but we note that another, second mode with two real phase speeds is in the process of having its phase speeds coalesce as the point M_{2} is approached from above. At this point the second mode becomes unstable. At longer waves, that is, smaller *y* wavenumber, between M_{2} and J_{m} there are now *two* unstable modes. They have the same real part of the phase speed (again equal to 0.5) with different growth rates. It is remarkable that in this range, although the beta effect is *O*(1), the phase speed of the unstable wave remains precisely equal to the mean flow velocity of the basic shear. However, as the point J_{m} is approached from above the growth rates of the two modes coalesce. For wavenumbers less than that critical coalescence value the growth rates of the two modes, remarkably, remain identical to each other but now the real part of the phase speed, *c*_{r}, symmetrically diverges from the mean flow velocity. One mode, with decreasing wavenumber has its real phase speed approach the minimum velocity of the mean flow, here zero, while the other mode has it value of *c*_{r} approach the maximum value of the mean flow velocity, here equal to unity. As these values are reached the growth rate for each mode goes to zero. In the interval between J_{m} and the long-wave cutoff the modes satisfy the invariance relation described in section 2 (i.e., *s* → −*s**). They have the same growth rate and values of *c*_{r} reflected about the value 0.5. We shall see below that the eigenfunctions satisfy the predicted reflection properties between the layers given by the invariance properties. We shall be able to predict the long- and short-wave cutoff point for each mode analytically using our knowledge of the real part of the phase speeds at those points. It is perhaps not obvious from the figures, but the short-wave cutoff corresponds to a branch point, square root behavior of the dispersion relation familiar from the zonal problem. The long-wave cutoff is more subtle, as we shall see below.

Figure 2 shows the eigenvalues, for the same parameters, of the first six channel modes. Only the unstable modes are now shown for the sake of clarity. Note the persistence of *c*_{r} at and around the value of 1/2. Also evident is the complicating pairing of modes as their growth rates coalesce to produce modes with identical growth rates but distinct, symmetrically placed phase speeds. The boxes on the real axis of the lower panel show the position of the long-wave cutoff for each mode. Although 6 modes are shown there are really an infinite number of cross-stream modes that are unstable (again in contrast to the zonal flow case). The low wavenumber end of the spectrum has an infinite number of modes with cross-channel structures that vary rapidly in *x* with very small growth rates. This accumulation point at long wavelength but increasing cross-channel mode structure is reminiscent of the structure of the Rossby normal modes in a meridional channel. There the frequency of the *n*th cross-channel mode for a given *l* is

which has a point of accumulation at zero frequency for large *n.* We speculate that the weakly growing modes in this part of the spectrum are destabilized Rossby modes when the meridional shear is added, although we have not been able to demonstrate this except indirectly as discussed below.

### c. Eigenfunctions

For the value of *β* in Figs. 1 and 2 we present in Fig. 3 the eigenfunctions in the channel at selected wavenumbers. The panel in the upper left-hand corner identifies the eigenvalue associated with each eigenfunction. Each panel in the rest of the figure is identified by a letter identifying the associated eigenvalue. Each such panel is doubled, the left-hand side representing the upper layer and the right-hand panel representing the lower layer.

It is important to bear in mind that the *y* wavenumber is *O*(1) and the channel width is *L* = 10, so for the purposes of spatial compression the aspect ratio has been artificially foreshortened. We see that unlike the zonal flow case the horizontal structure of the eigenfunctions strongly varies in parameter space as a function of wavenumber. We start by examining the solutions at high wavenumbers. The panels A1 and A2 show the neutral solutions to the right of the short-wave cutoff. By the time we reach point B the two modes have coalesced and become a growing mode with *c*_{r} = 0.5. There is a slight phase shift between the two layers. In this case, as predicted by the invariance properties of the equations, with the value of *s* = 0, the eigensolutions slope oppositely in the two layers. The point C represents a point in the wavenumber interval where a second mode has become unstable, again with *c*_{r} = 1/2. We see from the figure that it is indeed recognizable as a second cross-stream mode with, again, oppositely directed lines of constant phase in the two layers. The point D represents the coalescence of the solutions on the B and C branches. For longer wavenumbers there are two solutions with the same growth rates but different phase speeds. These are the solutions E1 and E2, and we note they obey the symmetry properties previously described. Each mode is a mirror of the other upon reflection between layers and a reversal of the direction of phase shift in the *x*–*y* plane. Following the solution curves to lower wavenumbers we note that by point F the amplitudes of the eigenfunctions have become quite different in the two layers. The amplitude in the second layer, *on this branch,* is much reduced. The opposite would be true on the upper branch, that is, it would be the amplitude in the upper layer that would be small. At the point G, very near the long-wave cutoff, the amplitude in the second layer has nearly vanished. The point H is actually in the stable, long-wave region and is not labeled on the upper-left panel but is at *l* = 0.45 and a strongly negative phase speed (*c*_{r} = −3.1). Note for this parameter setting, the most unstable wave corresponds to the *pair *E1, E2 with the same growth rate but differing phase speeds. A superposition of these two modes would yield a form of linear amplitude vacillation. We also note the tendency for each of these modes to be somewhat western intensified in at least one of the layers.

### d. The short-wave cutoff

We use the fact, observed above, that at the short-wave cutoff for the transition from stability to instability the real part of the phase speed is equal to the mean velocity of the meridional current, that is, *c*_{r} = 0.5. If we assume that value and set the imaginary part of *c* to zero the equations for the sum and difference of the eigenfunctions *A*_{n} can be obtained from (2.7) as

Eliminating *A*_{T} from the above equations yields a single equation for the barotropic field,

Solutions corresponding to the vanishing of the eigenfunction at *x* = 0 and *x* = *L* can be found in the form,

leading to the condition relating *l*_{o}, the short-wave cutoff and *β*; that is,

for all integer *m.*

We note that for the case of zero *β,* the short-wave cutoff corresponds to the value given by the Eady problem (in its two-layer version), namely, *l*_{o} = (2 − *m*^{2}*π*^{2}/*L*^{2})^{1/2} corresponding to *A*_{T} = 0. For small *β,* as Fig. 3 shows, the eigenfunctions are largely barotropic near the short-wave cutoff. The interesting consequence of (3.4) is that the position of the short-wave cutoff shifts to *larger* values of *y* wavenumber as *β* is increased. Figure 4 shows a plot of the critical value of *y* wavenumber, *l,* (on the abscissa) vs the value of *β* on the ordinate. In this particular setting *L* = 100 and *m* = 1. For *β* = 0 the solution asymptotes to the Eady value, very near 2^{1/2} for this broad channel and then increases as *β* increases. The results shown here agree remarkably well with our numerical results for the cutoff. The branch point character of the critical curve can be used to obtain a splitting of the roots for smaller wavenumbers than the critical value in much the same manner as the perturbation method in Pedlosky (1970). For the sake of brevity those details are not given here but a similar analysis is given in section 4 in our discussion of the asymptotic solution for a very broad channel.

### e. Longwave cutoff

We have noted that our results indicate that, as the wavenumber of the unstable wave is reduced, the phase speed of the mode tends to either the maximum or minimum value of the current, that is, either 0 or 1. At that point the mode becomes stable and the eigenfunction is limited to the layer whose current speed is not equal to *c*_{r}. We exploit that fact to first find an analytical expression for the long-wave cutoff. More importantly, the analysis and its perturbation to find the imaginary part of *c* near the cutoff emphasizes the relationship between the unstable mode and the Rossby normal mode in the channel.

Consider now the limit in which *c*_{r} → 0 at a wavenumber, *l*_{c}, to be determined.

In a neighborhood about that point we write,

where *μ* ≪ 1 and *δ* is an order one trace parameter measuring the deviation of the *y* wavenumber from its critical value and *μ* the order of magnitude of that difference.

On the basis of our numerical results the asymptotic solution is found in the form:

At lowest order the equation for the upper layer then yields

This is simply the wave equation for a Rossby normal mode in a channel of width *L* where the frequency of the wave is set at the advective frequency in the upper layer of the wave with wavenumber *l*_{c}. The solution satisfying zero boundary conditions at *x* = 0 and *x* = *L* is

In order to satisfy the conditions on the channel walls,

leading to a condition on the *y* wavenumber,

Figure 5 shows the eigenvalues as a function of wavenumber for the most unstable mode at various values of *β* in the left-hand panels, as well as showing the most unstable eigenfunction at that value of *β.* The long-wave cutoff is indicated by a small box in the figure on the wavenumber axis. The eigenvalue curves are generated by numerically finding the roots of the dispersion relation as described above, but the boxes are placed at the point predicted by (3.10). The agreement is very good. The important implication of the agreement is our interpretation that the instability springs from the destabilization, for small *l* of the Rossby normal mode in the channel whose frequency is the advective frequency in the other layer at *l*_{c}. The case we have described is for *c* = 0 but the symmetric case where *c* = 1 yields the same threshold with an interchange of layers. If the perturbation expansion is carried further to calculate the higher-order wave field in the lower layer,

The higher-order wave field in the upper layer satisfies

In order to find the eigenvalue *c*_{2} whose imaginary part will determine the growth rate, it is necessary to first find *A*^{(1)}_{2} and then use that in (3.12), whose solvability condition will yield *c*_{2}. The process is somewhat lengthy and almost straightforward. The one subtle feature is that the order of (3.11) as written is too low to satisfy both boundary conditions on *x* = 0 and *x* = *L.* To the general solution of (3.11) one must add a boundary layer solution of the problem:

whose solution is

For small *μ* this yields a boundary layer solution near *x* = 0 for the growing mode and this solution must be added to the solution of the interior problem given by (3.11) to obtain the full solution for *A*^{(1)}_{2}. When the solvability condition for (3.12) is used, an equation for *c*_{2} is obtained that may be put in the form:

where *T* is given in appendix A. It is not difficult to show that for instability Im(*T*) < 0 is needed and is obtained. This turns out to be the case *only if* the boundary layer is placed on *x* = 0 and if *δ* > 0, that is, only for waves shorter than the long-wave cutoff. If one looks for a solution with *δ* < 0, that is, on the long-wave part, it appears that *c*_{i} < 0. But then one must place the boundary layer on *x* = *L.* When that happens the sign of the imaginary part of *T* changes and a contradiction occurs. It follows that near the long-wave cutoff only for wavenumbers to the right of the cutoff are solutions possible with *c*_{i} ≠ 0 in agreement with the numerical results. Indeed the instability of waves shorter than the long-wave cutoff depends critically on the existence of the boundary layer solution (3.13b). We note that the solution is similar in structure to Stommel's boundary layer solution for the western intensification of the oceanic circulation except that here, instead of dissipation, it is the growth rate that acts to choose the boundary on which the intensification occurs. This interchange of role between growth and dissipation is quite common in wave problems. Figure 6 shows the predicted real and imaginary parts of *c* in the neighborhood of the long-wave cutoff for the case *β* = 0.9 and *L* = 10. The linearity of both the real and imaginary parts of *c* with respect to wavenumber is quite distinct from the behavior at the short-wave cutoff.

## 4. The broad channel approximation and the shortwave cutoff

A particularly simple but illuminating analytical approximation to the stability problem occurs when the channel width is large, that is, when *L* ≫ 1. In dimensional units the width is much greater than the deformation radius. A formal expansion can then be developed whose details will only be sketched here. The meridional variable, *y,* is unchanged, but it can be anticipated that, at least near the critical wavenumber for instability where the growth is slow, that the *x* scale will be large, of the order of the channel width so that the disturbance can avoid the stabilizing effect of *β* as much as possible. The *x* variable can then be rescaled

A new time variable *T* = ɛ*t* is also introduced. When these rescaled variables are inserted in (2.3), a series of equations is obtained by the usual methods of multiple time and space scales. The lowest-order problem retrieves the Eady solution for the infinite domain; that is,

where

and where the amplitude *A,* equal to lowest order in both layers is a function of *ξ* and *T.* Continuing the expansion in powers of ɛ eventually yields, by a solvability condition on the problem at *O*(ɛ)^{2}, the governing equation for *A,*

where Δ*l* is the small, *O*(ɛ^{2}) shift in the *y* wavenumber from its critical value of 2^{1/2}. The reader can check that, if the domain were truly infinite in *x,* the plane wave solution of (4.4) would match the standard dispersion relation for the two layer model (see, e.g., Pedlosky (1987, sec. 7.13) in the limit of small *x* wavenumber.

For the channel problem solutions of (4.4) can be sought in the form

so that *F*(*ξ*) satisfies

The solution of the form *F*(*ξ*) = sin*mπξ* that satisfies the boundary conditions at *ξ* = 0, 1 yields the dispersion relation,

so that instability will occur only for Δ*l* < 0, that is, on the long-wave side of the short-wave cutoff and arise as a square root branch point as predicted in section 3d. Indeed the critical condition arising from (4.7) when *ω* is set equal to zero is simply the limit of (3.4) when *L* is very large. What is more illuminating is the form of the eigenfunction. The exponential factor in(4.5) will be of the form

so that for the growing mode there is a western intensification of the disturbance proportional to the growth rate. As *β* is increased, this scale becomes ever shorter until finally the scale is short enough (of the order of the deformation radius in dimensional units) for our expansion based on slow *x* variations to fail. Nevertheless, the qualitative behavior agrees with the results of the numerical calculations (compare the solutions at points B and D in Fig. 3).

## 5. Results

### a. Dependence on β

Were the channel oriented in the zonal direction the flow would be stable if *β* were greater than unity. Figure 5 shows the dispersion curves for various values of *β* as well as the eigenfunction for the most unstable wave. We note that even for *β* = 1.5 the flow is still unstable. The growth rate is reduced for *β* > 1, but there is a tendency for the most unstable mode to shift to higher wavenumber and this somewhat compensates for the reduction in the value of *c*_{i} in determining the growth rate. For *β* < 1 the most unstable mode is the gravest mode in the channel; this is no longer true for *β* > 1 and for each value of *β* it is difficult to predict a priori either the most unstable wavenumber or the character of the cross-stream structure. For the case where *β* = 1.5, it is clear that the most unstable mode is strongly western intensified. There are, moreover, many unstable modes in the short-wave interval for *β* > 1 with each mode being unstable in a small interval of *y* wavenumber. Figure 7 shows in detail the curves for phase speed and growth rate for *β* = 1.5 for 1.2 < *l* < 1.5 while the panels to the right and also below show the eigenfunctions at the labeled points. This interval in wavenumber contains a large number of modes. As the wavenumber is altered, branches of the phase speed coalesce and become unstable. Moving from the points marked 4 and 3 to point 2 in the upper right of the phase speed window to the points C, F, and H a complex series of coalescence and bifurcations occur. Throughout the transitions the upper-layer eigenfunction appears to remain in a gravest mode while the lower-layer structure becomes increasingly oscillatory. The growing mode at point 4 occupies a very small range of *l* and occurs as the result of the coalescence of two neutral modes. It has a structure similar to that of point B of Fig. 3, namely a gravest mode in both layers with its real phase speed = 0.5. Where the mode becomes stable, the two neutral modes each have a *c*_{r} that diverges from 0.5. The branch with *c*_{r} < 0.5 then merges with a neutral mode with a second cross-channel mode structure and becomes a growing mode at point 3 with *c*_{r} < 0.5. This mode has a gravest modelike structure in the upper layer and a second mode structure in the lower layer. The same progression occurs for points C, F, and H as neutral modes with progressively higher cross-channel structures become unstable. After point H this mode does not stabilize until the long-wave cutoff is attained at *l* = 0.84. The obvious connection between the instability and the existence of Rossby neutral modes and the long-wave cutoff's relation to the layer-limited Rossby mode again suggests that the instability described here is a destabilization of the infinite set of Rossby normal modes by the shear. Even a small amount of vertical shear is apparently able to alter the vertical phase of the neutral mode to allow energy release. It would seem possible to exploit that fact to produce a perturbation theory at small shears showing the destabilization. So far, we have failed in that attempt but our conjecture of its validity is based on a broad suite of the numerical results that we have obtained. We have found instabilities at larger values of *β* but with ever smaller growth rates and increasingly oscillatory structures. These weak modes are not presented here.

### b. Channel width

For the case of a zonal channel there is a minimum width of the channel below which the flow is stable. This is the Eady short-wave cutoff and would occur at *L* = *π*/2^{1/2} = 2.222. When the flow is meridional and *β* is not zero, even a very narrow channel can still support instability for a narrow range of wavenumbers. Figure 8 shows the eigenfunctions and lists the growth rates for various channel widths. In the first panel the case where *L* = 0.5 is indicated. At these narrow channel widths we have found that instability occurs when two neutral modes coincide with the same frequency (or real phase speed at a given *l*). This is shown in Fig. 9. For the case *L* = 2, a channel width for which the zonally oriented flow would be stable, the range *l* = [0.35–0.52] is now unstable. The associated eigenfunctions are shown in Fig. 8. For very large *L,* for example, *L* = 50, the eigenfunction is similar to that at *L* = 10 except that the function occupies a smaller zone of the total channel in its western region.

## 6. Unequal layer depths

When the layer depths are unequal, as might be considered more pertinent for the ocean, the symmetries of the equal layer case are broken. The governing equations (2.7) become

where

We choose the scaling length to be the deformation radius of the baroclinic mode of the two-layer model

and this sets *F*_{1} = 2*H*_{2}/(*H*_{1} + *H*_{2}), *F*_{2} = 2*H*_{1}/(*H*_{1} + *H*_{2}) so that their sum is always 2 and each reduces to unity in the case of equal layer depths.

We have redone our calculations for the case where *H*_{1}/*H*_{2} = 0.25, representative of the thermocline and deep ocean. The equal layer limit might be better appropriate for an attempt to resolve instabilities *within* the thermocline itself. Figure 10 shows the numerically calculated dispersion relation for *L* = 10 and *β* = 0.5. The growth rates are not significantly different than the equal layer case. However the along-channel wavenumber of the most unstable mode has been reduced from *l* = 0.99 to *l* = 0.8. A similar reduction is seen at the higher value of *β* (1.5) shown in Fig. 11. The most unstable wave now lies on a branch with phase speed >0.5 and this produces eigenfunctions that slope in the *x*–*y* plane oppositely to those previously seen, which had phase speeds equal to or less than 0.5. The eigenfunctions are again western intensified. The long-wave cutoff is not altered by the unequal layer thickness since the problem for that cutoff is independent of the parameters *F*_{n}.

## 7. The initial value problem

In order to verify the eigenvalue calculations, especially in the case where *β* > 1, we have carried out integrations of Eqs. (2.3) when the initial condition consists of a disturbance of a single *l* in the *y* direction. The initial perturbation was a sine function in *x* satisfying the boundary conditions. The equations were finite differenced in *x* (Δ*x* = 0.1) and a second-order Runge–Kutta scheme was employed to march in time. Figure 12 shows the evolution in time of the maximum amplitude as a function of time in the panels to the left while the right-hand panels show the corresponding most unstable modes. After some time the growth is seen to approach the growth rate predicted by the normal mode theory, which is indicated by the dotted curve in the upper-left panel. In the first two panels the calculation is done for equal layer thickness. Due to the symmetry of the problem, there are two eigenfunctions sharing the same maximum growth rate but with distinct phase speeds and both emerge from the initial data. This leads to the linear vacillation of the amplitude around the predicted growth rate as the two modes alternately reinforce and interfere. On the other hand, in the last panel for the unequal layer case and for *β* = 1.5, there is only a single mode with the largest growth rate and the solution monotonically approaches the normal mode growth rate.

It is natural to suppose, and Spall (2000) has observed this behavior, that, if the channel is broad and the initial perturbation is also broad in the cross-stream direction compared to the deformation radius but is still confined initially to not feel the presence of the boundaries, that the initial growth that would be self-selected would be of a disturbance with very small cross-channel wavenumber. That is, locally the disturbance would look like an Eady mode and hardly sense the effect of *β* until the disturbance naturally broadens or propagates to interact with the boundary. Figure 13 show the results of a initial calculation in a channel 50 times wider than the deformation radius. The initial disturbance is a Gaussian in each layer, broader than the deformation radius but narrow enough to initially avoid interacting with the boundaries of the channel. The lower panel shows the growth of the disturbance at a value of *l* corresponding to the most unstable mode at this value of *β* = 0.5. We see that the initial growth matches the Eady growth rate as the disturbance first evolves with nearly no *x* variation within the Gaussian envelope. As time goes on and the disturbance senses the boundaries of the channel, the rate of growth approaches the normal mode growth rate that we have found from our analysis. Note the signature of the amplitude vacillation due to the presence of the superposition of the two normal modes with the same growth rate but diverse phase speeds. Figure 14 shows the evolution of the disturbance in the channel and the eventual emergence of the western intensified mode. Of course for smaller channel widths the eventual emergence of the normal mode will occur more rapidly.

## 8. Summary and discussion

The normal mode problem for the instability of meridional flows, even in this simple model where the flow has no horizontal shear, has a significantly different character than the zonal flow example first described by Phillips (1954). In the zonal flow case a number of a priori conditions for instability establish a threshold for the shear below which the flow must be stable. We have seen in the meridional flow case that there is no such threshold. It appears that for any value of *β* arbitrarily weak shears will become unstable although the associated growth rates will be small. We have not been able to prove the nonexistence of a threshold, but the combination of our failure to do so and the results of direct calculation strongly suggest its absence. We speculate that, from one point of view, this results from the destabilization of the Rossby normal modes that would exist in the channel in the absence of shear. This is further suggested by our analytical determination of the long-wave cutoff. This cutoff arises when a Rossby normal mode, limited to one layer, has a purely advective frequency (*Vl*) of the mean flow in one of the layers (either 0 or 1). In addition the appearance of instability for narrow channels reinforces that suggestion.

Unstable modes appear for an extended range of wavenumbers beyond the Eady cutoff when the flow is meridional and *β* is not zero in distinction to the zonal flow case, and again we believe this is due to a destabilization of the Rossby modes.

From another viewpoint it is worth pointing out that the potential vorticity gradients in the two layers are never colinear, so particle trajectories sensing potential vorticity gradients of opposite sign are always available. Although such opposition is necessary for zonal flows, in the absence of a stability theorem requiring such opposition in the meridional case that observation is only suggestive.

One of the fascinating features of the instability of the meridional flow is the strong dependence of the horizontal structure of the modes on the wavenumber and growth rate of the disturbance. The growing modes tend to be intensified in the western side of the domain. One can conjecture that this can be thought of as the result of the constant growth as a disturbance moves westward. The reflection of ever larger unstable disturbances from the western boundary would render the western boundary a source of perturbations larger than the arriving perturbations from the east leading to the observed western intensification. A somewhat similar result is seen in a nonmodal reflection process (Pedlosky 1993).

We have presented these results in an attempt to establish the basic properties of the normal modes for confined meridional flow. It is quite likely that for very broad channels or currents, the initial stage of growth in which the disturbance arranges itself to minimize the role of *β* will proceed to a stage of large enough amplitude such that secondary instabilities and nonlinear interactions will intervene before the normal modes of the present study may manifest themselves. For narrower regions or narrower currents that is less clearly the case.

We believe the subtle interactions between the nonzonal character of the flow and the *β* effect, illustrated by the normal modes described here, may be of general significance for the dynamics of nonzonal flows. It will be of special interest to extend the results of the present paper to include the effects of lateral structure of the basic current and the further nonlinear evolution and equilibration of the modes.

## Acknowledgments

This research was supported in part by NSF OCE 9901654. The authors wish to thank Mike Spall for several valuable conversations. A portion of this work is a part of the thesis of AW in fulfillment of a Master of Science degree of the MIT/WHOI Joint Program in Physical Oceanography.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### The Function *T* in (3.12)

The use of the solvability condition for Eq. (3.12) involves eliminating forcing terms that project on the adjoint operator for the solution of (3.7). This standard analysis yields

## Footnotes

*Corresponding author address:* Dr. Joseph Pedlosky, Woods Hole Oceanographic Institution, Clark 363 MS#21, Woods Hole, MA 02543. Email: jpedlosky@whoi.edu

^{*}

Woods Hole Oceanographic Institution Contribution Number 10472.