MULTIPLE EQUILIBRIA AND STABILITY OF THE NORTH-ATLANTIC WIND-DRIVEN OCEAN CIRCULATION M. JEROEN MOLEMAKER

AND HENK A. DIJKSTRA ∗

Abstract. The large scale ocean circulation is an important component of the global climate system and controls much of its low frequency variability. Evidence is growing that different mean circulation patterns in the North-Atlantic ocean are possible under the same forcing conditions. In this paper, the question of multiple equilibria is considered within an idealized finite element ocean model of the wind-driven ocean circulation in the North-Atlantic. Using pseudo-arclength continuation, branches of steady states are followed in parameter space and their stability is determined by solving the associated linear stability problem. An overview of the numerical methods to handle the large dimensional dynamical systems and their implementation (our code BAGELS) is given. Multiple circulation patterns are found and their existence is shown to be related to internal ocean dynamics rather than to the particular continental geometry. Temporal variability of the flows on intra seasonal time scales is shown to be related to only a small set of eigenmodes to which the steady flows are unstable.

1. Introduction. On the large scale, the ocean circulation is driven at the ocean-atmosphere interface by wind stress and by fluxes of heat and fresh water. The circulation driven by the latter component is called the thermohaline ocean circulation. Recently, much progress has been made in advancing our understanding of the role of the oceans in the present climate. Due to its large heat capacity, the oceans provide the main source of memory of the climate system. The surface, mainly wind-driven, circulation determines the structure of the sea surface temperature and is involved in variability up to decadal time scale. The thermohaline circulation, involving high latitude sinking and distributed upwelling and mixing at lower latitudes, is heavily involved in the longer time scales of climate variability. However, we remain largely ignorant of the processes controlling such variability and much remains to be done to achieve some skill in predictability. There are indications that the North-Atlantic ocean circulation has been different in the past with corresponding different climate states [4]. Modeling as well as observational studies have provided evidence of the existence of multiple circulation patterns under similar forcing conditions [5]. The stability of the present large scale ocean circulation has therefore become an important issue in climate research. Although many large scale numerical ocean models have been developed over the last decades, most are not very suited to address the stability question in detail. A systematic way to study possible patterns of the ocean circulation and internal ocean variability is to trace stationary solutions from the low ∗ Institute for Marine and Atmospheric Research, University of Utrecht, Princetonplein 5, 3584 CC Utrecht, the Netherlands

1

2

M. JEROEN MOLEMAKER AND HENK A. DIJKSTRA

forcing limit towards realistic forcing regimes with a large internal variability. Primary transitions to time-dependent behavior (e.g. Hopf bifurcations) often set the dominant time-scales of variability. The direct computation of equilibria and their linear stability of ocean models leads to large dimensional (typically O(105 ) degrees of freedom) systems of nonlinear algebraic equations and large generalized eigenvalue problems. A crucial step for performing these calculations is an efficient solution of large non-symmetric linear systems. In this paper we use our code BAGELS to study the structure of equilibria of the wind-driven circulation in the North Atlantic. This code is a highly improved version of that used in [12] which combines pseudoarclength continuation, Newton-Raphson iteration and the Simultaneous Iteration Technique. To solve the linear systems of equations, an efficient version of a preconditioned conjugate gradient method [30] is used. This leads to a very efficient method to compute steady non-parallel flows and their linear stability in parameter space [15] of which the main aspects are given below. A central feature in the North-Atlantic circulation is the Gulf Stream, which has fascinated physical oceanographers since its early description by Benjamin Franklin and Timothy Folger [24]. From the enormous amount of data now available, for example through altimeter satellite data, the timemean path of the Gulf Stream (Fig. 1) is quite well known [1, 22]. The southern part of the Stream (the Florida Current) flows almost parallel to the coastline. At Cape Hatteras, the Gulf Stream leaves the NorthAmerican continent and moves further eastward along 40 W. It is accompanied by recirculation regions to the north and the south [19]. It appears that the position of the separation point is quite stable, changing less than about 50 km over several years. Further to the east, the boundaries of the Gulf Stream display significant annual and inter annual variations, with about equal magnitude [1]. Further in the open ocean, the jet spreads out due to meandering and displays an enormous amount of variability due to instabilities of the main current. In this area, eddies are formed which move away from the mainstream generally in westward/south-westward direction. Their average wavelength is about 100 km and propagation speeds of these eddies are about 10 km/day. The separation behavior of the Gulf Stream is a serious problem in physical oceanography. Most ocean models have trouble simulating the correct mean path. Many mechanisms have been suggested as the dominant cause of the separation, for example the wind stress shape, stratification and the continental geometry. However, several of these processes seem to be involved and the detailed physics of the separation remains unclear. In constant density models with a flat bottom, Dengg [10] has stressed the importance of the shape of the wind stress and the no-slip conditions on the continents for obtaining a correct separation. Using a sharp convex corner, he showed that the modeled Gulf Stream may separate correctly

MULTIPLE EQUILIBRIA OF OCEAN CIRCULATION

3

Fig. 1. (a) Geography and bathymetry of Gulf Stream region (after [11]). (b) Sketch of the near surface circulation. Bold lines: Florida Current (FC) and Gulf Stream (GS), branching into the North Atlantic Current (NAC) and Azores Current (AC).

through inertial overshoot. Interestingly, he also finds two different separation states, albeit at different wind forcing amplitude. The latter result motivates us to explore the steady state solution structure within a shallow water ocean model with realistic continental geometry. The formulation of the model is given in section 2, together with a recapitulation of the numerical techniques that are used. If friction is small enough, multiple equilibria are found within this model with different separation behavior of the Gulf Stream (section 3). In the same section, oscillatory instabilities of the steady states are shown to lead to intra seasonal variability of the Gulf Stream. In the last section, the results are summarized and discussed. 2. Formulation. 2.1. The ocean model. Consider an ocean basin Ω with an arbitrary horizontal boundary Γ. A locally flat, Cartesian coordinate system is used in which the only effect of the earth’s sphericity is the variation of the Coriolis parameter f = f0 + β0 y around a central latitude set by f0 ; this is called the β-plane approximation [23]. Vertically, an active layer of mean depth H with density ρ is situated above a very deep layer having a slightly

4

M. JEROEN MOLEMAKER AND HENK A. DIJKSTRA

larger density ρ + ∆ρ and which is supposed to be motionless. The flow is driven by a wind stress (τ x , τ y ) with characteristic amplitude τ0 and spatial patterns (f x , f y ) in the zonal and meridional direction. The shallow water equations governing the flow in the active layer are given by [23]

(2.1)

∂h τx ∂u + u · ∇u − f v = −g ′ + A∇2 u − Ru + ∂t ∂x ρh y ∂v ∂h τ + u · ∇v + f u = −g ′ + A∇2 v − Rv + ∂t ∂y ρh ∂(hu) ∂(hv) ∂h =− − ∂t ∂x ∂y

The velocities in the eastward and northward directions are u and v respectively, h is the thickness of the upper layer and g ′ is the (reduced) gravity (g ′ = g∆ρ/ρ). Lateral and bottom friction [23] are represented by coefficients A and R, respectively. Let L be a typical horizontal spatial dimension and U a typical horizontal velocity of the flow. The governing equations are non-dimensionalized using scales L, H, U , L/U and τ0 for length, layer depth, velocity, time and wind stress, respectively. The system of equations (2.1) becomes   ∂h αf x ∂u + u · ∇u − (1 + ǫβy)v = −ǫF + E∇2 u − µu + ǫ ∂t ∂x h   ∂h αf y ∂v (2.2) ǫ + u · ∇v + (1 + ǫβy)u = −ǫF + E∇2 v − µv + ∂t ∂y h ∂(hu) ∂(hv) ∂h =− − ∂t ∂x ∂y Several dimensionless parameters appear in these equations: the Rossby number ǫ = U/(f0 L), the Froude number F = g ′ H/U 2 , the Ekman number E = A/(f0 L2 ), the Rayleigh friction coefficient µ = R/f0 and the wind stress coefficient α = τ0 /(f0 HU ρ). Standard values of the parameters are listed in Table 1. These values correspond to the following dimensional values: L = 106 m, H = 500 m, f0 = 0.5 × 10−4 s−1 , ρ = 103 kgm−3 , g ′ = 0.03 ms−2 , τ0 = 0.05 N m−2 , R = 5 10−8 s−1 , A = 200 m2 s−1 and U = 1 ms−1 . The Ekman number E is used as control parameter since its value is very poorly known. On the boundary of the domain Γ, no-slip conditions are prescribed, i.e. (2.3)

xεΓ:

u=0

To calculate a steady state solution of the equations (2.2 − 2.3) an extra condition for h is required to regularize the equations, because h is determined up to an additive constant. The correct condition to use is

MULTIPLE EQUILIBRIA OF OCEAN CIRCULATION

5

Table 1

Standard values of non-dimensional parameters

Parameter ǫ µ E β F α

Value 0.2 × 10−1 1.0 × 10−3 4.0 × 10−6 2.0 × 101 1.55 × 101 2.0 × 10−3

the integral mass balance over the upper layer, which becomes an integral condition for h over the domain Ω, i.e. (2.4)

Z

h dxdy = |Ω|

Ω

where |Ω| is the (dimensionless) area of the domain. In models that integrate the equations in time, this regularization problem is absent, since the integral of the layer depth is set by the initial conditions. 2.2. Numerical methods. 2.2.1. Discretisation. To obtain both an accurate representation of the geometry of the North-Atlantic basin and a high accuracy at positions of large gradients in velocity, a finite element (FEM) approach [8, 9] is chosen to discretize the partial differential equations (2.2). The FEM is widely used in engineering problems where complicated geometries are involved. In oceanography, however, most ocean models use finite differences in combination with a land-sea masking to incorporate continent geometries. Triangular Taylor-Hood elements are used that employ quadratic interpolation functions for the velocities and linear interpolation functions for the layer depth [29]. On such an element, six nodal points for the velocities and three nodal points for the layer depth are defined. This results in approximations of the velocities that are third order accurate with the element size whereas the layer depth is second order accurate. Since for adjacent elements, the shared side has the same nodal points and the same linear approximation, the layer depth is continuous on the boundaries between elements and therefore continuous over the whole flow domain. The necessity to have approximations for the velocities that are one order higher than the approximations for the layer depth arises from a solvability condition that is known as the Brezzi-Babuska condition [3]. A standard Galerkin approach and a Newton-Cotes integration rule over each element is employed for obtaining the discretized problem.

6

M. JEROEN MOLEMAKER AND HENK A. DIJKSTRA

2.2.2. Steady states. To obtain steady states, a (large dimensional) system of nonlinear coupled algebraic equations of the form (2.5)

F(x, p) = 0

must be solved. Here x is the solution vector and p the vector of parameters. Solutions of the equations (2.5) are traced in parameter space using a pseudo-arclength method [21]. With λ as control parameter, branches (x(s), λ(s)) are parametrized by an arclength parameter s. An additional equation follows from the normalization of the tangent along the branch, i.e. (with the dot denoting differentiation to s), (2.6)

x˙ T0 (x − x0 ) + λ˙ 0 (λ − λ0 ) − ∆s = 0

In the equation above, (x0 , λ0 ) is an analytically known starting solution or a previously computed point on a particular branch and ∆s is the step length. To solve the system of equations (2.5, 2.6) Euler-Newton continuation is used [15, 21]. To monitor simple bifurcation points on a particular branch, several indicator functions are used. For example, limit points are detected by ˙ For transcritical or pitchfork bifurcations the sign of the defollowing λ. terminant of the Jacobian matrix may be used. When this determinant is not readily available, a family of test functions τpq can be used following [26]. For the computation of each test function τpq , one additional linear system has to be solved. Other singularities, like Hopf bifurcation points, must be detected by solving the linear stability problem. 2.2.3. Linear stability. The linear stability analysis of a computed steady state amounts to solving a generalized eigenvalue problem of the form (2.7)

Ax = σBx

where A is a non-singular, non-symmetric matrix. Through boundary conditions and/or the continuity constraint, B may become singular. Traditional eigenvalue solvers (e.g. the QZ algorithm [17]) which determine all eigenvalues and, if desired, all eigenvectors are impossible to use since A has a very large dimension. However, in many hydrodynamic stability problems, the instability of a steady flow occurs only through a small number of modes. One is therefore only interested in computing the most dangerous eigenmodes, i.e. those with eigenvalues closest to the imaginary axis. This has motivated the development of several specific algorithms for this task (e.g. [7, 16]), of which we have implemented a variant. The first step of the eigenvalue solver is the application of a Cayley transformation (2.8)

σ =b+a

ν −1 ν +1

MULTIPLE EQUILIBRIA OF OCEAN CIRCULATION

7

where b ∈ R, a ∈ R+ . The parameter b introduces a shift of the spectrum over the real axis, whereas the parameter a stretches the spectrum. The left complex plane ℜ(σ − b) < 0 is mapped within the unit circle |ν| = 1 and ℜ(σ − b) = 0 is mapped onto this unit circle. The eigenvalue problem (2.7) transforms with (2.8) to (2.9)

(A + (a − b)B) x = −(A − (a + b)B)ν x

Let C = A + (a − b)B and D = −A + (a + b)B. Although B is singular, the matrices C and D are generically not singular and we therefore consider the problem (2.10)

D−1 Cx = νx

Through this transformation, most dangerous modes related to σ are transformed to dominant modes associated with eigenvalues ν. To determine a prescribed number of these dominant modes, the Simultaneous Iteration Technique [28] is used. This technique, which is a generalized power method, may converge slowly for some cases but turns out to be fairly robust for a wide range of problems. 2.2.4. Solution of the linear systems. For both the computation of the steady states as well as the determination of its linear stability, linear systems must be solved. The Bi-CGSTAB method [33] is used as an iterative method for the solution of these linear systems. This is a conjugate gradient type method which can be used for systems of linear equations in which the coefficient matrix is non-symmetric. Its convergence behavior is strongly influenced by the location of the eigenvalues of the coefficient matrix, and it appears to be very important that the spectral condition number of this matrix is small [2, 32]. For this reason, the linear system is preconditioned first before Bi-CGSTAB is applied. There is a wide choice of preconditioners, see for example [25]. An incomplete LU-decomposition is used as a preconditioner in which the sparsity pattern of L + U is based on a drop tolerance εp . The construction of the factors L, U and the choice of εp is described in detail in [30]. When several systems of linear equations have to be solved in which the coefficient matrices do not differ very much, as in continuation methods such as used here, it is possible to use the same preconditioner many times. In order to increase the efficiency of the incomplete decomposition, it is possible to perform a renumbering of the unknowns which is based on the same basic idea as in multi-grid methods. Many iterative methods can eliminate high-frequency errors very effectively, but they are inefficient at eliminating a long wavelength error. A couple of iteration steps result in an approximation with a smooth error. This error can therefore be well corrected on a coarser grid. Solving the equations on the coarse grid gives the two-grid method. Applying this idea recursively on successively coarser

8

M. JEROEN MOLEMAKER AND HENK A. DIJKSTRA

grids leads to the multi-grid method. The preconditioning technique uses a partition of the unknowns based on a similar sequence of grids. Renumbering the unknowns according to this partition enables the construction of an incomplete LU-decomposition which effectively eliminates both highand low-frequency errors [31]. 2.2.5. Summary of BAGELS. The total algorithm combining the continuation method with the eigenvalue solver and the iterative linear systems solver is summarized as follows. Suppose a point (x0 , λ0 ), the tangent (x˙ 0 , λ˙ 0 ) and n eigenvectors (and eigenvalues) are computed, then the algoritm consists of the following steps 1. Compute the Euler guess: x = x0 + ∆s x˙ 0 ; λ = λ0 + ∆s λ˙ 0 . 2. Compute the Jacobian of F and the preconditioning matrix. 3. Solve the system (2.5 − 2.6) using Newton-Raphson with an initial guess from step 1 in m iterations. This requires the solution of 2∗m systems of linear equations using a fixed preconditioning matrix. 4. Compute a desired number of test functions τpq . For each test function, one linear system must be solved. 5. Compute the matrices corresponding to the (Cayley-transformed) eigenvalue problem (2.9) and a preconditioning matrix. Start the Simultaneous Iteration Technique with n starting (eigen)vectors. With l iterations until convergence this requires the solution of n∗l linear systems (2.10) with fixed preconditioning matrix.

y

3. Results. In Fig. 2, the triangular grid is shown on which the calculations were done for a realistically shaped North-Atlantic basin. This grid consists of 2489 elements and amounts to a total of 11591 degrees of freedom.

x

Fig. 2. Grid used for computations on a realistically shaped North-Atlantic

basin.

Elimination of errors in this relatively complex code has been achieved by

MULTIPLE EQUILIBRIA OF OCEAN CIRCULATION

9

reproducing results from a similar study in a simple rectangular basin [27]. 3.1. Multiple equilibria. To represent a realistic wind forcing, the shape of (f x , f y ) is determined from the Hellerman and Rosenstein [18] data set. This data set consists of wind stress values on a 2◦ by 2◦ grid that is interpolated using a cubic spline fit to the finite element grid. At the southern part of the basin, the wind field is easterly (from east to west), reflecting the trade winds. At higher latitudes (larger y), westerly winds are particularly strong.

Fig. 3. Pattern of the Hellermann and Rosenstein [18] wind stress forcing field. Maximum amplitude (maximum length of vector) of the dimensional wind stress is scaled with 0.05 N m−2 .

Using this wind stress forcing pattern and the Ekman number E as control parameter, the bifurcation diagram is shown in Fig. 4a. On the vertical axis, the maximum dimensional northward volume transport φ (measured in Sverdrup: 1 Sv = 106 m3 s−1 ), over a section is shown. This quantity is calculated as   Z x ′ (3.1) φ = HLU maxy maxx vh dx 0

Two solution branches are found of which one exists over the whole E interval. For a solution on this branch, of which the h anomaly field is shown for E = 3.0 × 10−6 in Fig. 4b, the Gulf Stream passes along Cape Hatteras northward. This will be called the ’deflected’ Gulf Stream solution since it is deflected as compared to the real mean path of the Gulf stream. In this solution the Gulf Stream remains attached to the coast until is separates from the coast at a too northward position near the New England Seamount Chain (see Fig. 1). There is a very weak northern circulation region, a feature common in relatively low resolution ocean models of the North-Atlantic. At this value of E, its maximum transport φ is about 21.5 Sv. The other branch in Fig. 4a exists only for values of E smaller than 3.2 × 10−6 , which is the position of the limit point on this branch. The

10

M. JEROEN MOLEMAKER AND HENK A. DIJKSTRA

(a)

(b)

(c)

Fig. 4. Bifurcation diagram showing φ for varying Ekman number. Solid (dotted) curves indicate stable (unstable) steady states and the labels refer to locations where the solution is plotted in subsequent panels. H indicates the Hopf bifurcation point. (b) Layer thickness anomalies h of the ’deflected’ Gulf Stream solution at E = 3.0 × 10−6 . (c) Layer thickness anomalies h of the ’separated’ Gulf Stream solution at E = 3.0 × 10−6 . Contours of h roughly correspond to streamlines except at positions very close to the no-slip boundaries.

solution of the upper branch at E = 3.0 × 10−6 shows (Fig. 4c) a separated Gulf Stream which turns into the open ocean near Cape Hatteras as it does in reality. There is now a strong recirculation region north of the Gulf Stream, which is slightly too concentrated near the coast compared to reality. At this value of E, the maximum transport φ is about 26.2 Sv, which is substantially larger than that of the ’deflected’ Gulf Stream. The latitude and longitude at which the maximum in φ occurs for both solutions is very similar and hence it is the separation behavior which determines this different transport.

MULTIPLE EQUILIBRIA OF OCEAN CIRCULATION

11

3.2. Effect of the wind stress shape. One of the older theories on the Gulf Stream proposes that separation is mainly determined by the shape of the wind stress. Here, we investigate the change in solution structure when the realistic wind stress is deformed to a simple sinusoidal wind stress shape which is an idealization of the large scale features of the real wind stress shape, i.e. westerlies in the central part of the basin with easterlies to the north and south. This change is accomplished through a homotopy parameter ph , according to f x = ph f x + (1 − ph ) cos(πy) (3.2)

f y = ph f y

A value ph = 0.0 corresponds to a cosine-shaped wind stress whereas ph = 1.0 corresponds to the realistic wind forcing. The bifurcation diagram in the parameter ph at E = 7.6 × 10−6 is shown in Fig. 5a. The transport increases substantially with decreasing ph , although the maximum amplitude of the wind stress forcing remains the same. For the simple wind forcing (ph = 0), the northward transport has increased to about 80 Sv. The most interesting feature of Fig. 5a is that over an interval in ph , multiple steady states also appear, i.e. on the single branch two limit points appear. Two solutions for the same wind forcing (ph = 0.5) are shown in Figs. 5b and 5c. For the realistic forcing, only the deflected solution exists (see Fig. 5a). When ph is decreased, the deflection becomes less pronounced and the transport increases to 32.2 Sv at ph = 0.5 (Fig. 5b). Below ph = 0.3, the deflected solution does not exist anymore. The solution connects to the branch of the separated solutions via an unstable branch. The solution for ph = 0.5 (Fig. 5c) has a maximum northward transport of 41.5 Sv and is the unique solution which exists for the cosine wind forcing at this value of E (Fig. 5d). Hence, the shape of the wind stress is important for obtaining a unique separated solution. In this model, the separated solution does not exist for realistic wind stress forcing when the friction is too large (large E). 3.3. Oscillatory instabilities. The first Hopf bifurcation is found on the deflected Gulf Stream solution branch at E = 7.6 × 10−6 (Fig. 4a). At the Hopf bifurcation, a complex conjugate pair of eigenvalues of the problem (2.7), say σ = σr + iσi , crosses the imaginary axis. Here σr is the growth factor (which is zero exactly at the bifurcation) of the perturbation and σi is its angular frequency. The dimensional period of the oscillation T follows directly from the angular frequency, i.e. T = 2π/σi . The oscillation can be represented by (3.3)

Φ(t) = exp (σr t) [xr cos σi t − xi sin σi t]

where xr , xi are the real and imaginary parts of the eigenvector x.

12

M. JEROEN MOLEMAKER AND HENK A. DIJKSTRA

(a)

(b)

(c)

(d)

Fig. 5. Bifurcation diagram in the parameter ph . Labels refer to locations where layer thickness fields are shown in subsequent panels. (b) Deflected state at ph = 0.5 (φ = 32.2 Sv). (c) Separated state at ph = 0.5 (φ = 41.5 Sv). (d) Steady state at ph = 0.0 (φ = 80.3 Sv).

At E = 3.0 × 10−6 , the deflected solution in Fig. 4b is unstable to one oscillatory mode. This mode is represented in Figs. 6a-b by layer thickness perturbations corresponding to the real and imaginary part of the eigenvector. In this way, two phases of the oscillation (3.3) are presented. The propagation of the pattern can be determined by first looking at xi = Φ(−T /4) (Fig. 6b) and than at xr = Φ(0) (Fig. 6a). The period of the oscillation is about 5 months. It is clear that the center of action for the oscillation is located in the Gulf Stream. The scale of the perturbations is about 500 km and the axis connecting the maxima and minima of the thickness perturbations makes about a 45◦ angle with the x-axis. This orientation does not appear directly related to the orientation of the jet itself. The disturbances propagate south-westward against the flow direction of the mean current of the steady state. During propagation, the axis of orientation of the perturbations does not change. The separated solution (point (c) in Fig. 4a) is also unstable to an

MULTIPLE EQUILIBRIA OF OCEAN CIRCULATION

(a)

(b)

(c)

(d)

13

Fig. 6. Eigenfunctions at marked locations in Fig. 4a for (E = 3.0 × 10−6 ); shown are the layer thickness deviations. (a-b) Real and imaginary part of the unstable oscillatory mode on the basic state shown in Fig. 4b. (σr = 0.71 × 10−1 , σi = 0.51). (c-d) As (a-b) but for the steady state shown in Fig. 4c (σr = 0.15 × 10−1 , σi = 0.61).

oscillatory mode having a period of about 6 months. Layer thickness perturbations with a scale of about 400 km move around the jet (Figs. 6c-d) with slow westward (eastward) propagating motion south (north) of the jet. Again the centers of action are located in the jet and the response outside of the jet (e.g. for x > −1) is weak. Contrary to the unstable mode for the deflected solution now the orientation of the axis connecting extrema of the perturbation changes during one cycle of the oscillation; the perturbation moves around the jet in a clockwise manner. 4. Summary and Discussion. The numerical technique which is described in this paper is useful for investigating the structure of steady solutions of large dynamical systems in parameter space. Such a system arises here through discretization of the shallow water equations on a midlatitude β-plane. The advantage of this type of analysis is that branches of unstable solutions can be found, which can be done only under special circumstances using time-integration techniques. Although these unstable solutions might have no direct physical relevance they are often necessary for constructing possible transitions between stable patterns and to demon-

14

M. JEROEN MOLEMAKER AND HENK A. DIJKSTRA

strate how the stable states are related. The main result of the application of this technique to the wind-driven ocean circulation in the North-Atlantic is the existence of multiple mean paths of the Gulf Stream. For one of these solutions, the deflected solution, separation is very diffuse and occurs at the wrong location. The other solution, the separated solution, shows a better mean path of the Gulf Stream. For the model used here, the correct separated solution exists, but small friction is required to be able to reach this solution. Hence, large scale ocean models may have trouble obtaining the correct separated mean path of the Gulf Stream simply because they are in the wrong parameter regime. The wind stress shape was shown to have a substantial impact on the results. The separated solution becomes a unique solution at relatively large friction for the simple cosine wind stress, and the deflected solution is unique for realistic wind forcing. Of course, within the model context here, mechanisms other than continental geometry and wind stress cannot be tested and therefore cannot be ruled out as being important. The existence of the multiple equilibria is expected to be robust. The origin of these equilibria is related to the imperfect pitchfork bifurcation (Fig. 4a). This imperfect bifurcation was already found in a very simple geometry [20, 27] and is therefore unrelated to the shape of the basin. It is the internal ocean dynamics represented by the nonlinear shallow water system which is behind the imperfect pitchfork [6, 13, 14]. The solution structure carries over, qualitatively unmodified, to the realistic geometry. At small enough friction, the steady states become unstable to oscillatory modes which consists of propagating disturbances along the Gulf Stream. Their time scale (on the order of a few months) and pattern indicates that these are modifications of the modes found in a rectangular domain [27]. The modes are prototypes to explain the variability in the meandering intensity of the Gulf Stream. Using infrared images for the period April 1982 through December 1989, for example Lee and Cornillon [22] found a 9-month dominant periodicity. It is likely that internal ocean dynamics is the source of this variability. Acknowledgments. All computations were performed on the CRAY C90 at the Academic Computing Centre (SARA), Amsterdam, the Netherlands within the project SC498. Use of these computing facilities was sponsored by the National Computing Facilities Foundation (NCF) with financial support from the Netherlands Organization for Scientific Research (NWO). The work was supported by an NWO PIONIER grant (030-76187) to H.D. REFERENCES

MULTIPLE EQUILIBRIA OF OCEAN CIRCULATION

15

[1] S. J. Auer, Five-year climatological survey of the Gulf Stream system and its associated rings, J. Geophys. Res., 92 (1987), pp. 11709 – 11726. [2] D. Axelsson, Solution of linear systems of equations: Iterative methods, in Sparse Matrix Techniques: Copenhagen, V. A. Barker, ed., Springer-Verlag, Berlin, 1977, pp. 1–51. [3] A. Babuska and A. K. Aziz, Lectures on the Mathematical Foundations of the Finite Element Method, Academic Press, New York, 1972. [4] W. S. Broecker, The great ocean conveyor, Oceanography, 4 (1993), pp. 79–89. [5] F. Bryan, High-latitude salinity effects and interhemispheric thermohaline circulations, Nature, 323 (1986), pp. 301–304. [6] P. Cessi and G. R. Ierley, Nonlinear disturbances of western boundary currents, J. Phys. Ocean., 23 (1993), pp. 1727–1735. [7] K. N. Christodoulou and L. E. Scriven, Finding leading modes of a viscous free surface flow: An asymmetric generalized eigenproblem, J. Sci. Comput., 3 (1988), pp. 355–406. [8] T. J. Chung, Finite Element Analysis in Fluid Dynamics, McGraw-Hill, 1978. [9] C. Cuvelier, A. Segal, and A. van Steenhoven, Finite Element Methods and Navier-Stokes Equations, D. Reidel Publishing Company, 1986. [10] J. Dengg, The problem of Gulf Stream separation: A barotropic approach, J. Phys. Ocean., 23 (1993), pp. 2182–2200. [11] J. Dengg, Beckmann, and R. Gerdes, The Gulf Stream separation problem, in The warmwatersphere of the North Atlantic Ocean, W. A. Kraus, ed., Borntraeger, 1996, pp. 253–290. [12] H. A. Dijkstra, On the structure of cellular solutions in Rayleigh-B´ enardMarangoni flows in small-aspect-ratio containers, J. Fluid Mech., 243 (1992), pp. 73–102. [13] H. A. Dijkstra and C. Katsman, Temporal variability of the wind-driven quasigeostrofic double gyre ocean circulation: Basic bifurcation diagrams, Geophys. Astrophys. Fluid Dyn., 85 (1997), pp. 195–232. [14] H. A. Dijkstra and M. J. Molemaker, Imperfections of the North-Atlantic winddriven ocean circulation:continental geometry and windstress shape, J. Marine Research, Submitted (1998). [15] H. A. Dijkstra, M. J. Molemaker, A. J. Van der Ploeg, and E. F. F. Botta, An efficient code to compute nonparallel flows and their linear stability, Comp. Fluids, 24 (1995), pp. 415–434. [16] I. Goldhirsch, S. A. Orszag, and B. K. Maulik, An efficient method for computing leading eigenvalues and eigenvectors of large asymmetric matrices, J. Sci. Comput., 2 (1987), pp. 33–58. [17] G. Golub and C. F. Van Loan, Matrix Computations., The Johns Hopkins University Press, 1983. [18] S. Hellerman and M. Rosenstein, Normal monthly wind stress over the world ocean with error estimates, J. Phys. Ocean., 13 (1983), pp. 1093–1104. [19] N. Hogg, R. Pickart, R. Hendry, and W. Smethie, The northern recirculation gyre of the Gulf Stream, Deep-Sea Res., 33 (1986), pp. 1139–1165. [20] S. Jiang, F. Jin, and M. Ghil, Multiple equilibria and aperiodic solutions in a wind-driven double gyre, shallow water model, J. Phys. Ocean., 25 (1995), pp. 764–786. [21] H. B. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems, in Applications of bifurcation theory, P. H. Rabinowitz, ed., Academic Press, 1977. [22] D. Lee and P. Cornillon, Temporal variation of meandering intensity and domain-wide lateral oscillations of the Gulf Stream, J. Geophys. Res., 100 C7 (1995), pp. 13,603–13,613. [23] J. Pedlosky, Geophysical Fluid Dynamics, Springer-Verlag, New York, 1987. [24] P. Richardson, Benjamin Franklin and Timothy Folger’s first printed chart of the Gulf Stream, Science, 207 (1980), pp. 643–645.

16

M. JEROEN MOLEMAKER AND HENK A. DIJKSTRA

[25] Y. Saad, Krylov subspace methods on supercomputers, SIAM J. Sci. Stat. Comput., 10(6) (1989), pp. 200–232. [26] R. Seydel, Numerical computation of branch points in nonlinear equations, Num. Math., 33 (1979), pp. 339–352. [27] S. Speich, H. Dijkstra, and M. Ghil, Successive bifurcations in a shallow water model, applied to the wind-driven ocean circulation, Nonl. Proc. Geophys., 37 (1995), pp. 289–306. [28] W. J. Steward and A. Jennings, A simultaneous iteration algorithm for real matrices, ACM Trans. Math. Software, 7 (1981), pp. 184–198. [29] C. Taylor and P. Hood, A numerical solution of the Navier-Stokes equations using the finite element technique, Int. J. Computers and Fluids, 1 (1973), pp. 73–100. [30] A. v.d. Ploeg, Preconditioning techniques for non-symmetric matrices with application to temperature calculations of cooled concrete, Int. J. Num. Methods Eng., 35(6) (1992), pp. 1311–1328. [31] A. v.d. Ploeg, E. Botta, and F. Wubs, Grid-independent convergence based on preconditioning techniques, report W-9310, Department of Mathematics, Groningen, (1993). [32] A. v.d. Sluis and H. A. v.d. Vorst, The rate of convergence of conjugate gradients, Numer. Math., 48 (1986), pp. 543–560. [33] H. A. v.d. Vorst, Bi-CGSTAB: A fast and smoothly converging variant of bicg for the solution of nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 13 (2) (1992), pp. 631–644.