Elongation and Contraction of the Western Boundary Current Extension in a Shallow-Water Model: A Bifurcation Analysis

The double-gyre circulation, formulated in terms of the quasigeostrophic equations, has a symmetry about the basin midlatitude ( y → (cid:1) y , (cid:2) → (cid:1)(cid:2) ), which is absent in a formulation based on the shallow-water equations. As a result, the shallow-water model does not have the pitchfork bifurcation structures that, in the case of the quasigeostrophic model, connect together multiple solution branches with elongated and contracted recirculation gyres. For the shallow-water model, solution branches with elongated recirculation gyres are disconnected, and a one-parameter bifurcation analysis is unable to detect their existence. The deeply penetrating jet solution branches do, however, continue to exist, and can be found using a bifurcation analysis couched in terms of two parameters. An effective pair of parameters is the viscosity and a parameter controlling the symmetry of the wind stress profile. A bifurcation analysis with these parameters reveals the existence of new solution branches that were not found in previous bifurcation analyses of the shallow-water model. The new solutions have a jet extension that penetrates farther eastward and that is more stable than the jet-up and jet-down solutions found in previous studies. Furthermore, the origin of the low-frequency variability at low viscosities is associated with a sequence of bifurcations originating from one of the new steady-state solution branches. In particular, the eigenmode analysis of the new branch reveals that a so-called gyre mode is at the origin of the model’s low-frequency variability at decadal time scales.


Introduction
Satellite-based observations of sea surface height in the Gulf Stream and Kuroshio regions of strong air-sea heat flux are revealing large-scale patterns of variability on interannual to decadal time scales. Several studies of these regions (Vivier et al. 2002;Kelly 2004;Dong and Kelly 2004) have shown the importance of advective heat transport in the interannual to decadal variability of the upper-ocean heat budget. These observational results in regions where current systems are highly inertial and capable of manifesting intrinsic variability are possible examples where internal oceanic variability actively contributes to midlatitude climate variability. One intriguing feature of the Kuroshio and Gulf Stream extension regions is that the dominant mode of lowfrequency variability has been characterized as an oscillation between a state in which the jet extension and associated recirculation gyre are elongated, and one in which they are contracted (Kelly et al. 1996;Qiu 2000).
We show that equilibrium states with elongated and contracted jet extensions can coexist for identical windforcing and dissipation parameters in a shallow-water model of the midlatitude ocean circulation. Our results build on the previous work of Primeau (2002) and Simonnet (2005), who demonstrated the existence of similar multiple equilibria in a quasigeostrophic (QG) model. We investigate the stability of the different equilibria and show that the steady-state bifurcation tree is useful for organizing the behavior of timedependent solutions.
We are motivated to reexamine the bifurcation structure of the double-gyre model because of an important qualitative difference in the penetration scale of the intergyre jet in the steady-state solution branches for the quasigeostrophic equation (QGE) model (Primeau 2002) and those published so far for models based on the primitive equations (PEs) (Jiang et al. 1995;Speich et al.,1995;Schmeits and Dijkstra 2000, 2001Simonnet et al. 2003a,b).  and Dijkstra and Ghil (2005) provide an introduction and a review of bifurcation theory applied to the dynamics of the ocean circulation. In the PE system, increasing the nonlinearity by decreasing the friction or increasing the strength of the wind forcing reveals a perturbed pitchfork bifurcation structure with three solution branches. On these solution branches the intergyre jet-penetration scale increases very little when the forcing is increased beyond its value at the pitchfork. As a result, none of the steady-state solution branches in the above-mentioned PE studies have recirculation gyres that match the eastward extension seen in observations. For example, Schmeits and Dijkstra's (2001) solution for the Pacific basin shows recirculation gyres that are confined to the west of 145°E, while the observations show an inertial jet and recirculation gyre that extends eastward beyond 160°E or approximately 1800 km (Qiu 2000;Qiu and Chen 2005). In contrast, the QGE model forced by symmetric winds possesses a symmetric solution branch whose intergyre jet-penetration scale increases continuously as the nonlinearity is increased. The most energetic solution branches computed by Primeau (2002) have a penetration scale of 1800 km, which is similar to that observed for the Kuroshio Extension.
Because ocean models have a very large number of degrees of freedom, finding steady-state equilibriums in the nonlinear regime is possible only if one has a good starting guess for the equilibrium solutions. Thus, in the absence of any symmetry, we only obtain solutions by continuously varying some parameter. In the QG model formulation (Primeau 2002), the appearance of coexisting equilibrium states with contracted and elongated recirculation gyres happens through a sequence of pitchfork bifurcations as the model nonlinearity is increased. Because all of the solution branches are connected to the symmetric branch at easily found symmetry-breaking pitchfork points, it is possible to detect their existence and compute them. In a primitive equation model, however, this straightforward detection and computation of solution branches is not possible because the existence of the pitchfork bifurcations depends critically on the QG symmetry where is the streamfunction and y is the north-south coordinate (Cessi and Ierley 1995). When the symmetry is perturbed, as is the case for the shallow-water model, the pitchfork bifurcation points are destroyed (Jiang et al. 1995;Dijkstra and Molemaker 1999). The solution branches that were connected at the pitchfork in the symmetric case now appear as disconnected "isolated" branches. Thus, if solution branches with elongated recirculation cells exist in the shallow-water model, they will appear as isolated branches. Their existence would not be detectable in the one-parameter continuation studies of Speich et al. (1995), Schmeits and Dijkstra (2000, 2001, and Simonnet et al. (2003a,b). We nevertheless expect that for a sufficiently small Rossby number the other steady solution branches found in the QG model will also exist in the more accurate shallow-water model; the problem is finding them. Bifurcation theory tells us that solution branches that appear disconnected in a one-parameter setting may in fact be part of the same continuous solution surface when viewed in an appropriate two-parameter setting. Figure 1 shows a schematic diagram of a folded solution surface. A cross section through the surface shows the perturbed pitchfork structure with a disconnected branch. On the horizontal projection plane spanned by two parameters one sees a wedge-shaped curve with a spike known as a cusp point. The curve divides parameter space into two regions. Each point in the smaller wedge-shaped region corresponds to three distinct points on the folded surface. Each point corresponds to a different steady equilibrium of the model. Outside the wedge-shaped region the model has only one steady-state equilibrium solution for each parameter value. Bifurcation theory tells us that the cusp and folds depicted schematically in Fig. 1 are the only generic bifurcations in a two-parameter setting (e.g., Arnold 1992). We can thus expect that the sequence of pitchfork bifurcations found for the QG model should manifest themselves as a sequence of cusp points and folds in an ocean model without the QG symmetry.
Because two parameters need to be varied to locate cusp points, our computations will be couched in a twoparameter space. We introduce a parameter that controls the degree to which the model forcing departs from the QG symmetry, in addition to the parameter that controls the nonlinearity of the model. The pur- pose of this paper is to show by a direct computation using a two-parameter continuation algorithm that the steady-state solution branches with elongated recirculation gyres also exist for the PE model. In the process we will compute the first three cusp points and map out where in parameter space we find solutions with elongated recirculation gyres, and in particular where we find the coexistence of solutions with contracted and elongated recirculation gyres.
A further motivation for our work comes from the recent characterization of the Kuroshio Extension variability by Qiu and Chen (2005) in terms of stable and unstable states. Based on an analysis of sea surface height data, Qiu and Chen (2005) found that for the stable (unstable) state the jet extension is stronger (weaker), has smaller (larger) meanders, and exhibits weaker (stronger) kinetic energy variability. This behavior was considered to be counterintuitive by Qiu and Chen (2005) because conventional baroclinic instability theory for zonal flows would predict that the stronger jet would be more unstable and thus have larger kinetic energy variability. The suggestion that the flows with the weaker jet are less stable than those with the stronger jet is however qualitatively consistent with the linear stability results for the steady-state solution branches of the idealized QGE model (Primeau 2002). There it was found that the states with the elongated recirculation gyres and smaller jet meanders were more stable than the less energetic states with shorter recirculation gyres and larger jet meanders. We are thus encouraged by the qualitative consistency between the observed variability and the stability properties of the steady-state solutions of the QGE model. This consistency suggests that the idealized model behavior is physically relevant.
With the observational results as a backdrop, the goal of the present study is to extend the QGE results to the shallow-water equations. This study contributes an additional step in a hierarchy of models of increasing complexity and realism with which to investigate and interpret the low-frequency variability of the winddriven ocean circulation. In future work we hope to improve the model by considering other deviations from the QG symmetry, such as realistic basin geometry, wind stress, and bottom topography.
The paper is organized as follows: In section 2 we describe the model setup and review the formulation of the shallow-water equations. We also describe the numerical methods used in the paper. The results of the steady-state bifurcation analysis are presented in section 3, along with results from a suite of time-dependent simulations. Conclusions are presented in section 4.

a. Basin geometry and wind stress profile
We study an idealized model of the wind-driven ocean circulation in a rectangular basin with solid walls at x ϭ 0, x ϭ L x , and y ϭ ϪL y /2, y ϭ L y /2, where L x and L y are the east-west and north-south widths of the basin, and forced by a steady zonal wind stress in which 0 is the maximum wind stress, and the parameter A s controls the asymmetry of the wind stress profile and the relative input of vorticity in the subpolar and subtropical gyres. The form in Eq.
(2) is chosen to produce the same wind stress curl used in the QGE studies of Primeau (2002) and McCalpin and Haidvogel (1996). Figure 2 shows profiles of the wind stress and its curl, for several values of A s . For A s ϭ 0, the wind stress curl profile is antisymmetric about y ϭ 0, that is, it satisfies the QG symmetry (1). For A s Ͻ 0, the input of positive vorticity in the subpolar gyre is larger than the input of negative vorticity in the subtropical gyre. Conversely, when A s Ͼ 0, the input of negative vorticity in the subtropical gyre is larger than the input of positive vorticity in the subpolar gyre. As we will see, the parameter A s plays an important role in the two-parameter strategy used to find isolated solution branches.

b. Shallow-water equations on a ␤ plane
The model dynamics are governed by the reducedgravity shallow-water equations In (4), u and are the eastward and northward velocity components in a homogeneous fluid layer of depth H ϩ h, with mean depth H and density 0 . The layer of fluid overlies an infinitely deep lower layer, which is assumed to be at rest. The dynamically active upper layer feels a reduced gravity gЈ in which g, the usual gravity, is reduced by the fractional change in density. The Coriolis parameter f ϭ f o ϩ ␤y varies linearly with latitude. The friction is parameterized as an interfacial Rayleigh drag with coefficient r plus a horizontal bihar- In addition, because of the higher-order diffusion operator, we must impose one more boundary condition. For this we choose Ѩ 3 Ѩx 3 ϭ 0 at x ϭ 0, L x , and

c. Model parameters
The bifurcation analysis results are presented in terms of two dimensionless parameters-the horizontal Ekman number, which is a measure of the viscosity of the flow and the parameter A s in Eqs.
(2) and (3), which controls the asymmetry of the wind stress profile. Typical values for the fixed parameters are given in Table 1.

d. Numerical methods and continuation strategy
We obtain numerical solutions by spatially discretizing the shallow-water equations using second-order fi- nite differences. We denote the resulting discrete system of ordinary differential equations symbolically as where x ϵ [u, v, h] T is the solution vector defined at the mesh points and where ϭ [A s , E h ] are the two parameters we use for the bifurcation analysis. To obtain time-dependent solutions we start with the initial condition x(t 0 ) ϭ [0, 0, 0] T and integrate Eq. (8) forward in time using a third-order Adams-Bashforth time-stepping scheme. In practice, for sufficiently large E h , Eq. (8) is essentially linear and time-dependent trajectories converge to a stable steady-state equilibrium. However, when E h is decreased below some critical threshold, steady-state solutions become unstable and the model trajectory equilibrates to a time-dependent attractor, and a timestepping solver can no longer be used to find steadystate solutions. To continue tracking the steady-state solutions, we set dx/dt to zero in Eq. (8), and solve the resulting system nonlinear algebraic equations using Newton's method. Provided the initial iterate is sufficiently close to a fixed point of Eq. (9), Newton's method will converge quadratically. To make sure that we always have a good initial iterate and to compute solutions over a broad range of model parameters, we use a continuation method (e.g., Seydel 1994). Assuming for the moment that is a scalar (either A s or E h ), Eq. (9) defines x as a function of . We call x() a solution branch. The continuation method for computing solution branches consists of using a previously obtained solution at a point o on the branch as an initial iterate for the Newton's method solver applied at o ϩ d. Repeating the process allows the solver to march along the branch one step of length d at a time. Primeau (1998) or , for example, describe the implementation of the above algorithm. Our strategy to find the multiple solution branches is to start in a viscous regime (large E h ) where nonlinearity is weak and Newton's method, applied to (9), converges for the choice x 0 ϭ 0 as the initial iterate. We then gradually decrease E h using the arc length continuation approach to reach the nonlinear regime, and vary A s to look for turning points where dx/d becomes infinite. We then trace back the locus of such turning points to their organizing center at a cusp point by alternating between E h and A s as the continuation parameter and switching continuation parameter after rounding each turning point. Figure 3 illustrates the strategy in a schematic diagram. For example, the upper boundary, denoted by AC, is obtained by first decreasing the parameter E h until a turning point is  reached. After rounding the turning point we switch to the A s parameter for one continuation step and then switch back to E h until the branch turns again. The process is repeated until the solution surface unfolds at the turning point. A similar strategy can be used for the lower boundary CB in Fig. 3 to trace the locus of turning points away from the cusp.

STABILITY ANALYSIS
To determine the stability properties of equilibrium solutions we perform a linear modal stability analysis. To this end, we substitute the modal anzats x ϭ x s ϩ xЈe t into Eq. (8) and linearize by neglecting terms nonlinear in xЈ to obtain where F x is the Jacobian matrix. Equation (10) is in the form of a matrix eigenvalue problem in which is the eigenvalue and xЈ is the eigenvector. We solve Eq. (10) by using the sparse matrix eigenvalue package ARPACK (Lehoucq et al. 1998).
If the real part of all the eigenvalues is negative, then the equilibrium is stable. If one or more eigenvalues has a positive real part, then the equilibrium is unstable. The eigenvector xЈ, corresponding to the eigenvalue with the most positive real part, gives the spatial pattern of the fastest-growing mode. If the eigenvalue is complex, then the mode is oscillatory; the period of the oscillatory mode is given by 2 times the reciprocal of the imaginary part of the eigenvalue. In that case the pattern associated with the imaginary part of xЈ leads the pattern associated with the real part of xЈ by onequarter period. Parameter values at which an oscillatory mode switches from being stable to being unstable are known as Hopf bifurcation points.
To map out the locus of Hopf bifurcation points in a two-parameter space, we follow a zigzag continuation path on which we compute both the fixed point and its linear stability. The path alternates between continuation steps in the A s direction followed by continuation steps in the E h direction. We switch between E h and A s as continuation parameters every time we cross a Hopf bifurcation point.

Results
We present the results from the steady-state bifurcation analysis in sections 3a-3c, and in section 3d we present results from a suite of time-dependent simulations. We used the strategy outlined in section 2d to find steady-state solution branches for our model, and the results are first presented in terms of a two-param-eter regime chart in section 3a. The two-parameter regime diagram shows how solution branches, which would appear disconnected in a one-parameter setting, are in fact part of a continuous folded surface. A onedimensional bifurcation diagram is then presented in section 3b. The bifurcation diagram is the slice through the two-parameter regime chart at A s ϭ 0. It shows the "disconnected" solution branches and is used to help organize the behavior of a suite of time-dependent solutions presented in section 3d.

a. Two-parameter regime chart
Using our two-parameter continuation algorithm (section 2d), we traced six folds in the solution surface back to their organizing centers at three cusp points. Figure 4 shows the E h Ϫ A s plane regime chart, indicating the location of these three cusp points. Also shown in the figure are the locations of the pairs of folds (loci of turning points) that emanate from the cusp points in the direction of decreasing E h to form three expanding wedge-shaped regions. At the cusp points the solution surface folds and introduces two additional equilibrium solution surfaces within the wedge-shaped region.
The three cusp points shown in Fig. 4 organize solution branches according to their jet-penetration scale. They are the generic analogs for a model without the QG symmetry of the first three symmetry-breaking pitchfork bifurcation points found by Primeau (2002) for the QGE model. If one fixes A s to the value at a cusp point and decreases E h , the folding of the solution surface at a cusp point appears as the merging of three solution branches. Table 2 lists the position of the first three cusp points in the parameter space spanned by E h and A s . The table also lists the position of the first three pitchfork bifurcation points computed for the QGE model from Primeau (2002). Close to the cusp point the flow fields on each branch are nearly identical, and the merging branches approximately satisfy the QG symmetry [cf. Eq. (1)]. Figure 5 shows the equilibrium flow field for the solutions close to each cusp point. As we move away from the cusp point by decreasing E h , the distinct asymmetry pattern of the flow field associated with each solution surface becomes increasingly pronounced. This can be seen in Fig. 6, where we plot the steady-state solutions from each of the seven solution surfaces at E h ϭ 2.0 ϫ 10 Ϫ11 and A s ϭ 0. Because this point is situated within the three overlapping wedgeshaped regions shown in Fig. 4, there are seven equilibrium solutions. The naming convention for the solution branches is chosen to be in accord with the convention that Primeau (2002) used for the QG model: solutions are labeled in pairs starting with A and AЈ for the two solutions whose jet extension has the weakest eastward penetration. The primed solutions have the jet extension turning immediately south after separating from the coast, whereas in the unprimed solutions the jet first turns northward. The labels progress down the alphabet for pairs of branches that penetrate farther eastward and have an additional half meander in the jet extension.
Comparing the steady-state thickness field at the cusp point (Fig. 5) to the thickness field at a lowerviscosity value (Fig. 6), we see that the jet-penetration scale for the equilibrium solution pair A and AЈ changes little beyond its value at the C 1 cusp. Similarly, the penetration scale for the equilibrium solution pair B and BЈ does not change beyond its value at the C 2 cusp point. The three cusp points plotted in Fig. 4, together with their corresponding wedge-shaped regions, orga-nize parameter space into regions that admit circulation fields with jet extensions that have distinct penetration scales. Despite the lower viscosity at E h ϭ 2 ϫ 10 Ϫ11 compared to that at cusp points C 1 , C 2 , and C 3 , the jet is decreased. The loci of turning points are labeled T 1 to T 6 and the cusp points are C 1 , C 2 , and C 3 . The shaded areas in the regime diagram indicate regions of parameter space where at least one steady-state solution is stable. The thin lines indicate the position of the loci of Hopf bifurcation points. For comparison with the cusp points, the locations of the first three symmetry-breaking pitchfork bifurcations for the QG system are indicated with black squares. See text in section 3a for explanation of the arrows. The particular solution type that one obtains if E h is decreased with A s held fixed depends on the chosen starting value of A s . If one starts from the viscous regime (large E h ) and continuously decreases E h with fixed A s , then one will end up either on the A branch or the AЈ branch, depending on whether the fixed value of A s is smaller or greater than the A s value at the C 1 cusp point. We indicate this in Fig. 4 by the small arrows with the A and AЈ labels. Similarly, the pairs of labels (primed or unprimed letters) at the end of the small arrows pointing away from the locus of turning in Fig.  4 indicate the solution type one would obtain if E h was decreased with fixed A s starting on the new branches that come into existence at the turning points. These solution branches are generally not accessible by a oneparameter continuation algorithm.

b. Bifurcation plot
In the previous section we have shown how a twoparameter continuation method can be used to find multiple steady-state solution branches, some of which have elongated recirculation gyres and some of which have contracted recirculation gyres. Here we show how these solution branches appear as disconnected branches in a one-parameter bifurcation plot. In the next section we use the one-parameter bifurcation plot to organize results from a suite of time-dependent model simulations with decreasing E h .
The bifurcation plot in Fig. 7 shows the basinintegrated potential energy plotted as a function of E h for all the solution branches we found for the case of A s ϭ 0. For reference, the bifurcation plot corresponds to a one-dimensional cut in Fig. 4 along the line A s ϭ 0. The three turning points shown in Fig. 7 are located in Fig. 4 at the intersection of the line A s ϭ 0 and the three loci of turning points labeled T 1 , T 3 , and T 5 .
For sufficiently large E h there is a unique steady-state solution, but as E h is decreased, pairs of new steadystate equilibriums come into existence at the turning points. Three primary turning points labeled T 1 , T 3 , and T 5 are shown in the bifurcation plot. Solution branches C and D merge at T 5 , solution branches BЈ and CЈ merge at T 3 , and solution branches A and B merge at T 1 . To the right of T 1 only the AЈ solution branch exists.
We performed a linear stability analysis for each solution along the branches shown in Fig. 7 Fig. 7. The location of Hopf bifurcation points as a function of both A s and E h is shown in Fig. 4, in which the shaded region indicates regions of parameter space where we found at least one stable steady-state equilibrium. For A s ϭ 0, and for the range of parameters indicated in Fig. 7, we did not find any cases for which multiple stable steady-state solutions coexisted. We did not attempt to continuously track attractors other than steady-state solutions as a function of E h because this is technically challenging and beyond the scope of the present work. We did however find several examples where distinct attractors (neither of which are steady) coexisted for the same parameter values. These are discussed briefly in section 3d.

c. Summary of steady-state bifurcation analysis
We have shown that the two-parameter continuation strategy, with one parameter controlling nonlinearity and one controlling the asymmetry of the wind forcing, is an effective way of finding and computing solution branches that would appear disconnected in a oneparameter space. In the process, we computed the first three cusp points that play the role; in the case of the shallow-water equation (SWE) model, we computed the role played for the symmetric QGE model by the first three pitchfork bifurcation points. These cusp points organize in parameter space solutions whose jet extensions are successively more elongated.

d. Time-dependent solutions
To explore the time-dependent behavior of the model we computed a suite of time-dependent simulations with the wind profile asymmetry parameter fixed at A s ϭ 0 while varying the biharmonic viscosity from E h ϭ 3.5 ϫ 10 Ϫ10 to 7 ϫ 10 Ϫ12 . All runs were started from rest and spun up until the model trajectory was beyond a transient period and, depending on E h , was either a steady-state attractor, a limit cycle or torus, or a chaotic attractor. FIG. 7. Bifurcation plot as a function of the horizontal Ekman number, E h . The thin solid lines indicate steady-state solution branches that are linearly unstable. The thick solid lines indicate stable steady-state solutions. The small black dots indicate the time-average of the basin-integrated potential energy for a sequence of time-dependent model simulations that were initiated from a state of rest. The gray shading indicates the range of the potential energy time series excluding the spin-up period for each run. There is a critical value near E h Ϸ 5 ϫ 10 Ϫ11 that separates time-dependent flows whose time average is similar to the steady-state flow on the AЈ branch or to the CЈ branch.
In terms of their time-averaged flows, all of our simulations can be classified into two regimes separated by a critical value of E h Ϸ 5 ϫ 10 Ϫ11 . In the higherviscosity regime, the time-averaged flow has many similarities with the flow of steady-state solutions on the AЈ branch, while in the lower-viscosity regime, the timeaveraged flow is similar to the flow of the steady-state solutions on the CЈ branch. This is illustrated in Fig. 8 in which we show plots of the steady-state flows for the AЈ and CЈ branches at E h ϭ 8.38 ϫ 10 Ϫ11 and E h ϭ 2.15 ϫ 10 Ϫ11 , together with the time average of the corresponding time-dependent simulations. The upper panels show the steady-state solution and the lower panels show the time-averaged flow. The standard deviation of the interface height fluctuations from the timedependent simulation is indicated by the gray shading. To show the position of the eddy field relative to both the steady-state and time-averaged flow, we have overlaid the standard deviation pattern onto both the steady-state and time-averaged flow. At both parameter values the time-dependent flow has a vigorous eddy field, but time averaging reveals a flow pattern that is strikingly similar to the steady-state equilibrium. In the high-viscosity regime (E h Ͼ 5 ϫ 10 Ϫ11 ), the timeaveraged flow and the AЈ steady-state solutions have a short jet extension that turns southward immediately after separating from the coast. In the low-viscosity regime (E h Ͻ 5 ϫ 10 Ϫ11 ), the time-averaged flow and the CЈ steady-state solution both have a jet extension that penetrates deep into the basin interior. Furthermore, the meandering pattern of the jet extension is similar for both the steady-state and time-averaged flows: the jet extension turns slightly southward immediately after separating from the coast, and then executes one and a half meanders before fanning out to connect to the Sverdrup interior. At the jet exit region, the jet extension is oriented slightly southward for both the timeaveraged flow and the steady-state solution. The meandering of the jet can be understood as a stationary Rossby wave confined to the jet axis (Primeau 2002;Simonnet 2005).
To further summarize the time-dependent behavior of the model and put it in the context of our steadystate bifurcation analysis, we computed the mean and range of the globally integrated potential energy for each time-dependent model run and overlaid the result on the bifurcation plot in Fig. 7. The time mean potential energy is indicated by the small black dots, while the range of the potential energy time series is indicated by the gray shading. When computing the mean and range of the energy, we excluded the spinup period that, depending on the value of E h , lasted between 30 and 250 yr.
As a function of increasing E h we found two separate sequences of bifurcations from steady to chaotic flow: one at high viscosity that originates from the Hopf bifurcation of the AЈ branch, and a second at a lower viscosity that originates from the Hopf bifurcation of the CЈ branch. This can be seen in Fig. 7, where the critical value E h Ϸ 5 ϫ 10 Ϫ11 separates black dots that roughly track the AЈ at higher E h values and those that roughly track the CЈ branch at lower E h values.
Parameter values with a large potential energy range typically corresponded to chaotic flows with significant eddy variability. For reference, the flow shown in the left panel of Fig. 8, corresponding to E h Ϸ 8.38 ϫ 10 Ϫ11 , was highly chaotic and showed large eddy variability. The flow shown in the right panel of Fig. 8 corresponded to the early onset of decadal chaotic variability for the flow in the CЈ regime.
At large E h the time-dependent simulations converge to a steady state on the AЈ branch. Decreasing E h produced a sequence of bifurcations in which the model behavior transitioned from steady to highly chaotic flow. With a further decrease of E h beyond the critical value of E h Ϸ 5 ϫ 10 Ϫ11 , the model trajectory returns to a steady state, but on a new equilibrium branch. Without the guidance of the steady-state bifurcation analysis, the sudden transition from time-dependent flow back to steady, stable flow, when E h is decreased beyond E h Ϸ 5 ϫ 10 Ϫ11 , would be completely unexpected because decreasing the viscosity usually leads to more unstable flows. With the steady-state bifurcation results in hand, however, the return to a steady, stable flow can be anticipated by the existence of stable part of the CЈ branch. Eventually, the CЈ solution branch also becomes unstable as E h is decreased, and the model then undergoes a second series of bifurcations that culminate in chaotic turbulent flow.
In the CЈ regime, the flow first becomes time dependent through a Hopf bifurcation, labeled H CЈ in Fig. 7, at E h Ϸ 3.9 ϫ 10 Ϫ11 . The Hopf mode has a period of approximately 1.5 yr. Figure 9 shows the basic state together with the real and imaginary parts of layer thickness anomaly for the Hopf mode. The time evolution of the pattern is given by h(t) ϭ cos(⍀t)h R Ϫ sin(⍀t)h I , where ⍀ is the frequency of the mode, so that the imaginary pattern leads the real pattern by a quarter period. At E h Ϸ 2.5 ϫ 10 Ϫ11 , the model trajectory switches to a nearby periodic attractor. The period of the new attractor is 2/⍀ 1 ϭ 6.6 months. Figure 10 shows a sequence of power spectral density plots for the standardized time series of the basin-integrated potential energy for a sequence of runs with decreasing E h . At E h Ϸ 2.3 ϫ 10 Ϫ11 , a second oscillatory mode with a period 2/⍀ 2 ϭ 1.4 yr becomes excited. A third mode with a period 2/⍀ 3 ϭ 15 yr is visible in the power spectrum for E h ϭ 2.21 ϫ 10 Ϫ11 . This third mode is only weakly excited, and by E h Ϸ 2.15 ϫ 10 Ϫ11 it is no longer detectable when yet another mode with a period 2/⍀ 4 ϭ 10 yr becomes excited. With the excitation of the fourth mode the flow becomes chaotic, as shown by the broadening of the spectrum for E h Ϸ 2.15 ϫ 10 Ϫ11 . The filling in of the spectrum continues as E h is decreased further, but the peaks of the two primary modes with frequencies ⍀ 1 and ⍀ 2 remain visible.
The spatial patterns associated with the oscillatory modes are shown in Fig. 11 for the E h Ϸ 2.15 ϫ 10 Ϫ11 FIG. 9. Basic state together with the real and imaginary part of the first Hopf mode for the CЈ branch at E h Ϸ 3.9 ϫ 10 Ϫ11 . The period of the mode is 2/⍀ 1 ϭ 1.5 yr.
case. At this value of E h most of the interface height variability is associated with the mode with a period of 6.6 months. This mode is characterized by a train of eddies that propagate westward along the northern edge of the northern recirculation gyre. The mode with the period of 10 yr corresponds to a small amplitude meandering of the jet extension, which causes the easternmost tip of the jet to change orientation. A linear FIG. 10. Power spectral density for the standardized time series of the basin-integrated potential energy for a sequence of runs with decreasing E h in the region of parameter space where the flows of the CЈ regime make a transition from periodic to chaotic flow. The periods of the fundamental peaks are 2/⍀ 1 ϭ 6.6 months, 2/⍀ 2 ϭ 1.4 yr, 2/⍀ 3 ϭ 15 yr, and 2/⍀ 4 ϭ 10 yr. stability analysis of the nearby CЈ steady-state solution at the same E h value revealed an unstable mode with a period of 10.6 yr and a spatial pattern that is very similar to the pattern associated with the 10-yr oscillation in the time-dependent simulation (Fig. 12). The spatial pattern of the eigenmode is also very similar to the so-called one-bump gyre mode computed by Simonnet (2005) for the quasigeostrophic model. We therefore speculate that the origin of the peak in decadal variability is due to the destabilization of the gyre mode. FIG. 11. Oscillatory modes corresponding to the power spectral peaks at three different frequencies for the case of E h ϭ 2.15 ϫ 10 Ϫ11 overlaid on the time-averaged state. The contours indicated the interface height anomaly. The contour interval is 20 m for the mean state indicated by the thin contours. The contour interval for the modes indicated by the thick contours is 10 m for the mode with frequency ⍀ 1 ϭ 2/6.6 months, 0.5 m for the mode with ⍀ 2 ϭ 2/1.4 yr, and 2 m for the mode with ⍀ 4 ϭ 2/10 yr. The patterns in the left panels lead the patterns in the right panels by one quarter period.
Confirmation of this hypothesis would require Floquet analysis of the oscillatory solution, which is not at present computationally feasible. The origin of the lowfrequency gyre mode in the quasigeostrophic model is due to the merging of two purely real modes to produce a complex conjugate pair of oscillatory modes. We tracked the pair of oscillatory eigenmodes as a function of decreasing E h and did find that its period does increase. However, we were unable to find the actual merging point before reaching the turning point T 3 . We hypothesize that both parameters E h and A s need to be varied to find the merging point, but leave this for future work.
In summary, we have shown that the existence of a new steady-state solution branch (CЈ) at low viscosity with an elongated jet extension is relevant to understanding the intrinsic low-frequency variability of the wind-driven ocean circulation in the low-viscosity limit.

Summary and conclusions
We have carried out a bifurcation analysis for the SWE model in a two-parameter space and showed that the two-parameter bifurcation analysis allows us to find multiple steady-state solution branches with elongated and contracted recirculation gyres. The solution branches with flows that have a jet extension that penetrates far eastward into the interior of the basin appear as disconnected branches in a one-parameter setting. These solution branches are not accessible when a bifurcation analysis is performed with only one parameter, and thus had not been found in previous bifurcation studies of the shallow-water model.
We also found that one of the new disconnected steady-state solution branches is stable in a parameter regime where other branches are unstable. The model thus exhibits two transitions from steady to chaotic flow-one sequence occurs in the high-viscosity regime FIG. 12. (top) The oscillatory mode associated with the power-spectral peak with frequency ⍀ 4 ϭ 2/10 yr for the case E h ϭ 2.15 ϫ 10 Ϫ11 . (bottom left) The imaginary and (bottom right) real parts of the only unstable eigenmode for the steady state on the nearby CЈ branch at E h ϭ 2.15 ϫ 10 Ϫ11 . The unstable eigenmode has a period of 10.6 yr and a growth rate with an e-folding time scale of 5.1 yr. The contours indicate the interface height anomaly. The amplitudes of the top and bottom modes have been normalized. The patterns in the left panels lead the patterns in the right panels by one quarter period. and a second sequence occurs in the lower-viscosity regime.
The first sequence of bifurcations from steady to chaotic flow is associated with the higher-viscosity branch. This sequence has been discussed in several previous studies (Jiang et al. 1995;Speich et al. 1995;Schmeits and Dijkstra 2000, 2001Simonnet et al. 2003b). The second sequence of bifurcations from steady to chaotic flow associated with the lower-viscosity branch is completely new and has not been discussed in previous studies. We believe that it is this second sequence of bifurcations that is most relevant to understanding the intrinsic low-frequency dynamics of the wind-driven circulation, because it occurs in the physically more relevant low-viscosity regime and has a jet-extension penetration scale that better matches observations. The origin of low-frequency variability in our double-gyre model at sufficiently low viscosity can thus be traced back to a sequence of bifurcations originating from the stable part of the lower-viscosity steady-state solution branch. We wish to emphasize that the lower-viscosity branch corresponds to a flow with elongated recirculation gyres and is one of the new disconnected steadystate solution branches that exists only for sufficiently low viscosity. In particular, the origin of the decadal variability of the model appears to be due to the destabilization of the one-bump gyre mode that Simonnet (2005) found in a quasigeostrophic model. In future work it would be interesting to investigate if at an even lower viscosity there is any symmetry increasing bifurcations and crisis associated with the merging of different attractors associated with the different steady-state solution branches. We anticipate that our mapping of the steady-state bifurcation structures will provide a starting point for such investigations.
The two-parameter continuation method in which one of the parameters controls the symmetry of the forcing can be applied to models with realistic basin geometry and wind forcing. Our particular choice of parameter to control the symmetry of the wind stress as an artifice for finding the multiple solution branches is not unique. We have also tested the method by rotating the direction of the sinusoidal wind stress profile and easily obtained each solution branch with a different jet-penetration scale. Other choices are also possible. Once the multiple solution branches have been found for idealized wind forcing, one can introduce a homotopy parameter ␣ to continuously follow each solution branch from the idealized wind-forcing case to the realistic wind-forcing case, ͑ x , y ͒ ϭ ͑1 Ϫ ␣͒͑ x , y ͒ ideal ϩ ␣͑ x , y ͒ real , ͑11͒ as ␣ is increased from 0 to 1. Such an approach will allow the tools of bifurcation analysis to be applied to more realistic models. Another interesting aspect of our study is the existence of regions of parameter space where our model can, depending on the initial condition, exhibit either large-amplitude chaotic behavior or steady, stable flow. The more stable solutions were associated with flows with a stronger, more deeply penetrating jet extension with weaker meanders. This behavior is qualitatively similar to observations of the Kuroshio Extension system. Qiu and Chen (2005) characterized the Kuroshio Extension in terms of a stable and an unstable regime, with the stable regime corresponding to a stronger jet extension with small meanders and the unstable regime corresponding to the weaker jet with stronger meanders. We conjecture that the behavior of the Kuroshio Extension system can be explained in terms of the existence of multiple attractors with different barotropic stability properties. Our model was forced with steady winds so that the time-dependent trajectory of our model did not show spontaneous transitions between regimes. We speculate that with some time dependence in the wind forcing the model might transition between the stable and unstable states.