The governing equations

(0)2. 0 We assume an infinitely deep, or equivalently a dynamically passive, hydrostatic lower layer. This assumption leads to the relations,

where B is the undisturbed depth of the upper layer, b is the depth of the thermocline, and h is the height of the free surface, using the undisturbed free surface as the reference level, h=0. is the ratio of the density of the upper layer to that of the lower layer; we assume it to be 0.997, implying a density difference of 0.3% across the thermocline. It is clear from the above relations that most of the variation in the surface layer depth is accounted for by the movement of the themocline b since it is more heavily weighted than h.

The equations governing the evolution of the double gyre system in the classic ( layer) SW approximation are

where is the thickness of the surface layer, is the vertically averaged horizontal velocity vector, and is reduced gravity. The term is the usual Coriolis acceleration due to the vertical component of the Earth's rotation (f) and is a zonally uniform steady wind stress, described in more detail later. Here is supposed to represent the enhanced diffusivity due to eddy processes and is Rayleigh friction which models momentum losses at the interface. These two terms are parametrizations of sub-grid scale processes.

At the Green-Naghdi level (system 2), the following singular perturbation terms,

are added on the right-hand-side of the momentum equation, with b related to the layer depth according to (). Finally, other additional terms appear on the right-hand-side of the momentum equation in system 3 that account for the vertically-integrated effects of the leading order horizontal rotation terms. These terms are

where is the (horizontal) component of earth's rotation in the direction ,and u is the horizontal velocity component in the direction. Details of the derivations of the additional terms shown in equations () and () are given in Appendix A.

The numerical model we use is based on a forward-in-time (2 time level) semi-Lagrangian discretization of the above equations that is second-order accurate in both space and time. Gravity waves are treated implicitly and the resulting elliptic equations are formulated in terms of the horizontal velocity . The method of conjugate residuals is used to solve these elliptic equations [14].We impose no-slip boundary conditions on all four lateral boundaries of the basin. Details of the numerical scheme are given in Appendix B. All the computations were carried out on the Connection Machine CM-5 using either 32 or 64 nodes.

The timestep is limited to ensure that there is no crossing of the Lagrangian trajectories. This is the only limitation required to guarantee the stability of the computations. Thus, we are free to choose any value of and unlike some other numerical formulations. Mass is conserved in the computations during the course of 110 years to better than 0.3% and so the lack of conservation in the semi-Lagrangian formulation is not important.