## Abstract

A two-layer shallow-water model is used to investigate the transition of wind-driven double-gyre circulation from laminar flow to turbulence as the Reynolds number (Re) is systematically increased. Two distinctly different phases of turbulent double-gyre patterns and energy trajectories are exhibited before and after at Re = 95: deterministic and fully developed turbulent circulations. In the former phase, the inertial subgyres vary between an asymmetric solution and an antisymmetric solution and the double-gyre circulations reach the aperiodic solution mainly due to their barotropic instability. An integrated kinetic energy in the lower layer is slight and the generated mesoscale eddies are confined in the upper layer. The power spectrum of energies integrated over the whole domain at Re = 70 has peaks at the interannual periods (4–7 yr) and the interdecadal period (10–20 yr). The loops of the attractors take on one cycle at those periods and display the blue-sky catastrophe. At Re = 95, the double-gyre circulation reaches a metastable state and the attracters obtained from the three energies form a topological manifold. In the latter, as Re increases, the double-gyre varies from a metastable state to a chaotic state because of the barotropic instability of the eastward jet and the baroclinic instability of recirculation retrograde flow, and the eastward jet meanders significantly with interdecadal variability. The generated eddies cascade to the red side of the power spectrum as expected in the geostrophic turbulence. The main results in the simulation may indicate essential mechanisms for the appearance of multiple states of the Kuroshio and for low-frequency variations in the midlatitude ocean.

## 1. Introduction

Ocean prediction of the Kuroshio and the Gulf Stream will soon begin using an ocean general circulation model (OGCM) (Traon et al. 1999). Before this takes effect, it is important to investigate the predictability of the OGCM simulation and also to determine whether the interannual variability and/or the interdecadal variability (Berloff and McWilliams 1999; Meacham 2000) are caused solely by the ocean general circulation. Some types of intermediate models, such as quasigeostrophic (QG) and shallow-water (SW) models, are useful for studying the predictability of ocean currents and low-frequency variations of the ocean general circulation. In the present study, a two-layer SW model is used to investigate the nonlinear behavior of a double-gyre circulation that changes from laminar flow to turbulent flow as the Reynolds number increases.

To clarify the transition from laminar to turbulent flow, many researchers have studied the uniform flow through a right-circular cylinder based on the Reynolds number. In the case of a rotating fluid other nondimensional parameters, the Ekman number and *β* parameter are needed to investigate the transitions of the flow pattern (Matsuura and Yamagata 1985, 1986).

The flow through the right-circular cylinder in a rotating fluid is exterior flow around the solid boundary. Conversely, the oceanic double-gyre circulation is interior flow surrounded by the rectangular solid boundary. These two flow situations occupy the same fluid dynamic point: the interaction between the flows and the solid boundary in a rotating fluid. Therefore, the nondimensional parameters are related to the type of boundary layer, the distribution of vorticity, the location of boundary layer separation, and the inertial subgyre (Ierley 1990). The model used in this study employs Laplacian-type diffusivity as a horizontal diffusion and imposes nonslip boundary conditions in both zonal and meridional directions. The focus of this paper is the investigation of the Reynolds number dependence of double-gyre circulations.

Hidaka (1949) and Munk (1950) first introduced lateral friction in an ocean general circulation model. Munk’s frictional boundary layer was formed with a thickness of *δ _{M}* [=(

*A*)

_{H}/β^{1/3}, where

*A*is a horizontal frictional coefficient and

_{h}*β*is a latitudinal rate of change of the Corioliis coefficient]. His equation with nonlinear terms had a solution that indicated damped standing Rossby waves at the northwestern corner of the single gyre (Moore 1963). Bryan (1963) investigated numerically a single-gyre ocean general circulation based on the Reynolds number using a QG equation with the Laplacian-type horizontal diffusion. In the nonlinear model, the inertial sublayer formed and its thickness was denoted as

*δ*[=(

_{I}*U/β*)

^{1/2}]. The parameter

*δ*/

_{I}*δ*is crucial in simulating systematically the characteristics of an inertial subgyre.

_{M}OGCM simulations changed drastically with the development of an eddy-resolving ocean general circulation model (EGCM) (Holland and Lin 1975; Holland 1978). The EGCM simulation began with the observation of midocean eddies in order to research where and what instabilities occur in the ocean and also to investigate the energy and vorticity exchange between the upper and lower oceans. In general, biharmonic diffusion and a slip boundary condition were employed for the EGCM. In these simulations the recirculations extended zonally and longer than in non-eddy-resolving OGCM simulations (Holland et al. 1983).

Böning (1986) suggested that a barotropic single-ocean gyre using a non-eddy-resolving QG model with slip boundary conditions reached a steady state, even if the nonlinear effect was concentrated, and resulted in a pattern similar to that of average EGCM results. This solution was characterized by free advection dynamics of the inertial sublayer where the *β* term balances nonlinear advection terms. Based on this finding, the formation mechanism of an inertial subgyre was investigated and it was found that surplus vorticity generated at the western boundary advected to form the inertial subgyre (Cessi et al. 1987). Greatbatch (1988) investigated the relationship between the zonal length *λ _{M}*, the meridional width

*λ*, and the inertial length

_{p}*λ*for Böning’s steady inertial subgyre and EGCM recirculations. The former solution was defined by

_{I}*λ*. This result may be related to the difference between the laminar flow of ocean general circulations and the fully developed turbulence of ocean general circulations.

_{I}Recently, research on double-gyre circulations has been carried out from the viewpoint of nonlinear bifurcation dynamics (Dijkstra 2000). The solution of a barotropic double-gyre circulation in the QG equation displayed a pitchfork bifurcation, at Re ∼30. The antisymmetric steady solution of double-gyre circulation bifurcated to asymmetric steady solutions: the cyclonic inertial subgyre twins around the anticyclonic part and vice versa (Dijkstra and Katsman 1997). Another steady solution, layer antisymmetric model (Cessi and Ierley 1995) became unstable (Meacham 2000). The pitchfork bifurcation observed in symmetric QG models became an imperfect pitchfork in SW models (Jiang et al. 1995; Speich et al. 1995; Dijkstra and Molemaker 1999; Nauw and Dijkstra 2001; Nauw et al. 2004). In addition, as the Reynolds number increased, a Hopf bifurcation occurred and a phase trajectory depicted a limit cycle (Jiang et al. 1995).

The difference in behavior of the double-gyre between the QG and the SW models is attributed to the effect of a large amplitude interface height: the elevation for the cyclonic inertial subgyre and the depression for the anticyclonic inertial subgyre (Jiang et al. 1995; Shimokawa and Matsuura 1999). In the SW model, even if the wind stress profile is symmetrical meridionally, asymmetric perturbations occur in the double-gyre circulation and induce asymmetry of the subtropical and subpolar gyres. The degree of asymmetry depends on the parameterization of the lateral momentum friction (Nauw and Dijkstra 2001).

The generation mechanism of a low-frequency oscillation of the double gyre suggests two different scenarios: the time transition of a steady solution along a high-order stable manifold (Cessi and Ierley 1995; Primeau 1998) and the appearance of homoclinic and heteroclinic orbits (Berloff and McWilliams 1999; Meacham 2000; Chang et al. 2001; Simonnet et al. 2005). The periodic variation of the double-gyre circulation appears to be due to the merging of two nonoscillatory eigenmodes: a P mode, associated with multiple steady solutions, and L mode, associated with the strength of double-gyre energy (Simonnet and Dijkstra 2002). Merging of these two modes induces the gyre mode. Although there has been much research on double-gyre circulations in terms of nonlinear transition ordering to steady, periodic, quasiperiodic, and aperiodic states (Jiang et al. 1995; Berloff and McWilliams 1999; Meacham 2000; Chang et al. 2001; Simonnet et al. 2003a), the transition of flow regimes depending on a larger Reynolds number should also be studied in detail.

The Kuroshio and the Oyashio are the intensive currents of the subtropical circulation and the subpolar circulation in the North Pacific (Stommel and Yoshida 1972). In the 1980s, the two phases of the Kuroshio were discussed in connection with the multiple equilibria of nonlinear differential equations. Charney and Flierl (1981) constructed a model employing coastal geometry. Masuda (1982) drew attention to the two narrow straits near the Kyushu and the Izu Ridge as the bottom topography. These ideas are important to the boundary conditions (inflow and outflow) and limited the of degrees of freedom. Thus, if we are able to show multiple equilibria in a freer situation, such as a double-gyre ocean, the generality of the two phases of intensive currents increases (Schmeits and Dijkstra 2001).

Three mechanisms have been suggested for the appearance of interdecadal variability in the North Pacific Ocean: First, there is an atmospheric teleconnection between the Tropics and midlatitudes (Nitta and Yamada 1989); second, the interdecadal variability of tropical SSTs is due to ocean ventilation in midlatitudes (Gu and Philander 1997); and third, a decadal oscillation is advected into the Tropics by variations in the transport of the surface branch of the North Pacific subtropical cell (Kleeman et al. 1999). All three mechanisms play a crucial role in atmosphere–ocean interaction. However, if the ocean general circulation alone induces the interdecadal variability due to the nonlinear dynamics and if the low-frequency oscillation is generated from different time-scale nonlinear interactions, we propose a new mechanism for the interdecadal variability of the ocean system between the subtropical and subpolar circulations.

The purpose of this paper is to investigate the transition of the double-gyre pattern and to clarify the low-frequency mechanism of the ocean general circulation using a double-gyre, two-layer SW model. In section 2, the model configuration, the parameters used, the energy definition, and the numerical methods are explained. The double-gyre flow patterns dependent on the Reynolds number are described in section 3. In section 4, the transition from periodic to aperiodic solutions, the energy spectrum for the low-frequency oscillation, and the strange attractors for three energy phases—upper kinetic energy KE_{1}, lower kinetic energy KE_{2}, and available potential energy APE*—*are discussed. The study is summarized in section 5 and the Kuroshio with bimodal structures is discussed.

## 2. Model and methods

Studies of the double-gyre circulation are carried out using QG or SW models. Two-layer SW equations are used because it is important to consider the large amplitude effect of the density interface (Huang 1986) and the baroclinicity in the direction of greater realism in the investigation of the double-gyre ocean general circulation, particularly the difference in interfacial displacement between the cyclonic and anticyclonic inertial subgyres influence on their asymmetric evolution (Jiang et al. 1995; Shimokawa and Matsuura 1999). The pitchfork bifurcation from the steady antisymmetric to the steady asymmetric solution is different between the QG and SW models (Dijkstra and Molemaker 1999; Dijkstra 2000).

### a. Model

A two-layer SW model in transport form is used to study the nonlinear behavior of the double-gyre ocean circulation in midlatitudes. The model domain is confined to a rectangular basin given by 0 ≤ *x* ≤ *L* and 0 ≤ *y* ≤ *L* with upper and lower layers: *H*_{1} and *H*_{2}.

The two-layer SW equations with a *β* approximation are

Here *U _{i}*

**i**+

*V*

_{i}**j**=

*h*(

*u*

_{i}**i**+

*υ*

_{i}**j**) =

*h*

_{i}

*υ**, (*

_{i}*i*= 1, 2) are the two layer mass-flux vectors, where

*u*,

_{i}*υ*represent the eastward and northward components of the upper-layer and lower-layer thicknesses. The Coriolis parameter is represented by

_{i}*f*=

*f*

_{0}+

*βy*,

*g′*= (Δ

*ρ*/

*ρ*)

*g*is reduced gravity, and

*g*is the acceleration due to gravity. Horizontal eddy diffusion is represented by a harmonic operator with a coefficient

*A*. Using the quiescent upper- and lower-layer thicknesses

_{h}*H*(

_{i}*i*= 1, 2) and the interface elevation

*η*, the upper- and lower-layer thicknesses are related by

and pressures by

The upper layer of the ocean is driven by a zonal wind stress denoted by *τ*, which is constant in time but varies with latitude in a simple sinusoidal pattern:

where *τ*_{0} is the amplitude. This leads to the formation of a double-gyre circulation: anticyclonic subtropical and cyclonic subpolar. The parameters for the experiments performed are listed in Table 1.

In ocean general circulation modeling an Ekman diffusion and/or a harmonic-type horizontal diffusion or a biharmonic-type horizontal diffusion are generally used. The present research employs harmonic diffusion. The values of the dimensionless inertial and Munk’s viscous boundary layer thickness, *δ _{I}* = (2

*π*

*τ*

_{0}/

*ρ*

*β*

^{2}

*H*

_{1}

*L*)

^{1/2}and

*δ*= (

_{M}*A*/

_{h}*β*)

^{1/3}, are crucial in determining the whole midlatitude basin double-gyre circulation (Jiang et al. 1995; Cessi and Ierley 1995; Simonnet et al. 2003a, b; Berloff and McWilliams 1999).

It is assumed that the normal flow is zero at the lateral and horizontal boundaries and the tangential flow is also zero at both boundaries as follows:

To avoid the barotropic instability at the western boundary (Bryan 1963) and to reach the steady solution, Böning (1986) and Cessi and Ierley (1995) used the slip boundary condition. In other double-gyre simulations, to approximate actual oceanic boundary conditions, the experiments were carried out with partial slip boundary conditions (Haidvogel et al. 1992; Jiang et al. 1995; Ghil et al. 2002). The nonslip boundary conditions are used in this study because this approximation is more compatible with other studies (Berloff and McWilliams 1999; Meacham 2000, etc.).

The three representative parameters in rotating fluid are

and

where *R*_{0} denotes the Rossby number, *E _{L}* is the horizontal diffusion parameter, and

*F*is the rotating Froude number. The Reynolds number Re is substituted for

_{r}*R*

_{0}and

*E*as

_{L}Sixteen experiments are conducted to investigate the transition of the double-gyre circulations using the parameter values shown in Table 2; *F _{r}* is fixed at 0.021.

A convenient description of the globally oscillatory behavior of the system can be given in terms of box-integrated energies. The upper- and lower-layer kinetic energy KE* _{i}* (

*i*= 1, 2) and the available potential energy APE integrated over the entire box area are defined by

Variations of these energies are examined as the representative behavior of flow patterns for each Reynolds number dependence.

### b. Numerical methods

The primitive equations in Eqs. (1)–(6) are solved numerically using the finite difference scheme formulated by Holland and Lin (1975) for a two-layer system. The reduction method of Hockney (1971) for Poisson’s equation with the spectrum method is used. The finite difference in time is carried out using a second-order leapfrog technique and Euler’s backward scheme is employed to eliminate time splitting of the developing solutions. The frictional terms always lag by one time step to avoid linear numerical instability.

The integration time step is 20 min and spatial resolution is 20 km. Because this model uses a nonslip boundary condition on the sides of the basin, it is necessary to resolve Munk’s viscous boundary layer and/or the inertial boundary layer at the western boundary with numerical experiments. In all experiments the inertial boundary layer thickness is a constant 24.8 km and is resolved by two grid points. The range of values of Munk’s boundary layer thicknesses is from 39 to 17 km corresponding to a Reynolds number from 26 to 314 (cf. Table 2).

Until Re = 157 (expt 14), Munk’s boundary layer is also resolved by two grid points. In experiments 15 and 16, Munk’s boundary layer thicknesses are less than 20 km. In those cases, however, the inertial boundary layer is dominant and the results are qualitatively correct. When the fine-resolution simulation of 10 km is run in experiment 12 (Re = 95), nearly the same results are obtained between 20- and 10-km resolution.

## 3. Reynolds number dependence

The Reynolds number [Eq. (12)] changes in proportion to strength of the wind stress *τ*_{0} and in inverse proportion to the horizontal frictional coefficient *A _{h}* in the present numerical setting (

*β*= 2 × 10

^{−11}m

^{−1}s

^{−1},

*H*

_{1}= 1000 m,

*H*

_{2}= 4000 m,

*ρ*

_{0}= 1000 kg m

^{−3}). Only

*A*is changed, not

_{h}*τ*

_{0}(cf. Table 2). A chart classifying the double-gyre circulation based on the Reynolds number from Re = 26 to 314 is schematically shown in Fig. 1.

### a. *Steady solution (*Re ∼ *39)*

The pattern of steady double-gyre circulation appears until Re of approximately 39 in the upper layer. Antisymmetric small inertial subgyres appear in experiment 1 (Re = 26) and experiment 2 (Re = 31) (cf. Fig. 1b). In experiment 3 (Re = 39), the cyclonic and anticyclonic inertial subgyres evolve asymmetrically (the northward cyclonic gyre is somewhat stronger than the southward anticyclonic gyre) and the separation point of the double-gyre circulation at the western boundary shifts slightly southward (cf. Nauw et al. 2004). The different upper-layer thicknesses between cyclonic and anticyclonic inertial subgyres induce a steady asymmetric pattern: the cyclonic vortex with interface elevation intensifies unlike the anticyclonic vortex (Chassignet and Gent 1991). This steady-state asymmetry is different from the asymmetric solution obtained by Cessi and Ierley (1995) and is dependent on the lateral frictional form (Nauw et al. 2004). An imperfect pitchfork bifurcation in the present model takes the branch with the jet-down solutions connected to the branch with nearly antisymmetric solutions (cf. Fig. 1b in Nauw et al. 2004).

In experiment 1, the antisymmetric inertial subgyres in the lower layer circulate baroclinically in the opposite direction of the upper layer. In experiment 2, the wavenumber 2 circulation appears in the lower layer (not shown).

### b. *Periodic solution (*Re ∼ *52)*

As the Reynolds number increases to 52 (*δ _{I}*/

*δ*∼ 0.8), the energy trajectory forms a periodic attractor, a limit cycle due to Hopf bifurcation (Fig. 4a), as well as the reduced-gravity SW model case (Speich et al. 1995; Jiang et al. 1995). Figure 2a shows a time series of energy from experiment 4 (Re ∼ 52) in which the upper-layer kinetic energy KE

_{M}_{1}, the lower-layer kinetic energy KE

_{2}, and the available potential energy APE are integrated in the entire rectangular box for 250 years. A prominent spectral peak appears at year 5.937 for both KE

_{1}and KE

_{2}(cf. Fig. 3a). The periodic pattern of KE

_{2}indicates a double-peak variation: a larger amplitude and a smaller amplitude (see Fig. 2a). Another power spectrum peak appears at year 2.90. A reduced-gravity SW model simulation (Jiang et al. 1995) had a peak of 2.8 yr. They observed that the peak period was prolonged as the aspect ratio

*L*/

*D*increased. The model setting for

*L*/

*D*in the present study is twice as large as that of Jiang et al. (1995) and Speich et al. (1995) and our peak period is nearly two times longer. Spectral peaks in Fig. 3a appear at

*T*= 5.937,

*T*/2,

*T*/3,

*T*/4, and

*T*/9 because the variations of energy are not sinusoidal. The distribution of spectral peaks of KE

_{2}is nearly the same as that of KE

_{1}, although a power spectral peak at

*Τ*/2 is comparable with that at

*T*in KE

_{2}.

### c. Aperiodic solutions (*Re* ∼ 57 and *Re* ∼ 70)

As Re increases to 57, an aperiodic variation of antisymmetric and asymmetric patterns appears in the results of experiment 10. In this case, time series of energy show that there are aperiodic times (10–100 and 140–220 yr) in which interannual variability intermingles with decadal variability and periodic terms with approximately 4-yr oscillations (100–130 and 220–250 yr) over 250 yr (Fig. 2b). Spectral peaks appear in years 4.45, 8.91, and 12.47 (Fig. 3b) (8.91 yr is almost double the period of 4.45 yr). During intermingled periods with interannual and decadal variations (10–100 yr), the periods of low-energy level including high-frequency variations (2.4 months) continue for 4–5 yr. The inertial subgyres vary with asymmetric patterns in the low-energy periods.

In the transition from an antisymmetric period to an asymmetric period, KE_{2} increases rapidly due to energy exchanges with APE and KE_{1} (cf. Fig. 2b). The variations of APE and KE_{1} display similar time-dependent patterns. At Re ∼ 57, the asymmetric solution and the antisymmetric solution move in cycles and the asymmetric period continues longer than the antisymmetric period. This occurs because the asymmetric solution is stable in the local bifurcation and the antisymmetric solution is unstable (Dijkstra and Molemaker 1999).

A phase trajectory of energy in experiment 10 is shown in Fig. 4b. The trajectory consists of large decadal-period loops and small interannual-period (4–5 yr) loops. In experiment 10, the low-energy periods, denoted in red, continue longer than the high-energy periods.

In the case of Re = 70, experiment 11 (cf. Fig. 1), the time series of energy are more complex than those at Re = 57, experiment 10. The lower-layer energy, KE_{2}, becomes one order larger than at Re = 57 (cf. Figs. 2c, 2b). High-energy periods with nearly constant values of energy exist alternately with low-energy periods in which the amplitude varies significantly over 250 yr. In high-energy periods, KE_{1}, KE_{2}, and APE are in phase and, when APE and KE_{1} decrease rapidly, the lower-layer energy, KE_{2}*,* increases rapidly due to energy exchanges.

Power spectra of KE_{1} and KE_{2} show peaks at *T*_{1} = 5.937 yr ( *f*_{1} = 0.1684 yr^{−1}) and at its rational numbers, 7*T*_{1} = 41.56, (7/2)*T*_{1} = 20.78, (7/3)*T*_{1} = 13.85, and (7/11)*T*_{1} = 3.775. The maximum peak appears at year 7.335, as shown in Fig. 3c. This peak emerges from the mode coupling of a 5.937-yr periodic oscillation at Re = 52 and a 6.234-yr ( *f*_{2} = 0.1604 yr^{−1}) periodic oscillation at Re = 54 (cf. Fig. 10a) and the period is the inverse of the frequency, 4*f*_{2} − 3*f*_{1}.

During high-energy periods, antisymmetric inertial subgyres extend over the zonal distance of a quarter basin and the anticyclonic inertial subgyres remain steady (cf. Figs. 1d and 5b). Alternatively, during low-energy periods, an anticyclonic inertial subgyre shrinks significantly and the eastward current meanders widely from the straight eastward jet exhibited during high-energy periods (cf. Figs. 1c and 5a).

### d. Metastable state (*Re* ∼95)

When Re is increased to 95, the behavior of double-gyre circulation is similar to that in experiment 11 until approximately the 10th year. After that it becomes a metastable state and the recirculation, exhibiting a nearly antisymmetric pattern, stretches eastward toward the midposition of the zonal extension (Fig. 6a). At Re = 95 an extreme switch between 10 yr before and 10 yr after shows the transition from the Böning (1986) type solution of inertial subgyre, which is dependent on the inertial length scale, *λ _{I}*, to the EGCM solution, which is not dependent on

*λ*.

_{I}As can be seen in the time series of KE_{1}, KE_{2}, and APE in Fig. 2d, it takes approximately 2 years to spin up from the initial quiet state and then a variation of inertial subgyres become significant until about year 10. The energy variations decrease and become quiet from 10 yr until the simulation ends at 250 yr (see Fig. 2d). Enlarging the amplitude of energy KE_{1}, KE_{2}, and APE for the low variable time after 10 yr causes a high frequency oscillation to be superimposed on a regular low frequency oscillation (see Fig. 7). Spectrum peaks of high frequencies appear at year 0.157 (57.3 days) and year 0.154 (56.2 days). Two clusters of spectrum peaks are exhibited at about 1/2 (0.077 yr) and 1/3 (0.054 yr) in the cluster of maximum peaks (Fig. 3d). In experiment 12 the low-frequency power spectrum peak also appears at 7.793 yr. Period 0.157 yr (57.3 days) corresponds to the observed time scale in which mesoscale eddies occur because of baroclinic instability. The time series of energy in experiment 12 show resonance patterns, but the amplitude is very low (Fig. 7).

Figure 4d shows the trajectories of energies KE_{1}, KE_{2}, and APE for experiment 12. During the first 10 yr the trajectory drifts in the three-dimensional phase space of KE_{1}, KE_{2}, and APE. After approximately 10 yr it converges into a narrow region while oscillating. This trajectory is depicted by the time series of KE_{1}, KE_{2}, and APE in Fig. 2d. The loop of the trajectory is explained as a movement to an attractor. The expanded feature of the attractor shows a topological manifold (Fig. 7d). The critical Reynolds number is Re ∼ 95, that is, the transition from a behavior of recirculation in a reduced-gravity model type to that in a two-layer model type (cf. Figs. 17 and 20 in Berloff and McWilliams 1999).

### e. Transition to fully developed turbulence (*Re*\E ∼ 112 → 314)

When Re = 112, the variation due to an instability of return flow in the recirculation is more remarkable than that of the eastward jet formed by the subtropical gyre and the subpolar gyre. The leading edge of the eastward jet fluctuates intensely relative to the case of Re = 95. The return flow in the subpolar recirculation varies more intensely and forms more eddies than that in the subtropical recirculation (Fig. 8a). At Re = 112, spectral peaks appear at 55, 73, 116, and 226 days, which are mesoscale time scales in the ocean (not shown). The characteristics of the spectral peak appearance show that clusters form at about 55, 73, 116, and 226 days. Time series of energy show that the spin up is complete at 2 yr, the inertial subgyres fluctuate until 6 yr, the quiet state continues from 6 to 16 yr, and then modulation occurs with periods of about 10 yr. A trajectory forms a torus and converges becoming bloblike after 16 yr (not shown).

When Re increases to 157, an eastward jet extends prominently to the open ocean and varies intensely. However, the expected meandering of the eastward jet does not occur (cf. Fig. 8b). The amplitude of energy variations increases and, specifically, the value of KE_{2} increases about twofold when compared with the case of Re = 112 (cf. Fig. 12a, below). Irregular high-amplitude and low-amplitude variations appear repeatedly in the time series of energys with multiyear periods (Fig. 2e). Phase relations of KE_{1}, KE_{2}, and APE show that the variation between APE and KE_{2} is out of phase and that there is a phase shift of 90° between KE_{1} and KE_{2}. In the power spectrum (Fig. 3e), a beat with a period of 31.1714 yr occurs at 0.8540 and 0.8781 yr. With the beat period defined as *T _{1}*, spectral peaks appear at (2/5)

*T*, (4/13)

_{1}*T*, (1/4)

_{1}*T*, and (1/11)

_{1}*T*in Fig. 3e. As shown in Fig. 4e, the trajectory forms oval and three-dimensional shapes.

_{1}At Re = 209, the variation of the double-gyre circulation becomes more significant compared to the case of Re = 157, and the jet forms by the meeting of subpolar and subtropical circulations extending farther eastward from the zonal midpoint (Fig. 8c). Mesoscale eddies generate actively in the recirculation regions. The variation pattern in the lower layer is similar to that in the upper layer and the disturbances become more barotropic than in the case of smaller Re (cf. Figs. 8 and 9). The amplitude of the energy time series in experiment 15 increases more rapidly than those in experiments 13 (Re = 112) and 14 (Re = 157) (Fig. 12), and interdecadal variation starts to occur (not shown).

Although in experiment 16 with Re = 314 the intense variations of the streams in the box are similar to those in experiment 15, the extension length of the eastward jet and eddy generation continually increase and the jet meanders significantly (Figs. 8d and 9d). A more barotropic state is displayed and, hence, the values of KE_{2} reach those of KE_{1} (Figs. 2f and 12, below) and the double-gyre circulation generates fully developed turbulence. The interdecadal power spectra peaks appear at 41.51, 17.81, and 11.34 yr (Fig. 3f). In high-energy periods, the eastward jets extend more zonally and are longer, while in low-energy periods the eastward jets are shorter. The generated mesoscale eddies are barotropic and two-dimensional, and create a red cascade (Rhines 1977).

## 4. Discussion

The numerical experiments of double-gyre circulation using a two-layer shallow-water equation show that the double-gyre dynamics change drastically at Re = 95 (*δ _{I}*/

*δ*= 0.97) within the range Re = 26 to 314. Here Re = 95 is defined as the critical Reynolds number Re

_{M}*at which the double-gyre circulation is divided into two flow classes: a deterministic turbulence and a fully developed turbulence.*

_{c}At Re < Re* _{c}*, as Re increases from 26 to 70, steady solutions are obtained at first (section 3a), then periodic solutions (section 3b), and finally quasiperiodic solutions (Figs. 10b and 11b). After a two-frequency oscillation breakdown, aperiodic solutions lead to intermittency and the blue sky catastrophe (section 3c) (see Fig. 1). In reduced-gravity QG models, the solutions make a pitchfork bifurcation: a single antisymmetric solution and two, mirror image, asymmetric solutions (Dijkstra and Molemaker 1999; Chang et al. 2001). In this model, the upper-layer thicknesses between the subpolar and subtropical gyres are very different and hence the nonlinearity of the continuity equation cannot be ignored (Cushman-Roisin et al. 1992; Matsuura 1995). Because the symmetric solution in SW models is absent, the steady solution takes one side, and the bifurcation diagram becomes an imperfect pitchfork. When Re is increased further, limit cycles occur over the bifurcation line of the asymmetric solution due to a Hopf bifurcation. The solutions change from periodic to aperiodic through quasiperiodic as Re increases.

Hopf bifurcation first occurs at approximately Re = 52 and a periodic time of 5.937 yr ( *f*_{1} = 0.1684 yr^{−1}) (cf. Figs. 2a, 3a and 4a). As Re increases slightly, a second Hopf bifurcation occurs at Re = 54 and a periodic time of 6.234 yr ( *f*_{2} = 0.1604 yr^{−1}) (Figs. 10a and 11a). The ratio of frequencies between these two periods, *f*_{2}/*f*_{1}, is 20/21 and hence this value is a rational number. As Re increases to 55, a quasiperiodic solution appears (Fig. 10b). In this case, power spectral peaks of the upper-layer energy KE_{1} appear at *f*_{1}, *f*_{2}, *f*_{1} + *f*_{2}, and 2*f*_{1} + *f*_{2}, and those of the lower-layer energy KE_{2} appear at *f*_{1}, *f*_{1} + *f*_{2}, and 2*f*_{1} + *f*_{2}. The trajectory at Re = 55 is shown in Fig. 11b and a part of the torus evolves. By increasing Re, spectral peaks occur at (8/7) *f*_{1} (5.195 yr) (Fig. 10c) for Re = 55.8, (9/7) *f*_{1} (4.619 yr) (Fig. 10d) for Re = 56, and (4/3) *f*_{1} (4.453 yr) (Fig. 3b) for Re = 57. Thus, the spectral peaks move to the high-frequency side as Re increases and the integrated energies (Fig. 12) change from a quasiperiodic oscillation to an aperiodic oscillation. Basic frequencies *f*_{1} and *f*_{2} have the rational number, |*n*_{1}*f*_{1 }*+ n*_{2}*f*_{2}| (where *n*_{1} and *n*_{2} are integers), and the energy spectral peaks appear at their frequencies (cf. Fig. 10c): the spectral peaks appear at *nf* (where *f* = *f*_{1} − *f*_{2} = *f*_{1}/21 = *f*_{2}/20). This results in phase locking due to torus breakdown in the neighborhood of the quasiperiodic oscillation at Re = 55. The necessary conditions of frequency locking are the coupling of more than two independent frequencies and a dissipation system (Berge et al. 1984). Our numerical experiments satisfy the above conditions. Moreover, the results in the simulation show the devil’s staircase of frequency locking in the range Re = 55.8–57, at which point the frequency peak increases as Re increases.

Attractors of periodic (Re = 52 and 54), quasi-periodic (Re = 55), and additional aperiodic solutions depict a complicated looping of the limit cycle (Figs. 4a and 11a); the torus (Fig. 11b) and the breakdown of the torus (Figs. 11c and 11d); and the chaotic trajectory (Fig. 4b), respectively. The loops take on two-band structure at Re = 55.8 and three-band at Re = 56.1. At Re = 57, the trajectory reaches a chaotic pattern. In this case, a particularly interesting variation of the time series of energy reveals that aperiodic and periodic phases emerge alternately (see Fig. 2b). Both actual and model situations have shown that the behavior of aperiodic chaos and the regular oscillation are able to coexist (Nicolis and Prigogine 1989). The case of Re = 57 corresponds to the previous description. This results in an intermittency phenomenon switching to chaos (cf. Fig. 1).

Three theories have been proposed for the transition from a periodic oscillation to an aperiodic oscillation or chaos generation. The first theory indicates the transition from periodic solutions to aperiodic solutions through quasiperiodic solutions with two or three basic frequencies (*T* ^{2} or *T* ^{3} torus) (Ruelle and Takens 1971; Curry and York 1977). The second theory states that, when associated with the control parameter, a basic frequency, *f,* evolves to periodic doubling *f* /2, *f* /4 (corresponding periods 2*T,* 4*T*) and reaches chaos due to a periodic-doubling cascade (Berge et al. 1984). The third theory states that the periodic attractor changes the strange attractor directly due to instabilities (Berge et al. 1984), when associated with intermittency. Numerical experiments in our study show the root of the periodic attractor (Re = 52 and 54), the quasiperiodic attractor (Re = 55), phase locking (Re = 55.8 and 56), intermittency (Re = 57), and the chaotic attractor (Re = 70).

As Re → 70, high-energy period values remain nearly constant and low-energy periods in which the time-dependent energy amplitude varied intensely repeat alternately. During high-energy periods, the inertial subgyres fully develop and extend eastward (Fig. 5b). Conversely, during low-energy periods, the inertial subgyres shrink and the eastward jet meanders while entering the open ocean (Fig. 5a). Aperiodic variations of high and low energy periods in the double-gyre circulation were simulated by Chang et al. (2001, their Fig. 19) and Meacham (2000, his Figs. 13 and 14). The barotropic QG equation shows two exact solutions in the double-gyre circulation of the circular boundary (Yamagata and Greatbatch 1985): a high-amplitude antisymmetric solution and a low-amplitude antisymmetric solution. These antisymmetric solutions are unstable (Meacham 2000).

The variations of inertial subgyres and energy for experiment 11 are discussed from the viewpoint of nonlinear dynamics. As shown in Fig. 1, recirculation patterns alternately select a stable asymmetric solution in low-energy periods (stable fixed point; see Fig. 1c) and an unstable antisymmetric solution in high-energy periods (unstable fixed point; see Fig. 1d). A saddle-node bifurcation exists in the trajectories of phase spaces KE_{1}, KE_{2}, and APE and trajectories converge spirally to the saddle focus (B in Fig. 4c) during low-energy periods. Additional trajectories leave the saddle focus to connect to the saddle-node periodic orbit (A in Fig. 4c) with a homoclinic orbit while increasing total energy. When the total energy reaches a value in high-energy periods, the saddle and node collided and the long period of high energy continues with short oscillation periods. In this case, the Rossby basin mode and the two nonoscillatory modes, P and L modes, exist intensely. This Rossby mode is composed of baroclinically unstable Rossby waves (50–60-day period). The slow oscillation induced by the gyre mode, which originates from the merging of the P and L modes, switches between the antisymmetric and asymmetric solutions (cf. Fig. 1): they appear as the barotropic instability of the inertial subgyre with a period of 7–20 yr. This type of bifurcation is called the blue-sky catastrophe in nonlinear dynamics (Turaev and Shilnikov 1995).

The two-layer ocean in the SW equations shows a metastable state of recirculations at *H*_{1}/*H*_{2} = 0.25 and Re = 95, which differs from that of the reduced-gravity ocean. While experiment 12 reaches the metastable state after spinning up, inertial subgyres in the subpolar and subtropical gyres fluctuate and form mesoscale eddies. Berloff and McWilliams (1999, Fig. 20c) discussed the metastable state of recirculation for the two-layer QG equations and showed its result at *A _{h}* = 6 × 10

^{6}cm

^{2}s

^{−1}and

*H*

_{1}/

*H*

_{2}= 0.08. The parameter values at this situation are critical because of the division into deterministic turbulence and fully developed turbulence. In the metastable state, the barotropic instability for the recirculations is suppressed by the baroclinic instability (Berloff and McWilliams 1999). Although in experiment 11 (Re = 70) the spectral peaks appear prominently at 3.78, 7.335, and 20.78 yr, in a metastable state at Re = 95 (expt 12) the time series of energy exhibits long-term modulation with a period of 7.793 yr consisting of 0.157-yr and 0.154-yr high-frequency oscillations.

In APE and KE_{2}, high-frequency oscillations are superimposed over a low-frequency oscillation (see Figs. 7c and 7b). For KE_{1}, the beat appears from high-frequency variations (Fig. 7a). In this case, maximum envelop peaks for high-frequency oscillations shift from minimum peaks. The long-term variations of APE and KE_{2} are out of phase and the envelope peak of KE_{1} is slightly ahead of KE_{2} (Fig. 7).

A long-term period of 7.793-yr variation occurs in experiment 12. Sinusoidal oscillations with nearly the same periods, cos(2*π*/*T*_{1})*t* and cos(2*π*/*T*_{2})*t* (where *T*_{1} = 0.157 and *T*_{2} = 0.154 yr), are assumed. The energies KE_{1}, KE_{2}, and APE are written for linear and nonlinear combinations of sinusoidal functions as

The sum of two sinusoidal waves, the first and second terms on the right-hand side of Eqs. (15)–(17), produces a beat with a period of ∼8 yr. As shown in Fig. 6b, two baroclinic instability waves with nearly the same periods of 0.157 and 0.154 yr occur in the recirculation region and in the eastern half of the basin. The product of two sinusoidal waves with slightly different frequencies can be represented as a sum of waves that have both longer and shorter wavelengths than the original waves—a phenomenon that is crucial in nonlinear interactions. The third term on the right-hand side of Eqs. (15)–(17) reflects the aforementioned phenomenon. The time series of energy in Fig. 7 for experiment 12 corresponds to those of Eqs. (15)–(17) as shown in Fig. 13. Hence, it is found that the time series of energy in experiment 12 constitutes the linear and nonlinear combination of two slightly different sinusoidal waves.

For the relationship between an inertia size *λ _{I}* and a meridional size

*λ*or a zonal size

_{M}*λ*of recirculation, Greatbatch (1988) indicated that in Böning’s (1986) steady solution

_{p}*λ*and

_{M}*λ*depended on

_{p}*λ*and that

_{I}*λ*and

_{M}*λ*in the EGCM recirculations (Holland 1978; Schmitz and Holland 1982) were not dependent on

_{p}*λ*. In the numerical experiments based on Re, although solutions at 39 < Re < 95 are not steady, they are related to the Böning (1986) type solution where

_{I}*λ*and

_{M}*λ*depended on

_{p}*λ*. The numerical solutions at Re > 95 correspond to a type of nondependence on

_{I}*λ*. At Re = 95,

_{I}*δ*/

_{I}*δ*is nearly equal to 1, which is a critical value dividing dynamic types.

_{M}Mesoscale eddy generation is induced by the collapse of an inertial subgyre when Re < Re* _{c}*. Using two-layer SW equations, variation of the subpolar inertial subgyre is significant because the interface displacement of the two layers is large and the nonlinearity of continuity equations cannot be ignored (cf. Matsuura 1995; Shimokawa and Matsuura 1999). In experiment 11 the inertial subgyres are much larger than the radius of deformation, and the subpolar inertial subgyre collapses into several small eddies due to instabilities (cf. Fig. 5a).

The deterministic chaos (blue-sky catastrophe) at Re = 70 changes to the metastable state at Re = 95. When the total kinetic energy (TKE) and the available potential energy (APE) are stored in the rectangular box by the wind stress forcing, the ratio of APE to KE_{1} increases at Re = 95 more than at Re = 70 (cf. Fig. 12b). The ratio of APE to KE_{1} at Re = 70 is 6.67 and at Re = 95 it is 6.83. This implies that the buildup in APE (conversion of KE_{1} to APE) increases from the case at Re = 70 to the case at Re = 95. In the two-layer ocean, not only the forcing or the lateral friction, but also the energy exchange between kinetic energy and available potential energy play important roles in forming the circulation pattern.

When Re > Re* _{c}*, mesoscale eddies form as the geostrophic turbulence occurs at several genesis locations: the westward return flow of recirculation and the interior jet (Figs. 8 and 9). The return flow of recirculation shows a prominent wavy pattern due to baroclinic instability (Pedlosky 1987) and sheds eddies while amplifying the wavy pattern (Fig. 8) at Re = 112. The variation of the recirculation at Re = 112 radiates basin-mode Rossby waves more significantly than at Re = 95. The basin-mode Rossby waves with wavenumber 1 meridionally begin to divide into several wavenumbers (Fig. 9a). Mesoscale eddies in the open ocean emerging from the basin-mode Rossby waves evolve barotropically. When Re is increased to 112 (expt 13), the trajectories oscillate around a fixed point (not shown).

Mesoscale eddies in the open ocean strengthen and become larger as Re increases from 157 to 314 (cf. Fig. 9). This occurs because the scale of mesoscale eddies in the open ocean is estimated by *λ _{β}* = 2

*π*2

*U*/

*β*, where

*U*is a square mean velocity of turbulence (Rhines 1977). This result implies that the scale of mesoscale eddies in the lower ocean increases because

*U*increases as Re increases. Actually, kinetic energy in the lower-layer KE

_{2}increases rapidly over Re = 157 (Fig. 12).

Another genesis area of mesoscale eddies is the intense eastward jet where its length elongates from Re = 112 toward Re = 314. The amplitude of the meandering jet also becomes significant due to increasing shear instability (Fig. 8). At about Re = 157, detached eddies form at the leading edge of the eastward jet.

In the double-gyre numerical experiments using a two-layer SW model, when Re < Re* _{c}* stationary solutions lose their stability because of Hopf bifurcations that give rise to periodic solutions with two basic frequencies:

*f*

_{1}(

*T*

_{1}= 5.937 yr) and

*f*

_{2}(

*T*

_{2}= 6.234 yr). Physically, the gyre mode appears through the merging of two eigenmodes: the P mode, which is associated with the occurrence of multiple equilibria, and the L mode, which controls the intensity of the entire double-gyre flow (cf. Simonnet and Dijkstra 2002). Spectral peaks in experiment 11 appear at 3.778 and 7.335 yr as the interannual variability and at 13.85, 20.78, and 41.56 yr as the interdecadal variability. The interdecadal frequency corresponds to the switch from the antisymmetric to the asymmetric inertial subgyre through a homoclinic orbit, which has a period of 10–20 yr. High-energy states (antisymmetric gyres) alternate with low-energy states (asymmetric gyres) at 5-yr intervals.

When Re > Re* _{c}*, mesoscale eddies with a 50-day period are initially generated due to baroclinic instability and the interannual and interdecadal variability emerges as Re increases (cf. Fig. 3). In experiment 16 prominent interdecadal variability appears at 41.56 and 17.81 yr in the white noise spectrum profile (Fig. 3f) and may occur from the nonlinear stochastic resonance.

## 5. Summary and conclusions

Our numerical experiments show that there are significant differences in the double-gyre flow pattern and the time series of energy before and after the critical value of Re = 95 (*δ _{I}*/

*δ*∼0.97) (see Table 3). When Re < Re

_{M}*, variation of the double-gyre circulation leads to deterministic turbulence from the asymmetric and antisymmetric changes of inertial subgyres. This transition of double-gyre circulation is similar to the results obtained from the reduced-gravity models (Chang et al. 2001; Simonnet et al. 2005) where periodic to aperiodic solutions show the movement from a Hopf bifurcation to a homoclinic bifurcation. Instability phenomena exhibit an obvious gyre mode emerging from the merging of nonoscillatory P and L modes rather than from the Rossby wave basin mode (Simonnet and Dijkstra 2002; Simonnet et al. 2005). At Re = 57, the homoclinic orbit is shown as the connection of the limit cycle of the asymmetric solution with the unstable antisymmetric state. Strange attractors display two-dimensional-like patterns in the KE*

_{c}_{1}, KE

_{2}, and APE phase space and long-term variations of 3–4 and 7–8 yr appear in the power spectra of the KE

_{1}and KE

_{2}energy. At Re = 70, the blue-sky catastrophe and the interdecadal variability are found.

When Re > Re* _{c}*, fully developed geostrophic turbulence occurs through baroclinic instability in the recirculation return flow of the formation of mesoscale eddies in the open ocean, and the ring formation in the intensive eastward jet. In this case, long-term variability has an interdecadal time scale of 10–20 yr and/or 40 yr.

As Re increases from Re* _{c}*, three events occur: 1) the Rossby basin mode with meridional wavenumber 1; 2) isotropic eddies due to instability; and 3) mesoscale eddies evolve on a larger scale because of the red cascade. As described in Rhines (1977), in the evolution of mesoscale eddies changes the baroclinic eddy is much larger than the internal radius deformation, the baroclinic eddy is nearly as large as the internal radius deformation, and the barotropic eddy is larger than the radius deformation. According to Williams (1978) and Rhines (1977), initial geostrophic turbulence evolves to a nonhomogeneous zonal current in the

*β*plane. This tendency could not be duplicated in our numerical experiments.

There are two distinct paths of the Kuroshio western boundary current (Robinson and Taft 1972). The bimodality of the Kuroshio is reproduced in Fig. 14 using a 1/12° high-resolution ocean general circulation model (the Modular Ocean Model version 3, MOM3) forced by the wind stresses of National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis. Until the 1980s, the cause of bimodality was attributed to the bathemetry and/or topography south of the Japanese islands (Robinson and Taft 1972; White and McCreary 1976; Chao 1984; Charney and Flierl 1981; Masuda 1982; Yasuda et al. 1985). In the double-gyre numerical experiments discussed here, when Re = 70, the bimodality appears with high-energy periods corresponding to the nonlarge meandering path (Fig. 5b), and low-energy periods corresponding to the large meandering path (Fig. 5a). Spectral peaks of the Kuroshio exhibit interannual to interdecadal periods with the most dominant period ∼20 yr and the next prominent period at 7–8.5 yr (Kawabe 1987). The model shows interannual variability of 7 yr and interdecadal variability of 20 yr at Re = 70. Experiment 11 exhibits the dynamic essence corresponding to the Kuroshio bimodality (cf. Dijkstra and Ghil 2004). Therefore, we contend that the bimodality of strong currents is a gyre-scale phenomenon and not a limited, local current phenomenon (Qiu and Miao 2000; Schmeits and Dijkstra 2001).

In future investigations of the double-gyre circulation based solely on the Reynolds number, the Froude number dependency, implying the effect of stratification, will be extended. In addition, the effect of a time-dependent forcing to account for external forcing resonance will be investigated. Also, we will introduce the energy balance model of KE_{1}, KE_{2}, and APE and investigate the behavior of their solutions compared to the numerical experiments reported here.

## Acknowledgments

The authors thank W. Sasaki for simulating an ocean general circulation using MOM3. This research was supported by the “Study on Extreme Weather and Water-Related Disasters Due to Climate Changes” research project at the National Research Institute for Earth Science and Disaster Prevention of Japan.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Dr. Tomonori Matsuura, National Research Institute for Earth Science and Disaster Prevention, 3-1 Tennodai, Tsukuba, Ibaraki 305-0006, Japan. Email: matsuura@bosai.go.jp