## Abstract

Intrinsic low-frequency variability is studied in the idealized, quasigeostrophic, midlatitude, wind-driven ocean gyres operating at large Reynolds number. A robust decadal variability mode driven by the transient mesoscale eddies is found and analyzed. The variability is a turbulent phenomenon, which is driven by the competition between the eddy rectification process and the potential vorticity anomalies induced by changes of the intergyre transport.

## 1. Introduction

We pose the problem, discuss the background, and describe the ocean model, in this introduction. The low-frequency variability mode and the role of the ocean eddies in maintaining the mean flow are discussed in sections 2 and 3, respectively. Section 4 focuses on the dynamics of the mode, and the conclusions follow in section 5.

### a. Statement of the problem

The ocean variability may be classified into three major categories: 1) the directly *forced* variability resulting from the slaved response of the ocean to time-dependent fluctuations of the atmosphere, 2) the *intrinsic* ocean-only variability, and 3) the *coupled* variability of the combined ocean–atmosphere (or ocean–atmosphere–ice) system. The goals of this paper are to identify large-scale, low-frequency, intrinsic variability of the ocean circulation and to explain it in terms of the fluid dynamics. We use an ocean model that is rather idealized in terms of the basin configuration, stratification, and surface forcing, but its dynamics is highly nonlinear and even turbulent. An important point is that the low-frequency variability is driven by dynamically resolved, rather than parameterized, mesoscale oceanic eddies. The atmospheric forcing is approximated by a stationary wind, and buoyancy forcing is excluded. Thus, the novelty of this study is in the combination of the *longtime* ocean model integrations, necessary for capturing the low-frequency variability, and in the explicit representation of the *mesoscale eddies*.

The main steps of this study are the following. An eddy-resolving solution of the dynamical ocean model is obtained in the statistically equilibrated regime. The dominant low-frequency variability is diagnosed by decomposing fluctuations of the flow into empirical orthogonal functions with time-dependent principal components. The dynamical role of the eddies in maintaining the mean flow configuration is revealed. The low-frequency mode is studied by calculating relevant potential vorticity (PV) balances at different phases of the variability cycle. Variability of the balances is understood in terms of the eddy–mean flow interactions. Some arguments of the Eulerian analysis are supported by estimates of the Lagrangian transport.

### b. Background

Although significant evidence of the decadal and longer variability of the ocean is found in the midlatitude sea surface temperature (e.g., Kushnir 1994; Hansen and Bezdek 1996), hydrography (Qiu and Joyce 1992), and altimetry (Qiu 2003; Qiu and Chen 2005) data, the underlying dynamical mechanisms of this variability remain to be established. Present theoretical understanding is based on the null hypothesis: that the ocean integrates high-frequency atmospheric forcing and responds in terms of the red variability spectrum (Hasselman 1976). However, the possibility that intrinsic ocean dynamics may generate variability has not been ruled out, and the role of mesoscale oceanic eddies in driving such intrinsic variability remains poorly understood. This variability had been first noted by Bryan (1963), and the fundamental role of the mesoscale eddies in shaping the background flow was later emphasized by Holland (1978). Since then, different variability patterns and mechanisms have been suggested on the basis of idealized gyre dynamics.

One of the intensively studied ideas is that the ocean variability can be understood in terms of early bifurcations of the forced and dissipative dynamical systems (Cessi and Ierley 1995; Ierley and Sheremet 1995; Speich et al. 1995; Jiang et al. 1995; Sheremet et al. 1997; Dijkstra and Katsman 1997; Meacham and Berloff 1997a, b; Berloff and Meacham 1997, 1998; Meacham 2000; Nadiga and Luce 2001; Primeau 2002; Simonnet et al. 2003a). In this approach, it is argued that the ocean is steered by the underlying, unstable low-dimensional attractors such as fixed points (steady), limit cycles (periodic), tori (quasiperiodic), strange attractors (chaotic states), and homoclinic orbits. This argument is based on the analysis of early bifurcations of ocean dynamics in the vicinity of steady states. Bifurcations of the attractors are found with respect to varying control parameters, such as the eddy viscosity coefficient, of the chosen ocean model.

A significant mode of ocean variability that has been discovered using the dynamical systems approach has been called the *gyre mode* (Simonnet and Dijkstra 2002). It has been found in a number of different models with a double gyre circulation (e.g., Nauw and Dijkstra 2001; Simonnet and Dijkstra 2002; Simonnet et al. 2003a, b; Nauw et al. 2004; Simonnet 2005) and across a significant range in parameter space. The gyre mode has been identified as a homoclinic orbit arising from the merger of two stationary modes (Simonnet and Dijkstra 2002).

The most uncertain component of *all* of the models analyzed with the dynamical systems language is the parameterization of the mesoscale eddy effects by diffusion.^{1} There would be no need to resolve the eddies if the diffusive parameterization were accurate, but generally this is not the case. For example, it is well established that eddies act to sharpen (rather than diffuse) atmospheric jets through up-gradient eddy momentum flux (McIntyre 1970). There are numerous examples of such *negative viscosity*, or antidiffusive phenomena in geophysical turbulence (e.g., Starr 1968; McIntyre 2000) and, more specifically, in the ocean gyres (Berloff 2005a). Thus, the relevance of the low-dimensional attractors and, thus, the usefulness of the dynamical systems language has yet to be established in turbulent regimes with more realistic ocean stratification. To our knowledge, only one study (Berloff and McWilliams 1999a) has attempted to verify, without success, the relevance of low-dimensional attractors. This raises the question—still open—of whether the gyre mode, as well as other modes of variability found in simple ocean models, are relevant to the turbulent ocean circulation. Clearly, further study is needed. Among the most recent developments, on the one hand, it is argued that, despite some spatial pattern similarities between the gyre mode and the turbulent-ocean variability, the latter operates on a significantly longer time scale (Hogg et al. 2005). On the other hand, it is conjectured that eddy effects can be parameterized by Rayleigh friction in a barotropic model (Simonnet 2005) that improves representation of the eastward jet and increases the time scale of the gyre mode to the decadal scale.

Different types of variability modes may be excited if some idealizations of the gyre models are relaxed. In particular, some studies suggest that interaction of the deep western boundary current and the wind-driven gyres is a cause for the low-frequency variability (Spall 1996; Katsman et al. 2001). Others argue that the variability is that of the thermohaline circulation (e.g., Weaver et al. 1993). Last, on the basis of a comprehensive eddy-permitting simulation, Qiu and Miao (2000) suggest that interannual variability of the Kuroshio corresponds to the cycle of gradual strengthening of the adjacent inertial recirculation interrupted by its episodic instability and subsequent intrusions of the opposite-sign potential vorticity.

In summary, past studies have not yielded a universally accepted view as to the origins and dynamics of gyre-scale low-frequency variability.

### c. Ocean model

The idealized ocean model employed in this study is a fairly standard *double gyre* paradigm for the midlatitude circulation (e.g., Holland 1978). In this study the physical parameters of the model, which specify geometry of the ocean basin, stratification, planetary vorticity, and external forcing, are kept close to those used in the common practice. The exception is the eddy viscosity coefficient, which is specified to be relatively small. The eddy viscosity is the main free parameter of the problem,^{2} and is meant to represent effects of the dynamically unresolved oceanic eddies. When combined with the relatively fine spatial resolution of the model, the low eddy viscosity guarantees that the mesoscale eddies are dynamically resolved to an almost unprecedented level. Finally, the double-gyre model is used for very long, multicentury integrations that are required for the accurate description of the intrinsic low-frequency variability.

The dynamics of the model used in the present study represents that of the midlatitude ocean, with a prescribed density stratification, in a flat-bottom square basin with north–south and east–west boundaries and with size *L* = 3840 km. The quasigeostrophic PV equations (Pedlosky 1987) for three dynamically active isopycnal layers are

where the layer index starts from the top; *J*() is the Jacobian operator; *ρ*_{1} = 10^{3} kg m^{−3} is the upper-layer density; *β* = 2 × 10^{−11} m^{−1} s^{−1} is the planetary vorticity gradient; and *ν* = 100 m^{2} s^{−1} is the eddy diffusivity coefficient. The isopycnal layer depths are *H*_{1} = 250, *H*_{2} = 750, and *H*_{3} = 3000 m. The ocean model’s PV anomalies, *q _{i}*, are related to the velocity streamfunctions,

*ψ*, through the elliptic subproblem:

_{i}and

where the stratification parameters are

Here *g*′_{1} and *g*′_{2} are reduced gravities associated with the density jumps across the upper and lower internal interfaces between the isopycnal layers, and *f*_{0} is the Coriolis parameter. Stratification parameters are chosen so that the first, *Rd*_{1}, and the second, *Rd*_{2}, Rossby deformation radii are 46 and 27 km, respectively.

The wind stress curl *W* is given by

if *W*_{0}(*y*) is larger than *W*_{min} = −3.7*τ*_{0}/*L*, and otherwise: *W* = *W*_{min}. The forcing has a double-gyre structure and it is chosen to be asymmetric with respect to the middlelatitude of the basin in order to avoid the artificial symmetrization of the gyres in the turbulent regime, which impacts low-frequency variability (Berloff and McWilliams 1999a). The wind stress amplitude is *τ*_{0} = 1.0 N m^{−2}. The model operates at a large Reynolds number,

where *U* = *τ*_{0} (*ρ*_{1 }*H*_{1 }*L β*)^{−1} is the upper-ocean Sverdrup velocity scale (about 5.0 cm s^{−1}).

The partial-slip lateral boundary conditions,

where *α* = 15 km and the partial derivatives are taken in the normal direction to the boundary, are imposed for each isopycnal layer. If one Taylor expands *ψ _{i}* near the wall, then contributions of the linear and quadratic terms become comparable over the length scale

*α*, so it can be interpreted as a boundary sublayer. Solution sensitivity to

*α*is studied in section 4c, and a more detailed study at smaller Re is in Haidvogel et al. (1992).

The mass conservation constraints,

are enforced as well.

The governing equations (1)–(6) are discretized with the second-order finite differences, including PV-conserving formulation of the Jacobians. The prognostic equations (1)–(3) are then marched in time with the leapfrog scheme and 30-min time step, and the elliptic problem (4)–(6) is solved for the streamfunctions at each time step. The horizontal grid resolution is uniform, with 513 × 513 grid points and 7.5-km intervals between them, so that the Munk length scale, *δ _{M}* = (

*ν*/

*β*)

^{1/3}, is resolved by more than two grid points. This grid size is close to the

^{1}/

_{12}° resolution at which, it is argued, the mesoscale eddies are substantially resolved (e.g., Bleck et al. 1995). After the initial spinup from the state of rest, the reference solution equilibrates statistically. It is then computed for another 5 × 10

^{4}days and stored for the analysis.

## 2. Low-frequency variability of the gyres

The time-mean circulation has the classical double-gyre structure (Figs. 1a,b). Instantaneous flow patterns are characterized by coexisting types of eddies: meanders and vortices around the western boundary current (WBC) and its eastward jet (EJ) extension, and planetary waves in the rest of the basin (Figs. 1c,d). The wind stress pattern is shown in Fig. 1e.

Figures 1c,d show that the mesoscale eddy field is well developed. To discern low-frequency variability in the large-scale circulation we need to separate it from the eddy field: this is achieved through temporal filtering of the velocity streamfunction. For the purposes of this paper we define the large-scale flow, 〈*ψ*〉, as the filtered field, and the “eddy” field, *ψ*′, as the flow component eliminated by this filtering. We thus write,

where layer subscripts have been omitted for simplicity. The filters employed in this work are the time mean, a conditional ensemble mean, and a filter with a cutoff frequency based on eddy time scales (500 days in the standard case). The choice as to which one is used is predicated on the particular question at hand and will be indicated in the sections to follow.

To detect the large-scale, low-frequency variability of the gyres, we use the empirical orthogonal function (EOF) pattern recognition method (Preisendorfer 1988). The method is applied to the reference solution, which is low-pass filtered in time by the 500-day running-window averaging and sampled on the uniform 65 × 65 spatial grid (this grid is sufficient to resolve details of the large-scale flow). The 10 leading EOFs are found as eigenvectors of the covariance matrix of the upper-ocean streamfunction fluctuations, and the corresponding principal components (PCs) are obtained by projecting the reference solution time series on the eigenvectors. The leading EOF pair (Figs. 2a,b) contains about ¾ of the total variance of the low-passed flow and represents a coherent pattern with the sense of north–south propagation. The first EOF, which has a tripolar structure, captures coherent meridional shifts of the EJ, and the second EOF, which has a dipolar structure, describes variations of the intensity of the EJ. The tripole is asymmetric, which indicates that, when the EJ shifts to the most northern state, it also reaches the maximum of its intensity. The corresponding PCs clearly indicate that the associated variability is decadal (Figs. 2c,d). The time scale of propagation is pinned down by the time-lag correlation function of the first and second PCs (Fig. 2e). This function is roughly antisymmetric around zero, and the maximum correlation is about 0.7 at ±3 yr (i.e., at about one-quarter of the dominant period)—such a large correlation magnitude suggests that the signal is nearly periodic. The time-lag correlation function also indicates that the mode is weakly asymmetric in time: transition from the weak- to strong-jet regime is about 10% longer than the opposite transition. The mode’s time scale is consistent with the idea of the nonlinear adjustment of the combined EJ–eddies system and is unrelated to the basin size.^{3}

To summarize, the gyres exhibit robust large-scale, 12-yr variability concentrated around the EJ. We believe that a dynamically equivalent, but less pronounced, mode was reported in Berloff and McWilliams (1999a) for a somewhat different model setting and for a lower Reynolds number (Re). A similar mode is found for a different set of parameters and for somewhat different formulation of the model (Hogg et al. 2005). It is also plausible that the mode reported in Qiu and Miao (2000) is of a similar kind. All of the above suggests that the found mode is robust and generic. It should be noted that the gyre mode (Simonnet and Dijkstra 2002) possesses similar spatial characteristics, but different temporal characteristics, and does not require explicit turbulence for its existence. On the other hand, more recent results (Simonnet 2005) suggest that with the eddy effects parameterized by Rayleigh friction in a barotropic model the time scale of the gyre mode can become decadal. It is unclear from the literature whether the dynamics of the gyre mode and the present mode are equivalent; the examination of this issue is an objective of this paper.

To understand the underlying dynamics, we define and obtain the “key states” of the low-frequency mode. States A and C are given by the corresponding ensemble of flow states in which the first PC normalized by its variance (Fig. 2c) is larger than +1 or less than −1, respectively. States B and D are defined similarly but key on the amplitude of the normalized second PC (Fig. 2d). Overall, the flow is in one of these states 60% of the time, and it spends roughly equal intervals in each state. In each isopycnal layer and for each of the key states, the velocity streamfunction can be represented as a sum of the conditional ensemble-mean^{4} and fluctuation components: *ψ ^{K}* = 〈

*ψ*〉

*+*

^{K}*ψ*′

*, where 〈···〉*

^{K}*denotes ensemble averaging and*

^{K}*K*stands for the corresponding key state. The key states are shown in Fig. 3. Representing the low-frequency variability in terms of transitions between just a few states is very useful because it is extremely simple. More detailed, and also more complicated, representation of the variability would involve continuous separation into the large-scale and eddy components by a more elaborate filtering.

The differences between the key state patterns are significant. The amplitude of the meridional shift of the EJ is about 200 km, and the amplitude of upper-ocean transport variation in the eastern half of the EJ is more than 20 Sv (Sv ≡ 10^{6} m^{3} s^{−1}). Such variability, if it exists in the real ocean, is likely to have a significant impact on the midlatitude climate. Furthermore, Hogg et al. (2006) suggest that such a mode can couple with the atmospheric variability. Another important point is that the subtropical gyre is completely embedded in the negative wind forcing region, but over the cycle it occupies different fractions of it; thus net wind forcing of the gyre changes noticeably. The PV anomalies associated with the key states are predominantly negative in states A and B, and positive in states C and D (Fig. 4), suggesting that dynamics of the mode can be understood in terms of variability of different PV inputs.^{5}

The evolution of the low-frequency mode can be described in terms of the sequence of transitions: A → B → C → D → A. In state A, the EJ is both strong and shifted to the north (i.e., to high latitudes) relative to the time-mean state. Transition to state B is characterized by the southward (i.e., to low latitudes) migration of the EJ, which remains relatively strong. Transition to state C is characterized by farther southward migration and substantial weakening of the EJ. Transition to state D is characterized by the northward migration of the EJ to its most northern position, and the jet remains relatively weak. Strengthening and moderate southward retreat of the EJ characterize transition to state A, and then the cycle repeats itself.

## 3. Role of the mesoscale eddies in maintaining the mean flow

Since the pioneering work of Holland (1978), it is well known that the double-gyre models that dynamically resolve mesoscale eddies generate time-mean flows that are very different from the flows generated by the models that approximate the eddy effects in terms of the enhanced momentum diffusion. This is because the diffusion is often a poor representation of the actual large-scale–eddy interactions. The eddy forcing, *f* (*t*, **x**), diagnosed from the full model, represents history of nonlinear interactions of the eddies with the large-scale flow:

Here the angle brackets indicate that the large-scale flow is given by the conditional ensemble means, that is, by the key states. Similarly, *f* can be defined with respect to other choices of the large-scale flow, such as the infinite-time mean or a space–time filtered full flow. Given a large-scale–eddy flow decomposition, the corresponding *f* is defined uniquely. However, dynamical interpretation of *f* requires additional effort. A method to obtain such interpretation is based on applying the history of *f* (or one of its components) to force a non-eddy-resolving ocean model^{6} (Berloff 2005a; Berloff et al. 2007). In this framework, statistics of *f* can be related to the dynamical response.^{7}

We applied the method discussed above to the reference solution and arrived at the conclusion that the eddy forcing properties are similar to those reported in Berloff et al. (2007). Therefore, we just briefly list and discuss them: (i) eddy effects on the large-scale flow are dominated by the eddy-forcing fluctuations rather than the time-mean eddy forcing; (ii) in the EJ extension of the WBC, the fluctuations are dominated by their isopycnal-stretching (i.e., buoyancy) rather than the relative vorticity (i.e., Reynolds stress) component; (iii) the effect of the upper-layer eddy forcing dominates over the deep-ocean one. The main effect of the eddies, referred to as the flow *rectification*,^{8} is to enhance the upper-ocean EJ and its adjacent recirculation zones. This is a negative viscosity phenomenon associated with generation of the upgradient PV eddy fluxes across the eastern part of the EJ (McIntyre 1970). These fluxes involve strong correlations between the corresponding eddy forcing and the meridional velocity induced by it. This process resists the erosion of the EJ by the baroclinic instability acting upstream. Since the generic low-frequency mode appears only in the turbulent flow regime and actively involves the EJ (Berloff and McWilliams 1999a), it is likely that the rectification plays a key role not only in shaping the mean flow but also in driving the mode.

Along with space and time correlations, the variance is the key property of the eddy forcing. The variance is defined as

It is spatially inhomogeneous, with the maximum values located in the western part of the basin and around the EJ, and in the upper ocean (Fig. 5a). Also, *σ _{f}*(

**x**) substantially varies between the key states.

The flow transition from state A to B corresponds to the southward shift of the variance pattern (Fig. 5b). At the same time, the integral of the upper-ocean variance over the domain increases by about 10%. All of this is consistent with the fact that in state B, the EJ shifts to the south but remains enhanced (Fig. 3b). Below, in section 4 and in the appendix, we argue that these changes of the eddy forcing are associated with inhibited material and PV exchange between the gyres, and the transition from state A to B is induced by the associated potential vorticity anomaly that enhances potential vorticity contrast between the gyres.

The transition from state B to C is characterized by only a 2% increase of the total upper-ocean variance, by continuing meridional shifting of the EJ, and by substantial meridional rearrangement of the variance pattern (Fig. 5c). The rearrangement is characterized by the partial loss of the variance around the eastern end of the EJ, which is compensated by the gain of the variance in the western part of the EJ. The eastern-end variance is more efficient^{9} in maintaining the rectified EJ enhancement and the associated upgradient PV flux; therefore, its reduction induces a downgradient PV flux correction. The gain of the western-end variance has a relatively weak rectifying effect. In fact, the local intergyre PV flux even weakens due to the reduced PV contrast between the gyres. Implications of the above processes for the subtropical-gyre dynamical balance are discussed in section 4a.

The transition from state C to D is associated with a northward meridional shift of the variance. The variance remains relatively confined to the western boundary and, therefore, it is less efficient in terms of the induced rectification. This is consistent with the relatively weak EJ in state D. Finally, the transition from state D to A recovers the eastern-end variance and the corresponding rectified enhancement of the EJ reaches its maximum.

Although interpretation of the eddy kinetic energy in terms of the associated eddy effects on the large-scale flow is less straightforward than interpretation of the eddy forcing, it is a quantity that is relatively easy to diagnose from the solution. For comparison with *σ _{f}* , we analyzed the variability of the eddy kinetic energy and found a qualitatively similar picture, except that in terms of the basin integral, state B is more energetic than C. In particular, in the weak-EJ state the eddy kinetic energy along the EJ axis, except for its eastern end, is at maximum, consistent with observations of the Kuroshio Extension jet (Qiu and Chen 2005).

## 4. Dynamics of the low-frequency mode

We analyze dynamics of the low-frequency mode by calculating PV balances for different stages of the variability cycle. Let us refer to the process of adding/removing negative PV to/from the subtropical gyre, which is characterized by the negative-PV pool, as the local “heating”/“cooling.” Next, we focus on the upper-layer dynamics because here it plays a main role. The upper-ocean gyre balance, which is analyzed below, consists of PV inputs generated by the wind, created on the boundaries and fluxed to/from the subpolar gyre. In each key state, *K*, the ensemble-average eddy flux is given by

In the next section and further below, fluctuations are defined relative to the ensemble-mean key states.

However, transitions between the key states represent part of the low-frequency variability that drives its own eddy fluxes. Are these fluxes significant and important? We address this question directly by calculating the unconditional eddy PV flux of the low-passed reference solution (section 2) and by comparing it with the full-solution flux—the difference between those is associated with the mesoscale eddies. We find that the meridional flux, characterized by the positive values and well-defined maxima around the EJ and near the zonal boundaries, is completely dominated by the eddies (Fig. 6a). The zonal flux, characterized by cyclonic circulation around the EJ, is also dominated by the eddies but to a lesser extent (Fig. 6b). This is because of the low-frequency, coherent meridional shifts of the EJ. Thus, we find that the mesoscale contribution completely dominates the eddy PV flux: below we assess its role in driving the low-frequency variability mode.

### a. Dynamical evolution of the low-frequency mode: Upper ocean

#### 1) PV fluxes and main PV balance

For each key state *K*, the domain is divided into the subtropical and subpolar gyres with the intergyre boundary defined as the ensemble-mean streamline that connects the western and eastern basin boundaries, found by linear interpolation of 〈*ψ _{i}*〉

^{K}. The choice of a streamline as the intergyre boundary is very convenient because it guarantees that there is no ensemble-mean intergyre PV flux and, therefore, simplifies dynamical balances considered below.

We first calculate the PV inputs and sinks for the subtropical gyre. For accuracy, these calculations use the full solution, that is, with flow fields taken on each time step and at each grid point. The PV balance in the subtropical gyre is summarized in Table 1, showing four sources of PV. First, wind input integrated over the gyre area, a negative PV source, drives the system. Second, the diffusive flux of the relative vorticity, is given by

and the total boundary flux for the subtropical gyre is

where *C* is the contour of integration along the boundary and **n** is the corresponding normal vector. This diffusive flux has two components—the intergyre and boundary fluxes. Finally, there is the eddy flux across the intergyre boundary, calculated by finding normal vectors of unit length along the intergyre boundary, which are used to project eddy fluxes on the normal direction. The mean state of the gyre is thus negative PV, which is enhanced (reduced) by negative (positive) PV input. The gyre may respond to PV forcing by altering the mean flow speed, or by changing its shape and area.

We find that the main upper-ocean dynamical balance is always between the wind forcing and the diffusive boundary flux.^{10} Although diffusive boundary flux is large, its variance between the key states (i.e., rms value of the flux anomalies) is moderate. Intergyre eddy PV flux is only the third major player in the global balance, but its role in driving the low-frequency variability mode is crucial, as indicated by its large variance. The variance of the wind input is also large, reflecting variability of the gyre shape. The intergyre boundary is located substantially to the south of the zero wind curl line: therefore, when the gyre shrinks (expands), the corresponding wind input anomaly becomes positive (negative)—we refer to this as the geometric wind effect.

The total intergyre PV flux is always negative (Fig. 7a). Given that, in the upper-ocean subtropical (subpolar) gyres are characterized by negative (positive) pools of both PV anomaly, *q*_{1}, and absolute PV, *Q* = *q*_{1} + *βy*, implying downgradient mixing of *q*_{1} and *Q* between the gyres. Most of this exchange takes place near the western boundary where both the eddy kinetic energy and the intergyre PV gradient are maximal. In the interior of the basin, intergyre fluxes tend to be small but positive (i.e., upgradient): this is a characteristic feature associated with the eddy rectification process (section 3).

The diffusive boundary PV flux is largest along the western boundary, but it is also substantial along the southern boundary (Fig. 8a). Along the eastern boundary it is negligible. In the subtropical gyre, the total boundary flux varies due to two physical processes: first, there is meridional migration of the WBC separation point; second, the boundary layer flow dynamics changes. The total flux, *G*_{tot}, increases when either the separation point retreats to the south or the boundary current becomes more unstable. The former effect is simply the loss of the (positive) PV source due to a shorter contact base with the boundary. The latter effect is due to the Reynolds stresses generated by the flow instability: they sharpen the near-wall gradient of the mean relative vorticity, thus enhancing *G*_{diff} (Berloff and McWilliams 1999b). These two effects may be of opposite sign but can be discerned by the spatial location of the flux (Fig. 8).

In the following sections, the dynamics of the low-frequency mode is described, step by step, in terms of the transitions between the key states. The dynamical picture is summarized with the schematic diagram (Fig. 9), introduced here in order to help readers to navigate through the rest of section 4a.

#### 2) State A

An analysis of the differences between PV balances of the key states (Table 1) suggests the following dynamical picture. State A is characterized by intense heating of the gyre due to the negative wind forcing anomaly, inhibited intergyre eddy flux, and inhibited boundary flux. The heating is partially compensated by the diffusive intergyre flux that is at its maximum due to the strongest relative vorticity gradient across the EJ. The inhibited intergyre flux indicates that the fully developed EJ acts as a partial barrier to PV transport. This jet is also a partial transport barrier to the material transport—this argument is supported by calculating the Lagrangian exchange between the gyres (in the appendix). Inhibited permeability of the eastward jets has been a subject of different studies (e.g., Bower and Lozier 1994; Berloff et al. 2002). Here, we focus on the variations of this permeability over the low-frequency cycle and on implications of that for the dynamics. The lateral boundary PV flux is inhibited as well—this indicates that the WBC is relatively stable. According to the PV balance of state A, the subtropical gyre has to pick up negative PV: therefore the EJ has to strengthen, and PV contrast between the gyres and PV gradients in the WBC have to increase. However, there must be a physical process that terminates this process and reverses the low-frequency cycle. As argued below, this process is given by the change of the gyre shape induced by the PV input that increases PV contrast between the gyres.

#### 3) Dependence of the gyre shape on PV and eddy forcing anomalies

The following question requires clarification: why does the subtropical gyre shrink when it strengthens (i.e., heats up) and the intergyre PV contrast reaches its maximum in state A? Equivalently, why does it expand when it loses strength (i.e., cools down) in state C and the intergyre PV contrast reaches its minimum? These changes of the gyre shape are characterized by the corresponding meridional shifts of the EJ. We find that the low-frequency changes of the gyre shape are determined by the sign and magnitude of the PV input anomaly to the gyre. We have checked this statement directly by imposing positive or negative, upper-ocean PV forcing anomalies to the subtropical gyre, with a magnitude equal to that given in the last line of Table 1. These anomalies are exactly balanced by the imposed opposite-sign anomalies in the subpolar gyre, so they can be interpreted as a result of variability in the intergyre PV exchange. We find that PV anomalies that enhance PV contrast between the gyres (corresponding to a stronger intergyre barrier) induce contraction of the subtropical gyre and shift the EJ to the south. Anomalies reducing PV contrast work in the opposite direction. This qualitative result is found to be insensitive to the shapes of the imposed anomalies.

An alternative explanation of the EJ shift invokes changes of the eddy forcing variance around the EJ. We checked this hypothesis by imposing either a localized depression or enhancement of *σ _{f}*(

**x**) (14) around the eastern end of the EJ. This was achieved by multiplying interactively calculated eddy forcing fluctuations,

*f*′, with a smooth but highly localized function, which is either less or more than unity in the vicinity of the EJ but is equal to unity elsewhere. It is found that an enhancement of the eddy forcing tends to expand the subtropical gyre and, therefore, to shift the EJ to the north. On the other hand, a reduction of the eddy forcing works in the opposite way. Overall, the effect of the eddy forcing anomaly is exactly opposite to that of the PV anomaly.

The experiments described above confirm that PV anomalies drive changes in the gyre shape but do not indicate the mechanics of this phenomenon. This is an important issue, but one that is outside the scope of this paper and will be addressed in the future.

#### 4) State B

State B is a transitional state with a neutral PV balance (Table 1) in which heating by the moderately inhibited intergyre flux is balanced by the wind forcing and the boundary flux. This state is also characterized by some weakening of the transport barrier efficiency (the appendix): Flow transition from state A to B is characterized by the enhanced intergyre PV exchange in the western quarter of the basin, whereas in the interior of the basin the average flux remains the same (Fig. 7b). The flux enhancement can be attributed to an activated instability of the subtropical WBC and its EJ extension at some critical flow configuration and strength. The absence of such an enhancement farther to the east is a result of the relatively energetic eddies that continue to enhance the EJ, resisting the intergyre PV exchange.

State B is also characterized by the increase of the total boundary flux. This is so because enhancement of the flux along the western and southern boundaries, due to the instability, is not fully compensated by the flux loss because of the southward migration of the WBC separation point (Fig. 8b). Given the magnitude of this boundary effect, it seems problematic that the underlying instability can be accurately captured by linear stability analysis of simple geostrophic flows such as a zonal jet or near-boundary meridional current. We anticipate dynamical insight from stability analysis of a more relevant, two-dimensional flow configuration: however, such analysis is beyond the scope of this paper.

In state B, while the PV forcing is neutral, the intergyre PV contrast remains enhanced. This induces a farther shift of the EJ to the south [section 4a(3)], that is, toward state C. This behavior is analogous to a pendulum whose inertia causes it to swing past its stable equilibrium. In this analogy, the PV content of the gyre corresponds to the inertia.

#### 5) State C

Transition from state B to C is accompanied by increased permeability through the intergyre boundary, which is a consequence of the flow derectification. The main axis of the EJ migrates to its most southern position, according to the results of section 4a(3), but the associated loss of the subtropical gyre area is substantially compensated by the gain at the eastern end of the derectified EJ (Fig. 3).

State C is characterized by the most intense cooling, which is due to the minimal size of the gyre (and, hence, due to the weakest PV input from the wind forcing), the most intense intergyre flux, and the moderately enhanced boundary flux. The intergyre flux dominates, so the situation can be described in terms of the breached intergyre barrier and the associated accumulation of the PV anomaly that noticeably reduces PV contrast between the gyres. We find that transition from state B to C is characterized by the weakening of the downgradient intergyre PV flux in the western half of the basin. This flux reduction is not induced by reduction of the local intergyre mass exchange; on the contrary, the local mass exchange is at its maximum (see the appendix), driven by the energetic eddy field. This local weakening is due to the reduced intergyre PV contrast. On the other hand, this local intergyre flux reduction is overcompensated by the weakening of the upgradient fluxes at the eastern end of the EJ (Fig. 7c). This is a direct consequence of the derectification process, which, roughly speaking, opens a breach at the eastern end of the EJ.

In state C the boundary PV flux anomaly remains positive but reduces somewhat, thus indicating that instabilities in the WBC are still active. The combined PV input in state C is strongly positive, and the associated PV anomaly reduces PV contrast between the gyres. Thus, conditions become favorable for migration of the EJ to the north [section 4a(3)].

#### 6) State D

State D is a transitional state characterized by moderate heating, which is mostly due to the geometric wind effect. The gyre expands and is fully forced by the wind, not only because the EJ is shifted to the north but also because the EJ is relatively short, so that the intergyre boundary at its eastern end has its most northern position. In state D the boundary flux stays moderately inhibited, thus slightly heating the gyre. Due to the system’s inertia (weaker intergyre PV contrast: Fig. 4d) the EJ continues to migrate to the north [section 4a(3)], that is, toward state A. Rectification intensifies because the increased geometric wind effect translates into a more energetic eddy field and closes the intergyre transport breach at the eastern end of the EJ, thus, restoring efficiency of the barrier. The flow returns to state A and, then, the cycle repeats.

The above results demonstrate that the low-frequency variability depends on the geostrophic turbulence; hence we call it the *turbulent oscillator*. The dynamical picture of the turbulent oscillator is presented in terms of the simple diagram in Fig. 9. The essential dynamical ingredients involve eddy-driven flow rectification, the EJ acting as a transport barrier, and the response of the double-gyre flow configuration to PV forcing anomalies. The analogy between the low-frequency mode and the classical pendulum is useful: the latter is driven by competition between the gravity force and inertia, whereas the former is largely driven by competition between the gyre PV content and changes in PV forcing (both wind forcing and the intergyre transport barrier).

### b. Dynamical evolution of the low-frequency mode: Deep ocean

The analysis of PV balances is extended from the upper to deep isopycnal layers. We find that dynamics of the middle and bottom layers are qualitatively similar: therefore the presentation focuses only on the former. In the upper ocean, PV anomalies generated by the wind are balanced by PV fluxes through the lateral and intergyre boundaries. Qualitatively, the main difference between the upper and deep layers is the absence of the direct wind forcing in the latter. Quantitatively, all terms in the deep-ocean dynamical balance are significantly smaller than their upper-ocean counterparts (Table 2).

The deep-ocean gyres are forced by the eddy-driven downgradient mixing and homogenization of the background absolute PV, *Q* = *ζ* + *βy* (Figs. 10a,b). The intergyre flux is significant only in the western 1/5 of the basin, and its structure is similar to that in the upper ocean (Fig. 7a). This mixing process is the main process that maintains the pool of positive PV anomaly in the subtropical gyre. Intergyre diffusion of the relative vorticity adds to the positive PV anomaly by fluxing away negative relative vorticity from the anticyclonic recirculation located to the south of the EJ. The intergyre fluxes are balanced by the boundary fluxes that supply the gyre with negative PV. The deep WBC, similar to the upper-layer WBC, is characterized by positive PV flux through the boundary (Fig. 11a). Therefore, negative PV flux, which has to close the gyre balance, must come from a qualitative difference between the upper- and deep-ocean flow configurations, away from the western boundary. The difference consists of the elongated cyclonic (anticyclonic) recirculation cells appearing near the southern (northern) basin boundaries (Fig. 1b). The cells are responsible for the negative PV boundary fluxes (Fig. 11a) because, in comparison with the upper ocean, they have opposite-sign relative vorticity gradients on the adjacent boundaries. These features can be explained by weakly nonlinear interactions of the basin modes excited in the flow regimes characterized by the free- and partial-slip boundary conditions (Berloff 2005b).

The deep-ocean part of the low-frequency variability mode is largely driven by variations of the boundary PV flux; however their impact is substantially reduced by counteracting variations of the intergyre flux (Table 2). States A and C are characterized by heating and cooling of the subtropical gyre, respectively. In states B and D, different PV inputs are well balanced and their residuals are close to zero. The time evolution of the low-frequency mode can be described as the following (Figs. 10c–f). In state A both intergyre and western boundary PV fluxes are inhibited, as in the upper ocean, but the former is more inhibited. This balance is equivalent to the injection of positive PV through the western boundary (Fig. 11b). The injected PV partially annihilates the negative PV anomaly inherited from state D. This anomaly is advected to the west from the eastern basin. The rest of the injected negative PV is advected to the east by the EJ and its adjacent recirculations. In state B, the accumulated positive PV anomaly is advected into the eastern basin, where it is relatively well shielded from the intergyre and boundary fluxes. State C is characterized by the most active intergyre and western boundary PV fluxes, in accordance with the upper-ocean processes. The western boundary contribution dominates, and the outcome can be interpreted as an injection of negative PV anomaly through the western boundary. Relatively neutral state D mirrors state B. Overall, the deep-ocean low-frequency mode is characterized by mean-flow advection of PV anomalies generated on the western boundary, and its dynamics is slaved to the upper-ocean one. Neither eddy-driven flow rectification nor the intergyre transport barrier are the essential parts of the deep-ocean low-frequency dynamics.

### c. Parameter sensitivity study

How generic is the dominant large-scale, turbulent oscillator mode; that is, how sensitive is its space–time structure to parameters of the ocean model? Here, the most uncertain parameters are eddy viscosity *ν* and partial-slip coefficient *α* that for the reference solution are equal to 100 m^{2} s^{−1} and 15 km. We computed three series of solutions and several additional solutions, which differ from the reference solution in terms of *ν*, *α*, wind shape, and wind forcing amplitude. We analyzed the resulting time-mean states (Fig. 12) and the low-frequency variability patterns. The results are presented in Figs. 12 and 13 and are discussed below.

We start the parameter sensitivity study by varying the wind forcing. One solution is calculated with *τ*_{0} = 0.075 N m^{−2}, which is 75% of the reference solution wind amplitude. The turbulent oscillator mode is recovered, but it has a 17-yr period, which is 26% longer than the reference one. The asymmetry of the mode cycle increases: transition from state C to A is about 40% longer than the opposite transition. Sensitivity to the north–south asymmetry of the wind forcing is explored by considering an extreme modification of the reference wind curl in which all positive curl values are replaced by zero. In this case, given *α* = 0, wind drives a single subtropical gyre (Fig. 12o) and single EJ; the turbulent oscillator persists and is in the decadal frequency range; however, the spectral peak of the leading EOFs is not so well defined.

Next, we find that the turbulent oscillator appears only in the most turbulent regimes—having *ν* = 100 and 200 m^{2} s^{−1}—characterized by a well-developed EJ or a pair of EJs. A single EJ emerges for intermediate values of *α*, and in this regime the variability is characterized by a large fraction of the total variance and narrow spectral peak. Small-*α* and small-*ν* (i.e., large Re) regimes, characterized by the pair of EJs, exhibit a modified version of the generic variability. In this case, individual variability modes belonging to each EJ are spatially correlated with each other, thus forming EOF patterns that cover both EJs. Also, frequency spectra of the leading EOFs become substantially broader but still retain the decadal maxima. Last, we computed one solution with *ν* as small as 20 m^{2} s^{−1} (and all other parameters as in the reference solution) and found essentially the same low-frequency variability.

The transition to the small-Re regime occurs at around *ν* = 400 m^{2} s^{−1}, consistent with the earlier regime study by Berloff and McWilliams (1999a) using a two-layer quasigeostrophic (QG) model. Small-Re variability has been extensively studied in the past (section 1b) and, here, with one typical example in Fig. 13, we just want to emphasize how dramatically different the structure of this variability is from the generic one. The leading EOFs of the small-Re mode operate on much shorter, interannual time scale and exhibit weak time-lag correlations. Their spatial structure captures some mesoscale variability around the recirculation zones and some large-scale variability in the return westward currents of the gyres. The spatial patterns of the first and second large-Re EOFs are characterized by the asymmetric tripole and dipole (Figs. 2a,b) respectively; whereas it is hard to argue that similar patterns are present in the small-Re EOFs (Figs. 13a,b). The small-Re variability is so different because the actual eddy effects include negative viscosity phenomena (e.g., formation of the EJ through the rectification process: Berloff 2005a) or nondiffusive phenomena (e.g., EJ as a transport barrier: Berloff et al. 2002). These are too crudely parameterized by the enhanced eddy diffusivity.

By keeping the reference value of *ν* and varying *α*, we find strong sensitivity of the time-mean flow to the partial-slip parameter (as in Haidvogel et al. 1992). This is largely a result of the basin idealization, which assumes vertical lateral walls and no bottom topography (Leonov 2005). The turbulent oscillator mode of variability is found in both small-*α*/double-EJ (discussed above) and large-*α*/single-EJ regimes. However, the latter regime is physically unrealistic because it involves very regular generation of small and intense anticyclonic vortices from the subtropical WBC inflection near the separation point and their subsequent propagation to the northern boundary. (This process, in turn, substantially modifies the subpolar WBC and elevates the role of the intergyre PV fluxes in closing PV budget of the subpolar gyre.) Overall, the dominant period of the variability increases with *α:* for *α* = ∞ it reaches about 20 yr.

Last, we compare the turbulent oscillator with the gyre mode. To enable unambiguous comparison, first, we reproduced the three-layer QG simulations of Nauw et al. (2004). This simulation has asymmetric wind stress, *ν* = 900 m^{2} s^{−1}, larger Rossby radii, free-slip boundaries, and a 2000 km squared basin. We found a qualitatively similar solution, which can be compared to the small-Re solutions for our ocean configuration. In accord with Nauw et al. (2004), about 40% of the variance of this solution is captured by the gyre mode; and the gyre mode is recovered through the tripole–dipole pattern of the two primary EOFs. While one can argue that this spatial pattern is similar to that of the turbulent oscillator, the time scale of the gyre mode is very short—1 yr—compared to the turbulent oscillator time scale. This time scale is a factor of 4 shorter than that of the small-Re variability shown in Fig. 13. The temporal and spatial pattern differences between the gyre mode and the small-Re variability (Fig. 13) can be attributed to the factor of 2 difference in basin sizes. The smaller basin focuses on the EJ region and removes patterns outside of it: this makes the oscillation frequency spectrum more peaked but shifts it to higher frequencies. The gyre mode is created by generation of a single eddy in the EJ. By conducting a PV balance analysis of the gyre mode, we find that the intergyre PV exchange is dominated by diffusion, which is largest when the contrast of relative vorticity between the gyres is maximal. This state occurs when the eddy creates an EJ meander before it pinches off. This behavior is opposite to the turbulent oscillator, where the greatest intergyre PV flux occurs when the contrast of relative vorticity is weakest. Thus, although diffusion replaces the eddies as a major aspect of the variability mechanism, both the variability pattern and mechanism of the gyre mode are very different from those of the turbulent oscillator. This is so because diffusion does not correctly account for the eddy effects.

The gyre mode simulations shown here have a very short time scale, in keeping with previous simulations of this mode (Nauw and Dijkstra 2001; Simonnet and Dijkstra 2002; Simonnet et al. 2003a, b; Nauw et al. 2004). The only exception is the work of Simonnet (2005) who uses a barotropic model to show that the time scale of the gyre mode can be extended to decadal periods. However, this result requires the use of both Rayleigh friction to parameterize the effect of eddies and excessive wind stress; neither of these aspects is fully justifiable. Therefore, while the time scale of the gyre mode in Simonnet’s (2005) simulations is long, the processes that set it are not connected to the nonlinear and mesoscale dynamics, which are observed to dominate the turbulent oscillator.

To summarize the findings of this section, the turbulent oscillator studied in this paper seems to be a generic property of the double-gyre ocean with mesoscale turbulence. This mode is very different from the modes generated when the mesoscale turbulence is significantly suppressed. The spatial structure and the dynamical mechanism of the turbulent oscillator are remarkably robust, but the dominant period varies by a factor of 2 within the explored range of parameters. It is plausible that this mode is relevant to the ocean, but this hypothesis needs to be verified with more comprehensive ocean models and long-term observations.

## 5. Conclusions and discussion

Intrinsic large-scale, low-frequency variability is studied in the idealized, quasigeostrophic, midlatitude, wind-driven ocean gyres operating at large Reynolds number. The model solution is characterized by dynamically resolved, abundant and transient, mesoscale eddies as well as by the robust large-scale, decadal variability mode: the turbulent oscillator. This study finds that the mode is driven by the eddies and uncovers the underlying dynamical mechanism involving eddy–mean flow interactions. We hypothesize that the turbulent oscillator mode is dynamically equivalent to the 12-yr variability found in the more realistic model of the Kuroshio (Qiu and Miao 2000) and to the modes reported in the earlier studies using two-layer (Berloff and McWilliams 1999a) and three-layer (Hogg et al. 2005) QG models. It is also plausible that dynamically similar decadal variability can be observed in the western boundary currents of all midlatitude gyres. The mode of interest has the potential to influence the global climate because it can effectively couple with the atmosphere (Hogg et al. 2006). Finally, it is argued that, if the eddy effects are parameterized in terms of diffusion, the variability mode of interest degenerates and gives way to nongeneric modes associated with early bifurcations of the flow in the vicinity of the steady state (e.g., Berloff and McWilliams 1999a).

The turbulent oscillator is characterized by changes in meridional position and strength of the oceanic eastward jet and by the associated changes of the subtropical gyre. The stronger (weaker) jet state corresponds to the larger (smaller) subtropical gyre and to the EJ that is shifted to the north (south) [i.e., to the higher (lower) latitudes]. The mode is well represented by the leading pair of time-lag correlated empirical orthogonal functions. The first EOF captures meridional shifts of the EJ and the second EOF captures variability of the EJ strength. We find that the EJ is substantially enhanced by the nonlinear rectification process (aka, the turbulence backscatter process), which is driven by the eddy potential-vorticity flux divergence fluctuations. In the regime of interest, the main dynamical balance in the subtropical gyre is between PV inputs from the wind and boundaries, both intergyre and lateral. The intergyre PV exchange is smaller than the other two inputs in terms of the time-mean values, but its variance dominates the low-frequency cycle. Therefore, this exchange is fundamentally important for maintaining the variability. Another dynamical factor is the geometric wind effect: when the subtropical gyre expands (contracts) as a result of the negative (positive) net PV input, it covers a larger (smaller) area of the region forced by the negative wind curl; therefore, net wind forcing of the gyre changes. The third dynamical factor is the PV flux generated on the lateral boundaries.

On the one hand, the low-frequency variability is a fundamentally turbulent phenomenon, driven by the mesoscale eddies. On the other hand, the oceanic mode is largely driven by the competition between the eddy rectification process and the PV anomalies induced by changes in the intergyre transport barrier. This is analogous to the classical pendulum, which is driven by the competition between an external force and inertia.

The oscillation cycle has the following properties. In the strong-jet state, when the oceanic EJ is shifted to high latitudes, all dynamical inputs tend to enhance PV contrast between the gyres. The intergyre PV flux is inhibited because the strong EJ, maintained largely by the eddies, acts as a partial transport barrier—this property is confirmed by calculating Lagrangian particle transport. The mechanism limiting unbounded strengthening of the eastward jet in the northern state is represented by the low-frequency PV input anomaly that enhances the PV contrast across the barrier and pulls the jet away from the northern state. The jet migration is associated with shrinkage of the subtropical gyre. Transition to the opposite, low-latitude state occurs via partial derectification of the EJ, which breaches the intergyre transport barrier and reduces the intergyre PV contrast. The derectification of the EJ is a result of the geometric wind effect: the smaller subtropical gyre absorbs a lesser amount of PV from the wind. The reduction of the intergyre PV contrast creates a PV anomaly that reverses the cycle by driving the EJ away from the low-latitude, weak-jet state. As the subtropical gyre expands, the geometric wind effect starts to dominate, and the EJ rectifies back. This process restores the efficient intergyre transport barrier and the cycle repeats.

In the deep ocean, which is shielded from the direct effect of wind forcing, the roles of the other terms in the dynamical balance of the subtropical gyre are opposite to those in the upper ocean. The deep gyres are largely maintained by the intergyre PV fluxes, and this process is balanced by PV fluxes through the lateral boundaries. We find that the deep-ocean component of the low-frequency mode is slaved to the upper-ocean component and is characterized by mean-flow advection of PV anomalies generated on the western boundary.

We have shown that the turbulent oscillator is a different phenomenon to the gyre mode (e.g., Simonnet and Dijkstra 2002), although some features are similar. The gyre mode can be described mathematically, in terms of the dynamical systems language, as a pattern emerging from some successive bifurcations. We do not know the generic bifurcation sequence leading to the turbulent oscillator regime, but we can speculate—on the basis of some studies that trace the bifurcations toward the turbulent regime (Berloff and Meacham 1997, 1998)—that this sequence is incredibly complex and involves transitions to strange attractors, among other complicating issues. The later bifurcations, which are unknown, may be even more important than the earlier bifurcations. Given the above considerations, there are only two ways to proceed: The first way is to continue research deeper in the bifurcation sequences, relying on the dynamical system concepts. The second way is to introduce new concepts that are more compatible with the complex turbulent nature of the flow. In this manuscript we have followed the second path and investigated the statistics of the system, instead of attractors and bifurcations, keeping in mind that the diversity of methods is a scientific virtue.

The following developments of the results of this paper are anticipated. Given the importance of the geometric wind effect, the role of time-dependent wind forcing in driving the low-frequency mode has to be studied in ocean models, which are both uncoupled and coupled with the atmosphere. One has to study details of the wind-forcing pattern around the EJ. They can be important because their contribution constitutes the geometric wind effect. Also, an open question remains on how the low-frequency mode will be modified in dynamically more comprehensive models that account for more realistic coast lines, bottom topography, and thermohaline circulation. Finally, the mechanisms by which perturbations in the intergyre PV contrast induce meridional shifts of the EJ must be understood.

## Acknowledgments

Funding for Pavel Berloff was provided by NSF Grants OCE-0091836 and OCE-0344094, by the U.K. Royal Society Fellowship, and by the Newton Trust Award, A. M. Hogg was supported by an Australian Research Council Postdoctoral Fellowship (DP0449851) during this work, and William K. Dewar was supported by NSF Grants OCE-0424227 and OCE-0550139. Comments of the anonymous reviewers and G. I. Barenblatt led to improvement of the manuscript and, therefore, are gratefully acknowledged.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Lagrangian Intergyre Transport

Transport of material is represented by an ensemble of Lagrangian particles advected by the flow field. Lagrangian particle trajectories are obtained by solving the following equations:

The time integration of (A1) is performed by a fourth-order Runge–Kutta method with the right-hand side evaluated by cubic interpolation of the reference solution subsampled at every fourth point in each horizontal dimension.^{A1}

The method for estimating intergyre material flux is adopted from Berloff et al. (2002) and is applied here, both for the full solution and for the conditional key states. First, we release *m* = 10^{4} particles randomly distributed over the basin and follow them for *t* = 100 days. The ensemble of solutions is obtained for 200 uniformly distributed deployment times. A record is kept of locations where each particle that migrates from one gyre to the other first crosses the intergyre boundary. The ensemble-mean, intergyre flux in the *i*th layer, averaged over the time since deployment, is

where *V _{i}* =

*L*

^{2 }

*H*/

_{i}*m*is the fluid volume corresponding to each of

*m*deployed particles;

*N*(

_{i}*t*,

*x*) is the probability density of the first-time, intergyre boundary crossing; and the superscripts (

*n*,

*s*) indicate whether the crossing occurs from the northern or southern gyre.

^{A2}We focus on the upper layer (thus, the subscript is omitted), which possesses the key dynamics, and introduce

*N*= [

*N*

^{(}

^{n}^{)}+

*N*

^{(}

^{s}^{)}]/2. We find that the structure of

*N*(

*t*,

*x*) does not qualitatively depend on

*t*,

*m*and the number of deployment times as long as the latter two parameters are large.

On average only 2.6% of the total particle population migrates to the other gyre over *t* = 100 days. The corresponding *N*(*t*, *x*) indicates that most of the intergyre material exchange takes place in the western quarter of the basin (Fig. A1a). Next, we calculate *N* distributions for the key states, normalize them by the corresponding *N* distribution of the A state, and look at their differences (Figs. A1b–d).

Relative to state A, state B is characterized by the enhanced material flux near the western boundary. This difference is responsible for the enhancement of the corresponding PV flux. The material flux along the rest of the eastward jet remains largely unchanged. Thus, the Lagrangian diagnostics supports the idea that transition from state A to B involves an instability (and the associated burst of the eddy activity) of the combined western boundary current and EJ system rather than instability of the EJ extension alone. States B and C exhibit enhanced material fluxes in the locations of abrupt zonal changes of the intergyre boundaries (Figs. 4b,c); however, these fluxes do not translate into comparable peaks of the corresponding PV fluxes because of the weak local PV gradients. Relative to state B, state C exhibits not only enhanced flux near the western boundary, but also additional fluxes along the EJ extension—this indicates broad breaching of the EJ barrier. The latter flux enhancement is somewhat compensated by the flux reduction in the eastern part of the basin. Finally, transition from state D to A is associated with reduction of the material flux along the whole EJ, that is, with rebuilding of the strong intergyre barrier. However, this reduction is partially compensated by the enhanced flux in the east.

To summarize, state A, characterized by a relatively strong EJ, is more of the material intergyre transport barrier than the other states, and state C in particular. All of the above results are in accordance with patterns of the corresponding PV fluxes (Fig. 7). Therefore, descriptions of the variability cycle in terms of the intergyre material and PV transports are consistent with each other.

## Footnotes

*Corresponding author address:* P. Berloff, Clark Laboratory, Woods Hole Oceanographic Institution, MS 29, Woods Hole, MA 02543. Email: pberloff@whoi.edu

^{1}

Another common oversimplification is considering a space- and time-independent eddy diffusivity coefficient.

^{2}

The lateral boundary condition, which can be viewed as a simple parameterization of the unresolved boundary layer dynamics, has another free parameter of the problem.

^{3}

This has been originally demonstrated by Dewar (2003). We confirmed Dewar’s findings by applying a similar analysis to the ocean model and the circulation regime considered here.

^{4}

This can be viewed as a special filtering of the large-scale flow component.

^{5}

Qiu and Miao (2000) also find alternating PV anomalies in the EJ region of the Kuroshio.

^{6}

This requires coarse graining of *f* from the original fine grid to the coarse grid of the non-eddy-resolving model.

^{7}

Furthermore, if the *f*-forced non-eddy-resolving model acceptably approximates the original full flow, then an ocean modeler can think about approximating the actual *f* with some simple (stochastic) parameterization.

^{8}

Rectification is the generation of a nonzero time-mean flow response by a forcing with zero time mean. Another commonly used and analogous term is turbulence backscatter. Here, it requires nonlinearity of the dynamics and planetary vorticity gradient (Haidvogel and Rhines 1983).

^{9}

This has been verified by considering large-scale effects of localized eddy forcing patterns imposed at different locations within the double gyres.

^{A1}

Transport by unresolved fluctuations is not accounted for, but it is expected to be small (Armenio et al. 1999).

^{A2}

The *n* and *s* components of *M _{i}* are equal because of the mass conservation.