## Abstract

The Kuroshio Extension (KE) flow in the North Pacific Ocean displays a very distinctive decadal variability of bimodal character involving two completely different states (a large-meander “elongated” state and a small-meander “contracted” state) connected by very asymmetric temporal transitions. Although such a flow has been widely studied by means of a suite of mathematical models and by using several observational platforms, a satisfactory theoretical framework answering quite elementary questions is still lacking, the main question being whether such variability is induced by a time-varying wind forcing or, rather, by intrinsic oceanic mechanisms. In this context, the chaotic relaxation oscillation produced by a process-oriented model of the KE low-frequency variability, with steady climatological wind forcing, was recently recognized to be in substantial agreement with altimeter data. Here those model results are further compared with a comprehensive altimeter dataset. The positive result of such a comparison allows the conclusion that a minimal model for the KE bimodality has been identified and that, consequently, nonlinear intrinsic oceanic mechanisms are likely to be the main cause of the observed variability. By applying the methods of nonlinear dynamical systems theory, relevant dynamical features of the modeled flow are then explained, such as the origin of the relaxation oscillation as a consequence of a homoclinic bifurcation, the spatiotemporal character of the bimodal behavior, and the degree of predictability of the flow in the different stages of the oscillation (evaluated through a field of finite-time Lyapunov exponents and the corresponding Lagrangian time series).

## 1. Introduction

The bimodal behavior of the Kuroshio Extension (KE) in the North Pacific Ocean has fascinated physical oceanographers since indications of this phenomenon were found (Taft 1972). Why should a western boundary current extension switch between a large-meander and a small-meander state (denoted also as “elongated” and “contracted” states, respectively; e.g., Qiu 2002) and back in a few years time? Why does this bimodal behavior not appear in other western boundary currents, such as the Gulf Stream? With the analysis of satellite altimeter data (e.g., Qiu and Chen 2005) and those of in situ measurements (e.g., as obtained within the Kuroshio Extension System Study program: http://uskess.org/index.html), a quite detailed description of the different states of the KE and of their transition behavior is now available.

There is still, however, no consensus on what processes cause the observed bimodal low-frequency variability. A possible explanation relies on external causes, such as variable atmospheric forcing, that would act through an essentially linear time-dependent Sverdrup balance and with the intervention of oceanic Rossby wave propagation (e.g., Miller et al. 1998; Deser et al. 1999; Qiu 2002, 2003; Qiu and Chen 2005). In particular, Qiu and Chen (2005) showed that the weak, spatially broad wind-driven Rossby wave variability computed in the KE area appears to be in synchrony with the much stronger and spatially sharper KE bimodal variability observed in altimeter data, but they also showed (their Fig. 9) that the latter can by no means be explained by the former. The other possible explanation of the KE bimodality is the occurrence of highly nonlinear processes internal to the ocean system (e.g., Schmeits and Dijkstra 2001; Pierini 2006).

A theory of the bimodality of the KE should include a few essential ingredients. It should explain (i) why the KE can basically be in two different (transient) states and the origin of the spatial patterns of these states, (ii) the decadal time scale of the transition between the two states and the details of the transition itself, and (iii) why there is so much irregularity in the KE path and also in the temporal behavior during the transition from the small- to the large-meander state, while spatial and temporal variability are relatively low during the opposite transition.

An interesting issue concerns the identification of a minimal model that captures the KE bimodality in sufficient detail as to satisfactorily answer the questions above, while still being simple enough so that precise analyses of mechanisms can be made. On one end of the spectrum of models are eddy-resolving general circulation models (EGCMs), which seem to capture the bimodality in hindcast runs (e.g., Taguchi et al. 2005, 2007; Qiu et al. 2008). However, analysis of these results has not led to a satisfying theoretical framework, as indeed many phenomena act in concert to provide the variability, but it is not easy to figure out which processes are leading the melody and rhythm and which ones are just playing an arrangement role. On the other end of the spectrum are barotropic quasigeostrophic (QG) models of the double-gyre wind-driven flow in rectangular basins (e.g., Cessi and Ierley 1995; Dijkstra and Katsman 1997; see also Dijkstra 2005). Although these models display multiple states, oscillatory variability, and nonlinear mechanisms of transition between different states, their results cannot obviously be directly compared with any observational characteristics. So, working our way upward in the hierarchy of models (from the QG to the EGCMs), what will be a good candidate for a minimal model?

In a recent process-oriented model study, Pierini (2006) applied a reduced-gravity shallow-water model to an idealized North Pacific Ocean with a fairly realistic western coastline, and computed flows that are forced by a fairly realistic time-independent climatological double-gyre surface wind stress field. The results are striking in that the low-frequency variability thus produced by purely intrinsic oceanic mechanisms yields a relaxation oscillation between a large- and a small-meander state that is in surprisingly good agreement with that observed from satellite measurements (Qiu and Chen 2005). Could this model serve as a minimal model? Does it provide satisfactory answers to issues (i)–(iii)? If so, this would imply that baroclinic instability, bottom topography, complex time-dependent atmospheric forcing, as well as effects of the meridional overturning circulation somehow would not be essential to answer the main issues above. Although much analysis of the patterns and time evolution of the flow were already given in Pierini (2006), no definitive answers were provided to the issues (i)–(iii).

The main aim of this paper is to develop a nonlinear theory of the KE low-frequency variability that is based on the recognition that the model by Pierini (2006) is indeed an adequate minimal model for the problem under study, which in turn implies that internal oceanic mechanisms are responsible for the KE bimodal variability. First, in section 2, we confront the model results and observations in more detail to provide support that the model is indeed able to capture the most important characteristics in observations at least related to the issues (i)–(iii), and hence can serve as a minimal model. To provide answers to those questions we use the same approach followed in QG models in section 3 to investigate behavior of the model for different values of the lateral eddy viscosity parameter. We show that by going from high to low lateral friction, flow behavior and flow transitions are found that are already understood in the QG models. From a comparison with the results of those simpler models, we are then (section 4) able to provide answers to the issues (i)–(iii) in terms of physical mechanisms using the language of nonlinear dynamical systems theory. In the same section, the probability density function, the state space velocity, and a field of finite-time Lyapunov exponents, along with the corresponding Lagrangian time series, have been computed in an appropriate two-dimensional space in order to shed light on the structure of phase space and on the degree of predictability of the flow variability. Finally, in section 5 a summary is given and conclusions are drawn.

## 2. Comparison between model results and altimeter data

### a. The mathematical model

The model used here is the same as that in Pierini (2006, hereafter P06). This model is based on the reduced-gravity shallow-water equations forced by an analytic constant-in-time wind forcing derived from climatological winds. The domain of integration is shown in Fig. 1, which also contains a plot of the wind stress curl. The fundamental importance of the schematic coastline introduced in the western part of the ocean on the character of the low-frequency variability of the modeled KE is discussed by Pierini (2008). For technical details and parameter values used the reader should refer to P06; our only control parameter here is the lateral eddy viscosity, indicated by *A _{H}*.

The reduced-gravity (equivalent barotropic) approximation adopted in P06 is very appropriate for the KE [e.g., Qiu et al. (2006), see in particular their Figs. 4–5], and the chosen values of the depth of the surface layer and the reduced gravity fit well with observed values. Therefore, the sea surface height (SSH) field obtained from altimetry can in principle be compared with the same quantity obtained from a reduced-gravity model, as they both represent a good approximation for the upper-layer streamfunction.

### b. Model–data comparison

The issue we will address in this section is whether this model is able to provide the right transition time scale along with the correct characteristics of the path variability. To do this we compare the model results with altimeter data. A first encouraging answer in this respect was given in P06 (his section 4b), where a comparison between model results for a reference bimodal cycle and the merged Ocean Topography Experiment (TOPEX)/Poseidon, *Jason-1*, and *European Remote Sensing Satellites 1 and 2* (*ERS-1/2*) altimeter data presented by Qiu and Chen (2005, hereafter QC05) was discussed under the correspondence: *t* = 145 yr of model integration is equivalent to year 1993 of altimeter data. Here we extend that analysis to a further cycle and include the aid of visual comparison with appropriate time series and with altimetric maps, absent in P06.

We first compare the observed KE pathlength *L*_{KE} (Fig. 2a) and mean latitudinal position of the upstream KE axis (Fig. 2b) (both defined in QC05) with the corresponding model results (Figs. 2c,e and 2d,f, respectively) for *A _{H}* = 220 m

^{2}s

^{−1}. Let us begin by describing the observations (Figs. 2a,b). During the first two years the system is in a large-meander (elongated) state (see the color panels in Fig. 3), corresponding to a small

*L*

_{KE}and a nearly average , both being weakly variable. An abrupt transition occurs during the following year (1995), when

*L*

_{KE}suddenly increases, suddenly decreases, and both start showing a much larger variability. This marks the beginning of a transition period of about 7 yr (roughly from 1995 to 2002) called the unstable mode by QC05, and here denoted the transition phase. During this phase

*L*

_{KE}maintains virtually the same behavior it had at its beginning, while stops decreasing after about two years (1997); it then registers a positive trend that stops after about three more years (2000), and it eventually decreases up to the original value at around year 2002, which marks the end of the transition phase. During the remaining three years the situation is very similar to the initial one, corresponding again to a large-meander state (see the color panels in Fig. 3).

Let us now compare these observations with model results. The first cycle (Fig. 2c for *L*_{KE} and Fig. 2d for ) is the same as discussed in P06, with the same temporal synchronization with the altimetric signal. The successive cycle (Fig. 2e for *L*_{KE} and Fig. 2f for ) presents some quantitative difference but yields the same qualitative behavior. The model behavior agrees fairly well with the observed one, both quantitatively and qualitatively. Both Figs. 2a,c,e and 2b,d,f show two large-meander phases connected by a transition phase with the correct characteristics and timing. The time duration of the transition phase is 6–8 yr (for the first and second cycle, respectively) during which both *L*_{KE} and are highly variable. The shape of the time series of is remarkably similar to the observed one, especially for the first cycle. There is also much weaker temporal variability during the large-meander phase for both *L*_{KE} and , as is the case for the observations. In fact, the mean pathlength for both modeled cycles here agrees very well with the observed one during the large-meander phase (≈1500 km).

There are also features that are not well captured by the model. The observed temporal variability of *L*_{KE} during the northward KE migration is higher than the modeled one (by a factor ∼1.5), and on smaller (from monthly to weekly) time scales (cf. Fig. 2a with Figs. 2c,e; see also the KE paths shown in Fig. 3 of QC05 and in Fig. 4 of P06). This is not surprising since the real high-frequency variability is due to the shedding and merging of mesoscale eddies (as recognized by QC05), so it cannot be realistically reproduced by the present model, which apart from other causes, is only eddy permitting (the grid spacing is Δ*x* = 20 km while the internal Rossby deformation radius is *R _{i}* ≈ 50 km) and does not capture baroclinic instability. Nonetheless, it is important to notice that not only the low-frequency variability but also the high-frequency variability of the model results are

*qualitatively*similar to that observed. As far as is concerned, while the southernmost location of the jet is ≈ 33.5°N in both altimeter and model data, the most northerly position is ≈ 36°N in altimeter data, while it is ≈ 35°N in the model data.

We now analyze the spatial patterns. In P06 a preliminary model–data comparison of SSH maps was carried out for the first cycle (*t* = 145–157 yr), but without providing graphical evidence of the altimetric maps. Here we present the spatial patterns and extend the analysis in P06 with the second cycle (*t* = 157–169 yr). In Fig. 3, the annual-averaged SSH fields of QC05 (color images) are compared with the SSH fields obtained from model results (grayscale images) according to the same synchronization shown in Fig. 2 (note that the spatial scales of data and model maps are the same). It should be borne in mind that we are comparing the observed ocean variability with the results of an idealized model forced by a constant wind stress, so the several mismatches between data and model that are evident from Fig. 3 are the obvious consequence of such a daring comparison.

In 1993 the KE is in the large-meander, elongated state. The corresponding model solution (*t* = 157 yr) agrees well for the first large anticyclonic meander in position, shape, and strength. The second and third meanders are merging in the observations, while they are distinct in the model, but their positions and shapes are again in fairly good agreement. On the other hand, the cyclonic meander south of Japan is weaker in the observations. After 11 yr a similar situation is achieved (year 2004, *t* = 168 yr), and the agreement is now even better for the first and second anticyclonic meanders and for the cyclonic meander south of Japan (which is shifted slightly to the north in the model). For intermediate times, the variability found in model results is in surprising agreement with the observations. The disruption of the small-meander state, accompanied by the disappearance of the two main anticyclonic meanders, occurs in less than 1 yr (from year 1994, *t* = 158 yr, to year 1995, *t* = 159 yr). A further weakening of the KE leads to a very convoluted jet reaching its minimum energy at ∼ year 1997, *t* = 161 yr (see the dashed lines in Figs. 2e,f). From now on the jet is in its recharging phase (see P06, section 4a). The correspondence between the panels in the range years 1998–2004 and *t* = 162–168 yr is amazing, again bearing in mind the simplicity of the model.

One can argue whether the agreement between model and observations is sufficient to claim that P06 can serve as a minimal model: this will always be subjective. Regarding the issues (i) and (ii) mentioned in the introduction dealing with the low-frequency variability of the flow, the model behavior appears to be in significant agreement with observations. Also the irregular variability in the transition between large- and small-meander states (issue iii) is adequately simulated. So with respect to describing the dynamics of the SSH field, the results in this section provide support that all essential processes are contained in the reduced-gravity model (i.e., the latter can serve as a minimal model).

## 3. Analysis of the transition to the Kuroshio Extension relaxation oscillation

From analyses with reduced gravity or barotropic models in relatively small idealized basin geometries (Jiang et al. 1995; Speich et al. 1995) and with more realistic geometries (Dijkstra and Molemaker 1999; Schmeits and Dijkstra 2001), we know that for large values of *A _{H}* there exists a unique steady state. When

*A*is decreased, an imperfect pitchfork bifurcation occurs and a regime of multiple steady states exists. In the Pacific geometry used in Schmeits and Dijkstra (2001), the two stable steady states consist of a small-meander state and a large-meander state. With decreasing

_{H}*A*both steady states become unstable to transient disturbances through Hopf bifurcations, and periodic flow appears.

_{H}In this section, we study systematically the transition from these periodic equilibrium flows (still at relatively large values of *A _{H}*) in the present reduced-gravity model to the value of

*A*= 220 m

_{H}^{2}s

^{−1}, for which the comparison with observations was made in section 2. To do so, following P06, we use an empirical continuation method (e.g., Jiang et al. 1995; Chang et al. 2001), which requires a large number of time-dependent forward integrations for different values of a relevant parameter (here

*A*) and a successive monitoring of the changes in qualitative behavior (to this purpose the model code was parallelized and run on a high-performance computer facility). The transition behavior in similar systems has been systematically studied in the same way for quasigeostrophic models (Nadiga and Luce 2001; Simonnet et al. 2005) and shallow-water models (Speich et al. 1995; Simonnet et al. 2003) of rectangular basin flows.

_{H}### a. Bifurcation diagram

In P06 a preliminary bifurcation analysis was presented (his section 4c) based on a limited number of runs. The behavior was studied in the two-dimensional state space defined by the kinetic energy integrated in the sectors *A* and *B* (the Kuroshio Extension region and the region south of Japan, respectively) shown in Fig. 1. More precisely, the KE kinetic energy per unit mass *E*_{Ω} integrated in each sector Ω (Ω = *A*, *B*) is defined as

where *h* is the active upper-layer thickness and **u** is the velocity vector. Below we substantially extend the analysis of P06 by using many simulations at different values of *A _{H}*.

The time integrations have been performed within the relatively restricted range *A _{H}* = 200–400 m

^{2}s

^{−1}. Based on the solutions computed, the bifurcation diagram shown in Fig. 4 has been derived. For each value of

*A*the curves with labels min and max give the range within which

_{H}*E*varies after spinup. For

_{A}*A*< 230 m

_{H}^{2}s

^{−1}the two pairs of curves

*b*

_{min}–

*b*

_{max}and

*c*

_{min}–

*c*

_{max}are associated with the spontaneous switching of the orbit between two equilibrium flows. The power spectrum of each of the computed time series of

*E*is given as a function of

_{A}*A*and the period in the composite map of Fig. 5 (note the different scales in the three panels).

_{H}For values larger than *A _{H}* = 400 m

^{2}s

^{−1}the flow is steady. Over the whole range

*A*= 255–400 m

_{H}^{2}s

^{−1}, very weak periodic oscillations are present: we will refer to this as regime I. There is one major transition in the interval

*A*= 240–255 m

_{H}^{2}s

^{−1}, where the amplitude of the variability becomes much larger and a local transition to chaos occurs; we will refer to this as regime II. We next have a range

*A*= 235–240 m

_{H}^{2}s

^{−1}in which the most dramatic changes associated with a global bifurcation take place; we will refer to this as regime III. Finally, the large low-frequency variability characteristic for the model behavior as presented in section 2 occurs for viscosity values smaller or equal to

*A*= 235 m

_{H}^{2}s

^{−1}; we will call this regime IV (it is in this regime that the numerical solution shown in Figs. 2, 3 is obtained, corresponding to

*A*= 220 m

_{H}^{2}s

^{−1}). In the next subsections, we discuss each regime in detail.

### b. Regime I (A_{H} = 255–400 m^{2} s^{−1}): Weak periodic oscillations

Just after the first Hopf bifurcation, near *A _{H}* = 350 m

^{2}s

^{−1}, the spectrum of

*E*exhibits two low-frequency peaks: a stronger one at

_{A}*T*= 4 yr and a smaller one at

*T*= 2 yr, while

*E*displays the same peaks (but with smaller amplitude) plus a strong, very high-frequency peak at

_{B}*T*= 5.5 days. For lower viscosity (e.g.,

*A*= 280 m

_{H}^{2}s

^{−1}),

*E*exhibits the same two spectral peaks but with a higher amplitude and with a shift toward lower frequencies:

_{A}*T*= 7 yr and

*T*= 3.5 yr. In addition, at

*A*= 280 m

_{H}^{2}s

^{−1}, a high-frequency peak at 8.5 days is now present as well. For this value of

*A*,

_{H}*E*now has two high-frequency peaks (the same peak present for higher viscosity at 5.5 days plus the new one at 8.5 days) and two low-frequency peaks at the same periods of those of

_{B}*E*, but again with smaller amplitude.

_{A}In Fig. 6a, the spectral amplitudes of *E _{A}* and

*E*(thick and thin line, respectively) for

_{B}*A*= 255 m

_{H}^{2}s

^{−1}are plotted. For both

*E*and

_{A}*E*, the high-frequency peaks are basically unaltered with respect to those of

_{B}*A*= 280 m

_{H}^{2}s

^{−1}, so the associated dynamics appears independent of viscosity. In Fig. 6b the SSH anomaly corresponding to the 8.5-day peak is shown for

*A*= 255 m

_{H}^{2}s

^{−1}, which is obtained by subtracting the field at some

*t*

_{0}from the field at

*t*

_{0}+ Δ

*t*(in this case Δ

*t*=

*T*/2). The pattern is very localized in the western boundary current and south of Japan, where likely the horizontal shear is largest. Both the high-frequency of the variability as well as the pattern lead us to conclude that this is a so-called wall-trapped mode of variability that arises due to the viscous instability of the western boundary layer (Sheremet et al. 1997) and the anomaly tends to propagate with the mean flow, in this case northward. The reason why the high-frequency peak appears in

*E*at any

_{B}*A*while it appears in

_{H}*E*only for a sufficiently small value of

_{A}*A*(as noticed above) is now obvious: the wall-trapped mode (Fig. 6b) has a very limited penetration in region

_{H}*A*(which explains the higher amplitude in

*E*compared with

_{B}*E*), and for very high viscosity such penetration vanishes.

_{A}### c. Regime II (A_{H} = 240–255 m^{2} s^{−1}): Local transitions

For *A _{H}* values smaller than 255 m

^{2}s

^{−1}the range of

*E*experiences a very rapid increase, as shown in the diagram of Fig. 4. The flow behavior in the initial part of this interval (

_{A}*A*= 250–255 m

_{H}^{2}s

^{−1}) is illustrated by the time series of Fig. 7 (showing

*E*and

_{A}*E*in the left and right panels, respectively, for 100 yr of integration) and in the spectra of the two extreme values of this small

_{B}*A*interval in Fig. 6a and Fig. 8a, respectively. Apart from a broader energy distribution in the frequency domain, the most relevant feature is an increase by an order of magnitude of the spectral amplitudes of the two main low-frequency peaks, which undergo also a shift toward smaller frequencies (from 3.8 and 8 yr for

_{H}*A*= 255 m

_{H}^{2}s

^{−1}to 6 and 8.2 yr for

*A*= 250 m

_{H}^{2}s

^{−1}; see also Fig. 5).

The mean flow for *A _{H}* = 250 m

^{2}s

^{−1}is shown in Fig. 8b: it is very similar to that of

*t*= 163 yr of Fig. 3, corresponding to the minimum energy value of

*E*just at the beginning of the recharging phase of the relaxation oscillation (see the dashed line of Fig. 2e or 2f). This same feature is found for any value of

_{A}*A*in regimes II, III, and IV, and should be considered a very distinct character of our low-frequency variability: no matter what the nature of the oscillation is (whether an extremely small oscillation, a periodic oscillation, or a vigorous relaxation-type oscillation), the lowest energy state has always approximately the same structure and energy (for the latter compare lines

_{H}*b*

_{min}and line

*a*

_{min}for

*A*< 255 m

_{H}^{2}s

^{−1}). Figure 8c shows the SSH anomaly of the low-frequency variability about the mean flow of Fig. 8b. Unlike the high-frequency anomaly, a much broader scale pattern is evident in this case. Because of the large-scale pattern and its time scale, we interpret this oscillation as being caused by a gyre mode (Simonnet and Dijkstra 2002; Simonnet 2005).

The flow behavior in the interval *A _{H}* = 240–250 m

^{2}s

^{−1}displays a further remarkable increase of the energy range (Fig. 4) and a complete saturation of the periods of the two main low-frequency peaks (Fig. 5), which are finally located at 6.4 and 9 yr, respectively. It is interesting to notice (Fig. 5) that the periods of these two peaks remain relatively constant over a wide range of viscosity values (

*A*= 255 ∼ 280 m

_{H}^{2}s

^{−1}), while they increase and then saturate within the same range (

*A*= 240 ∼ 255 m

_{H}^{2}s

^{−1}) in which the variability of

*E*has a remarkable increase (Fig. 4).

_{A}### d. Regime III (A_{H} = 235–240 m^{2} s^{−1}): Global transitions

In this small interval of *A _{H}* values the system undergoes an impressive change in its behavior, as shown by the time series in Fig. 9. There is a transition from the weak amplitude irregular oscillations at

*A*= 240 m

_{H}^{2}s

^{−1}(corresponding to gyre-mode variability around an unstable steady state very similar to that shown in Fig. 8b) to a much larger-amplitude relaxation-type oscillation at

*A*= 235 m

_{H}^{2}s

^{−1}. The transition is extremely complex, and the six cases reported in Fig. 9 are indicative in this respect. Already for

*A*= 240 m

_{H}^{2}s

^{−1}a single relaxation oscillation appears just after

*t*= 20 yr, while the subsequent oscillations are irregular but do not represent bimodal flows. For

*A*= 239 m

_{H}^{2}s

^{−1}the situation is similar, but the local oscillations are quasi-periodic until

*t*∼ 120 yr, and then become irregular thereafter. For

*A*= 238 m

_{H}^{2}s

^{−1}and

*A*= 237 m

_{H}^{2}s

^{−1}, relaxation oscillations switch to local oscillations and vice versa in an irregular manner. For

*A*= 236 m

_{H}^{2}s

^{−1}a perfectly quasi-periodic local oscillation appears! The presence, in a single run, of both relaxation oscillations and local oscillations for

*A*= 236–240 m

_{H}^{2}s

^{−1}is represented in the bifurcation diagram of Fig. 4 by two different ranges (within the vertical dashed lines) delimited by

*a*

_{min}–

*a*

_{max}and

*b*

_{min}–

*b*

_{max}, respectively. This has its counterpart in the patterns of the spectral amplitude of Fig. 5.

Following P06, here we interpret the transition to the large-amplitude relaxation oscillation as a global bifurcation, whose consequence is to allow the trajectory in state space to leave the region close to the unstable steady state through its unstable manifold, and to move in a wide region of state space through what is usually denoted as a homoclinic orbit, which eventually goes back to the original region through the stable manifold of the steady state. We will analyze this transition in more detail in section 4.

### e. Regime IV (A_{H} = 200–235 m^{2} s^{−1}): Homoclinic and heteroclinic orbits

For flows in this last regime, time series of *E _{A}* and

*E*are plotted in Fig. 10. The relaxation oscillation presented in section 2 falls in this range and corresponds to the case

_{B}*A*= 220 m

_{H}^{2}s

^{−1}(a more detailed time series of

*E*for this case is also given by the dashed curves in Figs. 2c and 2e for two successive cycles). Figure 10 shows that this relaxation oscillation is a very robust feature, because it is present, with minor changes in its structure, for any viscosity value within the interval

_{A}*A*= 235–215 m

_{H}^{2}s

^{−1}.

However, for some values of *A _{H}* (e.g., for 230, 225, and 215 m

^{2}s

^{−1}in Fig. 10) the system is attracted by small-amplitude intermediate energy oscillations of the KE (see

*E*) accompanied by a strong anticyclonic recirculation gyre south of Japan (see

_{A}*E*). This happens, for instance, during the interval 130 yr <

_{B}*t*< 175 yr for

*A*= 230 m

_{H}^{2}s

^{−1}, after which the system is attracted back to the relaxation oscillation. The spatial pattern of the flow in this state resembles that at

*t*= 157 yr of Fig. 3, and, in fact, at that time both

*E*and

_{A}*E*for

_{B}*A*= 220 m

_{H}^{2}s

^{−1}show that the system has reached the region of state space associated with the small-amplitude oscillation, although it escapes from it very soon. Following P06 we interpret this behavior as a heteroclinic connection between those two attracting regions of phase space, which will be further analyzed in section 4.

The complexity of system behavior in this highly nonlinear regime is well represented by the case *A _{H}* = 210 m

^{2}s

^{−1}. After two typical KE relaxation oscillations and two further anomalously strong oscillations during the first 50 yr, the system is attracted to the small-amplitude oscillation state discussed in the previous paragraph. For this low-viscosity value, however, the system exhibits a perfectly repeating cycle (of 9.5 yr), and it consequently stays there forever. The same happens for smaller viscosity values, for example,

*A*= 200 m

_{H}^{2}s

^{−1}(see Fig. 7 of P06), in which case the final state is reached after

*t*∼ 150 yr. For even smaller viscosity the model (with the adopted spatial resolution) is numerically unstable.

## 4. Analysis of the relaxation oscillation

An analysis of the KE bimodal behavior as a self-sustained oscillation of intrinsic oceanic origin was proposed by P06 (his section 4a) on the basis of geophysical fluid dynamical considerations. In this section, following the same approach as used in section 3, an analysis in the language of nonlinear dynamical systems theory will be carried out. Although the time series discussed in the previous section provide an impression of the different flow regimes, they offer no direct insight into the geometric properties of the trajectories in state space (e.g., whether homoclinic and heteroclinic orbits actually occur and whether a particular flow is sensitive to initial conditions). We now develop a more geometric analysis by using projections of the trajectories in the state space represented by the *E _{B}*–

*E*plane, following the same choice of P06. In such a plane the system evolution is represented by the trajectories:

_{A}where *t _{k}* =

*k dt*(the sampling time of the time series is chosen as

*dt*= 1 day) and, for the sake of convenience, the energies are scaled by a factor 10

^{13}.

We first focus on the global bifurcation occurring between *A _{H}* = 235 m

^{2}s

^{−1}and

*A*= 240 m

_{H}^{2}s

^{−1}(regime III). In Fig. 11a, the abrupt transition from a small-amplitude oscillation to a large-amplitude relaxation oscillation is evidenced again in the time series. In the

*E*–

_{B}*E*phase plane (Fig. 11b), the view is even more dramatic; for

_{A}*A*= 240 m

_{H}^{2}s

^{−1}(blue curve) the trajectory occupies a relatively small area in state space while for

*A*= 235 m

_{H}^{2}s

^{−1}the trajectory suddenly explores a new high-energy area in state space. This change in state space behavior can be illustrated in more detail by considering the probability density function (PDF), which is computed as follows.

Let us subdivide a portion Γ of the plane (containing all trajectory points after spinup) in *n* × *m* rectangular cells *σ _{ij}* (

*i*= 1, … ,

*n*;

*j*= 1, … ,

*m*), each of area

*dσ*. The PDF

*P*of localization of the trajectory in the cell

_{ij}*σ*is defined as

_{ij}where *t*_{kmin} = *k*_{min }*dt* is the time at which the trajectory is first localized over the attractor, and where

Naturally, the normalization condition

is satisfied.

Figure 11c shows that the PDF (with *n* = *m* = 100) is confined in a restricted region around the unstable fixed point **X** for *A _{H}* = 240 m

^{2}s

^{−1}, and then collapses for

*A*= 235 m

_{H}^{2}s

^{−1}, embracing regions well beyond the original basin of attraction of

**X**. Such a small change in the viscosity has produced a small but extremely relevant change in the structure of state space:

**X**has now become an unstable saddle point so that the trajectory can now depart from it.

Figure 12a shows the PDF for the reference case *A _{H}* = 220 m

^{2}s

^{−1}(section 2b), again with

*n*=

*m*= 100, and Fig. 12b shows the corresponding contour plot limited to very low values of

*P*in order to evidence the main paths followed by the trajectory. In Fig. 13 the state space velocity field

**U**

*obtained from (2) in each cell*

_{ij}*σ*according to the formula

_{ij}is plotted. The point where the PDF is at a maximum corresponds to the one in which |**U**| attains its minimum value, and can therefore be identified with the unstable saddle point, **X** = (2.955, 0.356) (indicated by a dot in Figs. 12b and 13), whose corresponding spatial pattern is very similar to that for *t* = 161 yr of Fig. 3 (note that in **X** the KE energy *E _{A}* is minimum while the energy of the Kuroshio south of Japan

*E*has an intermediate value). A combined analysis of these figures shows that the trajectory has a long residence time near

_{B}**X**, it then leaves that region following the main direction indicated by the outgoing arrow in Fig. 12b and remains almost confined within the region denoted Π

_{out}, which is a portion of the unstable manifold of

**X**; the trajectory eventually approaches

**X**through the main direction indicated by the ingoing arrow, remaining almost confined within the region denoted Π

_{in}, which is a portion of the stable manifold of

**X**. It should be noted that Π

_{out}and Π

_{in}appear disconnected, but this is only because the homoclinic reconnection occurs in a wide region, roughly represented by the sector Π

_{rec1}∪ Π

_{rec2}in Fig. 13, in which very high speeds |

**U**| are accompanied by very small values of the PDF. In the same figure the sector Π

*delimited by the 2.0-isoline of*

_{m}*P*appears to be well representative of the region inside which the unstable and stable manifolds depart from the high residence time region surrounding

**X**.

Very useful information concerning the sensitivity of flow trajectories to initial conditions (and hence, concerning the degree of predictability of the flow) is provided by the rate of divergence of initially nearby points. We investigate this by computing a field of finite time Lyapunov exponents (FTLE) generalizing a technique typically used with experimental time series (e.g., Hilborn 2000, chapter 9). From the trajectory **X*** _{k}*,

*M*time instants

_{ij}*k*(

*i*,

*j*,

*ν*) (

*ν*= 1, … ,

*M*) are identified such that the points

_{ij}**X**

_{k(i,j,ν)}belong to distinct tracts of the trajectory and are contained in

*σ*(two different

_{ij}*k*are required to be separated by a time sufficiently large to ensure that the corresponding points belong to different relaxation oscillations). Then, the

*M*− 1 couples of nearby points (both falling initially inside

_{ij}*σ*)

_{ij}**X**

_{k(i,j,1)}and

**X**

_{k(i,j,1+μ)}(

*μ*= 1, … ,

*M*− 1) are chosen to compute the divergence of trajectories. Let us call

_{ij}*δ*

_{0}the initial distance (defined by the Euclidean metric) between two points, and

*δ*its evolution after time

_{τ}*τ*=

*dt k*′:

We define doubling time as the time *τ*_{2(i,j,μ)} (assuming it exists under the hypothesis of exponential divergence of nearby trajectories) after which *δ _{τ}* = 2

*δ*

_{0}. Then we can define an average doubling time relative to each cell:

It is customary to represent this information by introducing the Lyapunov exponent *λ*:

We can thus compute from (5) the finite-time (doubling time) Lyapunov exponent for each cell:

Finally, in order to characterize the trajectory **X*** _{k}* in terms of the degree of predictability of the state space regions it encounters, we propose to use a “Lagrangian” Lyapunov exponent following

**X**

*, defined as*

_{k}where, again, *t _{k}* =

*k dt*.

It should be borne in mind that *λ*_{2} and Λ** _{X}** are computed on the

*E*–

_{B}*E*phase plane, so the predictability implied by such Lyapunov exponents concerns only the energy of the system, not its actual state. In principle, cases might exist in which the evolution of the energy of different spatial sectors of the system could be very predictable and, at the same time, the flow evolution could be very sensitive to the initial conditions (see the discussion of Fig. 15 below), but the contrary would not be allowed. In practice, since the bimodality of the KE is described so well by

_{A}*E*and

_{A}*E*,

_{B}*λ*

_{2}and Λ

**are expected to provide significant information about the possibility of predicting the evolution of the KE state, at least in its gross features.**

_{X}It is also worth stressing that for the computation of the PDF, the state space velocity, the field of FTLE, and the Lagrangian FTLE defined by (3), (4), (6), and (7), respectively, trajectories have to be computed, and this is implied by the empirical continuation method adopted here. In this respect, the latter method presents, therefore, an important advantage over the more sophisticated analytical continuation methods, which, however, are not suitable for this purpose.

Figure 14 shows *λ*_{2} for *A _{H}* = 220 m

^{2}s

^{−1}, computed from a 400-yr-long time series. In Π

*,*

_{m}*λ*

_{2}is very small everywhere: this implies a slow divergence of trajectories and, consequently, a relatively high degree of predictability (in the sense stated above). On the other hand, it is in the reconnection regions Π

_{rec1}and Π

_{rec2}and in Π

*that*

_{a}*λ*

_{2}assumes very large values, implying a highly unpredictable state. Figure 15c shows Λ

**for the reference period**

_{X}*t*= 145–169 yr (corresponding to the two bimodal cycles discussed in section 2b), where

*E*,

_{A}*E*(Fig. 15a),

_{B}*L*

_{KE}, and (Fig. 15b) are also shown for the sake of comparison. Figure 15 allows us to investigate the variation of

*λ*

_{2}and to identify three main stages (S1, S2, and S3) of the relaxation oscillations.

In stage S1, during the intervals *t* ≈ 147–151 yr and *t* ≈ 159.5–165 yr (for the first and second cycles, respectively) Λ** _{X}** is very small, so

*E*and

_{A}*E*yield the highest predictability. Here we are in the first half of the transition phase (see section 2b), that is, just after the disruption of the large-meander phase and, subsequently, in the recharging phase of the relaxation oscillation (Fig. 15a, thick line) when the KE patterns are more variable (Fig. 3) and more convoluted (Fig. 15b, thick line) and is increasing (Fig. 15b, thin line). In other terms, the energy evolution proves less chaotic when the variability and convolution of flow patterns are maximum (this occurs in region Π

_{B}*of state space). Such property is somewhat unexpected, but is only apparently paradoxical [see the discussion just after Eq. (7)] and appears to be very distinctive of our KE theory. The concomitant evolution of*

_{m}*E*and

_{A}*E*in this stage of the RO differs little from cycle to cycle, but this hides the existence of highly variable mesoscale features (e.g., as shown by

_{B}*L*

_{KE}in Fig. 15b) that feed the KE jet and that, at the same time, are presumably very unpredictable.

In stage S2, during the intervals *t* ≈ 151–154 yr and *t* ≈ 165–168 yr (for the first and second cycles, respectively) Λ** _{X}** is highly variable, with peaks an order of magnitude larger than in the previous phase, so

*E*and

_{A}*E*yield a highly chaotic behavior and are, therefore, highly unpredictable (

_{B}*λ*

_{2}= 100 yr

^{−1}corresponds to a doubling time of just

*τ*

_{2}≈ 2.5 days). This occurs when the KE is reaching its maximum energy, the large-meander state is attained, is at its maximum, and

*L*

_{KE}has a weak variability. Moreover, the Kuroshio south of Japan yields a strong variability (Fig. 15a, thin line) that, according to P06 (see also Pierini 2008), is associated with the ejection of cyclonic mesoscale eddies detaching from a meridionally elongated cyclonic meander (a process well known to occur south of Japan; e.g., Qiu and Miao 2000). So, it is basically in connection with this process (confined outside the KE region and occurring mainly in Π

_{rec2}and in part of Π

_{rec1}) that the large values of Λ

**emerge.**

_{X}In stage S3, during the interval *t* ≈ 154–159.5 yr (referring to the first cycle; the corresponding interval for the second cycle is not included in the graph) Λ** _{X}** is less variable than in stage S2 but with several very large peaks (corresponding to

*τ*

_{2}≈ 1 day), implying a very unpredictable system. This occurs mainly in sector Π

*(Figs. 13, 14) when the KE is in the large-meander–elongated mode but is weakening, not being fed by the flux of negative relative vorticity coming from the south, which is dissipated within the region of the southern cyclonic meander (see P06, his section 4a).*

_{a}In conclusion, the concomitant evolution of *E _{A}* and

*E*is less chaotic when the KE jet is weak, more variable, and more convoluted, and in its recharging phase (stage S1). It is very chaotic when the KE jet is in its large-meander phase (stage S2), and also when it weakens, as long as it is in the large-meander state (stage S3), in this last case presenting distinct episodes of very large FTLE. These conclusions provide valuable insight into the predictability of the low-frequency variability of the KE system.

_{B}## 5. Summary and conclusions

In this paper we have proposed a theory for the KE bimodal low-frequency variability that is equivalent barotropic and intrinsically nonlinear. The theory is based on the transition behavior of double-gyre flows in QG models that is extended here for the Pierini (2006) model, and shows remarkably good agreement with observations. The control parameter in the transition behavior of the latter model is the lateral eddy viscosity parameter *A _{H}*. On the other hand, the amplitude of the wind stress forcing (another typical control parameter) was kept fixed. This choice is justified by the nature of our approach, which is based on an idealized process study but with essential elements of realism, among them a relatively realistic climatological wind stress field whose amplitude is therefore fixed. The time integrations have been performed for many values of

*A*within the relatively restricted range

_{H}*A*= 200–400 m

_{H}^{2}s

^{−1}. The investigation carried out in this paper on how the solution undergoes (even dramatic) changes in such a range has contributed to improving our understanding of the character of the KE relaxation oscillation.

For large values of *A _{H}* there is only a steady state [in agreement with results in Schmeits and Dijkstra (2001)], which is very similar to the small-meander KE state. Once

*A*is decreased, a first Hopf bifurcation occurs, leading to periodic behavior yielding high-frequency (on the order of days) and low-frequency (on the order of years) oscillations. Similar high-frequency oscillations have been found in QG models of the double-gyre flow in idealized basins and they could be attributed to either Rossby–basin modes or wall-trapped modes: here we have identified the latter. The low-frequency (interannual) variability has been attributed, like in simpler models, to a gyre-mode (shear) instability. For smaller

_{H}*A*, we have then provided evidence for the occurrence of a global bifurcation and the establishment of a quasi-homoclinic regime, in which the frequency of the variability saturates. We interpret this as a Shilnikov bifurcation (Simonnet et al. 2005), occurring in these models because the periodic behavior associated with gyre-mode variability becomes so large that it connects to a second steady state.

_{H}This transition scenario does not only provide insight into the origin of the complex behavior of the flow trajectories, it also provides a picture of the physical processes associated with time scales and patterns in the flow. Definitely, (at least) two steady state patterns of the KE play a crucial role in the bimodality of the jet; the two states differ considerably in their spatial structure and energy, with the small-meander state having kinetic energy that is about one-third of that of the large-meander state. The transition time scale between these patterns originates from a saturation of the gyre-mode instability (as it does in QG models) and hence from a saturation of the shear instability of the small-meander state. The high degree of irregularity that the flow has during the transition from the small- to the large-meander state is due to the fact that the high-energy state is only reached in a global bifurcation, which is accompanied with chaotic behavior.

In the analysis we have made use of mathematical tools whose application to the specific problem appears to be novel. In an appropriate two-dimensional phase space the probability density function, the phase space velocity, and a field of finite-time Lyapunov exponents and the corresponding Lagrangian time series have been computed from the long time series that have benefited from the empirical continuation method adopted here. We believe such tools could be applied to studies of the same nature in the future.

In view of these results we conclude that the KE bimodality is likely to arise from highly nonlinear interactions that are controlled by an equivalent-barotropic dynamics internal to the ocean system, and in which time-dependent atmospheric forcing, baroclinic instability, and the effects of bathymetry are of second-order importance.

## Acknowledgments

This research was supported by the “High Performance Computing - Europa” Project of the European Commission (Contract RII3-CT-2003-506079), under the VI Framework Program: “Support to Research Infrastructure Action.” A parallel version of the model code was run on a high-performance computing facility provided by SARA - Computing and Networking Services (Amsterdam). SP has conducted part of the research at the Institute for Marine and Atmospheric Research Utrecht (Utrecht University), whose hospitality is kindly acknowledged.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Stefano Pierini, Dipartimento di Scienze per l’Ambiente, Università di Napoli Parthenope, Centro Direzionale - Isola C4, 80143 Naples, Italy. Email: stefano.pierini@uniparthenope.it