Transitions of zonal flows in a two-layer quasi-geostrophic ocean model

We consider a 2-layer quasi-geostrophic ocean model where the upper layer is forced by a steady Kolmogorov wind stress in a periodic channel domain, which allows to mathematically study the nonlinear development of the resulting flow. The model supports a steady parallel shear flow as a response to the wind stress. As the maximal velocity of the shear flow (equivalently the maximal amplitude of the wind forcing) exceeds a critical threshold, the zonal jet destabilizes due to baroclinic instability and we numerically demonstrate that a first transition occurs. We obtain reduced equations of the system using the formalism of dynamic transition theory and establish two scenarios which completely describe this first transition. The generic scenario is that two modes become critical and a Hopf bifurcation occurs as a result. Under an appropriate set of parameters describing midlatitude oceanic flows, we show that this first transition is continuous: a supercritical Hopf bifurcation occurs and a stable time periodic solution bifurcates. We also investigate the case of double Hopf bifurcations which occur when four modes of the linear stability problem simultaneously destabilize the zonal jet. In this case we prove that, in the relevant parameter regime, the flow exhibits a continuous transition accompanied by a bifurcated attractor homeomorphic to $S^3$. The topological structure of this attractor is analyzed in detail and is shown to depend on the system parameters. In particular, this attractor contains (stable or unstable) time-periodic solutions and a quasi-periodic solution.


Introduction
Baroclinic instability is among the most important geophysical fluid dynamical instabilities playing a crucial role in the dynamics of atmospheres and oceans. In particular, this instability mechanism is the dominant process in atmospheric dynamics shaping the cyclones and anticyclones that dominate weather in mid-latitudes, as well as the mesoscale ocean eddies that play various roles in oceanic dynamics and the transport of heat and salt [36]. Much is known on the linear stability of zonal jets in a horizontally unbounded ocean in the quasigeostrophic (QG) flow regime. Classical models, such as the continuously stratified Eady model [14] and the Date: May 18, 2022. two-layer Phillips model [27], have lead to a detailed understanding of the mechanism of baroclinic instability of a zonal jet in the inviscid case. Long waves destabilize the zonal jet with maximum growth rates occurring for perturbations having wavelengths on the order of the Rossby deformation radius, typically 50 − 100 km for the mid-latitude ocean [35].
In case linear friction is included in the two-layer model, the neutral curve has a minimum at (k c , µ c ) where k c is the critical wavenumber and µ c the critical value of the control parameter (e.g. the maximum speed of the zonal jet). The nonlinear development of these perturbations has been extensively analyzed in the weakly nonlinear case [25,28,39,18]. In the regime |k − k c | = O( ) and |µ − µ c | = O( 2 ), [37] showed that on a long time scale T = 2 t and large spatial scale X = (x − c g t), where c g is the group velocity of the waves at criticality, the complex amplitude A of the wave packet destabilizing the jet satisfies a Ginzburg-Landau equation, written as where the γ i are complex constants. [37] also showed that the fixed point solution of this equation can become unstable to sideband instabilities. Subsequent analysis has shown [13] that upgradient momentum transport can occur due to the self-interaction of the instabilities leading to rectification of the zonal jet.
In reality, the ocean basins are zonally bounded by continents and the midlatitude zonal jets are part of the gyre system, for example the subpolar gyre and subtropical gyre in the North Atlantic, forced by the surface wind stress through Ekman pumping [26]. The problem of baroclinic instability of such non-parallel flows is much more complicated and has so far only been tackled numerically. When the wind-forced QG equations are discretized, the linear stability problem for the gyre flow results in a large-dimensional generalized eigenvalue problem, typically of dimension 10 4 . There are many results for the one-layer single-and double-gyre flows (for an overview, see [11,Chap. 5]), but in this case there is no baroclinic instability. There are relatively few results for the two-layer case. In [12], it was shown that in the two-layer case the double-gyre flow becomes unstable through a sequence of Hopf bifurcations. The perturbation flow patterns at criticality are "banana-shaped," locally resembling those of baroclinic instability in the Philips model. Stable periodic orbits result from these Hopf bifurcations, typically given rise to meandering motion of the gyre boundary.
As an intermediate, more analytically tractable case, we consider here the baroclinic instability of a zonal jet for a two-layer QG model in a zonally periodic channel. In this case, the properties of the bounded geometry are somehow represented, as the patterns of the unstable modes are restricted by the periodicity of the channel, so a sequence of Hopf bifurcations is expected just as in the more realistic gyre case. In addition, parallel flow solutions exist in the zonally periodic channel which simplifies the linear stability problem substantially such that a more detailed nonlinear analysis, akin to that in the horizontally unbounded case, can be performed. The parallel flow can also be connected to the surface wind stress, as in the full gyre case, but at the expense of adding an additional linear friction term to the upper layer vorticity equation; for more details, see Section 2 below. In this way, the situation studied is more relevant for the stability of the Antarctic Circumpolar Current, than for the midlatitude gyre circulations.
The case specifically studied in the paper is the circulation set up by a time-independent Kolmogorov windstress field (for k = 1, 2, . . .) where τ 0 is a characteristic mid-latitude wind-stress value. This wind stress forces an ocean in a periodic channel [0, R/(2L x Z)] × [−L y , L y ] on the β-plane. The case k = 1 and k = 2 are often referred to as the single-and double-gyre forcing. The stratification is modeled in terms of a two-layer system and the wind stress only directly forces the upper layer. As a response to this wind stress, the system supports a basic shear flow ψ s . The amplitude τ that controls the wind-stress curl, or equivalently the maximal velocity of the shear flow ψ s is chosen as the bifurcation parameter. We first perform a numerical linear stability analysis of this basic shear flow; for small values of τ , all associated eigenvalues have negative real parts such that the jet is stable. When the aspect ratio of the channel a = L y /L x is large, the eigenvalues remain in the left complex plane regardless of the value of τ . However, when the aspect ratio gets small, the basic shear flow loses stability at a critical τ in the form of a single conjugate pair or two conjugate pairs of eigenvalues crossing the imaginary axis, giving rise to either a Hopf or a double Hopf bifurcation. We next use the idea and method of the dynamic transition theory [20,21], which is aimed to determine all the local attractors near a transition. The approach allows for a classification near the instability onset of all transitions into three classes known as continuous, catastrophic and random types [21]. In this way, our study extends previous results using this approach on the single-layer barotropic case [32,10], the two-layer case for constant zonal jet velocities [3] and the barotropic Munk western boundary layer current profile case [15], to the zonally periodically bounded two-layer case.
Using the center manifold reduction, we obtain effective reduced (ordinary differential equation, ODE) models describing this transition. The dynamic transition theory identifies then transition numbers that qualifies the transition's type and are calculated from the reduced ODE's coefficients. The case of a Hopf bifurcation is generic while the case of a double Hopf bifurcation is degenerate and requires fine tuning of the aspect ratio to critical values where two conjugate pairs of eigenmodes with consecutive wavenumbers have their eigenvalues crossing the imaginary axis simultaneously. Using standard parameter values describing the midlatitude related oceanic flows, we perform numerical computations of the transition number for the forcing patterns corresponding to k = 1, 2, 3 and the aspect ratios a ≥ 3. We find that in the parameter regimes we are interested in, the Hopf bifurcation is supercritical and a stable limit cycle bifurcates. For the double Hopf bifurcation, we find that after the corresponding transition takes place, the system exhibits a bifurcated local attractor [20] near the basic shear flow which is homeomorphic to the 3D-sphere. The topological structure of this attractor is analyzed and depending on the parameters, it is found to contain a combination of limit cycles and a quasi periodic solution.
The paper is organized as follows. In Section 2, the quasi-geostrophic model is presented. This is followed by Section 3 where the theory and numerical results for the linear stability problem (Section 3.1), the Hopf bifurcation case (Section 3.2) and the double-Hopf bifurcation case (Section 3.3) are presented. These results are summarized and discussed in Section 4. Appendix A and Appendix B contains details regarding the proofs of the transition theorems and the explicit formulas of the corresponding reduced ODE's coefficients, in the Hopf and double-Hopf cases, respectively. Appendix C and Appendix D, as for them, provide details about for the numerical treatment of the linear stability analysis and the practical computation of the transition numbers. Finally the set of model parameters used in the numerical study of the problem is given in Appendix E.

The model
We consider two layers of homogeneous fluids, each with a different and constant density ρ 1 and ρ 2 and with equilibrium layer thicknesses H 1 and H 2 , on a mid-latitude β-plane with Coriolis parameter f = f 0 + β 0 y. The lighter fluid in layer 1 is assumed to lie on top of the heavier one in layer 2 so that the stratification is statically stable, i.e., ρ 1 < ρ 2 ; bottom topography is neglected.
This flow can be modeled by the two-layer QG model [25] using the geostrophic stream function ψ i and the vertical component of the relative vorticity ζ i in each layer (i = 1, 2). The quantities ψ i and ζ i are nondimensionalised by U L y and U/L y , respectively, wind stress with τ 0 , length with L y , and time with L y /U , where U is a characteristic horizontal velocity.
By choosing U = τ 0 /(ρ 0 β 0 L y H 1 ), where ρ 0 is a reference density, the dimensionless equations on the domain (0, 2/a) × (−1, 1) (with a = L y /L x ) become where {f, g} = f x g y − f y g x is the usual Jacobian operator and represents the damping of upper layer vorticity due to frictional processes. In the bottom layer, we include a linear (Ekman) friction term −r 2 ∆ψ 2 ; in both layers, Laplacian friction terms are neglected due to the absence of continental boundary layers making such terms much smaller than the other ones. The expressions for the dimensional and dimensionless parameters, with their standard values at a latitude 45 • N, are given in Table 2.
For the boundary conditions, we assume periodicity in the x-direction and free-slip boundaries in the ydirection. Hence, the conditions are In actual ocean basins, a steady zonal jet is generated by the applied wind stress through Ekman pumping, a Sverdrup balance and a western boundary layer flow [26]. Due to the periodic boundary conditions used here, such a flow cannot be captured in this model. However, due to the presence of the large upper layer friction, which corrects for the absence of the Sverdrup balance, the equations allow a steady state of the form (2.3) ψ s 1 = Ψ sin kπy, ψ s 2 = 0, which relates to the wind stress field, when F 1 − τ β sin kπy = 0. In this paper, we will assume that the wind-stress vorticity input is balanced by vorticity decay due to the linear friction term F 1 = −r 1 ∆ψ 1 , being aware that a larger friction coefficient r 1 is needed than can be justified from existing dissipative processes in the ocean. In this case, it follows that The parameter Ψ appearing in (2.3) can then be chosen as the control parameter as is the case in this study, instead of τ . It can be interpreted as the maximal (zonal) velocity of the shear flow (after taking the derivative in y of (2.3)) or equivalently the maximal amplitude of the wind forcing according to (2.4). By considering the perturbation ψ i = ψ i − ψ s i , i = 1, 2, we can write the system into the following operator form (after dropping the prime notation) where M and N are the linear operators defined as Lastly, the bilinear nonlinearity is given explicitly by The operators G and N are well-defined mappings, G : H 1 → H −1 and N : H 0 → H −1 , on the following functional spaces: Here, and H 4 (Ω), H 2 (Ω), L 2 (Ω) are the usual Sobolev and Lebesgue function spaces endowed with their natural inner products. These functional spaces account for spatial regularity of the solution ψ for which H p (Ω) denotes the space of square-integrable functions that possess pth-order derivatives (in the distribution sense) that are themselves square-integrable; see e.g. [2]. Using this functional framework, it can be shown that Eq. (2.5) is well-posed by application of e.g. a standard approximation method (such as Galerkin) [1].

Results
In this section, we first present the linear stability analysis of the basic shear flow and then we move on to describe the first transitions due to the instabilities, covering both the Hopf and double Hopf bifurcations. Although in realistic ocean basins the aspect ratio a would be small, we allow here the full range of a to also study the nonlinear interactions of localized instabilities.
3.1. Linear Stability Analysis. We first investigate the linear stability of the basic solution. For this purpose, we denote the eigenmodes of the linear problem by with eigenvalues σ m,j , i.e.
Since the linear operators M and N are real, we have This eigenvalue problem is solved numerically by means of a standard Legendre-Galerkin method; see Appendix C. A typical picture of the spectrum near the criticality is given in Figure 1. This figure shows that many eigenvalues are clustered near the imaginary axis at the critical parameter value Ψ = Ψ c as defined in (3.3) below.  We assume (as confirmed numerically for the parameter regimes considered below) that the eigenvalues are ordered so that for each integer m, σ m,1 has the largest real part among the σ m,j , for any nonnegative integer j.
For each (nonnegative) wavenumber m, we define Ψ m , when it exists, to be the value of Ψ for which the eigenvalue σ m,1 crosses the imaginary axis, that is Hence Ψ = Ψ m defines a neutral stability curve in the (a, Ψ)-plane. In Figure 2, these neutral curves are plotted for zonal wave numbers m = 1, 2, 3, 4.
Our numerical analysis suggests that for any m, Ψ m is well-defined only for aspect ratios of the basin characteristic lengths smaller than a threshold a m (depending on m), that is for a < a m . The threshold a m is defined by a vertical asymptote condition, namely Moreover, we numerically observe that the a m are ordered such that We define then the critical maximal shear flow's velocity, Ψ c , and the critical zonal wavenumber, m c , by and,

respectively.
A typical structure of the spectrum at the critical parameter value, Ψ = Ψ c , is shown in Figure 1 where a conjugate pair of eigenvalues is about to cross the imaginary axis. The eigenvalues on the real axis correspond to the wavenumber m = 0 and are always stable although they may be very close to zero as shown in Figure 1.
To describe the solutions near the onset of transition Ψ = Ψ c , we define the spatio-temporal function where σ m,1 is the first eigenvalue and ψ m,1 is its associated eigenfunction. The spatial structure of the eigenmodes ψ mc,1 is shown in Figure 3, revealing the well-known "banana-shaped" patterns characteristics of baroclinic instability.
The values of Ψ c with respect to the aspect ratio a for k = 1, 2, 3 is shown in Figure 4. By the previous remarks, By (3.6), for a > a 1 , the system is linearly stable. As is expected, the neutral stability curves (Figure 2) approach the asymptote Ψ c → ∞ as a converges to the critical aspect ratio a 1 over which the system is linearly stable for all Ψ. Our numerical results in Figure 4 show that for aspect ratios a in the range 3 ≤ a ≤ 20, the critical maximal shear velocity Ψ c ≈ 0.1, 0.04, 0.03 for k = 1, 2, 3 respectively. By (2.4), this value of critical maximal shear velocity corresponds to an upper layer friction that is approximately, which is indeed much larger than can be justified from dissipative processes in the ocean but, as explained in Section 2, is needed here to connect the background zonal jet to the applied wind stress as a Sverdrup balance is absent in the original model formulation.
The friction term in the lower layer however is physical (Ekman friction) and for this study, it is fixed at r 2 = 5.0, see Table 2. Also, from Figure 2, we see that for small aspect ratios, many modes become unstable as Ψ c is exceeded.
For a < a 1 , the system has a first transition at Ψ = Ψ c and exactly one of the following two principal of exchange of stability (PES) condition holds:  Table 1. The double Hopf transition numbers. Here k = 1, 2 is the wavenumber of the forcing, m c and m c + 1 are the wavenumbers of the first two critical modes which become unstable simultaneously at the critical aspect ratio a DH . A, B, C, D are the coefficients of the double Hopf transition normalized with respect to the maximum absolute value of those coefficients. The parameters δ and θ are defined by (3.14).
According to (3.7) and (3.8) either, two or four eigenvalues become unstable as Ψ crosses Ψ c . The case of two critical eigenvalues is generic and results into a Hopf bifurcation. The case of four critical eigenvalues results in a double Hopf bifurcation and requires the fine-tuning of the aspect ratio a = a DH so that Ψ c = Ψ mc = Ψ mc+1 . The values of a DH where double Hopf transition occurs is given in Table 1.
Although the double Hopf transition is not generic, its analysis gives an insight into regimes of moderate values of Ψ where multiple eigenvalues are unstable. In Figure 5, we show the dominant part of the spectrum of the linearized operator at a critical aspect ratio a DH .
We recall that the neutral stability curves (shown in Figure 2) are found by identifying the set of parameter values (here in the (a,Ψ)-plane) for which the real part of the critical eigenvalue is zero while the rest of the eigenvalues have a negative real part, i.e. are associated with a stable mode. By a continuity argument, we can thus infer that in a neighborhood of such critical curves the PES condition is satisfied, for instance when Ψ crosses the critical value Ψ c (depending on a and k) shown in Figure 2. The PES condition (3.7) has been rigorously verified for Kolmogorov flows in [22] via a continued fraction method. This method has later been extended for the single layer QG model for the k = 1 case in [32] and for k ≥ 2 in [19] where k is the forcing frequency in (2.1). It is still an open problem to rigorously verify the PES condition for the current problem.   The analysis of Lemma 3.1 yields the following theorem.
Theorem 3.1. Assume that the first critical eigenvalue is purely complex with a simple multiplicity, that is the PES condition (3.7) is satisfied. Then the following assertions hold true.
(1) If Re(P ) < 0, then Eq. (2.5) in H 1 (see (2.9)) undergoes a continuous transition accompanied by a supercritical Hopf bifurcation on Ψ > Ψ c , with Ψ c defined in (3.3). In particular, the steady-state solution bifurcates to a stable periodic solution ψ on Ψ > Ψ c , satisfying ψ → 0 as Ψ → Ψ c and has the following approximation (3.10) ψ(x, y, t) = − Re(σ mc,1 ) Re(P ) The spatial structure of the time periodic solution ψ is shown at t = 0 for different aspect ratios a in Figure 3. (2) If Re(P ) > 0, then Eq. (2.5) in H 1 (see (2.9)) undergoes a jump transition on Ψ < Ψ c accompanied by a subcritical Hopf bifurcation. In particular, an unstable periodic solution ψ given by (3.10) bifurcates on Ψ < Ψ c and there is no periodic solution bifurcating from 0 on Ψ > Ψ c . Moreover, there is a singularity separation at some Ψ s < Ψ c generating an attractor and an unstable periodic solution ψ.
When the PES condition (3.7) holds, the system exhibits a Hopf bifurcation as described by Theorem 3.1. The type of transition boils down to the determination of the transition number P given in (A.8) of Appendix A. For the practical calculation of this number, we refer to Appendix D. A numerical evaluation of P as reported in Figure 6 shows that for a low-frequency forcing (k = 1, 2, 3), generally Re(P ) < 0 and as a result only continuous transition (supercritical Hopf bifurcation) is possible for the parameter regime we have selected. In Figure 6, we also display the critical wavenumber m c . In the range 4 ≤ a ≤ 20, the critical wavenumber m c is found to range from 1 to 4.  Table 2. There are discontinuities in P vs a plot in Figure 6 of the transition number which are due to the changes in the critical zonal wavenumber m c . These discontinuities take place at double Hopf bifurcation aspect ratios where two consecutive zonal wavenumbers become critical simultaneously which is investigated in the next section. However, there are also discontinuities in Figure 6, k = 1 case (for example near a = 16) whose origin is mysterious.
As detailed in Appendix A, the transition number P accounts for two types of nonlinear interactions between the eigenmodes, and is written P = P 0 + P 2 , where P 0 accounts for nonlinear interactions between the critical modes and the zonally homogeneous (stable) modes m = 0 (see (A.9) below), while P 2 accounts for the interactions between the critical modes and the modes having a wavenumber twice that of the critical modes (see (A.10) below). A comparison of typical numerical values of P 0 and P 2 shows that P 2 is several orders of magnitudes smaller than P 0 ; see Figure 7. We refer to [4, Theorem III.1] for a similar transition number diagnosing also the type of Hopf bifurcations arising, generically, in delay differential equations, and whose nature is also characterized by the interactions of linearized modes through the model's non-linear terms.
In our case, where P 0,j measures the contribution to P of the j-th mode with zero-wavenumber (m = 0) interacting with the critical mode; see (A.9) again. Figure 8 shows that the dependence of P 0,j on j is essentially linear. We believe the results in Figure 7 and Figure 8 may help when choosing the modes to include in a simulation when the maximal shear velocity is well above the criticality.  Ly U compared to the (dimensional) length of the channel L dim = 2L y /a of the both stable and unstable bifurcated time periodic solution (3.10) in the Hopf bifurcation case. Here L y and U are the characteristic scales defined in Table 2. The jumps in the derivative of the time period of the bifurcated solution is due to the change of the imaginary part Im σ mc,1 of the critical eigenvalue at the double Hopf aspect ratios a DH .
We also compare the dimensional time period of the bifurcated solution (3.10) to the (dimensional) length of the channel in Figure 9. With the default parameters as chosen in Table 2, our simulations yield a solution with time period of 180-380 days depending on the channel length of 100-700kms.
3.3. Double Hopf Bifurcation. In this section we are interested in the transitions that take place at the critical aspect ratios a DH where four modes with consecutive wavenumbers m c , m c + 1 become unstable as given by the PES condition (3.8).
We first present the reduced equations in this case (for proofs, see Appendix B).
where z 1 (t) and z 2 (t) denote the complex amplitudes aiming at approximating the projection of the model's solution onto the critical modes ψ mc,1 and ψ mc+1,1 , respectively. Figure 10. The regions in the λ 1 -λ 2 plane with different dynamical behaviors. In region V, the basic steady state is locally asymptotically stable. In regions IV and VI, the system undergoes a supercritical Hopf bifurcation. The dynamics in regions I, II and III is the double Hopf bifurcation scenario and the details are given in Figure 11. The lines T 1 and T 2 in the figure have slopes 1/θ and δ (as defined in (3.14)) respectively.
The transition numbers A, B, C, D are determined by the nonlinear interactions of these critical modes with higher modes given by (B.2). More precisely, the terms A and D account for the self-interactions among the critical modes, while the terms B and C account for the cross-interactions between the critical modes with the higher modes.
It is known that the equation (3.11) exhibit a zoo of dynamical behaviors. We refer to [16] for a detailed analysis of all possible cases. Here, we restrict our attention to the case which is the only case we observe in our numerical experiments, see Table 1. Under these conditions it is known that the transition is continuous (see Theorem 2.3 in [17]). For the next theorem, let us define the numbers (3.13) λ 1 = Re(σ mc,1 ), λ 2 = Re(σ mc+1,1 ).
Recalling f mc defined in (3.5) (with m = m c ), we define the following spatio-temporal profiles Theorem 3.2. Assume that the conditions of Lemma 3.2 as well as the condition (3.12) hold. Then the equations (2.5) undergo a continuous transition at Ψ = Ψ c , and an S 3 attractor Σ bifurcates on Ψ > Ψ c , which converges to 0 as Ψ approaches to Ψ c from right. Depending on the values of θ and δ, there are three transition scenarios as shown in Figure 10. In each scenario, near the onset of transition (λ 1 , λ 2 ) = (0, 0), the λ 1 − λ 2 plane is dissected into several regions with distinct topological structures for the attractor Σ as given in Figure 11.  Figure 11. The dynamics in the regions given in the first quadrant of Figure 10. ψ mc p , ψ mc+1 p are time-periodic with zonal wavenumbers m c and m c + 1 respectively. ψ qp is the quasiperiodic solution given in (3.15).
(1) If z 2 = 0 or z 1 = 0, the equations (3.11) reduce to the equation (3.9) with A = P or D = P respectively. Thus Lemma 3.1 and Theorem 3.1 are special cases of Lemma 3.2 and Theorem 3.2 respectively.
(2) We note that the features of the spatial structures of upper vs lower layer of the bifurcated periodic solutions in Figure 3 and the quasi-periodic solution in Figure 12 do not alter much. We expect that the situation would be different if bottom topography is included.
The transition scenario of double Hopf transition is given by Theorem 3.2 by Figure 10 and Figure 11. We find that near the onset of transition, depending on the fluctuations the basic state experiences a transition either to a time periodic solution or to a quasi periodic solution. Our results in Table 1 show that the three scenarios sketched in Figure 10 are realizable.
In particular, near a double Hopf transition point, one of the following three possibilities must occur post transition, Ψ > Ψ c : (i) there is only a single stable limit cycle, (ii) there are two distinct stable limit cycles, and an unstable quasi-periodic solution (iii) there is a stable quasi-periodic solution and either one or two unstable limit cycles. For the double Hopf transition, Theorem 3.2 basically tells that all of the above local structures, the time-periodic solution and the 2D torus, if they exist, reside in a local attractor homeomorphic to the three dimensional sphere. The existence of this attracting 3D sphere is guaranteed by the attractor bifurcation theorem; see [20, Theorem 6.1].

Summary and Discussion
In this paper, we investigated the stability of a parallel zonal jet forced by a Kolmogorov-type wind stress in a periodic zonal channel, using a two-layer quasi-geostrophic (QG) model. This problem is, in terms of complexity, situated between the horizontally unbounded problem [27,37] and the fully bounded gyre problem [12]. More precisely, the effect of boundaries is captured by the interactions of only a few modes (solutions of the linear stability problem), while still keeping a parallel flow which is connected to the wind-stress field.
Our results show, in the continuity of earlier works [15,3,10,19], that the zonal shear flow is linearly stable if its maximal amplitude Ψ, or equivalently the maximal amplitude τ of the wind-stress curl (see (2. Our approach allows for computing the coefficients of the Hopf reduced equation (3.9) by means of explicit formulas involving the interactions between the critical mode and the (zonal) stable ones through the nonlinear terms; see (A.8)-(A.10). A numerical examination of these coefficients reveals that instead of infinitely many (stable) modes which would affect the type of transition (supercritical vs subcritical Hopf), the transition is in fact determined by the interactions between the conjugate pair of critical modes losing stability and only the first few zonally homogeneous (m = 0) modes. That the interactions with the zonal-mean modes are the determining factors governing the flow patterns emerging after the onset of instability may be viewed as consistent with other works that have shown (in more turbulent regimes) that baroclinic turbulence statistics can be well recovered for certain geophysical flows by omitting the effects of the interactions among the higher-order modes (eddy-eddy interactions); see e.g. [24,23,33,9].
We also investigated the double Hopf bifurcation scenario which takes place at critical length scales where four modes with consecutive wavenumbers become critical. By a rigorous center manifold analysis, we obtain the coefficients of the 4D-ODE system. Our results show that for the parameters we have considered, there exists a quasi-periodic solution that bifurcates. This quasi-periodic solution is a linear combination of two periodic solutions and may be stable depending on the parameters; see Section 3.3. From a transition point of view, in the double Hopf transition, an attractor homeomorphic to 3D sphere bifurcates; see Theorem 3.2. This attractor contains stable/unstable limit cycles and an stable/unstable invariant torus (supporting a quasiperiodic solution).
Overall, our results add more details to the nonlinear development of baroclinic instabilities on a non-constant parallel zonal jet, in that the periodic solutions can become unstable to torus bifurcations and give rise to quasiperiodic behavior. Such a scenario was also found for the barotropic double gyre flow [38], but only in a weakly nonlinear framework using a set of (reduced) amplitude equations. The transition scenario found for the zonally periodic zonal jet is likely to be more relevant for the ocean circulation than the sideband instabilities in the zonally unbounded zonal jet case which require a nearby band of wavenumbers to be unstable.
In our set-up, the linear friction coefficient in the upper layer is relatively large and as explained needed to balance the vorticity input by the wind stress for generating the zonal jet. When this friction is decreased, more modes will become unstable near the critical point and their interactions are expected to give a detailed view on the development of the unstable modes into ocean vortices, or eddies, due to baroclinic instabilities. In that respect, the dynamic transition theory adopted here, along with extensions based on the variational approach of [6,5] and its stochastic rectification [7] to handle regimes beyond the onset of instability, provides a way forward to develop a mathematical theory of such ocean-eddy formation processes. The results derived here combined with recent advances in the rigorous analysis of transitions arising in stochastically driven flows [8] open new prospects for the study of regime transitions in stochastically forced ocean models [34,29].
Appendix A. Proof of Lemma 3.1 and Theorem 3.1 We first proceed with the proof of Lemma 3.1 For this, we denote the adjoint modes by ψ * m,j = e iαmx Y * m (y). We denote the critical eigenmode and the critical eigenvalue by ψ c = ψ mc,1 , σ c = σ mc,1 .
We denote the bilinear operator G as The center part of the solution is where c.c. stands for complex conjugate of the terms before. The evolution of z(t) near the onset of transition is obtained by the projection onto the critical mode ψ c .
where Φ is the center manifold function. We will obtain its quadratic approximation Φ 2 given by Here o(n) = o (|(z, z)| n ) denotes higher than n-th order terms in z, z.
Using the notation (A.1), the reduced equation (A.3) can be written To obtain a closed system, we need to approximate the center manifold function. The approximation of the center manifold in this case reads, see [30], where L = Π s M −1 N and Π s is the projection on the stable space. Using the formula (A.5), we obtain the following expansion of the center manifold (see also [5,Theorem 2]) (A.6) Φ 2 = z 2 j≥1 g 2mc,j ψ 2mc,j + |z| 2 j≥1 g 0,j ψ 0,j + c.c.
Here (A.7) are the coefficients of the center manifold function. We write (A.4) as (3.9), that isż which finishes the proof of Lemma 3.1.
Recalling the definition of G s given in (A.1), the transition number P can then be written as denotes the contribution of the zero-wavenumber (stable) modes ψ 0,j while (A.10) denotes the contribution of the modes ψ 2mc,j on the transition number respectively. The transition type depends on the real part of the transition number P . The proof of Theorem 3.1 follows from the standard Hopf bifurcation analysis of the reduced equation.
Appendix B. Proof of Lemma 3.2 and Theorem 3.2 As the reduction in the case of (3.8) is similar to the case of (3.7) given in the previous section, we will only mention the differences between these two cases. Under the assumption (3.8), we write the center part of the solution as where the first two critical modes are ψ 1 = ψ mc,1 , ψ 2 = ψ mc+1,1 , ψ −1 = ψ −mc,1 , ψ −2 = ψ −mc−1,1 .

Appendix C. Numerical treatment of the linear stability problem
To solve the eigenvalue problem numerically, we first plugin the ansatz (C. 1) ψ(x, y) = e iαmx Y j (y), j ∈ N, m ∈ Z, α m := amπ.
into the eigenvalue problem σMψ(x, y) = N (y)ψ(x, y), to obtain (C.2) σ MY (y) = N (y)Y (y), Y (y) = (Y 1 (y), Y 2 (y)) Here the linear operators M and N (y) are defined as We use Legendre-Galerkin method to discretize and solve the (C.2) with boundary conditions (C.5). We refer to [31] for the details of the Legendre-Galerkin method and to [10] for its use in dynamical transition problems.
Let into (C.2) to obtain Here (C.9) A 4 = (cos kπyD 2 f j , f k ), A 5 = (cos kπyf j , f k ), with (f, g) = 1 −1 f (y)g(y) dy. The explicit expression of the matrices A i , i = 1, . . . , 5 can be found in [10].  Table 2. Set of model parameters used in the numerical study of the problem.