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

Mickael D. Chekroun  Department of Earth and Planetary Sciences, Weizmann
Institute, Rehovot 76100, Israel<EMAIL_ADDRESS>,
Henk Dijkstra  Centre for Complex System Studies, Department of Physics,
Utrecht University, Utrecht, The Netherlands<EMAIL_ADDRESS>,  Taylan
Şengül  Department of Mathematics, Marmara University, 34722 Istanbul, Turkey
<EMAIL_ADDRESS>and  Shouhong Wang  Department of Mathematics,
Indiana University, Bloomington, IN 47405<EMAIL_ADDRESS>http://www.indiana.edu/~fluid

###### Abstract.

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
a conjugate pair of modes loses stability and a Hopf bifurcation occurs as a
result. Under an appropriate set of parameters describing related 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.

###### Contents

  1. 1 Introduction
  2. 2 The model
  3. 3 Results
    1. 3.1 Linear Stability Analysis
    2. 3.2 Hopf Bifurcation
    3. 3.3 Double Hopf Bifurcation
  4. 4 Summary and Discussion
  5. A Proof of Lemma 3.1 and Theorem 3.1
  6. B Proof of Lemma 3.2 and Theorem 3.2
  7. C Numerical treatment of the linear stability problem
  8. D  Practical aspects for the calculation of the transition number
  9. E Model Parameters

##  1\. 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 quasi-geostrophic (QG) flow regime. Classical models,
such as the continuously stratified Eady model [14] and the 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},\mu_{c})$ where $k_{c}$ is the critical wavenumber
and $\mu_{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}|={\mathcal{O}}(\epsilon)$ and
$|\mu-\mu_{c}|={\mathcal{O}}(\epsilon^{2})$, [37] showed that on a long time
scale $T=\epsilon^{2}t$ and large spatial scale $X=\epsilon(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

$\frac{\partial A}{\partial
T}=\gamma_{1}A+\gamma_{2}\frac{\partial^{2}A}{\partial
X^{2}}-\gamma_{3}A|A|^{2},$

where the $\gamma_{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 wind-stress field (for $k=1,2,\ldots$)

$\tau^{x}(y)=-\tau\frac{\tau_{0}}{k\pi}\cos
k\pi\frac{y}{L_{y}}~{};~{}\tau^{y}=0$

where $\tau_{0}$ is a characteristic mid-latitude wind-stress value. This wind
stress forces an ocean in a periodic channel
$[0,\mathbb{R}/(2L_{x}\mathbb{Z})]\times[-L_{y},L_{y}]$ on the $\beta$-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 $\psi^{s}$. The amplitude
$\tau$ that controls the wind-stress curl, or equivalently the maximal
velocity of the shear flow $\psi^{s}$ is chosen as the bifurcation parameter.

We first perform a numerical linear stability analysis of this basic shear
flow; for small values of $\tau$, 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 $\tau$. However, when the aspect ratio gets small,
the basic shear flow loses stability at a critical $\tau$ 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\geq 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.

##  2\. The model

We consider two layers of homogeneous fluids, each with a different and
constant density $\rho_{1}$ and $\rho_{2}$ and with equilibrium layer
thicknesses $H_{1}$ and $H_{2}$, on a mid-latitude $\beta$-plane with Coriolis
parameter $f=f_{0}+\beta_{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., $\rho_{1}<\rho_{2}$; bottom topography is neglected.

This flow can be modeled by the two-layer QG model [25] using the geostrophic
stream function $\psi_{i}$ and the vertical component of the relative
vorticity $\zeta_{i}$ in each layer ($i=1,2$). The quantities $\psi_{i}$ and
$\zeta_{i}$ are non-dimensionalised by $UL_{y}$ and $U/L_{y}$, respectively,
wind stress with $\tau_{0}$, length with $L_{y}$, and time with $L_{y}/U$,
where $U$ is a characteristic horizontal velocity. By choosing
$U=\tau_{0}/(\rho_{0}\beta_{0}L_{y}H_{1})$, where $\rho_{0}$ is a reference
density, the dimensionless equations on the domain $(0,2/a)\times(-1,1)$ (with
$a=L_{y}/L_{x}$) become

(2.1) $\displaystyle\left[\frac{\partial}{\partial
t}+\\{\psi_{1},\cdot\\}\right]\left(\Delta\psi_{1}+F_{1}(\psi_{2}-\psi_{1})+\beta
y\right)=-r_{1}\Delta\psi_{1}-\tau\beta\sin k\pi y$
$\displaystyle\left[\frac{\partial}{\partial
t}+\\{\psi_{2},\cdot\\}\right]\left(\Delta\psi_{2}+F_{2}(\psi_{1}-\psi_{2})+\beta
y\right)=-r_{2}\Delta\psi_{2},$

where $\\{f,g\\}=f_{x}g_{y}-f_{y}g_{x}$ is the usual Jacobian operator and

${\mathcal{F}}_{1}=-r_{1}\Delta\psi_{1},$

represents the damping of upper layer vorticity due to frictional processes.
In the bottom layer, we include a linear (Ekman) friction term
$-r_{2}\Delta\psi_{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^{\circ}$N, are given in Table 2.

For the boundary conditions, we assume periodicity in the $x$-direction and
free-slip boundaries in the $y$-direction. Hence, the conditions are

(2.2) $\displaystyle\psi_{i}\mid_{x=0}=\psi_{i}\mid_{x=2/a},\qquad i=1,2.$
$\displaystyle\psi_{i}\mid_{y=\pm 1}=\frac{\partial^{2}\psi_{i}}{\partial
y^{2}}\mid_{y=\pm 1}=0,\qquad i=1,2.$

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) $\psi_{1}^{s}=\Psi\sin k\pi y,\qquad\psi_{2}^{s}=0,$

which relates to the wind stress field, when ${\mathcal{F}}_{1}-\tau\beta\sin
k\pi 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
${\mathcal{F}}_{1}=-r_{1}\Delta\psi_{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

(2.4) $\Psi=\frac{\tau\beta}{(k\pi)^{2}r_{1}}.$

The parameter $\Psi$ appearing in (2.3) can then be chosen as the control
parameter as is the case in this study, instead of $\tau$. 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 $\psi_{i}^{\prime}=\psi_{i}-\psi^{s}_{i}$,
$i=1,2$, we can write the system into the following operator form (after
dropping the prime notation)

(2.5)
$\mathcal{M}\partial_{t}\psi=\mathcal{N}\psi+\mathcal{G}(\psi),\qquad\psi=(\psi_{1},\psi_{2}),$

where $\mathcal{M}$ and $\mathcal{N}$ are the linear operators defined as

(2.6)
$\mathcal{M}\psi=\begin{bmatrix}\Delta\psi_{1}+F_{1}(\psi_{2}-\psi_{1})\\\
\Delta\psi_{2}+F_{2}(\psi_{1}-\psi_{2})\end{bmatrix},$

and

(2.7) $\mathcal{N}\psi=\begin{bmatrix}\Psi k\pi\cos k\pi
y\left((k\pi)^{2}\frac{\partial\psi_{1}}{\partial
x}+F_{1}\frac{\partial\psi_{2}}{\partial
x}+\frac{\partial\Delta\psi_{1}}{\partial
x}\right)-\beta\frac{\partial\psi_{1}}{\partial x}-r_{1}\Delta\psi_{1}\\\
-\Psi k\pi F_{2}{\cos k\pi y}\frac{\partial\psi_{2}}{\partial
x}-\beta\frac{\partial\psi_{2}}{\partial x}-r_{2}\Delta\psi_{2}\\\
\end{bmatrix}.$

Lastly, the bilinear nonlinearity is given explicitly by

(2.8)
$\mathcal{G}(\psi)=\begin{bmatrix}&-\\{\psi_{1},\Delta\psi_{1}+F_{1}(\psi_{2}-\psi_{1})\\}\\\
&-\\{\psi_{2},\Delta\psi_{2}+F_{2}(\psi_{1}-\psi_{2})\\}\end{bmatrix}.$

The operators $\mathcal{G}$ and $\mathcal{N}$ are well-defined mappings,
$\mathcal{G}:H_{1}\to H_{-1}$ and $\mathcal{N}:H_{0}\to H_{-1}$, on the
following functional spaces:

(2.9) $\displaystyle H_{1}=\\{\psi=(\psi_{1},\psi_{2})\in H^{4}(\Omega)\times
H^{4}(\Omega)|~{}\psi~{}\text{satisfies}~{}\eqref{BC}\\},$ $\displaystyle
H_{0}=\\{\psi=(\psi_{1},\psi_{2})\in H^{2}(\Omega)\times
H^{2}(\Omega)|~{}\psi~{}\text{satisfies}~{}\eqref{BC}\\},$ $\displaystyle
H_{-1}=L^{2}(\Omega)^{2}.$

Here,

(2.10) ${\Omega=(0,2/a)\times(-1,1),\textrm{ with }a=L_{y}/L_{x},}$

and $H^{4}(\Omega)$, $H^{2}(\Omega)$, $L^{2}(\Omega)$ are the usual Sobolev
and Lebesgue function spaces endowed with their natural inner products. These
functional spaces account for spatial regularity of the solution $\psi$ for
which $H^{p}(\Omega)$ denotes the space of square-integrable functions that
possess $p$th-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].

##  3\. 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

$\psi_{m,j}(x,y)=e^{i\alpha_{m}x}Y_{j}(y),\qquad
j\in\mathbb{N},m\in\mathbb{Z},\quad\alpha_{m}=am\pi.$

with eigenvalues $\sigma_{m,j}$, i.e.

(3.1) $\sigma_{m,j}\mathcal{M}\psi_{m,j}=\mathcal{N}\psi_{m,j}.$

Since the linear operators $\mathcal{M}$ and $\mathcal{N}$ are real, we have

$\overline{\sigma_{m,j}}=\sigma_{-m,j},\qquad\overline{\psi_{m,j}}=\psi_{-m,j},\qquad\forall
m\in\mathbb{N}.$

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
$\Psi=\Psi_{c}$ as defined in (3.3) below.

$-1$$0$$-40$$-20$$0$$20$$40$$\operatorname{Re}(\sigma_{m,j})$$\operatorname{Im}(\sigma_{m,j})$$m=-2m_{c}$$m=-m_{c}$$m=0$$m=m_{c}$$m=2m_{c}$
Figure 1. The first 240 eigenvalues at the critical parameter when $m_{c}=2$,
$a=10$, $k=1$ and $N_{y}=240$.

We assume (as confirmed numerically for the parameter regimes considered
below) that the eigenvalues are ordered so that for each integer $m$,
$\sigma_{m,1}$ has the largest real part among the $\sigma_{m,j}$, for any
nonnegative integer $j$.

For each (nonnegative) wavenumber $m$, we define $\Psi_{m}$, when it exists,
to be the value of $\Psi$ for which the eigenvalue $\sigma_{m,1}$ crosses the
imaginary axis, that is

(3.2) $Re(\sigma_{m,1})=Re(\sigma_{-m,1})=\begin{cases}<0,\text{ if
}\Psi<\Psi_{m}\\\ =0,\text{ if }\Psi=\Psi_{m}\\\ >0,\text{ if
}\Psi>\Psi_{m}\end{cases}.$

Hence $\Psi=\Psi_{m}$ defines a neutral stability curve in the
$(a,\Psi)$-plane. In Figure 2, these neutral curves are plotted for zonal wave
numbers $m=1,2,3,4$.

$0$$5$$10$$15$$20$$25$$5\cdot
10^{-2}$$0.1$$0.15$$0.2$$\Psi_{1}$$\Psi_{2}$$\Psi_{3}$$\Psi_{4}$$a$$\Psi$$k=1$

$0$$5$$10$$15$$20$$25$$3$$4$$5$$6$$7$$\cdot
10^{-2}$$\Psi_{1}$$\Psi_{2}$$\Psi_{3}$$\Psi_{4}$$a$$\Psi$$k=2$

$0$$5$$10$$15$$20$$25$$2$$2.5$$3$$3.5$$4$$\cdot
10^{-2}$$\Psi_{1}$$\Psi_{2}$$\Psi_{3}$$\Psi_{4}$$a$$\Psi$$k=3$

Figure 2. Neutral stability curves $\Psi_{m},m=1,2,3,4$ defined in (3.2), for
wind-stress profiles defined by $k=1,2,3$. For values of $\Psi>\Psi_{m}$, the
shear flow $(\psi_{1}^{s},\psi_{2}^{s})$ in (2.3) becomes unstable to a
perturbation pattern with wavenumber $\alpha_{m}$.

Our numerical analysis suggests that for any $m$, $\Psi_{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

$\lim_{a\to a_{m}-}\Psi_{m}=\infty.$

Moreover, we numerically observe that the $a_{m}$ are ordered such that

$\infty>a_{1}>a_{2}>\cdots$

We define then the critical maximal shear flow’s velocity, $\Psi_{c}$, and the
critical zonal wavenumber, $m_{c}$, by

(3.3) $\Psi_{c}=\min_{m\in\mathbb{N}}\Psi_{m},$

and,

(3.4) $m_{c}=\underset{m\in\mathbb{N}}{\mathrm{argmin}}\;\Psi_{m},$

respectively.

A typical structure of the spectrum at the critical parameter value,
$\Psi=\Psi_{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 $\Psi=\Psi_{c}$, we
define the spatio-temporal function

(3.5)
$f_{m}(x,y,t)=2\operatorname{Re}\left(e^{i\operatorname{Im}(\sigma_{m,1})t}\psi_{m,1}(x,y)\right),\quad
m\in\mathbb{Z},$

where $\sigma_{m,1}$ is the first eigenvalue and $\psi_{m,1}$ is its
associated eigenfunction. The spatial structure of the eigenmodes
$\psi_{m_{c},1}$ is shown in Figure 3, revealing the well-known “banana-
shaped” patterns characteristics of baroclinic instability.

Figure 3. Real part of the upper and lower layers of the time-periodic
solution $f_{m_{c}}$ (or equivalently of the dominant eigenmode
$\psi_{m_{c},1}$) at $t=0$ for $k=1$ and $a=10$ and $a=5$, respectively.

The values of $\Psi_{c}$ with respect to the aspect ratio $a$ for $k=1,2,3$ is
shown in Figure 4. By the previous remarks,

(3.6) $\lim_{a\to a_{1}-}\Psi_{c}=\lim_{a\to a_{1}-}\Psi_{1}=\infty.$

By (3.6), _for $a>a_{1}$, the system is linearly stable._ As is expected, the
neutral stability curves (Figure 2) approach the asymptote $\Psi_{c}\to\infty$
as $a$ converges to the critical aspect ratio $a_{1}$ over which the system is
linearly stable for all $\Psi$. Our numerical results in Figure 4 show that
for aspect ratios $a$ in the range $3\leq a\leq 20$, the critical maximal
shear velocity $\Psi_{c}\approx 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,

$r_{1}\approx\frac{\tau\beta}{k^{2}\pi^{2}\Psi_{c}}\approx\frac{{10^{3}}}{k},$

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.

$4$$6$$8$$10$$12$$14$$16$$18$$20$$1\cdot 10^{-2}$$5\cdot
10^{-2}$$0.1$$k=1$$k=2$$k=3$$a$$\Psi_{c}$ Figure 4. The critical maximal shear
flow’s velocity $\Psi_{c}$ dependence on the channel’s aspect ratio $a$ for
$k=1,2$ and $k=3$.

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 $\Psi_{c}$
is exceeded.

For $a<a_{1}$, the system has a first transition at $\Psi=\Psi_{c}$ and
exactly one of the following two principal of exchange of stability (PES)
condition holds:

(3.7) $\displaystyle Re(\sigma_{m,1})=Re(\sigma_{-m,1})\begin{cases}<0,\text{
if }\Psi<\Psi_{c}\\\ =0,\text{ if }\Psi=\Psi_{c}\\\ >0,\text{ if
}\Psi>\Psi_{c}\end{cases}$ $\displaystyle\text{ if }m=m_{c}$ $\displaystyle
Re(\sigma_{m,1})=Re(\sigma_{-m,1})<0$ $\displaystyle\text{ if }m\neq m_{c}$

(3.8) $\displaystyle Re(\sigma_{m,1})=Re(\sigma_{-m,1})\begin{cases}<0,\text{
if }\Psi<\Psi_{c}\\\ =0,\text{ if }\Psi=\Psi_{c}\\\ >0,\text{ if
}\Psi>\Psi_{c}\end{cases}$ $\displaystyle\text{ if }m\in\\{m_{c},m_{c}+1\\}$
$\displaystyle Re(\sigma_{m,1})=Re(\sigma_{-m,1})<0$ $\displaystyle\text{ if
}m\notin\\{m_{c},m_{c}+1\\}.$

According to (3.7) and (3.8) either, two or four eigenvalues become unstable
as $\Psi$ crosses $\Psi_{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 $\Psi_{c}=\Psi_{m_{c}}=\Psi_{m_{c}+1}$. The
values of $a_{DH}$ where double Hopf transition occurs is given in Table 1.

$k$ | $m_{c}$ | $a_{DH}$ | $\Psi_{c}$ | $Re(A)$ | $Re(B)$ | $Re(C)$ | $Re(D)$ | $\delta$ | $\theta$ | $\delta\theta$
---|---|---|---|---|---|---|---|---|---|---
1 | 1 | 11.26 | 0.10 | -0.25 | -1.00 | 0.01 | -0.33 | -0.05 | 3.08 | -0.15
1 | 2 | 7.30 | 0.10 | -0.37 | -1.00 | -0.13 | -0.40 | 0.36 | 2.51 | 0.89
1 | 3 | 5.38 | 0.10 | -0.95 | -1.00 | -0.22 | -0.61 | 0.23 | 1.65 | 0.37
1 | 4 | 4.23 | 0.09 | -0.81 | -1.00 | -0.29 | -0.61 | 0.36 | 1.63 | 0.58
2 | 1 | 11.18 | 0.05 | -0.15 | -1.00 | 0.12 | -0.23 | -0.77 | 4.39 | -3.38
2 | 2 | 7.20 | 0.05 | -0.34 | -1.00 | -0.04 | -0.46 | 0.13 | 2.20 | 0.28
2 | 3 | 5.27 | 0.05 | -0.37 | -1.00 | -0.14 | -0.45 | 0.38 | 2.24 | 0.86
2 | 4 | 4.14 | 0.05 | -0.38 | -1.00 | -0.21 | -0.44 | 0.55 | 2.25 | 1.25
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 $\delta$ and $\theta$ are defined by (3.14).

Although the double Hopf transition is not generic, its analysis gives an
insight into regimes of moderate values of $\Psi$ 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}$.

$-1$$-0.5$$0$$\cdot
10^{-2}$$-20$$0$$20$$\operatorname{Re}(\sigma_{m,j})$$\operatorname{Im}(\sigma_{m,j})$$m=-2$$m=-1$$m=0$$m=1$$m=2$
Figure 5. The spectrum near the double Hopf aspect ratio $a_{DH}=11.263$ for
$k=1$. The first four critical eigenvalues $\sigma_{m,1}$, $m$ in
$\left\\{-2,-1,1,2\right\\}$, can be seen on the imaginary axis.

We recall that the neutral stability curves (shown in Figure 2) are found by
identifying the set of parameter values (here in the ($a$,$\Psi$)-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 $\Psi$
crosses the critical value $\Psi_{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\geq 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.

###  3.2. Hopf Bifurcation

We first investigate the generic Hopf transition scenario based on the
attractor bifurcation theorem [20, Theorem 5.2] and the dynamical transition
theorem [21, Theorem 2.1.3]. For proofs of the following lemma and theorem,
see Appendix A).

######  Lemma 3.1.

Assume that the PES condition (3.7) holds. Then the transition and stability
of the steady state solution (2.3) to Eq. (2.5), in the vicinity of the
critical maximal shear flow’s velocity $\Psi=\Psi_{c}$ and for any
sufficiently small initial data, are equivalent to the stability of the zero
solution of the following reduced equation

(3.9) $\frac{dz}{dt}=\sigma_{m_{c},1}z+Pz\lvert z\rvert^{2}+o(\lvert
z\rvert^{3}),$

where $z(t)$ denotes the complex amplitude aiming at approximating the
projection of the model’s solution onto the critical mode $\psi_{m_{c},1}$.
The coefficient $P$ is defined in Eq. (A.8) of Appendix A and is called the
transition number.

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. (1)

If $\operatorname{Re}(P)<0$, then Eq. (2.5) in $H_{1}$ (see (2.9)) undergoes a
continuous transition accompanied by a supercritical Hopf bifurcation on
$\Psi>\Psi_{c}$, with $\Psi_{c}$ defined in (3.3). In particular, the steady-
state solution bifurcates to a stable periodic solution $\psi$ on
$\Psi>\Psi_{c}$, satisfying $\psi\to\mathbf{0}$ as $\Psi\to\Psi_{c}$ and has
the following approximation

(3.10)
$\psi(x,y,t)=\left(\frac{-\operatorname{Re}(\sigma_{m_{c},1})}{\operatorname{Re}(P)}\right)^{\frac{1}{2}}f_{m_{c}}(x,y,t)+o\left(\lvert\Psi-\Psi_{c}\rvert^{\frac{1}{2}}\right).$

The spatial structure of the time periodic solution $\psi$ is shown at $t=0$
for different aspect ratios $a$ in Figure 3.

  2. (2)

If $\operatorname{Re}(P)>0$, then Eq. (2.5) in $H_{1}$ (see (2.9)) undergoes a
jump transition on $\Psi<\Psi_{c}$ accompanied by a subcritical Hopf
bifurcation. In particular, an unstable periodic solution $\psi$ given by
(3.10) bifurcates on $\Psi<\Psi_{c}$ and there is no periodic solution
bifurcating from $\mathbf{0}$ on $\Psi>\Psi_{c}$. Moreover, there is a
singularity separation at some $\Psi_{s}<\Psi_{c}$ generating an attractor and
an unstable periodic solution $\psi$.

$4$$6$$8$$10$$12$$14$$16$$18$$20$$-1$$0$$a$$\operatorname{Re}P$$k=1$

$4$$6$$8$$10$$12$$14$$16$$18$$20$$-1$$0$$a$$\operatorname{Re}P$$k=2$

$4$$6$$8$$10$$12$$14$$16$$18$$20$$-1$$0$$a$$\operatorname{Re}P$$k=3$

$m_{c}=1$$m_{c}=2$$m_{c}=3$$m_{c}=4$

Figure 6. The real part of the transition number $P$ compared to the channel
aspect ratio $a$ normalized by the largest absolute value of
$\operatorname{Re}P$. The parameters are as set in Table 2.

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\leq a\leq 20$, the critical wavenumber $m_{c}$ is
found to range from $1$ to $4$.

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.

$5$$10$$15$$20$$10^{0}$$10^{2}$$10^{4}$$a$$\lvert P0/P2\rvert$$k=1$

$5$$10$$15$$20$$10^{0}$$10^{2}$$10^{4}$$a$$k=2$

$5$$10$$15$$20$$10^{0}$$10^{2}$$10^{4}$$a$$k=3$

Figure 7. The effect of $P_{0}$ (of $m=0$ modes) compared to the effect
$P_{2}$ (of $2m_{c}$ modes) on the transition number.

In our case,

$P_{0}=\sum_{j=0}^{\infty}P_{0,j}$

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.

$-2$$0$$2$$4$$6$$8$$10$$12$$14$$-1$$-0.8$$-0.6$$-0.4$$j=2$$j=4$$j=6$$j=8$$j=10$$a$$k=1$

$-2$$0$$2$$4$$6$$8$$10$$12$$14$$-1$$-0.95$$-0.9$$-0.85$$-0.8$$j=2$$j=4$$j=6$$j=8$$j=10$$a$$k=2$

Figure 8. $P_{0,j}/\lvert P_{0,2}\rvert$ for even values of $j$. $P_{0,j}=0$
for odd values of $j$.

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.

$100$$200$$300$$400$$500$$600$$700$$200$$250$$300$$350$$L_{\text{dim}}$(in
kilometers)$T_{\text{dim}}$(in days)$k=1$$k=2$$k=3$ Figure 9. The
(dimensional) time period
$T_{\text{dim}}=\frac{2\pi}{\operatorname{Im}(\sigma_{m_{c},1})}\frac{L_{y}}{U}$
compared to the (dimensional) length of the channel $L_{\text{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
$\operatorname{Im}\sigma_{m_{c},1}$ of the critical eigenvalue at the double
Hopf aspect ratios $a_{DH}$.

###  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).

######  Lemma 3.2.

Assume that the PES condition (3.8) holds. Then, the transition and stability
of the steady state solution (2.3) of Eq. (2.5) in $H_{1}$ (see (2.9)), in the
vicinity of the critical zonal shear flow’s velocity, $\Psi=\Psi_{c}$, and for
any sufficiently small initial data are equivalent to the stability of the
zero solution of the following system

(3.11) $\displaystyle\frac{dz_{1}}{dt}$
$\displaystyle=\sigma_{m_{c},1}z_{1}+z_{1}(A\lvert z_{1}\rvert^{2}+B\lvert
z_{2}\rvert^{2})+o(\lvert z\rvert)^{3}),$ $\displaystyle\frac{dz_{2}}{dt}$
$\displaystyle=\sigma_{m_{c}+1,1}z_{2}+z_{2}(C\lvert z_{1}\rvert^{2}+D\lvert
z_{2}\rvert^{2})+o(\lvert z\rvert)^{3}),$

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
$\psi_{m_{c},1}$ and $\psi_{m_{c}+1,1}$, 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

(3.12)
$\operatorname{Re}(A)<0,\,\operatorname{Re}(B)<0,\operatorname{Re}(B)+\operatorname{Re}(C)<0,\operatorname{Re}(D)<0,$

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)
$\lambda_{1}=Re(\sigma_{m_{c},1}),\qquad\lambda_{2}=Re(\sigma_{m_{c}+1,1}).$
(3.14) $\delta=\frac{Re(C)}{Re(A)},\qquad\theta=\frac{Re(B)}{Re(D)}$
$\eta_{1}=\left(\frac{\lambda_{1}-\theta\lambda_{2}}{(\theta\delta-1)\operatorname{Re}(A)}\right)^{\frac{1}{2}},\qquad\eta_{2}=\left(\frac{\lambda_{2}-\delta\lambda_{1}}{(\theta\delta-1)\operatorname{Re}(D)}\right)^{\frac{1}{2}}.$

Recalling $f_{m_{c}}$ defined in (3.5) (with $m=m_{c}$), we define the
following spatio-temporal profiles

(3.15)
$\displaystyle\psi_{p}^{m_{c}}(x,y,t)=\left(-\frac{\lambda_{1}}{\operatorname{Re}(A)}\right)^{\frac{1}{2}}f_{m_{c}}(x,y,t)+o\left(\lvert\Psi-\Psi_{c}\rvert^{\frac{1}{2}}\right),$
$\displaystyle\psi_{p}^{m_{c}+1}(x,y,t)=\left(\frac{\lambda_{2}}{\operatorname{Re}(D)}\right)^{\frac{1}{2}}f_{m_{c}+1}(x,y,t)+o\left(\lvert\Psi-\Psi_{c}\rvert^{\frac{1}{2}}\right),$
$\displaystyle\psi_{qp}(x,y,t)=\eta_{1}f_{m_{c}}(x,y,t)+\eta_{2}f_{m_{c}+1}(x,y,t)+o\left(\lvert\Psi-\Psi_{c}\rvert^{\frac{1}{2}}\right)$

$T_{1}$$T_{2}$$\lambda_{2}$$\lambda_{1}$IIIIII IVVVI (a) $\theta>0$,
$\delta>0$, $\theta\delta>1$

$T_{2}$$T_{1}$$\lambda_{2}$$\lambda_{1}$IVIIIII IVVVI (b) $\theta>0$,
$\delta>0$, $\theta\delta<1$

$T_{2}$$T_{1}$$\lambda_{2}$$\lambda_{1}$VIIVIIIVIIX IVV (c) $\theta>0$,
$\delta<0$

Figure 10.  The regions in the $\lambda_{1}$–$\lambda_{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/\theta$ and $\delta$
(as defined in (3.14)) respectively.

$\psi_{p}^{m_{c}}$$\psi_{p}^{m_{c}+1}$$\lvert z_{1}\rvert$$\lvert
z_{2}\rvert$O (a) Region I

$\psi_{p}^{m_{c}}$$\psi_{p}^{m_{c}+1}$$\psi_{qp}$$\lvert z_{1}\rvert$$\lvert
z_{2}\rvert$$O$ (b) Region II

$\psi_{p}^{m_{c}}$$\psi_{p}^{m_{c}+1}$$\lvert z_{1}\rvert$$\lvert
z_{2}\rvert$O (c) Region III

$\psi_{p}^{m_{c}}$$\psi_{p}^{m_{c}+1}$$\psi_{qp}$$\lvert z_{1}\rvert$$\lvert
z_{2}\rvert$$O$ (d) Region VII

$\psi_{p}^{m_{c}+1}$$\psi_{qp}$$\lvert z_{1}\rvert$$\lvert z_{2}\rvert$$O$ (e)
Region VIII

$\psi_{p}^{m_{c}}$$\psi_{qp}$$\lvert z_{1}\rvert$$\lvert z_{2}\rvert$$O$ (f)
Region IX

Figure 11. The dynamics in the regions given in the first quadrant of Figure
10. $\psi_{p}^{m_{c}}$, $\psi_{p}^{m_{c}+1}$ are time-periodic with zonal
wavenumbers $m_{c}$ and $m_{c}+1$ respectively. $\psi_{qp}$ is the
quasiperiodic solution given in (3.15).

######  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 $\Psi=\Psi_{c}$,
and an $S^{3}$ attractor $\Sigma$ bifurcates on $\Psi>\Psi_{c}$, which
converges to $\mathbf{0}$ as $\Psi$ approaches to $\Psi_{c}$ from right.
Depending on the values of $\theta$ and $\delta$, there are three transition
scenarios as shown in Figure 10. In each scenario, near the onset of
transition $(\lambda_{1},\lambda_{2})=(0,0)$, the $\lambda_{1}-\lambda_{2}$
plane is dissected into several regions with distinct topological structures
for the attractor $\Sigma$ as given in Figure 11.

######  Remark 3.1.

  1. (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. (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, $\Psi>\Psi_{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].

Figure 12. Real part of the upper and lower layers of the time quasi-periodic
solution $\psi_{qp}$ given in Theorem 3.2 at $t=0$ for $k=1$ and $a=5.722$
where the first two critical wavenumbers are $m_{c}=1$ and $m_{c}=2$.

##  4\. 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 $\Psi$, or
equivalently the maximal amplitude $\tau$ of the wind-stress curl (see (2.4)),
is below a critical threshold $\Psi_{c}$. Moreover, as this critical threshold
is exceeded, generically a Hopf bifurcation occurs; see Theorem 3.1. To derive
this result, we approached the problem by means of dynamic transition theory
determining all the attractors near the onset of instability. These attractors
then describe how the onset of linear instability translates into the
emergence of nonlinear patterns for the actual flow; see (3.10).

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 quasi-
periodic 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
quasi-periodic 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

$\psi^{*}_{m,j}=e^{i\alpha_{m}x}Y_{m}^{*}(y).$

We denote the critical eigenmode and the critical eigenvalue by

$\psi_{c}=\psi_{m_{c},1},\qquad\sigma_{c}=\sigma_{m_{c},1}.$

We denote the bilinear operator $\mathcal{G}$ as

$\mathcal{G}(u)=G_{2}(u,u)$

where $G_{2}(u,v)$ is linear in each component. Let us define now

(A.1) $G_{s}(u,v)=G_{2}(u,v)+G_{2}(v,u).$

The center part of the solution is

(A.2) $u_{c}=z(t)\psi_{c}+c.c.$

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 $\psi_{c}$.

(A.3)
$\dot{z}=\sigma_{c}z+\frac{1}{\langle\mathcal{M}\psi_{c},\psi_{c}^{*}\rangle}\langle\mathcal{G}(u_{c}+\Phi),\psi_{c}^{*}\rangle.$

where $\Phi$ is the center manifold function. We will obtain its quadratic
approximation $\Phi_{2}$ given by

$\Phi=\Phi_{2}(z,\overline{z})+o(2)\\\ $

Here

$o(n)=o\left(\lvert(z,\overline{z})\rvert^{n}\right)$

denotes higher than $n$-th order terms in $z$, $\overline{z}$.

Using the notation (A.1), the reduced equation (A.3) can be written

(A.4)
$\dot{z}=\sigma_{c}z+\frac{1}{\langle\mathcal{M}\psi_{c},\psi_{c}^{*}\rangle}\langle
G_{s}(u_{c},\Phi_{2}),\psi_{c}^{*}\rangle+o(3).$

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],

(A.5)
$\displaystyle\Phi_{2}=(2\sigma_{c}-\mathcal{L})^{-1}\Pi_{s}G_{2}(z\psi_{c},z\psi_{c})+(\sigma_{c}+\overline{\sigma}_{c}-\mathcal{L})^{-1}\Pi_{s}G_{2}(z\psi_{c},\overline{z\psi_{c}})+c.c.$

where $\mathcal{L}=\Pi_{s}\mathcal{M}^{-1}\mathcal{N}$ and $\Pi_{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) $\Phi_{2}=z^{2}\sum_{j\geq 1}g_{2m_{c},j}\psi_{2m_{c},j}+\lvert
z\rvert^{2}\sum_{j\geq 1}g_{0,j}\psi_{0,j}+c.c.$

Here

(A.7) $\displaystyle
g_{0,j}=\frac{1}{(\sigma_{c}+\overline{\sigma}_{c}-\sigma_{0,j})\langle\mathcal{M}\psi_{0,j},\psi_{0,j}^{*}\rangle}\langle
G_{2}(\psi_{c},\overline{\psi_{c}}),\psi_{0,j}^{*}\rangle$ $\displaystyle
g_{2m_{c},j}=\frac{1}{(2\sigma_{c}-\sigma_{2m_{c},j})\langle\mathcal{M}\psi_{2m_{c},j},\psi_{2m_{c},j}^{*}\rangle}\langle
G_{2}(\psi_{c},\psi_{c}),\psi_{2m_{c},j}^{*}\rangle,$

are the coefficients of the center manifold function.

We write (A.4) as (3.9), that is

$\dot{z}=\sigma_{c}z+Pz\lvert z\rvert^{2}+o(3).$

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

(A.8) $P=P_{0}+P_{2},$

where

(A.9) $P_{0}=\sum_{j\geq 1}P_{0,j},\qquad
P_{0,j}=\frac{1}{\langle\mathcal{M}\psi_{c},\psi_{c}^{*}\rangle}(g_{0,j}+c.c.)\langle
G_{s}(\psi_{c},\psi_{0,j}),\psi_{c}^{*}\rangle,$

denotes the contribution of the zero-wavenumber (stable) modes $\psi_{0,j}$
while

(A.10) $P_{2}=\sum_{j\geq 1}P_{2,j},\qquad
P_{2,j}=\frac{1}{\langle\mathcal{M}\psi_{c},\psi_{c}^{*}\rangle}g_{2m_{c},j}\langle
G_{s}(\overline{\psi_{c}},\psi_{2m_{c},j}),\psi_{c}^{*}\rangle,$

denotes the contribution of the modes $\psi_{2m_{c},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

$u_{c}=z_{1}(t)\psi_{1}+z_{2}(t)\psi_{2}+c.c.$

where the first two critical modes are

$\psi_{1}=\psi_{m_{c},1},\quad\psi_{2}=\psi_{m_{c}+1,1},\quad\psi_{-1}=\psi_{-m_{c},1},\quad\psi_{-2}=\psi_{-m_{c}-1,1}.$

with corresponding eigenvalues

$\sigma_{1}=\sigma_{m_{c},1},\quad\sigma_{2}=\sigma_{m_{c}+1,1},\quad\sigma_{-1}=\sigma_{-m_{c},1},\quad\sigma_{-2}=\sigma_{-m_{c}-1,1}$

The equation (A.4) becomes the system

(B.1)
$\dot{z}_{j}=\sigma_{j}z_{j}+\frac{1}{\langle\mathcal{M}\psi_{j},\psi_{j}^{*}\rangle}\langle\mathcal{G}(u_{c}+\Phi_{2}),\psi_{j}^{*}\rangle+o(3),\qquad
j=1,2.$

and the center manifold function (A.5) is replaced by

$\Phi_{2}=\sum_{\lvert j\rvert,\lvert
k\rvert=1}^{2}z_{j}z_{k}\Psi_{j,k}+o(2),\qquad\Psi_{j,k}=(\sigma_{j}+\sigma_{k}-\mathcal{L})^{-1}\Pi_{s}G_{2}(\psi_{j},\psi_{k}),$

where $\Pi_{s}$ denotes the projector onto the stable subspace.

Now, the equations (B.1) become (3.11) with the coefficients defined below:

(B.2) $\displaystyle A=g_{1,1,-1,1}+g_{1,-1,1,1}+g_{-1,1,1,1}$ $\displaystyle
B=g_{1,2,-2,1}+g_{1,-2,2,1}+g_{2,1,-2,1}+g_{2,-2,1,1}+g_{-2,1,2,1}+g_{-2,2,1,1}$
$\displaystyle
C=g_{2,1,-1,2}+g_{2,-1,1,2}+g_{1,2,-1,2}+g_{1,-1,2,2}+g_{-1,2,1,2}+g_{-1,1,2,2}$
$\displaystyle D=g_{2,2,-2,2}+g_{2,-2,2,2}+g_{-2,2,2,2}$ $\displaystyle
g_{i,j,k,l}=\frac{1}{\langle\mathcal{M}\psi_{l},\psi_{l}^{*}\rangle}\langle
G_{s}(\psi_{i},\Psi_{j,k},\psi_{l}^{*})\rangle,$
$\displaystyle\Psi_{j,k}=(\sigma_{j}+\sigma_{k}-\mathcal{L})^{-1}\Pi_{s}G_{2}(\psi_{j},\psi_{k})$

We note that the above coefficients contain only $g_{i,j,k,l}$ for which
$i+j+k=l$. The expansion of the center manifold coefficients can be written
more explicitly as (see also [5, Theorem 2])

$\Psi_{j,k}=\sum_{i=1}^{\infty}\frac{\langle
G_{2}(\psi_{m_{j}},\psi_{m_{k}}),\psi^{*}_{m_{j}+m_{k},i}\rangle}{\langle\mathcal{M}\psi_{m_{j}+m_{k},i},\psi^{*}_{m_{j}+m_{k},i}\rangle}(\sigma_{j}+\sigma_{k}-\sigma_{m_{j}+m_{k}})^{-1}\psi_{m_{j}+m_{k},i}$

Now we analyze the equations (3.11) by first putting them in polar form

$z_{j}=\rho_{i}e^{i\gamma_{j}},\quad j=1,2$

which yields

(B.3)
$\displaystyle\dot{\rho_{1}}=Re(\sigma_{1})\rho_{1}+\rho_{1}(Re(A)\rho_{1}^{2}+Re(B)\rho_{2}^{2})+\text{
h.o.t. }$
$\displaystyle\dot{\rho_{2}}=Re(\sigma_{2})\rho_{2}+\rho_{2}(Re(C)\rho_{1}^{2}+Re(D)\rho_{2}^{2})+\text{
h.o.t. }$

and

$\displaystyle\dot{\gamma_{1}}=Im(\sigma_{1})+\text{ h.o.t. }$
$\displaystyle\dot{\gamma_{2}}=Im(\sigma_{2})+\text{ h.o.t. }$

For the specific case of (3.12), the equations (B.3) always admit the
solutions which represent the periodic solutions

$\displaystyle(\rho_{1},\rho_{2})=\left(-\frac{Re(\sigma_{1})}{Re(A)},0\right),$
$\displaystyle(\rho_{1},\rho_{2})=\left(0,-\frac{Re(\sigma_{2})}{Re(D)}\right),$

with respective eigenvalues

$\displaystyle\kappa_{1}=-2\sigma_{1},\qquad\kappa_{2}=\sigma_{2}-\delta\sigma_{1}$
$\displaystyle\kappa_{1}=-2\sigma_{2},\qquad\kappa_{2}=\sigma_{1}-\theta\sigma_{2}$

Also, the equations (B.3) admit the following solution which represents a
quasi-periodic solution

$(\rho_{1},\rho_{2})=\left(\frac{\sigma_{1}-\theta\sigma_{2}}{Re(A)(\theta\delta-1)},\frac{\sigma_{2}-\delta\sigma_{1}}{Re(D)(\theta\delta-1)}\right).$

Since the Jacobian matrix of the right hand side of (B.3) at the quasi-
periodic solution has determinant

$\frac{-4(\sigma_{2}-\delta\sigma_{1})(\sigma_{1}-\theta\sigma_{2})}{\theta\delta-1}.$

With this information, the transition scenarios summarized in Figure 10 and
Figure 11 can be obtained by a standard analysis. To prove the claim on the
bifurcation of an $S^{3}$-attractor, we need to prove that
$(\rho_{1},\rho_{2})=(0,0)$ is locally stable equilibrium of (B.3) at
$\Psi=\Psi_{c}$, that is when
$\operatorname{Re}(\sigma_{1})=\operatorname{Re}(\sigma_{2})=0$. In this case
by assumption (3.12), from (B.3), we can obtain

$\frac{d}{dt}(\rho_{1}^{2}+\rho_{2}^{2})=Re(A)\rho_{1}^{4}+(\operatorname{Re}(B)+\operatorname{Re}(C))\rho_{1}^{2}\rho_{2}^{2}+\operatorname{Re}(D)\rho_{2}^{4}<0$

which proves the claim.

##  Appendix C Numerical treatment of the linear stability problem

To solve the eigenvalue problem numerically, we first plugin the ansatz

(C.1) $\psi(x,y)=e^{i\alpha_{m}x}Y_{j}(y),\qquad
j\in\mathbb{N},m\in\mathbb{Z},\quad\alpha_{m}:=am\pi.$

into the eigenvalue problem

$\sigma\mathcal{M}\psi(x,y)=\mathcal{N}(y)\psi(x,y),$

to obtain

(C.2) $\sigma\widetilde{\mathcal{M}}Y(y)=\widetilde{\mathcal{N}}(y)Y(y),\qquad
Y(y)=(Y_{1}(y),Y_{2}(y))$

Here the linear operators $\widetilde{\mathcal{M}}$ and
$\widetilde{\mathcal{N}}(y)$ are defined as

(C.3) $\widetilde{\mathcal{M}}=\begin{bmatrix}\Delta_{m}-F_{1}&F_{1}\\\
F_{2}&\Delta_{m}-F_{2}\end{bmatrix},\qquad\widetilde{\mathcal{N}}(y)=\begin{bmatrix}N_{11}&N_{12}\\\
N_{21}&N_{22}\end{bmatrix}$

where

(C.4) $\displaystyle N_{11}=c_{1}\cos k\pi
y\left((k\pi)^{2}+\Delta_{m}\right)-\beta i\alpha_{m}-r_{1}\Delta_{m}$
$\displaystyle N_{12}=c_{1}F_{1}\cos k\pi y$ $\displaystyle N_{21}=0$
$\displaystyle N_{22}=-c_{1}F_{2}\cos k\pi y-\beta
i\alpha_{m}-r_{2}\Delta_{m}$ $\Delta_{m}=D^{2}-\alpha_{m}^{2},\qquad
D=\frac{\partial}{\partial y}$ $c_{1}=\Psi k\pi i\alpha_{m}.$

The eigenvalue problem (C.2) is supplemented with the following boundary
conditions

(C.5) $Y_{i}(\pm 1)=D^{2}Y_{i}(\pm 1)=0,\qquad i=1,2.$

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 $\\{L_{j}\\}$ be the Legendre polynomials and consider compact
combinations of the Legendre polynomials

$f_{j}(y)=L_{j}(y)+\sum_{k=1}^{4}c_{jk}L_{j+k}(y)$

with $c_{jk}$ chosen so that $f_{j}$ satisfy the boundary conditions (C.5),
i.e.

$f_{j}(\pm 1)=D^{2}f_{j}(\pm 1)=0.$

To discretize the eigenvalue problem, we plug

(C.6)
$Y_{i}^{N_{y}}(y)=\sum_{j=0}^{N_{y}-1}y_{j}^{(i)}f_{j}(y),\qquad\widehat{Y}_{i}=[y_{0}^{(i)},\dots,y_{N_{y}-1}^{(i)}]^{T},\qquad
i=1,2.$

into (C.2) to obtain

(C.7) $\sigma\begin{bmatrix}\widehat{\Delta}_{m}-F_{1}A_{3}&F_{1}A_{3}\\\
F_{2}A_{3}&\widehat{\Delta}_{m}-F_{2}A_{3}\end{bmatrix}\begin{bmatrix}\widehat{Y}_{1}\\\
\widehat{Y}_{2}\end{bmatrix}=\begin{bmatrix}\widehat{N}_{11}&\widehat{N}_{12}\\\
\widehat{N}_{21}&\widehat{N}_{22}\end{bmatrix}\begin{bmatrix}\widehat{Y}_{1}\\\
\widehat{Y}_{2}\end{bmatrix}$

(C.8)
$\displaystyle\widehat{N}_{11}=c_{1}(k\pi)^{2}A_{5}+c_{1}\left(A_{4}^{T}-\alpha_{m}^{2}A_{5}\right)-\beta
i\alpha_{m}A_{3}-r_{1}\widehat{\Delta}_{m}$
$\displaystyle\widehat{N}_{12}=c_{1}F_{1}A_{5}$
$\displaystyle\widehat{N}_{21}=0$
$\displaystyle\widehat{N}_{22}=-c_{1}F_{2}A_{5}-\beta
i\alpha_{m}A_{3}-r_{2}\widehat{\Delta}_{m}$

Here

(C.9) $\displaystyle A_{1}=(D^{4}f_{j},f_{k}),\qquad
A_{2}=(D^{2}f_{j},f_{k}),\qquad A_{3}=(f_{j},f_{k}),$ $\displaystyle
A_{4}=(\cos k\pi yD^{2}f_{j},f_{k}),\qquad A_{5}=(\cos k\pi yf_{j},f_{k}),$
$\displaystyle\widehat{\Delta}_{m}=A_{2}-\alpha_{m}^{2}A_{3},$

with $(f,g)=\int_{-1}^{1}f(y)g(y)\,\mathrm{d}y$. The explicit expression of
the matrices $A_{i}$, $i=1,\dots,5$ can be found in [10].

##  Appendix D  Practical aspects for the calculation of the transition number

The practical calculation of the $P_{0}$-term in (A.9) and the $P_{2}$-term in
(A.10), boils down to the efficient calculation of the inner and trilinear
products involved therein. In that respect, we provide here explicit
expressions of the latter. They are given by

$\displaystyle\langle\mathcal{M}\psi_{m,j},\psi_{m,j}^{*}\rangle=i\alpha_{m}$
$\displaystyle\int_{-1}^{1}\left((D^{2}-\alpha_{m}^{2}-F_{1})Y_{m,j}^{1}+F_{1}DY_{m,j}^{2}\right)\overline{Y_{m,j}^{*1}}\,\mathrm{d}y$
$\displaystyle+i\alpha_{m}$
$\displaystyle\int_{-1}^{1}\left((D^{2}-\alpha_{m}^{2}-F_{2})Y_{m,j}^{2}+F_{2}DY_{m,j}^{1}\right)\overline{Y_{m,j}^{*2}}\,\mathrm{d}y,$

and

$\langle
G_{2}(\psi_{m,j},\psi_{n,k}),\psi_{p,l}^{*}\rangle=-\delta_{m+n-p}\left(\int_{-1}^{1}G_{2}^{1}\overline{Y_{p,l}^{*,1}}+G_{2}^{2}\overline{Y_{p,l}^{*,2}}\right)\,\mathrm{d}y,$

$\displaystyle G_{2}^{1}=i\alpha_{m}$ $\displaystyle
Y_{m,j}^{1}\left(D(D^{2}-\alpha_{n}^{2}-F_{1})Y_{n,k}^{1}+F_{1}DY_{n,k}^{2}\right)$
$\displaystyle-i\alpha_{n}\left((D^{2}-\alpha_{n}^{2}-F_{1})Y_{n,k}^{1}+F_{1}Y_{n,k}^{2}\right)DY_{m,j}^{1}$
$\displaystyle G_{2}^{2}=i\alpha_{m}$ $\displaystyle
Y_{m,j}^{2}\left(D(D^{2}-\alpha_{n}^{2}-F_{2})Y_{n,k}^{2}+F_{2}DY_{n,k}^{1}\right)$
$\displaystyle-i\alpha_{n}\left((D^{2}-\alpha_{n}^{2}-F_{2})Y_{n,k}^{2}+F_{2}Y_{n,k}^{1}\right)DY_{m,j}^{2}.$

In practice, the integrals can be evaluated by any commonly used quadrature
rules in which the values of the integrand are evaluated at quadrature points.
In our calculations, we use

$\int_{-1}^{1}f(y)\,\mathrm{d}y=\sum_{n=0}^{N_{y}}f(y_{n})\omega_{n},$

where $y_{n}$ and $\omega_{n}$ are Legendre-Gauss-Lobatto quadrature points
and weights respectively.

##  Appendix E Model Parameters

$L_{y}$ | meridional length scale |  $10^{6}$ m
---|---|---
$H_{1}$ | upper layer depth |  $250$m
$H_{2}$ | upper layer depth |  $750$m
$\Delta\rho$ | density difference |  $\rho_{2}-\rho_{1}=1$ kgm-3
$U$ | characteristic velocity |  $\tau_{0}/(\rho_{0}H_{1}\beta_{0}L_{y})$ = 0.02 m s-1
$\tau_{0}$ | characteristic zonal wind stress | 0.1 Pa
$\beta_{0}$ | planetary vorticity gradient |  $2\times 10^{-11}$ (ms)-1
$f_{0}$ | reference Coriolis parameter |  $10^{-4}$ s-1
$\epsilon_{0}$ | bottom friction coefficient |  $10^{-7}$ s-1
$g$ | gravitational acceleration |  $9.8$ ms-2
$g^{\prime}$ | reduced gravity |  $g\Delta\rho/\rho_{0}=0.01$ ms-2
$\rho_{0}$ | reference density |  $10^{3}$ kgm-3
$F_{1}$ | upper layer Froude number | $f_{0}^{2}L_{y}^{2}/(g^{\prime}H_{1})=4000$
$F_{2}$ | lower layer Froude number | $f_{0}^{2}L_{y}^{2}/(g^{\prime}H_{2})=4000/3$
$\beta$ | planetary vorticity factor | $\beta_{0}L_{y}^{2}/U=1000$
$\tau$ | wind-stress parameter | $1.0$
$r_{2}$ | lower layer linear friction coefficient | $r_{2}=\epsilon_{0}L_{y}/U=5.0$
Table 2. Set of model parameters used in the numerical study of the problem.

### Acknowledgments:

The authors are grateful for two anonymous reviewers for their insightful
comments.

This work has been partially supported by the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation program (grant
agreement No. 810370). This study was also supported by a Ben May Center grant
for theoretical and/or computational research and by the Office of Naval
Research Multidisciplinary University Research Initiative Grant N00014-20-
1-2023.

## References

  * [1] Bernier, C.: Existence of attractor for the quasi-geostrophic approximation of the Navier-Stokes equations and estimate of its dimension.  Adv. Math. Sci. Appl. (4), 465–489 (1994)
  * [2] Brézis, H.: Functional Analysis, Sobolev Spaces and Partial Differential Equations.  Springer (2010)
  * [3] Cai, M., Hernandez, M., Ong, K., Wang, S.: Baroclinic instability and transitions in a two-layer quasigeostrophic channel model.  arXiv:1705.07989 [physics] (2017)
  * [4] Chekroun, M.D., Koren, I., Liu, H.: Efficient reduction for diagnosing Hopf bifurcation in delay differential systems: Applications to cloud-rain models.  Chaos (Woodbury, N.Y.) 30(5), 053130 (2020)
  * [5] Chekroun, M.D., Liu, H., McWilliams, J.: Variational approach to closure of nonlinear dynamical systems: Autonomous case.  Journal of Statistical Physics 179(5), 1073–1160 (2020)
  * [6] Chekroun, M.D., Liu, H., McWilliams, J.C.: The emergence of fast oscillations in a reduced primitive equation model and its implications for closure theories.  Comp. & Fluids 151, 3–22 (2017).
  * [7] Chekroun, M.D., Liu, H., McWilliams, J.C.: Stochastic rectification of fast oscillations on slow manifold closures.  Proc. Natl. Acad. Sci. USA 118(48), e2113650118 (2021)
  * [8] Chekroun, M.D., Liu, H., McWilliams, J.C., Wang, S.: Transitions in stochastic non-equilibrium systems: Efficient reduction and analysis.  arXiv preprint arXiv:2202.07031 (2022)
  * [9] Chemke, R., Kaspi, Y.: The effect of eddy–eddy interactions on jet formation and macroturbulent scales.  J. Atmos. Sci. 73(5), 2049–2059 (2016)
  * [10] Dijkstra, H., Sengul, T., Shen, J., Wang, S.: Dynamic transitions of quasi-geostrophic channel flow.  SIAM Journal on Applied Mathematics 75(5), 2361–2378 (2015).
  * [11] Dijkstra, H.A.: Nonlinear Physical Oceanography: A Dynamical Systems Approach to the Large Scale Ocean Circulation and El Niño, second edn.  Springer, New York, 532 pp. (2005)
  * [12] Dijkstra, H.A., Katsman, C.A.: Temporal variability of the wind-driven quasi-geostrophic double gyre ocean circulation: Basic bifurcation diagrams.  Geophys. Astrophys. Fluid Dyn. 85, 195–232 (1997)
  * [13] Dijkstra, H.A., Van der Vaart, P. C. F.: On the physics of upgradient momentum transport in unstable eastward jets.  Geophys. Astrophys. Fluid Dyn. 88, 295–323 (1998)
  * [14] Eady, E.T.: Long waves and cyclone waves.  Tellus 1, 33–52 (1949)
  * [15] Kieu, C., Sengul, T., Wang, Q., Yan, D.: On the Hopf (double Hopf) bifurcations and transitions of two-layer western boundary currents.  Communications in Nonlinear Science and Numerical Simulation 65, 196–215 (2018).
  * [16] Kuznetsov, Y.: Elements of Applied Bifurcation Theory, third edn.  Applied Mathematical Sciences. Springer-Verlag, New York (2004).
  * [17] Li, K.: Dynamic transitions of the Swift-Hohenberg equation with third-order dispersion.  Discrete & Continuous Dynamical Systems - B (2020).
  * [18] Lovegrove, A.F., Moroz, I.M., Read, P.L.: Bifurcations and instabilities in rotating, two-layer fluids: II. $\beta$-Plane.  Nonlinear Processes in Geophysics 9(3/4), 289–309 (2002).
  * [19] Lu, C., Mao, Y., Sengul, T., Wang, Q.: On the spectral instability and bifurcation of the 2d-quasi-geostrophic potential vorticity equation with a generalized Kolmogorov forcing.  Physica D: Nonlinear Phenomena 403, 132296 (2020).
  * [20] Ma, T., Wang, S.: Bifurcation Theory and Applications, _World Scientific Series on Nonlinear Science Series A_ , vol. 53.  World Scientific (2005).
  * [21] Ma, T., Wang, S.: Phase Transition Dynamics.  Springer International Publishing, Cham (2019).
  * [22] Meshalkin, L.D., Sinai, I.G.: Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid.  Journal of Applied Mathematics and Mechanics 25(6), 1700–1705 (1961).
  * [23] O’Gorman, P.A., Schneider, T.: Recovery of atmospheric flow statistics in a general circulation model without nonlinear eddy‐eddy interactions.  Geophysical Research Letters 34(22) (2007).
  * [24] Panetta, R., Held, I.: Baroclinic eddy fluxes in a one-dimensional model of quasi-geostrophic turbulence.  J. Atmos. Sci. 45(22), 3354–3365 (1988)
  * [25] Pedlosky, J.: Finite-amplitude baroclinic waves.  Journal of the Atmospheric Sciences 27(1), 15–30 (1970).
  * [26] Pedlosky, J.: Geophysical Fluid Dynamics, second edn.  Springer-Verlag, New York (1987)
  * [27] Phillips, N.A.: A simple three-dimensional model for the study of large-scale extratropical flow patterns.  J. Meteor. 8, 381–394 (1951)
  * [28] Romea, R.D.: The effects of friction and $\beta$ on finite amplitude baroclinic waves.  Journal of Atmospheric Science 34, 1689–1695 (1977)
  * [29] Sapsis, T.P., Dijkstra, H.A.: Interaction of additive noise and nonlinear dynamics in the double-gyre wind-driven ocean circulation.  Journal of physical oceanography 43(2), 366–381 (2013)
  * [30] Şengül, T., Wang, S.: Dynamic transitions and baroclinic instability for 3d continuously stratified Boussinesq flows.  Journal of Mathematical Fluid Mechanics 20(3), 1173–1193 (2018).
  * [31] Shen, J., Tang, T., Wang, L.L.: Spectral Methods, _Springer Series in Computational Mathematics_ , vol. 41.  Springer Berlin Heidelberg, Berlin, Heidelberg (2011).
  * [32] Simonnet, E., Ghil, M., Wang, S., Chen, Z.M.: Hopf bifurcation in quasi-geostrophic channel flow.  SIAM Journal on Applied Mathematics 64(1), 343–368 (2003)
  * [33] Srinivasan, K., Young, W.: Zonostrophic instability.  J. Atmos. Sci. 69(5), 1633–1656 (2012)
  * [34] Sura, P., Fraedrich, K., Lunkeit, F.: Regime transitions in a stochastically forced double-gyre model.  Journal of Physical Oceanography 31(2), 411–426 (2001)
  * [35] Talley, L.: Descriptive Physical Oceanography.  Academic Presss, New York (2011)
  * [36] Vallis, G.: Atmospheric and Oceanic Fluid Dynamics.  Cambridge University Press, Cambridge, UK (2006)
  * [37] Van der Vaart, P.C.F., Dijkstra, H.A.: Sideband instabilities of mixed barotropic/baroclinic waves growing on a midlatitude zonal jet.  Physics of Fluids 9, 615–631 (1997)
  * [38] Van der Vaart, P. C. F., Schuttelaars, H.M., Calvete, D., Dijkstra, H.A.: Instability of time-dependent wind-driven ocean gyres.  Physics of Fluids 14, 3601–3615 (2002)
  * [39] Welander, P.: A two-layer frictional model of wind-driven motion in a rectangular oceanic basin.  Tellus 18(1), 54–62 (1966).