Bifurcation structure of a wind-driven shallow water model with layer-outcropping

The steady state bifurcation structure of the double-gyre wind-driven ocean circulation is examined in a shallow water model where the upper layer is allowed to outcrop at the sea surface. In addition to the classical jet-up and jet-down multiple equilibria, we ﬁnd a new regime in which one of the equilibrium solutions has a large outcropping region in the subpolar gyre. Time dependent simulations show that the outcropping solution equilibrates to a stable periodic orbit with a period of 8 months. Co-existing with the periodic solution is a stable steady state solution without outcropping. A numerical scheme that has the unique advantage of being diﬀerentiable while still allowing layers to outcrop at the sea surface is used for the analysis. In contrast, standard schemes for solving layered models with outcropping are non-diﬀer-entiable and have an ill-deﬁned Jacobian making them unsuitable for solution using Newton’s method. As such, our new scheme expands the applicability of numerical bifurcation techniques to an important class of ocean models whose bifurcation structure had hitherto remained unexplored.


Introduction
The wind-driven double-gyre ocean circulation is responsible for the transport of heat, mass and fresh water between low and high latitudes, and is therefore important to climate variability. An important feature of the mid-latitude ocean circulation -found both in observations and numerical simulations -is its lowfrequency variability. For example, in observations decadal time scale variability has been documented by Qiu and Chen (2005) for the Kuroshio extension system and in models intrinsic decadal variability due to non-linear dynamics is common (e.g. McCalpin and Haidvogel, 1996;Hogg et al., 2005). The mechanisms responsible for this variability however are not fully understood. The motivation for this paper is to continue a line of study that applies bifurcation theory and dynamical-systems theory to a hierarchy of increasingly more realistic ocean models with the goal of better understanding the dynamics of this flow.
Absent from the above hierarchy are models based on a layered formulation that allow isopycnal layers to outcrop at the sea surface or along sloping bottom topography. This gap might appear surprising given the historical and continuing importance of models with outcropping layers in the development of theories of the ocean circulation. An analytical solution to a layered formulation with outcropping was the key to the development of the theory of the ventilated thermocline (Luyten et al., 1983;Pedlosky, 1996), models with outcropping are also central to many theories explaining the separation of the Gulf Stream from the coast (Parsons, 1969;Veronis, 1973;Chassignet and Bleck, 1993), and layered model formulations remain an important class of general circulation models (Chassignet et al., 2000). However, the gap in the hierarchy is easy to explain if we consider the computational difficulties involved. Numerical bifurcation analysis relies on Newton's method for computing both stable and unstable steady state solutions, but Newton's method is not applicable to outcropping numerical models because the numerical schemes that allow outcropping are generally non-differentiable and do not have a properly defined Jacobian. Without a properly defined Jacobian, the convergence of Newton's method is unreliable making numerical bifurcation analysis impossible.
In the present article, we show how to overcome the numerical difficulties by modifying the numerical scheme so that it does have a properly defined Jacobian. We focus on a 1.5-layer wind-driven shallow water (SW) model because it is the simplest possible model with layer outcropping. We concentrate on finding steady-state solution in the presence of outcropping to demonstrate the feasibility of the scheme for bifurcation and stability analyses.

Model formulation and solution
Consider an ocean in a rectangular basin forced by wind stress. We represent this with the reduced-gravity wind-forced shallow water model where (u, v) is the velocity of a layer of fluid of depth h and density q 0 , overlying an infinitely deep layer assumed at rest and of slightly higher density such that the dynamics of the upper layer feel a reduced gravity g 0 . The Coriolis parameter is assumed to vary linearly with latitude, f = f 0 + by. The wind-stress forcing is (s x , s y ) and the friction is parameterized as an interfacial Rayleigh drag with parameter r plus a horizontal biharmonic diffusion of momentum with coefficient A. We will also refer to (1) as the 1.5-layer shallow water model. Computing the time-evolution of these shallow water equations, such as (1), is relatively straightforward. A variety of numerical schemes are available, such as Arakawa and Lamb's C-grid scheme (Arakawa and Lamb, 1977). But even with non-dimensionalization, there are many parameters that affect the dynamics of this winddriven flow. To get a more general and broader understanding of the overall dynamics of this model, we compute bifurcation diagrams and parameter charts. These diagrams and charts show, as a function of parameter value, where there are qualitative changes in the flow. The simplest class of solutions to (1) are equilibrium solutions, that is (u, v, h) that satisfy the equations with the time derivatives set to zero. Our goal is to trace out the locus of these equilibrium solutions, otherwise known as the branching structure, as a function of one or more parameter values.
To compute these equilibrium solutions, we discretize the steady state form of the above shallow water model using the standard Arakawa C-grid, producing a set of non-linear algebraic equations of the form where u = (u, v, h) is a n-dimensional vector of unknown velocities and depths at the grid points, and p is an m-dimensional vector of parameters. The parameter vector p could include parameters that tend to vary or that are not precisely known, for example the Rayleigh drag parameter r, the horizontal diffusion coefficient A, and the wind-stress intensity s. We compute branches of solutions by varying any one of the m parameters in p. Using pseudo-arclength continuation (Keller, 1977), we map out the branch structure and compute equilibrium solutions that are either stable or unstable. Primeau and Newman (2006), discuss how the ability to vary two parameters is key to mapping out the different parts of the solution surface that would not be accessible by varying only one parameter. In this study we also use this approach of varying two parameters. The pseudo-arclength continuation method incorporates Newton's method to solve the system of equation (2), for a fixed parameter vector p. While Newton's method is a well-known technique (e.g. Kelly, 2003), the ability to compute the Jacobian J ¼ oF ou is critical.

Keeping layer thickness non-negative
Numerical solutions to the multi-layer shallow water equations must ensure that the thickness of each layer never becomes negative. Preventing negative layer thickness (like preventing negative tracer concentrations), is a known problem with many solution approaches. Multiple-layer ocean (atmosphere) models often use an isopycnal (isentropic) vertical coordinate system, for example the MICOM (Bleck et al., 1992) or HIM (Hallberg and Rhines, 1996) ocean models. Layer outcropping can then be directly modeled by having coordinate surfaces appear to intersect with the boundary, which is numerically represented by massless layers extending to the edge of the domain, as shown in Fig. 1.
layer-2 layer-1 Fig. 1. Schematic of a multiple-layer isopycnal ocean model with layer-2 outcropping at the surface. Numerically, layer-1 extends to the left hand edge of the domain as a layer of negligible thickness. In this study we consider a 1.5-layer model where layer-2 is infinitely deep. The sharp front created by the layer-1/layer-2 interface rapidly meeting the surface needs to be accurately and stably handled by the numerical method.
There are a variety of numerical schemes for keeping layer thickness non-negative, most of which rely on an upwind scheme. In an upwind scheme the discretization of the flux terms depends on the sign of the flux. For example, a one-dimensional thickness advection equation discretized in flux-divergence form using a firstorder upwind scheme, approximates thickness flux term as follows: Such a scheme is a simple way to prevent negative values but has the disadvantage of being highly diffusive. A less diffusive scheme that is also an upwind scheme is the flux-corrected transport method (Boris and Book, 1973). Hsu and Arakawa (1990) developed a positive-definite advection scheme to discretize the continuity and momentum equations for a layered-atmosphere model. Even though, by design, this scheme is meant to maintain positive thickness while minimizing dispersion errors, the authors explain that due to truncation or roundoff errors a further modification is required to ensure that layer thicknesses do not become less than a prescribed lower limit. This modification redistributes mass vertically to restore layer thickness to at least this lower limit. However, all these differencing schemes are non-differentiable, primarily because they all include an upwind or upwind-like term. By non-differentiable, we mean that the Jacobian J is not well-defined. This can be quickly seen by considering a function f(u) that includes absolute values juj, which is the case in upwinding. In this case, the Jacobian of ou is not defined at u = 0. This lack of differentiability is not an issue for time integrations, but is crippling for solving the discretized steady state equation (2). Convergence of Newton's method is at best not helped by non-smoothness of F(u, p), and at worst not guaranteed. Consequently we need to find an alternative numerical scheme for solving the shallow water equation with outcropping. Salmon (2002) proposed a new numerical scheme for solving the two-layer shallow water equations that can -with appropriate modifications -avoid non-differentiable terms. To prevent layer thicknesses from becoming negative, he modified the dynamics by inserting an artificial form of potential energy that becomes very large as the layer becomes very thin. Salmon demonstrated the accuracy of the results, showing that it was a useful scheme for studying wind-driven ocean circulation with bottom topography.
Salmon proposed modifying the pressure gradient term in the momentum equations from Àg 0 $h to Àg 0 (1 + P(h))$h. The P(h) term represents a modification to the potential energy with the properties that P(h) is negligible when h ) h 0 , where h 0 is a small positive constant, and P(h) ? 1 as h ? 0. Salmon suggested the function P(h) = (h 0 /h) n and numerically determined that n = 4 worked well without incurring a heavy stiffness penalty due to the smaller time step required by higher wave speed. We will refer to h 0 as the Salmon thickness. The key feature of this modification is that it is differentiable, and can therefore be used to compute equilibrium solutions using Newton's method. In this study, we adopted this method for preventing negative layer thickness. We denote this transformation as the following modification to the reduced gravity: Fig. 2 shows the modified potential energy versus the layer thickness (dashed line) compared to the unmodified potential energy (heavy solid line). We point out that this is an artificial modification that may change the dynamics of the shallow water model, so our goal is to be able to select h 0 sufficiently small so as to minimize the non-physical nature of this term, and show that solutions are relatively insensitive to changes in h 0 .

Vertical distribution of wind-stress forcing
The forcing in the momentum equations in our 1.5-layer model (1) comes from the s x /q 0 h term (for zero north-south wind stress). The magnitude of this forcing term becomes infinite as the layer outcrops and h ? 0. Although the wind stress only acts at the interface, the 1.5-layer model does not realistically represent what happens physically. A more likely scenario is that there is significant vertical mixing occurring at the surface. A simple way to represent this mixing is to vertically distribute the wind stress input. We introduce a modification to this forcing term that will keep the momentum flux finite as the layer thickness approaches zero. If we multiply the forcing term by a factor that approaches unity for large h and zero as h ? 0 we obtain the following modified forcing: where h * represents a length scale over which the wind-stress is distributed. In this paper we will refer to h * as the mixing thickness. The above modification to the wind-stress forcing could be considered to be a parameterization of vertical mixing with length scale h * . Using the above expression, when h = 0 the forcing in the momentum equation approaches F x ¼ s x q 0 h Ã . It is more realistic to distribute wind stress through the layers, particularly if the top layer becomes very thin. Multi-layer models such as HIM and MICOM in fact distribute wind stress over one or more layers, depending on the thickness of the top layer.

Final model configuration
Incorporating the modifications in Eqs. (3) and (4) into our shallow water equations (1), and setting the time derivatives to zero gives where h 0 is the Salmon thickness, h * is the mixing thickness, s is the magnitude of the wind-stress forcing, and /(x) the wind-stress pattern characterized by a mid-latitude sinusoid where L y = 2000 km is the north-south extent of the basin. For this study, we consider simulations on the b plane in a rectangular domain of (1000 Â 2000) km 2 . The nominal upper layer thickness is h = 100 m (so the volume of fluid is constant 2 Â 10 14 m 3 ), with a density of q 0 = 1000 kg m À3 . The reduced gravity g 0 = 0.1 m s À2 , and the Coriolis parameter and the planetary PV gradient are set to f 0 = 1 Â 10 À4 s À1 and b = 2 Â 10 À11 m À1 s À1 . The Rayleigh drag is r = 5.3 Â 10 À7 s À1 and the horizontal diffusion is A = 5.2 Â 10 10 m 4 s À1 . At the boundaries we apply free-slip conditions. We consider two spatial resolutions, the high resolution case at Dx = Dy = 7.8 km and the low resolution case Dx = Dy = 15.6 km. This resolution is sufficiently high to resolve diffusive length scales controlled by our choice of hyper-diffusion parameter. With these parameters set, we compute equilibrium solutions for Salmon thickness h 0 = 16 m and h 0 = 32 m on the two grids. Our two continuation parameters are the wind-stress and the mixing thickness p = (s, h * ).
Details of the arclength continuation methodology can be found in Primeau and Newman (2006). To invert the sparse Jacobian matrix we use a multi-frontal approach as implemented in the MUMPS package (Amestoy et al., 2001).

Results
We have introduced two modifications to our original shallow water equations. The first modification (4), which is an artificial potential to prevent layer thickness from becoming negative, is not physical, so we need to show that solutions are independent of h 0 , for relatively small h 0 . The second modification, described by (5), is a reasonable representation of the physical process of vertical mixing, so we treat h * as a physical parameter. We compute parameter charts as a function of s and h * for different values of h 0 . We will also show that the results are relatively insensitive to h 0 and resolution.

Bifurcation structure and parameter chart
Using the arclength continuation approach in a two parameter setting as in Primeau and Newman (2006), we have traced out the locus of saddle node bifurcation (turning points) where the solution surface folds to produce multiple equilibria. The two parameters are the wind-stress s and the mixing thickness h * . The parameter chart in Fig. 3 shows the result. We computed the parameter chart for two values of the Salmon thickness, the solid lines with h 0 = 16 m and the dashed lines with h 0 = 32 m. These show that the boundary is relatively insensitive to the choice of the Salmon thickness parameter h 0 . This insensitivity is important, because it confirms that, given sufficiently high resolution, the modification of the potential energy has little effect, and can be tolerated.
On the parameter chart we see two boundaries, the upper one displaying a classic cusp type bifurcation, (e.g. Arnold, 1992)  The cusp point at s = 0.044 N/m 2 and h * = 62 m indicated by the filled circle in Fig. 3 shows where the solution surface first folds to produce multiple equilibria. Cuts through parameter space with fixed h * that intersect the wedge shaped region would yield an imperfect pitchfork bifurcation structure like the one first found by Jiang et al. (1995). The orientation of this cusp in parameter space, and the characteristics of the multiple equilibria in the wedge shaped region that emanates from the cusp are similar to those found by Primeau and Newman (2006), for a shallow water model without layer outcropping. The two-parameters used by Primeau and Newman (2006), included the horizontal Ekman number and a parameter that controls the relative input of vorticity into the subpolar and subtropical gyres. The correspondence between those parameters and the ones used here can be explained by the fact that increasing the wind-stress has a similar effect to decreasing the Ekman number -both increase non-linearity -and that varying h * has an effect similar to varying the wind-stress asymmetry parameter used in Primeau and Newman (2006). The reason why varying h * has the effect of varying the wind-stress asymmetry is due to the fact that when the upper layer becomes sufficiently thin the mixing thickness h * limits the fraction of momentum input into the upper layer by spreading the momentum input vertically into the deeper layer. Increasing h * thus decreases the relative input of vorticity in the subpolar gyre where the active layer is thin in comparison to the subtropical gyre. The similarity between our present results and those of Primeau and Newman (2006) is also due to the fact that for equilibria a, b and c, the top layer does not outcrop and without outcropping both models are essentially the same. Unlike the situation near the cusp point, the broad rounded locus of turning points in Fig. 3 is a feature that was not found in previous studies. Fig. 7 shows the solutions for s = 0.07 N/m 2 and h * = 6 m. These are indicated by d and e in Figs. 3 and 5. The right panel (e) of Fig. 7 shows a solution where the subpolar recirculation gyre has outcropped. In this solution the inter-gyre jet turns northward following the outcrop boundary. The left panel shows a solution (d) where the subtropical recirculation gyre is stronger than the subpolar one, and in which the inter-gyre jet turns southward. We have thus found a case where we have multiple equilibria, one of which has outcropped. Apart from the non-linearity due to the inertial terms it appears that the non-linearity associated with layer-thickness change plays an important role in the existence of these new multiple equilibria.
To see more clearly the nature of the outcropped front, we take a section at x = 250 km, and show the depth profile in Fig. 8. We can see that the outcropping produces a sharp front with a relatively small depth of h = 8 m in the outcrop region. Note that this minimum depth is related to the Salmon thickness h 0 which is set to h 0 = 32 m for this case. Recall that the artificial potential scales as (h 0 /h) 4 so that h must become somewhat smaller than h 0 before the artificial pressure force is strong enough to prevent further thinning of the layer.

Sensitivity results
To show that the results of the bifurcation analysis are not numerical artifacts we recomputed all the steady state solutions for both h 0 = 16 m and h 0 = 32 m at a different grid resolution. Fig. 9 shows the parameter chart for the coarser resolution of Dx = 15.6 km. The same qualitative features seen in Fig. 3 are recovered. Importantly, we note that for the higher resolution case (Fig. 3) the solution is less sensitive to the precise value of the Salmon thickness than for the low resolution case (Fig. 9). This confirms that the Salmon potential is a consistent scheme for doing numerical bifurcation analysis in models with layer outcropping.

Stability analysis
In this section we demonstrate that because our model has a properly defined Jacobian it is amenable to a modal stability analysis. To determine the stability we linearize the dynamical equations about the steady state solution and look for solutions of the form u exp(Àrt). This yields an eigenvalue problem for the Jacobian matrix J, which we solve using ARPACK (Lehoucq et al., 1998). Of the three solutions shown in Fig. 6 we find that the jet-up and the jet-down solutions (a and c) are stable while the intermediate solution (c) is unstable. The stability characteristics of these solutions near the cusp point is qualitatively similar to what was found in previous studies with either quasi-geostrophic or shallow water models (e.g. Cessi and Ierley, 1995;Jiang et al., 1995).
The modal stability analysis reveals that solution (d) is stable to small perturbations and that solution (e) is unstable with an oscillatory mode with a period of approximately 7.7 months. The amplitude and phase of the  Fig. 3, the chart shows the locus of turning points for Salmon thickness h 0 = 16 m (solid curves) and h 0 = 32 m (dashed curves). We see broad qualitative agreement between the features and structure in this parameter chart with those in Fig. 3. Note that the Hopf boundary is not plotted on this chart. unstable mode for solution (e) is shown in Fig. 11. The mode amplitude is focused around the front where the layer outcrops at the surface. Lines of constant phase indicate the propagation of wave disturbances around the outcrop in the counter-clockwise direction. We also verified the results of our linear stability analysis with fully non-linear time-dependent runs. Simulations started from equilibrium (d) perturbed with small amplitude noise confirmed that the solution is indeed stable. Simulations started from equilibrium (e) on the other hand moved away from the equilibrium point until the model trajectory converged to a finite amplitude limit cycle with a period of 8 months. The time history of the solution's basin integrated kinetic energy is shown in Fig. 10. A sequence of contour plots of the layer thickness at four phases during one period of the time-dependent solution is shown in Fig. 12. The figure shows how the outcrop region wobbles as the northward-curving jet contracts and sheds a 'warm' eddy into the subpolar gyre before extending eastward to repeat the cycle. The variability is consistent with the unstable mode of the steady state solution (e). The unstable mode's period of 7.7 is close to the 8 month period of the finite amplitude limit cycle. The amplitude and phase propagation of the mode along the outcrop front is also consistent. We also computed the Hopf bifurcation boundary separating the stable from the unstable part of the solution branch on which solution (e) lies. The boundary is shown in Fig. 13. To the left of the boundary the solution is stable. To the right of the boundary the steady state solution is unstable. Close to the marginal stability boundary the flow equilibrates to a finite amplitude limit cycle. At the point (e) in Fig. 13 the limit cycle solution co-exists with the stable steady state solution (d) on another branch. The initial condition determines which regime the time-dependent trajectory will ultimately tend to.

Conclusion
By applying a numerical scheme first proposed by Salmon (2002) to solve the shallow water equations, we have performed the first numerical bifurcation analysis of the double-gyre circulation in the case where layers outcrop at the surface. In one part of parameter space we recovered the multiple equilibria with the jet-up and jet-down solutions first found by Jiang et al. (1995). In other parts of parameter space we have found new situations in which the same forcing parameters can produce multiple equilibria, one of which has a large outcrop region in the subpolar gyre.
We have tested that the results of the steady-state bifurcation are not sensitive to an artificial potential energy term added to prevent layer thicknesses from becoming negative. The new scheme makes it possible to apply the techniques of numerical bifurcation analysis and dynamical systems theory to an important class of ocean models where layers outcrop either at the sea surface or along sloping boundaries.
Although a model with multiple active layers would be a better representation of the thermocline we have made the important step of addressing the outcropping problem in a 1.5-layer model. This step will allow us, in further studies, to continue this research with a multi-layer shallow water model. This is particularly important because outcropping of the thermocline drives seasonal changes in the ocean-atmosphere instability, which may contribute to the maintenance of persistent winter-time sea surface temperature anomalies that can produce mid-latitude climate variability.