Multiple Equilibria and Low-Frequency Variability of the Wind-Driven Ocean Circulation

The link between low-frequency time-dependent variability and the existence of multiple unstable steady- state solutions in a reduced gravity quasigeostrophic ocean model for the midlatitude wind-driven circulation is investigated. It is shown that a sequence of successive symmetry-breaking pitchfork bifurcations lead to multiple equilibria that differ from each other primarily in the elongation of the recirculation cell, in the amount of meandering present in the intergyre jet, and in a north–south shift in the eastward jet. The elongation of the recirculation cells and the meandering of the jet play compensating roles in the establishment of the global energy and vorticity balance. The solutions also have distinct energy levels, but general agreement between them and the bumps in a histogram of the total energy obtained from a 1200-yr time-dependent simulation is not found. Nevertheless, a substantial fraction of the variance (30%) can be accounted for by four coherent structures that capture the subspace spanned by four vectors that point from the mean state to four selected ﬁxed points. The steady-state solution with the most elongated recirculation cells acts most strongly in steering the trajectory of the time-dependent model in phase space and sets a rough upper bound on the energy of the ﬂow.


Introduction
Analyzing sea surface height in regions of the midlatitude western boundary currents, Kelly et al. (1996) found that the dominant mode of variation of the surface currents was a change in the structure of the recirculation regions. In the 2.5-yr record analyzed, they found elongation and contraction of the recirculation gyres in both the Kuroshio Extension and the Gulf Stream, with a trend toward a shorter recirculation gyre in the Atlantic. More recently, Qiu (2000) also found similar variability for the Kuroshio extension system in 7 years of altimetry data. Both these studies show that in the state with elongated recirculation cells, the jet extension has a greater zonal penetration and a more northerly zonal-mean path and that in the contracted state, the jet extension follows a more southerly mean path. They also found that the periods with elongated recirculation cells have weaker eddy kinetic energy in the upstream region while periods with more contracted recirculation cells have higher eddy kinetic energy.
Based on a study of Rhines and Schopp (1991), Kelly et al. (1996) speculated that the tilt of the zero wind stress curl line might be responsible for the change in the elongation of the recirculation cells but found that on timescales of months, there was no correlation between the gyre fluctuations and the nonseasonal curl line tilt, while on longer timescales the observational record was too short to draw any firm conclusions. Another possibility as suggested by Qiu (2000) is that the variability in the elongation/contraction of the recirculation cells is a manifestation of the intrinsic variability of the recirculation gyre dynamics and occurs independently of changes in the wind stress.
A mode of variability with a spatial structure very similar to the one identified by Kelly et al. (1996) and Qiu (2000) exists in nonlinear wind-driven ocean models forced by steady winds. For example, McCalpin and Haidvogel (1996) used a reduced gravity quasigeostrophic model forced by a steady wind stress to study the intrinsic variability of the double-gyre wind-driven ocean circulation. They found that the low-frequency variability of the model was associated with irregular transitions between different regimes that could be characterized as having either recirculation cells that were elongated or contracted. Spall (1996) also found variability associated with the transition between regimes with either long or short recirculation gyres, using a three-layer primitive equation model that was forced by P R I M E A U a steady deep western boundary inflow in addition to the steady wind stress.
In the present study we postulate that the existence of the multiple regimes found by McCalpin and Haidvogel and by Spall is due to the existence of multiple equilibria characterized by flow patterns having recirculation gyres that are either elongated or contracted. Even though the mechanism for the transition between states might be different between the models of Mc-Calpin and Haidvogel and of Spall-the first involved only barotropic instability and the second involved baroclinic instability-the similar nature of the preferred regimes suggest that the existence of multiple regimes (as opposed to the transition timescales and frequency of the associated variability) can be studied with a simple one-layer quasigeostrophic (QG) model.
Our strategy is to apply Newton's method to a onelayer QG ocean model with the same configuration as used by McCalpin and Haidvogel and investigate the possible connection between the steady-state solutions and the multiple regime behavior they identified. Unlike other bifurcation studies of the double-gyre model, which have focused on the first bifurcations away from an antisymmetric stable steady state (e.g., Katsman et al. 2001, Dijkstra andKatsman 1997), we find it necessary to extend the bifurcation tree past several bifurcations in order to find the connection between the fixed points and the time-dependent behavior of the model. In general, one cannot expect steady solutions that are far from the marginal stability curve to have any influence on the model trajectory, but we show that for the model configuration studied by McCalpin and Haidvogel (1996), solutions that bifurcate away from the antisymmetric state far down the bifurcation tree have the strongest influence on the time-dependent flow. Furthermore, there are some important similarities between the structure of the elongated and contracted recirculation cells observed in the recirculation systems of the Gulf Stream and Kuroshio Extension systems and the model's multiple equilibria. This allows us to study the global balances of energy and vorticity of the different flow states in the simple context of a QG model without the complications of time-dependence.
The plan of the paper is as follows. In section 2 we present the model formulation and the method of solution as well as the arclength continuation strategy. In section 3, we present the multiple equilibria and the global balances of energy and vorticity that allow the very different steady-state solutions to exist with the same steady forcing and dissipation. In section 4 we present the time-dependent solution and investigate the connection between the fixed point and the behavior of the model in phase space. In section 5, we present a linear stability analysis for the steady-state solutions. Finally in section 6 we present a discussion of the results.

Model formulation and method of solution
We use the same model formulation as McCalpin and Haidvogel (1996)-the reduced-gravity quasigeostrophic vorticity equation. In terms of the interface height anomaly, h, it is written as follows gЈH 0 (refer to Table 1 for a description of the symbols). The equation describes the time evolution of the interface anomaly between two immiscible, homogeneous layers of fluid of slightly different densities. The upper layer has thickness H ϩ h(x, y, t) and the lower layer is assumed to be infinitely deep and at rest.
The domain of integration is a rectangular basin with solid walls at x ϭ 0, x ϭ L x and y ϭ 0, y ϭ L y . The flow is forced by a zonal wind stress curl profile, where the parameter A s controls the north-south asymmetry of the wind stress curl; for A s ϭ 0, the wind stress curl profile is antisymmetric about the center of the basin at y ϭ L y /2. To conserve mass, we impose the integral constraint and no normal flow where t is a unit vector tangent to the basin walls. We also impose no stress, along the basin boundaries. Finally, because of the highorder viscosity term, one additional boundary condition is required. For this, we choose

a. Nondimensionalization
Most of the results we will present in this paper are in dimensional variables. We will nevertheless convert the governing equation into nondimensional form so that it will be clear how many parameters we can vary independently, and to make it easier for the reader to compare the results with those of other studies.
We put Eq. (1) in nondimensional form by introducing the following scales where time has been nondimensionalized by the time for long Rossby waves to cross the basin and the thickness h has been nondimensionalized by the east-west thickness difference obtained from Sverdrup balance. The nondimensional form of the Eq. (1) is (after dropping the asterisks) with the nondimensional wind stress curl profile given by We have introduced the following nondimensional parameters the basin aspect ratio ␣, the width of the inertial boundary layer ␦ I , the width of the viscous boundary layer ␦ H , and the width of the frictional Stommel boundary layer ␦ S . We have also introduced the nondimensional Laplacian operator 2 2 2 2 ٌ ϵ ␣ ‫ץ‬ ϩ ‫ץ‬ .
x y The reference numerical values for the parameters are given in Table 1.

b. Spatial discretization
The spatial discretization of the model is achieved via second-order finite difference approximations. Except for the stability analysis presented in section 5 a grid of 181 ϫ 141 points, corresponding to a horizontal grid spacing of 20 km is used. For the stability analysis, memory constraints on the available computer allowed a mesh of only 141 ϫ 121 grid points. To compensate for the decreased resolution, we introduced a stretched grid to concentrate the grid points close to the western boundary and to the jet axis along y ϭ L y /2. We also varied the stretching to verify that the stability results are not too sensitively dependent on the resolution. A cubic mapping of the form form x ϭ ai 3 ϩ bi ϩ c was used in both the x and y directions to map the grid points from computational space to nondimensional physical space.

c. Steady-state solutions: continuation strategy
To obtain the steady-state solutions, we use an arclength continuation algorithm and Newton's method. See Seydel (1994) for a practical description of the method.
To find multiple steady-state solutions at the reference parameter values, we make use of the symmetry property of the quasigeostrophic governing equation, h(x, y) ϭ Ϫh(x, Ϫy ϩ 1).
Initially, we keep the wind stress profile exactly antisymmetric by setting A s ϭ 0 in Eq.
(2) so that any pitchfork bifurcation structures leading to multiple equilibria are not destroyed. [For a discussion of symmetry breaking pitchfork bifurcations applied to the doublegyre problem see Jiang et al. (1995) and Cessi and Ierley (1995)]. Note that for A s ϭ 0, the prescribed wind stress curl profile also satisfies condition (12). For A s 0, the nonsymmetric solution branches that would be connected to the antisymmetric branches at pitchfork bifurcation points become disconnected and cannot all be found by continuously varying one of the parameters.
Plot of the maximum transport across the intergyre jet for the antisymmetric solutions as a function of the biharmonic viscosity coefficient A b . Also shown is the location of the pitchfork bifurcation points, labeled PF with subscripts A through H, as well as the location of two saddle node bifurcation points NPL and NPH. For reference the Sverdrup balance transport is approximately 40 Sv.
To map out the solution branches into the desired nonlinear regime, we begin by computing an antisymmetric solution branch in the viscous regime (␦ H k ␦ I ) where the solution is essentially linear and therefore unique. We then gradually decrease the viscosity so that the solution becomes progressively more nonlinear. Only once a solution on each of the distinct nonsymmetric solution branches has been found and followed to the desired forcing and dissipation parameters, do we gradually increase A s to make the wind-forcing nonsymmetric.
In summary, our continuation strategy is the following.
1) Hold A s ϭ 0 fixed and vary A b to compute the antisymmetric branch from the viscous regime to A b ϭ 8 ϫ 10 10 m 4 s Ϫ1 . Bifurcation points are detected along the way by monitoring the sign of the determinant of the Jacobian matrix. 2) Hold A s ϭ 0 fixed and vary A b to continue each of the pitchfork branches from the value of A b at their respective bifurcation point to A b ϭ 8 ϫ 10 10 m 4 s Ϫ1 . 3) Hold A b ϭ 8 ϫ 10 10 m 4 s Ϫ1 fixed and vary the wind profile assymetry parameter to continue each branch to the reference value used by McCalpin and Haidvogel (1996), that is, from A s ϭ 0 to A s ϭ 0.05.

a. Antisymmetric solutions
If the wind stress curl profile is antisymmetric [i.e., A s ϭ 0 in Eq. (2)], at least one antisymmetric solution branch exists for all parameter values. Figure 1 shows the maximum transport in the jet for the antisymmetric solution as a function of the viscosity parameter. All other parameters except for A s are those given in Table  1; that is, we have the following nondimensional parameters Ϫ3 ␣ ϭ 7/9, ␦ ϭ 5.29 ϫ 10 , S For reference, the transport across the jet corresponding to Sverdrup balance is 4 0 /(␤ 0 ␣) ഠ 40 Sv (Sv ϵ 10 6 m 3 s Ϫ1 ). Figure 1 also shows the location of the bifurcation points leading to multiple equilibria. The antisymmetric solutions at each of these pitchfork bifurcation points are shown in Fig. 2. As such, Fig. 2 shows how the antisymmetric solution changes as the viscosity is decreased. In this sequence, the recirculation cells, which are at first confined to a region near the western wall, expand progressively farther eastward as the viscosity is reduced. There is also an increase in the in- tensity of the recirculation cells as the viscosity is reduced from 1.5 ϫ 10 13 m 4 s Ϫ1 , (␦ H ϭ 1.66 ϫ 10 Ϫ2 ), (corresponding to the location of the transition to the rapidly increasing transport in Fig. 1) to A b ϭ 1.3 ϫ 10 12 m 4 s Ϫ1 , (␦ H ϭ 1.02 ϫ 10 Ϫ2 ), where the maximum transport is reached (84 Sv). When this maximum in the transport is reached, the recirculation cells extend 600 km eastward into the basin interior. A subsequent reduction in the viscosity causes the intensity of the recirculation cells to decrease but their eastward extent continues to increase in a monotonic fashion until a saddle node bifurcation point is reached at A b ϭ 6.0 ϫ 10 9 m 4 s Ϫ1 , (␦ H ϭ 2.47 ϫ 10 Ϫ3 ). The low nose point associated with this bifurcation is labeled NPL in Fig.  1. At this point the recirculation cells extend up to 180 km west of the eastern wall. To continuously follow the antisymmetric solution, the viscosity must be increased from the low nose point, NPL, up to the high nose point, NPH, at A b ϭ 7.6 ϫ 10 9 m 4 s Ϫ1 , (␦ H ϭ 3.64 ϫ 10 Ϫ3 ). The north-south extent of the recirculation cells does not change much until the jet extends across the entire basin. Once the jet can no longer expand in the eastwest direction because of the eastern wall, the recirculation cells expand in the north-south direction. However, due to the finite value of the interfacial drag parameter (␦ S ϭ 1.41 ϫ 10 Ϫ3 ), the limit ␦ H → 0, does not tend toward basin filling gyres (Ierley and Sheremet 1995;Cessi and Ierley 1995;Primeau 1998a). Instead, the solution tends to one in which the subtropical gyre resembles one of the intermediate single gyre solutions in the sequence computed by Veronis (1966).

b. Nonsymmetric solutions
As the viscosity is decreased beyond each of the pitchfork bifurcation points labeled PF A through PF H in Fig. 1, an additional pair of nonsymmetric equilibria comes into existence. The members of each new pair are mirror images of each other and are related to each other through the symmetry condition given in Eq. (12). The pair of nonsymmetric solutions bifurcating at PF A are labeled A and AЈ. The pair bifurcating at PF B are labeled B and BЈ, and so on for the other pitchfork bifurcation points. At each bifurcation point, one of the eigenmodes of the linearized system has a zero eigenvalue. The eigenmode corresponding to this null eigenvalue captures the essential difference between the bifurcating branches. Figure 3 shows contour plots of the eigenmodes that have a zero eigenvalue at each of the bifurcation points. The structure of these modes is confined to the region of the basin occupied by the recirculation cells (cf. Fig. 2 showing the basic state and Fig. 3 showing the null eigenmode). Since these modes are symmetric about the line y ϭ ½, they destroy the antisymmetry of the flow field. The solutions on the nonsymmetric branches can be distinguished from the antisymmetric solutions by the meandering of the jet separating the recirculation cells and by a northward or southward shift in the mean position of the jet. The solutions on the successive nonsymmetric branches can, in turn, be distinguished from each other by the number of meanders in the jet with the first pitchfork bifurcation giving rise to a solution with one meander, the second with two, and so on. The nonsymmetric solutions are similar to stationary waves. In the intergyre jet, the eastward flow is strong enough to arrest the westward propagation of the Rossby wave. As the viscosity is decreased the eastward velocity in the jet does not change much, but the recirculation cells become progressively more elongated. Once the recirculation gyres become sufficiently elongated, additional meanders of the stationary wave can fit between the western wall and the jet exit region. Each new pitchfork bifurcation corresponds to the destabilization or stabilization of the basic state to a stationary wave mode with one additional meander. The stationary wave nature of the resulting nonsymmetric branch can be seen by superimposing the null eigenmodes plotted in Fig. 3 onto the corresponding steady-state solution in Fig. 2. Each eigenmode consists of an elongated cell situated over the jet axis and two weaker counterrotating cells on the north an south flanks. The amplitude of the cell overlaying the jet axis varies in the downstream direction to produce meanders in the jet when superimposed on the fixed-point flow field.
The sequence of bifurcations also has an alternating pattern in the structure of the null mode relative to the underlying basic state. The cell overlaying the jet axis for the first mode at PF A extends past the eastern most extent of the recirculation cells in the basic state. For the second mode at PF B , the cell does not quite reach the eastern end of the recirculation cells. Instead, the weaker counter rotating cells wrap around the eastern end so that the cell on the northern flank overlays the eastern tip of the northern recirculation cell and one on the southern flank overlays the recirculation cell to the south. The next mode at PF C has a pattern similar to the one at PF A and the one at PF D has a pattern similar to the mode at PF B and so on down the sequence. The result is that, when the sequence of modes are superimposed on the basic state, the nonsymmetric branches have flow fields for which the northern and southern recirculation cells alternately wrap around the southern and northern recirculation cell at their eastern end. In this way, as the viscosity is decreased and the recirculation cells expand eastward, each new bifurcating branch has a flow field with one additional stationary meander in the jet extension.

c. Multiple equilibria for the reference parameter values
In order to make a direct comparison to the simulations of McCalpin and Haidvogel (1996), each solution branch was traced out using the continuation method VOLUME 32 from A s ϭ 0 (symmetric wind profile) to A s ϭ 0.05 (assymetric wind profile). For A s ϭ 0, the solutions bifurcating at bifurcation point PF A are labeled A and AЈ, those bifurcating at PF B are labeled B and BЈ and so on for the other bifurcation points. The nonsymmetric solutions in which the jet turns southwards after separating from the western wall are denoted by the primed letters. Those with the jet turning northwards at first are denoted by unprimed letters. For A s ϭ 0, the primed and unprimed solutions are mirror images of each other. For A s ϭ 0.05, this is no longer the case, but the solutions are still qualitatively mirror images of each other.
In order to avoid a proliferation of symbols, we retain the symbol names, A and AЈ, B, and BЈ, etc. for the continuation of the branches from A s ϭ 0 to A s ϭ 0.05. Unless stated explicitly in the remainder of the article, a reference to the fixed point AЈ, for example, will refer to the fixed point for the reference parameter values used by McCalpin and Haidvogel (1996) on branch AЈ.
In Fig. 4, the primed member from each pair of nonsymmetric solutions are contoured. In Fig. 4 we show only one of the three solutions from branch AЈ since the differences among the solutions are slight rearrangements of the many closed circulation cells. The solution labeled EЈ is the continuation of the antisymmetric branch from A s ϭ 0 to A s ϭ 0.05. For A s Ͼ 0, the assymetry of the wind stress introduces imperfections in the pitchfork bifurcations such that the part of the branch which is antisymmetric between PF D and PF E for A s ϭ 0 connects the DЈ and EЈ branches without going through a pitchfork bifurcation point.
To illustrate how the assymetry of the wind profile destroys the pitchfork bifurcation, Fig. 5 shows a continuation of each branch as a function of A b with A s ϭ 0.05 held fixed. The ordinate in Fig. 5 is the sum of the VOLUME 32 5. Plot of solution branches as a function of A b for the asymmetric wind stress profile (A s ϭ 0.05). The abscissa is the biharmonic vicosity coefficient A b and the ordinate is the sum of the interface height anomaly at two points situated 200 km to the north and south of the zero wind stress curl line, and 160 km to the east of the western wall. The plot shows how the branches that were connected at the pitchfork bifurcation points for A s ϭ 0.0 become disconnected for A s ϭ 0.05. interface height anomaly, h, at two points 200 km to the north and south of the zero wind stress curl line and 160 km east of the western wall. Antisymmetric solutions would plot along the zero line. Note that the branches which were connected at pitchfork bifurcation points for A s ϭ 0.0 become disconnected for A s ϭ 0.05.
The major differences between the solutions shown in Fig. 4 is in the eastward extent of the recirculation cells, and in the amplitude of the meanders in the jet and in the jet exit region. The solutions, which bifurcate at progressively smaller values of viscosity, have progressively more elongated recirculation cells and progressively less meandering of the streamlines. The differences in the degree to which a solution has long recirculation cells with weak meandering or vice versa indicates the differences in the way the solutions achieve global balances of energy and vorticity.

d. Global energy balance
In this section we discuss the differences in the energy balances of the different steady state solutions. Scott and Straub (1998) give a discussion of the global energy balance for symmetric and nonsymmetric steady state solutions. In addition to having different flow fields, each equilibrium state has a different energy level. Figure 6 shows a plot of the total energy, TE, as a function of the biharmonic viscosity parameter A b for the case with antisymmetric wind stress curl. The total energy is given by the sum of the potential and kinetic energy,  10 13 m 4 s Ϫ1 (␦ H ϭ 1.53 ϫ 10 Ϫ2 ) recirculation cells form in the region where the western boundary currents from the subtropical and subpolar gyres meet. As the biharmonic diffusivity is further decreased, the total energy of the antisymmetric branch increases rapidly. This rapid increase in energy is accompanied by a rapid increase in the zonal extent of the recirculation cells. As the recirculation cells continuously expand in the zonal direction they allow the successive pitchfork bifurcations to occur. Each pair of new equilibria has an additional meaner in the part of the jet separating the counter rotating recirculation gyres. In contrast to the antisym-metric branch, the energy level for the nonsymmetric branches A and AЈ, B and BЈ, C and CЈ, and D and DЈ decreases or remains nearly constant (Fig. 6) as A b is decreased.
The difference in the energy level maintained by each state is due to the fact that both the energy dissipation and the energy input by the wind stress are functions of the flow field. The energy input by the wind stress is given by the correlation between the curl of the wind stress and the streamfunction field 2 1 gЈ P ϭ Ϫh١ ϫ dx dy.
The energy dissipation due to interfacial drag is given by r 0 ͵͵ 2 2 f 0 and the energy dissipation due to lateral diffusion is given by Table 2 lists the basin integrated energy balance for each solution. There is a 37% difference in the energy level of equilibrium B, which has the lowest energy and equilibrium EЈ, which has the highest energy. The input of VOLUME 32 energy by the wind stress varies by 13% between these two equilibria. The energy dissipation by interfacial friction varies by 15% between equilibrium B and equilibrium EЈ while the energy dissipation by lateral diffusion varies only by 4.7%. The larger relative difference between the interfacial dissipation for equilibria B and EЈ reflects the fact that equilibrium EЈ has a much higher energy level than equilibrium B. The difference in interfacial friction, however, is not so large as the difference in the total energy level. Most of the difference in the energy levels can in fact be attributed to differences in the potential energy while interfacial dissipation is proportional to the kinetic energy. For comparison, Table 2 also gives the energy balance for the linearized model. For this solution, lateral diffusion accounts for more than half the energy dissipation. Since there are no inertial effects for the linearized model, all streamlines pass through the frictional boundary layer. The absence of inertial effects also prevents recirculation cells from forming, thereby eliminating important regions where interfacial friction dissipates energy.

GLOBAL VORTICITY BALANCE
In this section we discuss the global vorticity balance for the different steady-state solutions. Primeau (1998a) presents a discussion of the global vorticity balance for symmetric and nonsymmetric solutions. If the streamline separating the subpolar from the subtropical gyre is not coincident with the line of zero wind stress curl, the circulation can advect negative vorticity into a region of positive wind stress curl and vice versa. We can think of this advection of vorticity as an intergyre flux of vorticity provided we define the gyres to be the regions occupied by the subtropical and subpolar gyres of the linear Munk-like solution. From this point of view, the region occupied by the gyres is fixed and, consequently, the vorticity input by the wind stress curl is also fixed. Tables 3 and 4 give the gyre integrated vorticity budget for the subtropical and subpolar gyres respectively. The advection terms in the vorticity equation cannot generate any vorticity; they act to only redistribute it. Thus any net basin-integrated input of vorticity by the wind must be removed by the explicit fric-tion terms. For A s ϭ 0.05, the net input of vorticity by the wind in the subtropical gyre is Ϫ2.23 ϫ 10 Ϫ3 s Ϫ2 , and for the subpolar gyre it is 2.0177 ϫ 10 Ϫ3 s Ϫ2 . The subtropical gyre receives 5% more vorticity from the wind than the subpolar gyre.
Tables 3 and 4 show that the intergyre flux of vorticity is crucial for equilibria A, B, and BЈ, which are the first to bifurcate. Since these equilibria are farthest in parameter space from their bifurcation points, they are the least antisymmetric. Also, the solutions BЈ, CЈ, and DЈ, which have a jet that first turns south after separating from the western wall, have weaker intergyre fluxes of vorticity than their nearly mirror image counterparts B, C, and D, which have a jet that first turns north. This asymmetry is due to the weaker/stronger vorticity input in the subpolar/subtropical gyre.
The dominant explicit dissipation term in the vorticity equation is the biharmonic viscosity. It generally becomes more important for the more antisymmetric solution, although the relative differences are small compared with the changes in the advection and interfacial friction terms. The sink of vorticity through lateral diffusion is generally stronger for the unprimed solutions. This, along with the weaker intergyre flux of vorticity for the unprimed solutions, increases the importance of interfacial drag for removing the excess vorticity. To compensate for the weaker intergyre flux of vorticity, the more antisymmetric solutions dissipate much more vorticity through interfacial friction than do the more nonsymmetric solutions. For example, interfacial friction is 44% more important for solution EЈ than it is for solution CЈ. It can also be noticed that for the primed solutions, interfacial friction is always stronger than for the unprimed counter parts. This is consistent with the weaker intergyre vorticity flux and weaker lateral diffusion.

Fixed points and time-dependent simulations
In the previous section we have shown the existence of multiple equilibria whose main difference is the degree of elongation of the recirculation cells. Some of these equilibria are remarkably similar to the flows averaged within the high, medium, and low energy regimes identified by McCalpin and Haidvogel (1996) (see their Fig. 4). Compare, for example, their contour plots for the streamfunction averaged within the low, medium, and high energy levels to the fixed points BЈ, CЈ, and EЈ, respectively, shown in Fig. 4. These show similar zonal jet penetration as well as a similar meandering structure. The stability analysis to be presented in section 5 shows that all the fixed points are unstable to small perturbations. Nevertheless, these fixed points might still act to ''steer'' the model trajectory in phase space (Legras and Ghil 1985), in the sense that the model trajectory will at times follow orbits that are close to the stable manifolds of the fixed points before being expelled on orbits that lie close to the fixed point's unstable manifold. If this is the case, the flow field associated with the fixed point, along with its spectrum of unstable modes, will be useful in characterizing the state of the system during certain regimes.
To further investigate this possibility we compute time-dependent solutions by time-stepping the model. The no-normal flow boundary condition, Eq. (4), requires that h ϭ c(t) on the boundary. To obtain c(t), we impose the conservation of mass condition, Eq. (3), at each time step.
If the idea that the model's fixed points act to steer the time-dependent trajectory in phase space is correct, we would expect to see modes of variability associated with structures in phase space that point away from the time-mean state and toward the fixed points. To test this hypothesis, we have conducted a simulation of the corresponding time-dependent model with the standard parameter set given in Table 1 on the unstreched grid with 20-km resolution. The interface height anomaly was saved at 5-day intervals. In Fig. 7 the time-mean interface height anomaly, , is contoured. It is obtained by h averaging the field saved over a period of 1200 years excluding the spinup period. The amount of variability away from this mean state and toward the fixed points EЈ, DЈ, CЈ, and BЈ (those similar to the time-averaged flows within the high, medium, and low energy regimes identified by McCalpin and Haidvogel) was evaluated by projecting the variability onto a set of four orthonormal vectors spanning the directions in phase space that point away from the mean state and toward the four primed fixed points. Approximately 30% of the total variability is captured by the four modes. The amount is significant considering that the system has 24 882 degrees of freedom. Furthermore, most of the variance in the interface height anomaly captured by the four modes, is at low frequencies. Figure 8 shows a plot of the frequency times power density spectrum for the basin integrated variance of the interface height anomaly for the full field and for the field in which the projection onto the span of the four modes has been removed. The plot shows a significant part of the variance associated with periods longer than 1 year project onto the four modes.
A simple comparison between the energy histogram and the energy levels of the fixed points does not agree with the simple idea that existence of fixed points will be associated with peaks in the energy histogram. Figure  9 shows a histogram of the total energy for the 1500yr time series as well as the energy level for each steadystate solution. Except for the general agreement between the order of magnitude of the energy levels of the fixed points and the time dependent trajectory, there is no clear agreement between the energy levels of the steady state solutions and the peaks in the histogram near 3.55 ϫ 10 ϩ17 J (low energy) and 3.95 ϫ 10 ϩ17 J (medium energy). Note however, that very little of the distribution density spreads to energy levels higher than the level of equilibrium EЈ. As we will show below most of the density of high-energy realizations can be attributed to trajectories which tend toward equilibrium EЈ from lowenergy levels.
The high, medium, and low energy regimes define high-dimensional spherical shells centered on the origin [h(x, y) ϭ 0] in phase space. Two points in phase space which have the same energy will lie on the same shell, but they need not be close to each other. To determine if the model trajectory tends toward fixed point EЈ during high energy events and to quantify in a more objective manner the similarity between the fixed points and the time averaged flows within each of the high, medium and low-energy regimes, we computed the distance in phase space between a fixed point and the instantaneous model state. The distance d X (t) between a fixed point X and the model state at time t is given by in which potential energy is used as the norm. Note that the total energy norm would produce essentially the same result since the kinetic energy is an order of magnitude smaller than the potential energy. In Fig. 10, a typical segment of the time series of the VOLUME 32 J O U R N A L O F P H Y S I C A L O C E A N O G R A P H Y FIG. 8. Frequency times power density spectrum for the variance of the interface height anomaly of the total field (upper curve), and for the field in which the part of the variance, which projects onto four vectors that point from the time-averaged state to the fixed points BЈ, CЈ, DЈ, and EЈ has been removed (lower curve). The dashed lines are 95% confidence intervals. TABLE 5. Number of events in which the model stayed in a particular regime (as defined by the proximity to a fixed point) for periods of time between 3 and 5 yr, between 5 and 10 yr, and more than 10 yr.

Regime
3-5 yr 5-10 yr Ͼ10 distances between the model trajectory and the primed fixed points is shown. During persistent high energy events, the model trajectory is close to the fixed point EЈ, as indicated by the broad minimums of d EЈ which coincide with periods of high total energy. In Fig. 11 we plot the histograms of the distributions of the distance to the various fixed points. The relative proximity of the model trajectory to fixed point EЈ can be evaluated by comparing it to the spread of model state in phase space. The smallness of d EЈ during high energy events compared to the spread of the distribution of d EЈ confirms that the model trajectory does, in fact, come close to fixed point EЈ.
As another characterization we used the distance diagnostic to determined the fixed point nearest to the model trajectory for each time. We partition points, p, in phase space into regions defined as follows such that points in R X are closer to fixed point X than to any other fixed point. From this partitioning we found that, in general, the trajectory is nearest to fixed point EЈ during high and medium energy levels, nearest to fixed point DЈ during low and medium energy levels, nearest to fixed point CЈ during the medium and low energy levels, and nearest to fixed point BЈ during low energy levels. There is no one-to-one relationship between the regimes defined by McCalpin and Haidvogel and the proximity of the model trajectory to the fixed point. However, if the cross tabulations are restricted to the most persistent regimes, a simpler picture emerges for the role of fixed point EЈ. In Table 5 we give the number of events for which the model trajectory stayed  Table 6. within R X for a duration of 3-5 yr, 5-10 yr, and longer than 10 yr. From this table we can see that only regime R EЈ persists for periods of time greater than 10 yr. For comparison, Table 6 shows the number of occurrences of the high, medium and low energy regimes that persist for lengths of time between 3 and 5 yr, 5 and 10 yr, and longer than 10 yr. There is an exact correspondence between events where the model state is in the highenergy regime and the events when the model is closest to equilibria EЈ.
In Fig. 12 we plot the square of the minimum distance to the fixed point EЈ for the 23 persistent high-energy events listed in Table 6 as a function of the duration time of the corresponding events. Note that the squared distance plotted on the ordinate is substantially smaller than the spread of energy levels in the histogram of Fig.  9. This is not surprising since the high-energy regimes have weak eddy activity, and it is easy to see that the streamfunction pattern during those events is very similar to the streamfunction pattern for fixed point EЈ. The higher levels of eddy activity during the medium and low energy regimes makes it more difficult to see any connection between the model trajectory and the lower energy fixed points. Furthermore, from Fig. 12, we see that the closer the model trajectory gets to the fixed point the longer the high-energy event persists. This behavior is very similar to that of the Lorenz model (Lorenz 1963) where the amount of time the model spends spiraling around a particular lobe depends on how close the trajectory started from the unstable fixed point at the center of the lobe (Primeau 1998b).

Stability analysis
The stability properties of the fixed point solutions also sheds light on the structure of the model attractor.
To determine the stability of the fixed points, we linearize the governing equation about each of the fixed points and look for modal solutions of the form hЈ(x, y, t) ϭ ⌽(x, y) exp(t). The real part of gives the growth rate of the mode and the imaginary part gives the frequency.
The result of the stability analysis is that all the fixed points we have found are unstable at the standard parameter set given in Table 1. The spectrum of unstable eigenmodes is plotted in Fig. 13 for each fixed point. The higher energy fixed points are generally more stable than the lower energy ones. Most notably the high-en-ergy fixed point EЈ is unstable to only one oscillatory mode with a period of 1.8 yr and a long e-folding time of 2.9 yr. In contrast, the lowest energy regime AЈ is unstable to 12 oscillatory modes and 2 stationary modes. Its most unstable mode has an e-folding time of only 1.4 months. The modal structures for the low-frequency modes (period Ͼ 1.5 yr) are plotted in Fig. 14 The stability properties of the fixed points is qualitatively consistent with the time-dependent behavior described by McCalpin and Haidvogel (1996) in several ways. The low-energy fixed points have more unstable modes and higher growth rates consistent with the fact that the low-energy regime in the time-dependent simulation has the most eddy variability. In contrast, the high-energy fixed point EЈ has only one unstable ei- genmode with a weak growth rate and is thus the least unstable. This is consistent with the time-dependent simulation that has the least eddy variability during periods when the state of the system is in the high-energy regime. The weak growth rate of E's unstable mode might also explain the observed persistence of the highenergy regime.
The patterns (not all shown) associated with the unstable eigenmodes of each fixed point are also broadly consistent with the patterns of time-averaged eddy kinetic energy described by McCalpin and Haidvogel (1996). The unstable eigenmodes associated with AЈ and BЈ (not shown) have their amplitudes concentrated in roughly the same region as that of the strong time-averaged eddy kinetic energy identified for the low-energy regime. Similarly, the amplitudes of the unstable eigenmodes associated with CЈ and DЈ (not shown) are concentrated in roughly the same regions where the medium-energy regime has its strong eddy kinetic energy. The unstable eigenmode associated with EЈ has its amplitude concentrated along the jet axis in the same way that the time-averaged eddy kinetic energy is concentrated along the jet axis during high energy events. Dijkstra and Katsman (1997) and Katsman et al. (2001) have shown in their studies that the low-frequency variability of their double-gyre models can be traced back to a Hopf bifurcation of the steady flow on the first branch to bifurcate at a pitchfork bifurcation (the steady flow equivalent to AЈ in the present study). The mode they identify has an interannual period and they show that this mode is at the origin of the lowfrequency variability. In the present study, we also find several unstable oscillatory modes for the fixed point AЈ, which have interannual periods. Figure 14 shows the modal structures for these modes. For the parameter values used by McCalpin and Haidvogel (1996), the low-frequency unstable eigenmodes of fixed point AЈ have complicated structures with many closed recirculation cells. It would be difficult from these structures alone to predict that the low-frequency variability would involve periods with an elongated jet and weak eddy activity. On the other hand, other fixed points further down the bifurcation tree have an elongated jet. The fixed point EЈ in particular is unstable to only one weakly growing mode (Fig. 15). It thus appears that the high energy regime described by McCalpin and Haidvogel (1996) is associated with the existence of fixed point EЈ. The transition from low to high energy regimes might involve the low-frequency modes of fixed point AЈ, but establishing exactly how this happens is beyond the scope of the present study.

Discussion
The major point of this paper is that unstable steady solutions can be useful in describing and understanding the state of ocean models during different dynamical regimes. We have shown that a reduced-gravity quasigeostrophic model admits steady-state solutions that are very similar to the regimes visited in the time-dependent model. What is most remarkable is that even fixed points that appear relatively far down the bifurcation tree capture the essence of relavent dynamical regimes of the time-dependent trajectory (e.g., fixed point EЈ). As was shown in the stability analysis, the fixed point EЈ that captures the essential global vorticity and energy balances during the high-energy state is unstable to only one oscillatory mode with a relatively weak growth rate despite occuring relatively far down the bifurcation tree. This explains in part the persistent nature of the weakly meandering high energy state. In contrast, the low-energy fixed points which were the first to bifurcate are unstable to many more modes with larger growth rates.
The multiple equilibria that we find are the result of a sequence of symmetry-breaking pitchfork bifurcations the first of which was first identified by Jiang et al. (1995) and Cessi and Ierley (1995). We suggest that the symmetry-breaking pitchfork bifurcations can be interpreted as a stationary Rossby wave. A westward propagating Rossby wave superimposed on the intergyre jet can become stationary if its wavelength is such that its phase speed is equal and opposite to the flow speed. As one reduces the dissipation parameters in the model the recirculation cells expand eastward in order to dissipate the excess vorticity that is no longer dissipated in the western boundary layer. As the recirculation elongates sufficiently to allow an additional meander of the stationary wave to fit within the region with strong eastward flow, a new pitchfork bifurcation occurs and the new pair of nonsymmetric equilibria has one additional meander. Previous bifurcation analysis of the double gyre model used a basin configuration with a larger north-south than east-west extent and only captured one symmetry breaking pitchfork bifurcation before the jet extended completely across the basin (Jiang et al. 1995;Speich et al. 1995;Cessi and Ierley 1995;Dijkstra and Katsman 1997;Primeau 1998a). Because of this the stationary wave nature of the pitchfork bifurcation was not apparent in these studies. Dijkstra and Katsman (1997), for example, explain the physical mechanism that allows the stationary mode to become destabilized and thus lead to the first pitchfork bifurcation. They do not address the balances that allow the mode to be stationary in the presence of advection and background potential vorticity gradients. For the first few pitchfork bifurcations, the stationary wave picture is not so clear, but, as the recirculation cells expand, the balance between westward wave propagation and eastward advection becomes more apparent and the pitchfork bifurcations identified in previous studies are seen to be the first in a sequence of stationary waves with progressively more meanders. We hope to investigate in future work if one can identify the coalescence of oscillatory modes to form stationary modes that then cross the imaginary axis to give rise to the pitchfork bifurcations.
The wider basin also shows that the saddle-node bifurcation leading to multiple antisymmetric equilibria first identified by Ierley and Sheremet (1995) and Cessi and Ierley (1995) occurs at the same parameter value for which the jet reaches the eastern wall. In the narrower basin studies, this coincidence is not so striking, and the possibility that the interaction of the jet with the western wall is responsible for the saddle node bifurcation is not made apparent.
Another way in which the wider basin clarifies the dynamics is that it allows for a transition region with damped Rossby waves in between the recirculation cells and the Sverdrup interior. This Rossby wave field is similar to the solution proposed by Moore (1963) for the structure of the inertial western boundary layer where the Sverdrup flow is eastward. As discussed by Pedlosky (1996), the Moore solution cannot be regarded as a model of the western boundary layer. It should be viewed as a distinct dynamical regime for the region separating the recirculation cells and the Sverdrup interior. Before the present study, the stationary wave field had only been observed in a model with no-slip boundary conditions, because only with such boundary conditions does the recirculation cell remain limited enough to allow for a zonal Sverdrup interior that can support stationary Rossby waves. With a wider basin, the stationary Rossby wave field can exist even with free-slip boundary conditions.
Finally, we point out some interesting similarities between the multiple equilibria found in the simple QG model and the different regimes with elongated and contracted recirculation cells identified by Kelly et al. (1996) and Qiu (2000) from the altimetry observations of the Gulf Stream and Kuroshio Extension systems. We found that the fixed points with elongated recirculation cells, EЈ and DЈ, (Fig. 4) have weak meandering and a more deeply penetrating jet extension. Furthermore, in the time-dependent simulation, we found that when the model trajectory was closest to these equilibria in phase space the eddy kinetic energy was generally lower than at other times in the simulation. A similar type of behavior was found in the above mentioned observations for the elongated state with the more deeply penetrating jet extension. In contrast the fixed points with the more contracted recirculation cells (CЈ, BЈ, and AЈ), had stronger meanders and a more weakly penetrating jet extension. In the time-dependent simulation, flow fields most similar to these fixed points generally had higher levels of eddy kinetic energy. This again is similar to the observations for the state with the contracted recirculation cell. We also found that the primed fixed points, (AЈ, BЈ, CЈ, and DЈ) had a zonal mean jet position that moved progressively southward as the recirculation cells became more contracted. For the case with nonsymmetric wind stress (A s ϭ 0.05) the time-dependent flow field generally remained closer to these fixed points than to their nearly mirror image counterparts. Consistently we found that during periods, when the simulated flow had elongated recirculation cells, the jet extension followed a more northerly path and conversely, when the recirculation cells were contracted, the jet followed a more southerly path. This also is consistent with the observations. Presumably, the nonsymmetries of the real flow (sphericity, ageostrophic effects, slanting coastlines, and nonsymmetric wind stress patterns, etc.) act to give an effective asymmetry parameter that is positive. A negative A s would have given the opposite effect, with the more northerly path associated with the state with the contracted recirculation cells.
The present study shows clearly how very different regimes with either elongated or contracted recirculation cells can achieve energy and vorticity balance without any change in the forcing wind stress. This suggests that intrinsic nonlinear dynamics is a viable candidate to explain the observed low-frequency variability. If the flow fields associated with this variability can produce sustained interannual SST anomalies through their effects on the heat transport divergence and storage, intrinsic ocean variability might be a contributing cause to interannual and decadal climate fluctuations. tinuation code written in FORTRAN. His code was used to compute all the high-resolution, steady-state solutions. Numerous suggestions by Glenn Flierl, Joe Pedlosky, Paola Rizzoli, and Roger Samelson are gratefully acknowledged. The useful coments from the anonymous reviewers are also gratefully acknowledged. Funding for this research was provided by the Office of Naval Research (ONR) under Grant N00014-95-1-0226.