Nonlinear Processes in Geophysics Conditional nonlinear optimal perturbations of the double-gyre ocean circulation

In this paper, we study the development of finite amplitude perturbations on linearly stable steady barotropic double-gyre flows in a rectangular basin using the concept of Conditional Nonlinear Optimal Perturbation (CNOP). The CNOPs depend on a time scale of evolution te and an initial perturbation threshold δ. Under symmetric wind forcing, a perfect pitchfork perturbation occurs in the model. The CNOPs are determined for all linearly stable states and the time evolution of the CNOPs is studied. It is found that the patterns of the CNOPs are similar to those of the non-normal modes for small te and approach those of the normal modes for larger te. With slightly asymmetric winds, an imperfect pitchfork occurs in the model. Indications are found that the time evolution of the CNOPs is related to the value of the dissipation function of the underlying steady state.


Introduction
The so-called quasi-geostrophic double-gyre flow has been recognized as one of the characteristic problems to study the nonlinear dynamics of the wind-driven ocean circulation (Jiang et al., 1995;Dijkstra, 2005).Usually such a flow is considered in an idealized geometry, such as a rectangular ocean basin, on a midlatitude β-plane.The linear problem, neglecting inertia, is the basis for the Sverdrup-Stommel-Munk theory of the wind-driven ocean circulation.In this case, the Sverdrup balance holds over most of the basin and viscosity only affects the flow in thin boundary layers at the eastern and western boundaries.The Sverdrup transport is compensated only in the western boundary layer and hence the western boundary flow is much stronger than the eastern one.
Correspondence to: A. D. Terwisscha van Scheltinga (arjen.terwisschavanscheltinga@ualberta.ca)When the flow is forced by a symmetric wind stress (with respect to the mid axis of the basin) the quasi-geostrophic equations have a reflection symmetry.For the one-layer (barotropic) case much is known about stability bounds and bifurcation behavior of nonlinear flows as the viscosity is decreased.When only lateral friction is considered as a dissipation mechanism, there is basically only one control parameter, the Reynolds number Re.The finite amplitude stability of the linear (Munk) solution was studied using analytical methods (Crisciani and Mosetti, 1990;Crisciani et al., 1994Crisciani et al., , 1995)).The energy stability limit Re E of the unique antisymmetric solution, existing at high viscosity, was calculated numerically in Dijkstra and De Ruijter (1996) and guarantees a monotonic decay of the kinetic energy of any perturbation; this stability limit hence provides sufficient conditions for stability (Joseph, 1976).
The linear stability limits Re L of this barotropic doublegyre flow in a relatively small ocean basin were presented in Dijkstra and Katsman (1997).The first bifurcation is a symmetry-breaking pitchfork bifurcation where the antisymmetric solution becomes unstable and two asymmetric solutions appear.These asymmetric solutions become unstable at several Hopf bifurcations where periodic orbits appear.Eventually chaotic behavior occurs due to a homoclinic bifurcation, which can be either of Lorenz or Shilnikov type, depending on the parameters of the system (Nadiga and Luce, 2001;Simonnet et al., 2005).
In recent years, tools of generalized stability theory (Farrell and Ioannou, 1996;Moore, 1999;Moore et al., 2002) have also been applied to this problem.With these tools, one is interested in determining growth of perturbations on a particular reference state due the non-normality of the Jacobian of that state; the latter state may be either a linearly stable steady state or a time-mean state of a very irregular flow.These tools enable one to determine the response of the flow to stochastic perturbations and hence are interesting with respect to predictability issues.
A. D. Terwisscha van Scheltinga and H. A. Dijkstra: CNOPs of double-gyre flows In Moore (1999), a basin of size 1000×2000 km was considered for both the asymptotically stable as well as unstable regimes and the effects of stochastic wind forcing on the flow was studied.Focus was on the stochastic forcing patterns that account for the largest fraction of noise-induced variability, the so-called stochastic optimals.The structure of the Ekman pump velocity of the gravest stochastic optimal on time scales of 2-3 weeks corresponds to a single-gyre basin wide flow.This particular wind perturbation induces changes in vorticity which project strongly on the fastest linear singular vectors.The noise forcing is most effective in the asymptotically stable regime but otherwise not sensitivity to the chosen norm, basic state flow and geometry.In Moore et al. (2002), it was shown that the variability is maintained by Rossby waves that interact with the western boundary current.The perturbations that maintain the stochastically induced variance in the linearly stable regime have a large projection on some of the non-normal, least-damped eigenmodes.
In generalized stability theory, it is assumed that the initial perturbation is so small that its evolution can be described by a linearized system (the tangent linear model) and optimal growth is determined through the largest singular value of the forward propagator of the linearized system.A generalization of linear singular vectors is the concept of Conditional Nonlinear Optimal Perturbations (CNOPs) as introduced by Mu et al. (2003).The CNOP is the initial (finite amplitude) perturbation whose nonlinear evolution attains a maximum growth rate at a chosen end time t e given an initial bound δ on the norm of the initial condition.
The CNOPs of a steady state determines the dominant time-dependent nonlinear behavior of finite amplitude perturbations.On one hand, such behavior bridges the gap between the behavior below the energy stability boundary (monotonic decay of all perturbations) and above the linear stability boundary (exponential growth of infinitesimally small perturbations).On the other hand, when compared to non-normal modes, the CNOP displays how much nonlinearity affects the evolution of finite amplitude perturbations.In the case of linearly stable multiple equilibria, the CNOPs also provide a way to compute finite amplitude stability boundaries of each of the equilibria (Mu et al., 2004).It is thus important to be able to compute CNOPs for flows modeled by systems of partial differential equations.
The computation of CNOPs has so far been only accomplished in models having a small number of degrees of freedom such as ocean box models (Mu et al., 2004) and relatively simple atmospheric (Mu and Zhang, 2006) and ENSO models (Mu et al., 2003).As far as we know the CNOPs for the barotropic double-gyre ocean flow problem are here calculated for the first time.We use the implicit 4D-Var methodology (Terwisscha van Scheltinga and Dijkstra, 2005) which is relatively easily extended to compute CNOPs.We determine the CNOPs for the double-gyre problem both under symmetric and asymmetric wind stress forcing.

Model and methods
In this section, we will shortly recall the model used (Sect.2.1) and then briefly describe the CNOP methodology (Sect.2.2).

Model
Consider a flow domain V consisting of a rectangular ocean basin of size L×L having a constant depth D. The basin is situated on a midlatitude β-plane with a central latitude θ 0 =45 • N and Coriolis parameter f 0 =2 sin θ 0 , where is the rotation rate of the Earth.The meridional variation of the Coriolis parameter at the latitude θ 0 is indicated by β 0 .The density ρ of the water is constant and the flow is forced at the surface through a wind-stress vector T=τ 0 [τ x (x, y), τ y (x, y)].The governing equations are nondimensionalized using a horizontal length scale L, a vertical length scale D, a horizontal velocity scale U , the advective time scale L/U and a characteristic amplitude of the windstress vector, τ 0 .The effect of deformations of the oceanatmosphere interface on the flow is neglected.
The dimensionless barotropic quasi-geostrophic model of the flow for the vorticity ζ and the geostrophic streamfunction ψ is (Pedlosky, 1987) where the horizontal velocities are given by u=−∂ψ/∂y and v=∂ψ/∂x.The parameters in Eq. ( 1a) are the Reynolds number Re, the planetary vorticity gradient parameter β and the wind-stress forcing strength α τ .These parameters are defined as: where g is the gravitational acceleration and A H is the lateral friction coefficient.When the horizontal velocity scale is based on a Sverdrup balance of the flow, i.e., it follows that α τ =β.Consequently, there are only two free parameters, for example the dimensionless boundary layer thicknesses δ I and δ M defined by δ 2 I =1/β and δ 3 M =1/(βRe), respectively.
We assume no-slip conditions on the east-west boundaries and slip on the north-south boundaries.The boundary conditions are therefore given by The wind-stress profile considered is where the dimensionless parameter σ controls the shape of the zonal wind stress and τ 0 is a typical amplitude.For σ =1 (σ =0), the wind stress induces a single-gyre (double-gyre) flow.
The governing equations are discretized on a 60×40 grid with central spatial differences.The resolution in the eastwest direction is slightly higher because the flows are westerly intensified.An implicit time-integration scheme (Terwisscha van Scheltinga and Dijkstra, 2005) is used with a time step of t=1 day.Standard parameter values of the model are shown in Table 1.After discretization, the state vector x ∈ R d (of dimension d=2×60×40=4800) consists of the values of ψ and ζ at the grid points.

Calculation of CNOPs
The discretized equations governing the evolution of perturbations x on a particular state x can be written as: where t is time, x(t)=(x 1 (t), x 2 (t), ..., x n (t)) is the perturbation state vector and F is a nonlinear differentiable operator.Furthermore, x 0 is the initial perturbation, x is the basic state which we take here as a linearly stable steady state, (x, t)∈R d ×[0, t e ] and t e <+∞.Suppose the initial value problem (6) is well-posed and the nonlinear propagator M is defined as the evolution operator of Eq. ( 6) which determines a trajectory from the initial time t=0 to time t e .Hence, for fixed t e >0, the state is the result of the time-evolution at t=t e of the initial perturbation x 0 at t=0.For a chosen norm • measuring x, the perturbation x 0δ is called the Conditional Nonlinear Optimal Perturbation (CNOP) with constraint condition C(x 0 )= x 0 ≤δ, if and only if where The CNOP is the initial perturbation x 0 whose nonlinear evolution attains the maximal value of the functional J at time t e with the constraint condition x 0 ≤δ; in this sense it is  16), against the control parameter Re=U L/A H .The energy stability boundary Re E is about Re E ≈10 (Dijkstra and De Ruijter, 1996).(b) Streamfunction ψ of the anti-symmetric steady state for Re=25, (c) the jet-down steady state for Re=50 and (d) the jet-up steady state for Re=50.The contour values are scaled with respect to a maximum of ψ=2.2 for (b), which represents a transport of 5.5 Sv; and a maximum of ψ=1.1 for (c, d), which represents a transport of 10.9 Sv.The contour interval is 0.2.Table 1.Standard values of the parameters for the barotropic quasigeostrophic ocean model in the steady flow regime.For these values of the parameters, the dimensional parameters have values α τ =β=2.8×10 3 .

Parameter
Value called "optimal" (Mu et al., 2003).The CNOP can be regarded as the most (nonlinearly) unstable initial perturbation superposed on the basic state.
To numerically calculate the CNOP for the double-gyre problem, the kinetic energy norm Let L be a linear operator that maps x 0 to the velocity vector which is calculated from the values of ψ at four neighbouring grid points.The energy norm is now evaluated numerically as: where • 2 the L 2 -norm.Hence, the energy norm is calculated by multiplying the perturbation x 0 with the matrix L and calculating the norm of the result.Using this numerical approximation and Eq. ( 7) we find the following implementation J num of J : It is easy to derive the gradient: where M is the tangent linear model and the subscript T indicates the transpose.The constraint function is numerically implemented likewise as: The cost function is evaluated using forward integration.Here, the perturbation x 0 is added to the basic steady state x.This state is then propagated forward.At t=t e the basic steady state is subtracted and the kinetic energy norm of the perturbations is calculated.The gradient is evaluated by backward integration with the adjoint model M T .The same techniques developed for the implicit data assimilation (Terwisscha van Scheltinga and Dijkstra, 2005) are used here.We use the same tangent linear model, which is stored during the forward integration.For the evaluation of the gradient the stored tangent linear model is transposed and used as adjoint.
The constraint optimization problem ( 8) is implemented as: and is solved using the NAG routine E04UCF.This routine uses a method that is essentially identical to the method discussed in Gill et al. (1986).The basic structure uses an Sequential Quadratic Programming method (Gill et al., 1981) to solve a quadratic subproblem along the search direction and uses Lagrangian multipliers for the constraints.

Results
In the results below, we first consider the case of symmetric wind forcing (σ =0) and subsequently the case of a slightly asymmetric wind forcing (σ =0.05).

Symmetric case
For the case σ =0, the structure of the steady solutions is shown through the bifurcation diagram in Fig. 1a, where the value of the asymmetry of streamfunction , defined as is plotted against Re=U L/A H .At large values of A H (small Re), the anti-symmetric double-gyre flow (Fig. 1b) is a unique state.When lateral friction is decreased, this flow becomes unstable at the pitchfork bifurcation P 1 and two branches of stable asymmetric states appear for smaller values of A H (larger Re).The solutions on these branches (Fig. 1c, d) have the jet displaced either northward (negative ) or southward (positive ) and are exactly symmetrically related for the same value of Re.
We first focus on the case Re=25 which is in the asymptotically stable regime of the anti-symmetric state.For δ=0.1, the streamfunction patterns of the CNOPs are shown for four different values of t e =7, 14, 21, 28 days in Fig. 2. The pattern of the CNOP in Fig. 2a is similar to the pattern of the most energetic disturbance at a similar value of β (Fig. 7b of Dijkstra and De Ruijter, 1996) as determined through the energy stability analysis.As the energy stability boundary Re E is located near Re E ≈10 (Fig. 5 of Dijkstra and De Ruijter, 1996), this shows that in the conditional stability regime (Joseph, 1976)   still relevant for short times t e .When the CNOP is computed for larger times t e , the pattern keeps the same symmetric three cell structure, but the cells at top and bottom extend in size (Fig. 2b-d).
The streamfunction patterns at t e days when the CNOPs in Fig. 2 are taken as an initial condition are shown in Fig. 3.The pattern of Fig. 3a results after 7 days when the steady state at Re=25 is perturbed with the pattern in Fig. 2a with δ=0.1.It has not changed much in shape but it becomes more localized into the western boundary region.The patterns of Fig. 3b-d arise at 14, 21 and 28 days when the steady state is perturbed with the pattern in Fig. 2b-d, respectively, with δ=0.1.The final pattern in Fig. 3d is recognized as the least stable normal mode, the P-mode in Simonnet and Dijkstra (2002).Hence, the CNOPs for larger times t e induce a response into the direction of the normal mode, which indeed controls the long time evolution behavior.
The energy norm of the perturbation at t e for the steady state at Re=25 and several values of δ=0.1, 0.25 and 0.50 is plotted in Fig. 4 as a function of t e .For δ=0.1, the values correspond to the amplitudes of the patterns in Fig. 3.For each value of δ, the energy norm flattens for larger times t e and the value increases with increasing δ.
We next consider the case Re=50 on the asymmetric branches for which two asymmetric steady solutions (jet-up state and the jet-down state) are linearly stable.For δ=0.1, the streamfunction patterns of the CNOP of the jet-down state is shown for t e =7 days in Fig. 5a.The pattern is no longer symmetric because of the asymmetry of the back- ground state and it already quite localized near the western boundary current region.When the steady state is perturbed with the CNOP and with δ=0.1, the deviation from the steady state after 7 days (Fig. 5b) shows a bipolar pattern resembling one phase of a Rossby basin mode.This is an oscillatory normal mode to which the steady state becomes unstable at slightly larger Re (Dijkstra and Katsman, 1997).Figure 5c-d shows the CNOP and its evolution after 7 days for the jet-up steady state at Re=50.The patterns are simply related to those in Fig. 5a-b by the reflection symmetry.There are very minor differences due to accuracy set in the minimization algorithm and the implicit time-stepping schemes.For both asymmetric steady states the curves of the final amplitude of the energy norm versus t e are the same due to the reflection symmetry.For Re=50 the time-scale of flow changes in the system is set by the gyre advection which is  a few years.For each δ the energy norm (not shown) shows a monotonic increase in the range of t e contrary to the case of Re=25 (Fig. 4) where saturation occurs over a period of a month.

Asymmetric case
When the wind stress is taken slightly asymmetric, an imperfect pitchfork bifurcation results as can be seen in the the bifurcation diagram for σ =0.05 in Fig. 6a.The jet-up solution is now continuously connected with the near anti-symmetric solution at small values of Re.On the other hand, the jetdown solution becomes an isolated branch.The asymmetric wind-stress forcing gives a preference for the jet-up solution since the easterlies in the northern part of the domain are slightly weaker than those in the southern part of the domain.The position of the saddle-node bifurcation (at Re≈54 for σ =0.05 Fig. 6a) shifts to larger values of Re when σ increases.Along the branches for σ =0.05 the value of the dimensionless viscous dissipation function is plotted in Fig. 6b.As can be seen, values of differ between the jet-up and jet-down solutions for similar values of Re, with the (stable) jet-down steady state having a lower viscous dissipation.For σ =0.05 and Re=60, CNOPs for fixed δ=0.1 and t e =7 days are plotted for both jet-up and jet-down solutions in Fig. 7. Patterns now slightly differ between the two cases.The pattern for the jet-down state seems less deformed from the symmetric case (compare Fig. 7a with Fig. 5a).The CNOP pattern for the jet-up solution on the contrary has deformed substantially (compare Fig. 7c with Fig. 5c).The evolution of both CNOPs eventually leads to anomalies which have a pattern resembling a Rossby-basin mode, just as in the symmetric case.
For both states the final amplitude of the energy norm is plotted against t e for several values of δ in Fig. 8.The solid curves are those for the jet-up steady state while the dashed ones are those for the jet-down solution.For all values of δ there is a clear difference between the CNOP evolution from both states.For t e <16 days, the final amplitude of the energy norm is the largest for the jet-down state, i.e. the state with the lower viscous dissipation; for t e >16 days the opposite occurs.In these results, equilibration of the amplitude of the perturbations occurs on a longer (advective) time scale.

Conclusions
In this paper, we have explored the development of finite amplitude perturbations of linearly stable steady states of the double-gyre flow in the barotropic quasi-geostrophic model, by determining the Conditional Nonlinear Optimal Perturbations (CNOPs).These are the perturbations to the flow which have an optimal nonlinear evolution at a time t e (in a chosen norm) under the condition of a bound δ on the norm of the initial perturbation.The i4D-Var methodology as presented in Terwisscha van Scheltinga and Dijkstra (2005) was easily adapted to compute these CNOPs efficiently and hence provides a technique to determine CNOPs for fairly general systems of partial differential equations.By calculating the CNOPs for the symmetric (σ =0) double-gyre flow, we have added another detail of the behavior of this flow system when Re is changed.Up to the energy stability boundary Re E ≈10 (as determined in Dijkstra and De Ruijter (1996)) the anti-symmetric flow is monotonically stable, i.e., the kinetic energy of every finite amplitude perturbation decays monotonically to zero.Just above Re E , there exist perturbation patterns of which the kinetic energy grows in time and the CNOPs are the ones with optimal growth under the conditions of chosen t e and δ.The patterns of these CNOPs are basin wide and their spatial structure correspond to the ones of the non-normal modes as found in Moore et al. (2002).For small δ, the growth of the CNOPs is similar to that of the non-normal modes but for large δ it may be larger.For Re<Re L (the first pitchfork bifurcation) these CNOPs evolve in time to patterns resembling the least stable normal modes.Certainly, as soon as Re>Re L , the anti-symmetric state becomes linearly unstable and the perturbations with the largest growth rates are the normal modes.Fig. 8.The energy norm of the perturbation at t e , J (x 0δ ) as defined by Eq. ( 12), against t e , for different δ for the steady state at Re=60 and σ =0.05.The solid curves are for the jet-up solution while the dashed curves are for the jet-down solution.
For Re>Re L , two asymmetric linearly stable steady states exists which are (by their simultaneous existence) unstable to finite amplitude perturbations.The CNOPs for both solutions (for the same Re) are symmetry related and these patterns project during evolution on the normal mode patterns (associated with the first Hopf bifurcation on the asymmetric branches) which is most clearly seen at large evolution times t e .For the slightly asymmetric case, we showed that the growth of finite amplitude perturbations is different for the jet-up and jet-down steady states at similar Re.The physics of this difference is likely related to differences in the value of the viscous dissipation function of each steady state.
The separatrices (of attraction basins) are very difficult to calculate for the double-gyre flows; for the 60×40 grid used here a system with 4800 degrees of freedom results.However, the CNOPs may provide information on the finite amplitude stability boundaries in multiple equilibrium regimes.One can vary δ at fixed t e and determine for which critical δ the time evolution of the CNOP will not return to the original steady state.Such finite amplitude stability boundaries were determined in Mu et al. (2004) for a simple box-model (with 2-degrees of freedom).As this is not an easy computation for the double-gyre flow, with a very large CPU time needed for the minimization process, it is outside the scope of this paper.
Fig. 1.(a) Bifurcation diagram for the double-gyre (σ =0) for a square basin with the asymmetry of the streamfunction, defined as Eq.(16), against the control parameter Re=U L/A H .The energy stability boundary Re E is about Re E ≈10(Dijkstra and De Ruijter, 1996).(b) Streamfunction ψ of the anti-symmetric steady state for Re=25, (c) the jet-down steady state for Re=50 and (d) the jet-up steady state for Re=50.The contour values are scaled with respect to a maximum of ψ=2.2 for (b), which represents a transport of 5.5 Sv; and a maximum of ψ=1.1 for (c, d), which represents a transport of 10.9 Sv.The contour interval is 0.2.

Fig. 2 .
Fig. 2. Patterns of the barotropic streamfunction ψ for the CNOPs of the steady state at Re=25 for δ=0.1 and (a) t e =7 days, with an absolute maximum of 0.20; (b) t e =14, with an absolute maximum of 0.17; (c) t e =21 days, with an absolute maximum of 0.13 and (d) t e =28 days, with an absolute maximum of 0.11.

Fig. 5 .
Fig. 5.For Re=50, σ =0.0, δ=0.1 and t e =7 days: (a) CNOP for the jet-down steady state; (b) deviation of the flow from the jetdown steady state at t=t e ; (c) CNOP for the jet-up steady state; (d) deviation of the flow from the jet-up steady state at t=t e .

Fig. 6 .
Fig. 6.Bifurcation diagram for σ =0.05 where the asymmetry of the streamfunction ψ plotted against control parameter Re.(b) Dimensionless viscous dissipation along the branches in (a).

Fig. 7 .
Fig. 7.For the case Re=60, σ =0.05, δ=0.1 and t e =7 days: (a) CNOP for the jet-down steady state; (b) deviation of the flow from the jet-down steady state at t=t e ; (c) CNOP for the jet-up steady state and (d) deviation of the flow from the jet-up steady state at t=t e .