## Abstract

Material spreading and mixing by oceanic mesoscale eddies are analyzed in an idealized, numerical model of the wind-driven, midlatitude oceanic circulation. The main focus is on the regime with large Reynolds number, Re, and vigorous mesoscale eddies, although brief comparisons are made with less turbulent solutions at smaller Re. The analyses are based on ensembles of Lagrangian particle trajectories. The authors find that tracer transport by mesoscale eddies differs in many ways from the commonly used model of homogeneous, isotropic eddy diffusion.

The single-particle dispersion, which describes the spreading process, is generally anisotropic and inhomogeneous and in most places it is not diffusive (i.e., not linear in time) during intermediate-time intervals after tracer release. In most of the basin and especially in the deep layers, subdiffusive single-particle dispersion occurs due to long-time trapping of material by coherent structures such as vortices near the strong currents and planetary waves in the eastern part of the gyres. Superdiffusive dispersion behavior is found in the western part of the subtropical gyre and near the boundaries in fluctuating jetlike flows. Sub- and superdiffusion are associated with a strong first negative and second positive lobe, respectively, in the Lagrangian velocity autocorrelation function. The two-particle dispersion, which describes the mixing process, is characterized by initial exponential growth, and its exponent has strong geographical inhomogeneities, with faster rates in the upper western gyres. The finite-time Lyapunov exponents and recurrence times—other descriptors of the mixing—indicate a similar geographical partitioning, but are better alternatives to the two-particle dispersion on short and long time intervals, respectively.

Over large time intervals, due to inhomogeneity of the transport properties, material spreads and mixes much faster within the subtropical than the subpolar gyre. The eastward jet extension of the weaker subpolar western boundary current behaves as a strong barrier to the intergyre flux, but the eastward jet of the stronger subtropical gyre behaves as a weaker transport barrier. It is shown that net Lagrangian fluxes across the eastward jets are less intensive but more efficient in the deep ocean. The permeability of the barrier and the main pathways across it are measured with distributions of crossing particles, with finite-time Lyapunov exponents, and, at small Re, with the geometry of turnstyle lobes.

## 1. Introduction

The paper considers some aspects of material transport in midlatitude, wind-driven gyres that occur in all major oceanic basins. In addition to advection by the mean circulation, material tracers spread and mix by turbulent eddies, most rapidly by mesoscale eddies within isopycnal layers (beneath the surface turbulent boundary layer). The *spreading* (i.e., propagation) of material from its initial location is described by single-particle dispersion and the evolution of its distribution [i.e., the probability distribution function (PDF): *P*(**x**, *t* | **x**_{0}, *t*_{0}) for a source at **x** = **x**_{0} and times *t* ≥ *t*_{0}]. The *mixing* is related to the patchiness of individual realizations of tracer distributions, and it is described by finite-time Lyapunov exponent, two-particle dispersion, and homogenization rate.

Our approach is diagnostic: we analyze the transport properties exhibited by ensembles of Lagrangian particles deployed in the velocity field obtained as a solution of an eddy-resolving, quasigeostrophic (QG) model. The analyses are based on calculating ensembles of Lagrangian particle trajectories. We consider three solutions—at large, intermediate, and small Reynolds numbers (Re)—with most attention given to the more realistic, large-Re solution. In the introduction we provide some background and pose the problem. Most of the transport analysis deals with the large- and intermediate-Re solutions, and the results are presented in section 2. Section 3 focuses on the small-Re solution. Section 2a discusses general aspects of tracer evolution from localized releases, sections 2b and 2c focus on tracer spreading and mixing processes, respectively, and section 2d analyzes the large-scale meridional transports. Conclusions and discussion follow in section 4.

### a. Statement of the problem

We focus on statistical equilibrium solutions of an eddy-resolving, double-gyre, wind-driven, QG model of the midlatitude oceanic circulation. For the last two decades the dynamically simple double-gyre circulation has been a popular paradigm of the ocean, but it is still unknown to what extent it simulates real oceanic processes. Beyond the present model setup, we envision similar transport analysis of a hierarchy of idealized ocean models with growing physical complexity and, therefore, with related changes of the transport properties (see also section 4).

This paper analyzes the transport properties of the midlatitude oceanic circulation model by asking the following questions:

What is the geographical partitioning of the basin in terms of the local single- and two-particle dispersion behavior at intermediate times and the local value of finite-time Lyapunov exponent?

What are the principal long-time transport barriers and mixing zones? In particular, what are the most permeable regions and main pathways across the intergyre boundary and the subtropical eastward jet extension of the western boundary currents (WBCs)?

What are the roles of the time-mean circulation and coherent structures such as rings, eddies, and planetary waves? How do these behaviors vary depending upon the Reynolds number?

An important motivation for answering these questions is to provide a body of evidence for devising better transport parameterizations for large-scale, oceanic models. A companion paper (Berloff and McWilliams 2002) tests an hierarchy of stochastic transport models using the analyses presented here.

The complete model description is in Berloff and McWilliams (1999a). Here we state the particular parameters of the solutions analyzed below. The wind stress is steady and zonal,

where *τ*_{0} = 0.04 N m^{−2}; *L* = 3840 km is the size of the square basin (with 0 ≤ *y* ≤ *L*); and the intergyre asymmetry parameter^{1} is *ε* = 2. The bottom is flat. The vertical stratification is represented by three layers with mean thicknesses *H*_{1} = 200 m, *H*_{2} = 1200 m, and *H*_{3} = 2600 m. The ratio of the density jumps across the internal interfaces is *γ* = (*ρ*_{2} − *ρ*_{1})/(*ρ*_{3} − *ρ*_{2}) = 2, and the first internal Rossby deformation radius is 52 km. The only dissipative process is lateral eddy momentum diffusion with eddy viscosity coefficients *ν* = 100, 800, and 1680 m^{2} s^{−1} in the three solutions analyzed here. The large-scale Reynolds number is defined by

where *U* = *τ*_{0}(*ρ*_{1}*H*_{1}*Lβ*)^{−1} is the velocity scale derived from the Sverdrup balance (with *β* = 2 × 10^{−11} m^{−1} s^{−1}). Thus, the solutions here are characterized by Re ≈ 100 (large), Re ≈ 12 (intermediate), and Re ≈ 6 (small). The small Re regime is near the first Hopf bifurcation of the flow solution and it is periodic in time. The dynamical model is solved numerically on uniform horizontal grid. The resolution requirements are such that the viscous length scale (Munk scale), *δ*_{M} = (*ν*/*β*)^{1/3}, is resolved by more than two grid points (Berloff and McWilliams 1999a). In the large-, intermediate-, and small-Re regimes, the grid resolutions are 7.5 (513 × 513), 15 (257 × 257), and 30 km (129 × 129 grid points), respectively.

The transport of tracer is represented by ensembles of material particles advected by the velocity field of the dynamical solution. Lagrangian trajectories of the particles are obtained by solving the nonautonomous system,

where **x**(*t, ***X**) is the position of the particle initialized at (*t*_{0}, **X**), and **u** = (*u, **υ*) is the nondivergent, geostrophic, horizontal velocity, which can be decomposed into time-mean and fluctuating components. Equation (3) may be written for *i*th layer as

where *ψ*_{i} is the corresponding geostrophic streamfunction. The time integration of (3) is performed by a fourth-order Runge–Kutta method with the right-hand side evaluated by cubic interpolation of *ψ*_{i}(*t, **x, **y*) archived on the uniform 129 × 129 grid every day for 10 000 days. Sensitivity tests^{2} show that, provided spatiotemporal cubic interpolation the trajectories are calculated accurately but at Re larger than presented here, the data have to be sampled more densely both in space and time. Calculating material transport of the eddy-resolving solutions, we neglect transporting effects of unresolved, that is, subgrid-scale, fluctuations. It is expected that such fluctuations have small influence on the motion of the Lagrangian particles (Armenio et al. 1999).

In 2D the dynamics of an individual trajectory governed by (3) is a deterministic and time-reversible process that expresses the transport as stretching and folding of material lines (Ottino 1989). On the other hand, the evolution of the ensemble-mean PDF for particle positions at later times after a local initial release is an irreversible process, whose outcome can be identified with a spreading ensemble-mean tracer concentration field. The PDF can be estimated by solving (3) with an ensemble of release times, but it is commonly modeled as a stochastic process (Rodean 1996).

### b. Background

Hydrographic measurements of oceanic materials have provided the primary empirical basis for describing the general circulation through its transport rates (e.g., Reid 1997). Direct observations of transport with neutrally buoyant, subsurface floats are increasingly used (e.g., Freeland et al. 1975; Schmitz et al. 1981; Colin de Verdiere 1983; Rossby et al. 1983; Riser and Rossby 1983; Davis 1991; Owens 1991; Spall et al. 1993; Carr et al. 1997; Davis 1998; Boebel et al. 1999), though their sampling requirements are formidable for the global ocean. Float trajectories are used to estimate local velocity, **u** = **u** + **u**′, which is the sum of time-mean and fluctuation components; Lagrangian velocity autocorrelation function,^{3}

and its associated frequency spectrum; Lagrangian integral timescale,

Lagrangian velocity variance,

[where **u**′(*t, ***x**) is velocity fluctuation]; single-particle dispersion [see (8)], *D*(*t*); lateral eddy diffusivity, *K*_{ij}; and the time- and length scales of mesoscale eddies (e.g., Colin de Verdiere 1983; Krauss and Boning 1987; Sundermeyer and Price 1998; Rupolo et al. 1996). Ocean measurements show that *R*_{ii}(*τ*) decays, changes sign, oscillates, and often has a decay rate so slow that it is difficult to estimate *T*_{L}. In the real ocean, accuracy in estimating *R*_{ii}(*τ*) and σ_{ij}(*x, **y, **z*) is limited by the sparseness of available Lagrangian float datasets: the statistical error bars are large, and it is impossible to accurately remove the mean-flow contribution.

Single-particle dispersion functions (Taylor 1921) show a common behavior only in the limit *t* → 0 when the ballistic regime, *D* ∼ *t*^{2}, is recovered. For intermediate-time ranges (e.g., from 20 to 100 days), single-particle dispersion rates are not universal. Power-law fits, *D* ∼ *t*^{α}, typically yield 0.5 ≤ *α* ≤ 1.5 [cf. Fig. 10 in Krauss and Boning (1987)], straddling the eddy-diffusion law (i.e., *α* = 1). In some cases *D*(*t*) tends to saturate after about 100 days. The widespread occurrence of nondiffusive dispersion behavior, effects of boundaries and inhomogeneities, and long-time velocity correlations suggest that modeling the eddy-induced transport as a diffusion process (or, equivalently, as a random walk), may not be reliable. However, this is the transport parameterization commonly used in ocean general circulation models.

Single-particle dispersion describes spreading of a tracer, whereas two-particle dispersion describes its mixing, that is, disappearance of the patchiness in individual realizations, distortion, and homogenization (Batchelor 1952). The dispersion of pairs of passive particles has been addressed empirically by Richardson (1926) and advanced theoretically by Obukhov (1941) and Batchelor (1952). Theory predicts that in the inertial range of 3D, isotropic, homogeneous, stationary turbulence, the two-particle dispersion grows as *t*^{3}, as confirmed in laboratory experiments (Jullien et al. 1999) and atmospheric balloon measurements (Er-El and Peskin 1981). A two-particle analysis of oceanic float trajectories is in the very preliminary stage (LaCasce and Bower 2000) where it is found that mixing rates are faster in the western rather than eastern part of the North Atlantic; and it is confirmed that the observed mixing scenario is consistent with that of 2D turbulence where energy cascades toward large scales. Alternatively to the two-particle dispersion, the mixing process can be analyzed by the means of finite-time Lyapunov exponents, as it is shown in kinematic (e.g., Pierrehumbert 1991) and dynamic (e.g., Pierrehumbert and Yang 1993) models of atmospheric flows.

In addition to single- and two-particle descriptions of the ocean, a relatively large number of studies focus on transport across eastward jet extensions of the WBCs such as the Gulf Stream. Upper-ocean floats deployed near the Gulf Stream separation point at Cape Hatteras show inhibited crossings of the instantaneous Gulf Stream, but deep in the thermocline, where potential vorticity gradients are weaker, floats cross the Gulf Stream more frequently (Bower and Lozier 1994). Most of the floats that cross the current are not trapped in Gulf Stream rings, but they are advected by other velocity fluctuations. The hypothesis that jets are barriers to transport has been confirmed in a laboratory experiment (Sommeria et al. 1989); in simple kinematic models (i.e., with velocity that is not a solution of dynamical equations) where the flow consists of propagating disturbances superimposed on a steady zonal jet (Bower 1991; Samelson 1992; Liu and Yang 1994); and in dynamic, local, barotropic models of the zonal jet in linear (Pratt et al. 1995) and weakly oscillating, nonlinear (Rogerson et al. 1999) regimes. However, transport properties of eastward jets in more realistic and nonlinear flows at large Re remain unexplored.

Dynamics of idealized, eddy-resolving models of wind-driven gyres has been intensively studied (e.g., Holland 1978; Haidvogel et al. 1992). On the other hand, transport phenomenology of such models has been rarely looked at, and most of such studies focused on estimating either intergyre and cross-jet fluxes or eddy diffusivity coefficient. In a symmetrically forced, double-gyre, dynamical model with weak eddy variability, particles deployed near the WBC separation point and tracked for several weeks rarely cross the eastward jet extension of the WBC (Bower and Lozier 1994), in agreement with analogous observations; however, it is not clear what fraction of the cross-jet transport is diagnosed in this way. Figueroa (1994) shows that the cross-jet flux occurs because of meandering streamlines near the eastward jet core, rings detaching from the jet, and other types of fluctuations. The inhomogeneous eddy-diffusivity tensor estimated from a double-gyre, eddy-resolving solution is substantially anisotropic with a dominant zonal component (Rhines and Schopp 1991; Figueroa and Olson 1994); however, in these experiments *R*(*τ*) does not always decay to zero on inhomogeneity scales, thus precluding a diffusion regime. A coarse-grid simulation of advective–diffusive tracer transport using the fitted eddy-diffusivity tensor of Figueroa and Olson (1994) displays a substantially different evolution from an eddy-resolving solution (Figueroa 1994). This important result supports the idea that transport by mesoscale eddies requires a better parameterization than eddy diffusion. Recently we used a particular class of double-gyre solutions to investigate low-frequency intrinsic variability associated with large-scale changes in strong currents (Berloff and McWilliams 1999a), instabilities in the WBC (Berloff and McWilliams 1999b), and the emergence of mesoscale coherent vortices at very large Re (Siegel et al. 2001). Here we use these solutions to analyze the transport phenomenology.

The transport by well-structured events, such as ring detachment and merger or by other simple flows, can be analyzed geometrically using dynamical system techniques based on the calculation of the finite-time invariant manifolds of the hyperbolic fixed points of the flow (Poje and Haller 1999). The mathematical theory of invariant manifolds and their partitioning of the flow has been developed over the last decade (e.g., Beigie et al. 1994). In a weakly fluctuating, double-gyre model, the manifolds can be used to calculate the intergyre flux (Coulliette and Wiggins 2001; also see section 4).

## 2. Transport at large Re

The upper-ocean time-mean circulation (Fig. 1a) is calculated at each grid point. It consists of two gyres, filling about two-thirds (southern or subtropical) and one-third (northern or subpolar gyre) of the basin. The boundary between the gyres is defined as the time-mean streamline that originates at the separation point of the subpolar eastward jet and ends at the eastern boundary; by this definition the stronger subtropical eastward jet lies well south of the intergyre boundary. Instantaneous flows (Figs. 1b,d) contain a rich variety of coherent mesoscale fluctuations: meandering swift currents, intense vortices, eddies, and planetary waves. The circulation in the deep ocean is characterized by a weak time-mean flow in the eastern basin, a cyclonic recirculation zone under the subtropical WBC, and widespread eddies and waves (Figs. 1c,d). The middle and bottom layers have similar dynamic and transport properties, so in the text we will refer mostly to analyses for the middle layer. At large Re the WBCs become unstable and the resulting fluctuations contribute substantially to the local variability (Ierley and Young 1991; Berloff and McWilliams 1999b; Siegel et al. 2001). Given the dynamical simplicity of the model, an open question remains to what degree such variability represents real oceanic processes.

The geographical distributions of Eulerian and Lagrangian velocity statistics allow for a phenomenological partitioning of the domain into an eastern part dominated by planetary waves at all depths and a western part dominated in the upper ocean by swift, meandering currents and intense, sparse mesoscale vortices and at depth by weaker but more densely packed mesoscale eddies. In each layer the velocity variance tensor (7) has nonzero, nondiagonal components, and it is locally characterized by the pair of eigenvalues, *σ̂*_{ii}, and the principal angle that shows deviation of the direction of largest velocity variance eigenvector from the closest coordinate direction. The fluctuation amplitude, *σ̂*_{11} + *σ̂*_{22}, (Figs. 2a,b) is largest in the WBCs and their eastward jet extensions, in the upper ocean, and in the subtropical gyre. The principal angle (Figs. 2c,d) is mostly aligned with the rectangular boundaries, except in the strong currents where its diagonal orientation is related to Reynolds stress by eddies (e.g., Pedlosky 1987). The anisotropy ratio, *σ̂*_{11}/*σ̂*_{22}, (Figs. 2e,f) is large along the boundaries. The swift currents are nearly isotropic (i.e., ratio ∼ 1), and the interior of the basin is characterized by enhanced meridional fluctuations (ratio ∼ 2).

### a. Tracer evolution from localized releases

We diagnose the Lagrangian transport by solving (3) for ensembles of 10^{3} particles in each layer locally released within 50-km squares at a sequence of 40 deployment times uniformly distributed over the solution record. We illustrate different transport behaviors with particles initially at six (*x, **y*) locations (in km): 1 (35, 1050); 2 (115, 1050); 3 (1250, 1050); 4 (35, 3250); 5 (115, 3250); 6 (1250, 3250). Locations (sites) 1 and 4 are to the west and 2 and 5 are to the east of the WBC velocity maxima, and 3 and 6 are in the central parts of the gyres (Fig. 1a). The evolving PDFs of particle positions (i.e., ensemble-mean tracer concentrations) are found by binning the basin, counting the number of particles in each bin, and normalizing the particle density so that basin integral of the PDF is unity.

In the upper ocean, the particles from sites 1–3 spread all over the subtropical gyre within about 1000 days (Fig. 3). This is faster by an order of magnitude than in the small-Re case (section 3), even though the Sverdrup gyre circulation is the same for all Re. Over the same time interval, the subpolar gyre is only weakly invaded because the subpolar eastward jet acts as a strong transport barrier (section 2d). There are rather small differences between sites 1 and 2, although in the latter case the particles spread over a larger area before exiting the eastward jet so the concentration front propagating into the eastern basin is less sharp. The similarity between sites 1 and 2 is because near-boundary vortices and meanders effectively mix fluid within the WBC (which is not true at small Re; see section 3). The PDF evolution from site 3 (not shown) is similar to site 2, but it is delayed by the time required for particles to arrive at the subtropical WBC. In the eastern basin some particles remain concentrated in the concentration front for a long time. This is an indication of the slow single-particle dispersion rates in the eastern basin (cf. sections 2b and 2d), where planetary (Rossby) waves are dominant. Only a small percentage of particles from site 4 escapes to the subtropical gyre (Fig. 4) due to the intergyre barrier [section 2d(1)]. Most of these transiting particles escape from the near-boundary part of the subpolar WBC, in agreement with the patterns of intergyre fluxes and recurrence times [sections 2d(1) and 2d(3)]. The particles that are advected from sites 5 and 6 by the time-mean flow to the eastern and northern parts of the gyre remain grouped in filaments for thousands of days [sections 2c(2) and 2c(3)] and spread only within the inner part of the subpolar gyre (Fig. 4).

In the deep ocean, the particles from sites 1–3 slowly fill the western part of the subtropical gyre and only after that start to penetrate into the subpolar gyre (Fig. 5). Eastward spreading of the particles is inhibited because of very weak time-mean flow away from the WBC, and also because of slow single-particle dispersion rates (section 2b). Even slower spreading occurs in the middle of the subpolar gyre from site 6 (Figs. 6e,f). The particles released in the subpolar WBC at sites 4 and 5 propagate in two primary directions: along the northern boundary, and across the subpolar eastward jet and farther into the subtropical gyre (Figs. 6a–d). Both routes of spreading are dominated by velocity fluctuations rather than time-mean advection, and they are associated with fast single-particle dispersion rates (section 2b).

There is an overall qualitative similarity between the evolving PDF patterns in the different solutions at high-Re and intermediate-Re (not shown), consistent with both having mesoscale variability with broadband space and time spectra. However, the spreading and mixing rates in the latter case are uniformly smaller. This difference is particularly noticeable in a smaller intergyre flux and a sharper propagating front in the subtropical eastern basin. In the deep ocean, the relatively fast-filled western part shrinks to about one-fourth of the basin, and there is no enhanced spreading along the zonal boundaries.

### b. Tracer spreading

This section is the central part of the paper. Here we introduce the idea of intermediate-time, single-particle dispersion power laws and provide geographical partitioning of the basin in terms of the power-law exponents.

The single-particle (or absolute) dispersion tensor is

where the overbar denotes an ensemble average over particles released at the same initial position, **x**(0). Both fluctuating and time-mean components of the flow contribute to *D*_{ij}. The time-mean contribution is initially ballistic [i.e., *D*(*t*) ∼ *t*^{2}], because over small distances the mean flow is essentially unchanged. At long times, comparable to gyre circuits, the mean flow structure causes *D*(*t*) to oscillate. At very long times *D*(*t*) converges to a constant value that corresponds to a uniform particle distribution between the maximum and minimum streamlines in the initial conditions. The contribution from the fluctuating velocities is also initially ballistic, but only over Lagrangian velocity correlation time. To quantify the contribution from the fluctuating velocities, we calculate

where **x**′(*t*) evolves in (3) only by the local velocity fluctuation, **u**(*t, ***x**) − **u**(**x**). Then we fit to a power-law form,

dividing each layer in 40 × 40 spatial bins. In each bin we make 1000 individual realizations (25 particle trajectories for 40 different deployment times). We test the sensitivity of the estimate for *α*_{ij} by using a larger ensemble size and find that the geographical map of *α*_{ij} is robust. We further test the robustness of our estimates of *α* by using an alternative definition of the single-particle dispersion that differently separates the effect of the fluctuations from the time-mean flow; namely,

where **x**(*t*) evolves by the full velocity field, and **u** is calculated along individual trajectories. We find that *D̂*(*t*) and *D*′(*t*) yield similar distributions of *α* everywhere except in the subtropical WBC, where the former is slightly larger. The discrepancy is due to the fact that in swift currents the alternative definition is more limited by the local homogeneity assumption formulated below. Given the similarity, we base our analyses on *D*′(*t*).

The intermediate-time fitting interval for (10) has to be well separated from the ballistic timescale. We estimate an upper bound for this scale as the Lagrangian correlation time, *T*_{ii}, calculated as the time at which *R*_{ii}(*τ*) approaches its first zero. For a given range of *τ, **R*_{ii}(*τ*) may not reach zero; hence for the following we define local *T*_{ii}(*x, **y*) as the time at which the local ensemble-averaged *R*_{ii}(*τ*) reaches *R*_{0} = 0.3 for the first time. Since the flow is anisotropic, elements of the diagonalized *T*_{ii} are found from *R*_{ii}(*τ*) calculated along principal directions of the velocity variance tensor (7). The Cartesian-coordinate form of *T*_{ij}(*x, **y*) is then obtained by rotation over the principal angle.

Here *T*_{11} and *T*_{22} vary smoothly from about 3 days in the upper-ocean subtropical WBC to 11–13 days in the eastern basin and in the deep ocean. Usually *T*_{12} is small, and we neglect it for simplicity. The power law (10) is fitted for a time interval starting well after *T*_{ii}(*x, **y*), but stopping before the single-particle dispersion distance, *L*_{D′ii}(*t*) = *D*^{′}_{ii}(*t*), becomes much larger than the bin size in any direction. In each bin we define 〈*T*〉 = (*T*_{11} + *T*_{22})/2 (Fig. 7) and make a least squares fit to (10) from 3〈*T*〉 to the smaller of 200 days or 10〈*T*〉. The former upper limit holds in the eastern and northern regions characterized by relatively weak fluctuations and slow growth in *D*^{′}_{ii}(*t*). The latter upper limit is a compromise between having a sufficiently long temporal fitting interval and a minimal spatial bin size; in the worst case it leads to *L*_{D′ii} ≈ 200 km near the subtropical WBC separation point. We interpolate *α*_{ii}(*x*_{bin}, *y*_{bin}) (calculated at each bin) to obtain *α*_{ii}(*x, **y*). The fitting algorithm is based on the assumption of *local homogeneity*: during the fitting interval *L*_{D′ii} remains smaller than the local inhomogeneity scale, *L*_{Ii}, for Lagrangian statistical properties. This assumption is not fully satisfied near the subtropical WBC in the direction perpendicular to the boundary, where *L*_{It} is relatively small and velocity fluctuations are large.

The anisotropy and spatial variations of the *α*_{ii}(*x, **y*) are large, especially in the upper ocean (Fig. 8). We define a location as subdiffusive if *α*_{ii}(*x, **y*) < 0.8, as superdiffusive if *α*_{ii}(*x, **y*) > 1.2, and as diffusive if 0.8 < *α*_{ii}(*x, **y*) < 1.2. Most of the basin is subdiffusive in at least one direction. There are two main subdiffusive regions: the eastern and northern parts of the gyres and the subtropical WBC with its eastward jet extension. In Fig. 9 we show *D*^{′}_{ii}(*t*) averaged over four neighboring bins for several characteristic upper-ocean locations marked in Fig. 8. Subdiffusion is always associated with an oscillatory *R*(*τ*) with a pronounced first negative lobe, and superdiffusion is associated with a relatively small first negative and big second positive lobe of *R*(*τ*) (Fig. 10). The off-diagonal elements of *R*_{ij}(*τ*) are not small in some places [i.e., where the principal angle in Figs. 2c,d is not nearly aligned with the *x* or *y* axis, for *R*(0)], but they do not exhibit any different timescales beyond those in *R*_{ii}(*τ*) in Figs. 10 and 11. The oscillations of *R*(*τ*) indicate rotational motion of the particles trapped in the vortices or oscillatory displacements due to planetary waves. The asymmetry of *R*(*τ*) in the superdiffusive regions indicates that the oscillatory motion is combined with a sustained drift. A global autocorrelation function, 〈*R*(*τ*)〉, is calculated from 400 randomly initialized trajectories integrated for 10 000 days (Fig. 11). It has the same characteristic feature as the local *R*(*τ*) away from the swift currents (Fig. 10): oscillation on a timescale of about 40 days. Because 〈*R*(*τ*)〉 is averaged over regions with different autocorrelation properties, its anisotropy is relatively weak, and its negative lobe is relatively less pronounced.

In the upper ocean the single-particle dispersion is subdiffusive in the subtropical WBC and its eastward jet extension (Fig. 9a) due to material trapping in strong vortices (similar to Gulf Stream rings) [Fig. 1b; also see section 2c(1) for estimating the transport barrier between the fluid inside and outside the rings by finite-time Lyapunov exponents]. Analogous trapping behavior occurs for coherent vortices in 2D turbulence (e.g., Provenzale 1999). Superdiffusive regions occur, especially in the upper ocean, to the east of the subtropical WBC and south of its eastward jet extension, in the middle of the subtropical gyre, and between the two eastward jets near the western boundary (Figs. 8, 9c). These regions are affected by the low-frequency variability of the eastward jet, and they include the main population of strong vortices. In 2D turbulence results (Elhmaidi et al. 1993), the intermediate-time single-particle dispersion behavior of Lagrangian particles outside strong vortices is weakly superdiffusive due to the large shear strain rate. Large-scale, low-frequency oscillations (Berloff and McWilliams 1999a) also may contribute to superdiffusivity near the subtropical eastward jet and its local recirculation, but we have not quantitatively estimated that contribution. The subdiffusive regions in the eastern and northern parts of the basin (except near the boundaries in the tangential direction) are due to coherent oscillations by planetary waves. The zonal subdiffusion in the eastern basin (N.B., *D*_{11}(*t*) in Fig. 9d) is partly responsible for particle distributions showing a sharp front (Fig. 3), but the velocity anisotropy also contributes to the sharpness of the front (i.e., *σ*_{22}*σ*^{−1}_{11} ≈ 2). The subpolar gyre and its eastward jet are predominantly subdiffusive in the meridional direction (Fig. 9b). This property is consistent with the inhibited tracer penetration across the intergyre boundary [sections 2a and 2c(1)]. Approximately diffusive single-particle dispersion is found in transition zones (green color in Fig. 8) covering about 15% of the upper ocean. Even there, however, the single-particle dispersion often shows an oscillatory component (Fig. 9b) that cannot be described by a classical diffusion process. Along all the boundaries except the western one (where the local homogeneity assumption does not permit resolving near the boundary), we see subdiffusive behavior in the normal direction (Figs. 8, 9d), and superdiffusive behavior in the tangential direction indicative of alongshore fluctuating jetlike currents. Typically the tangential component of the single-particle dispersion decreases after about 150 days, marking the end of superdiffusive growth.

In the deep ocean there is a clear separation between the zonally subdiffusive eastern part of the basin and the diffusive and superdiffusive regions in the western basin and near the northern and southern boundaries (Fig. 8c). The subdiffusive region acts as a transport barrier for particles arriving from the west (cf. Fig. 5). The meridional single-particle dispersion is subdiffusive almost everywhere (Fig. 8d). The zonally superdiffusive strip along the northern boundary matches the enhanced eastward tracer spreading in Figs. 6c,d. The eastern-jet pathway for tracer penetration from the subpolar to the subtropical gyre (Figs. 6c,d) matches the meridionally and zonally superdiffusive strips in the northern and southern sectors of the subpolar WBC (Figs. 8c,d).

The intermediate-time, single-particle dispersion properties can be seen in individual trajectories (Fig. 12). In the subdiffusive subtropical WBC and its eastward jet extension, trajectories have loops due to particle recirculations within strong vortices. In the superdiffusive regions, trajectories exhibit relatively frequent and long intervals of “flight.” In the subdiffusive eastern basin, trajectories show small-amplitude wiggles around the time-mean streamlines, particularly where the fluctuations are weak.

The intermediate-Re solution is characterized by weaker fluctuations and much longer correlation times: *T*_{ii} is 100–150 days in the eastern basin and deep ocean, 5–10 days in the upper-ocean subtropical WBC and its eastward extension, and 30–50 days near the eastern boundary. The intermediate-time, single-particle dispersion regime lasts longer (e.g., we use 1000 instead of 200 days as the upper limit for intermediate-time interval), but the general pattern of *α*_{ii}(*x, **y*) is similar. Only two substantial differences occur: the primary zonally superdiffusive region is now centered around the WBC eastward jet extension, and the WBC itself is weakly superdiffusive in the meridional direction. Both effects are due to the absence of strong meanders and rings that enhance negative lobes of *R*(*τ*). Thus, in the single-particle dispersion, as in other measures discussed below, the transport rates do vary with Re but the qualitative transport behavior does not, as long as the eddy field is sufficiently active. This is contrasted in section 3 with the small-Re regime.

### c. Tracer mixing

In this section we focus on the process of mixing, which is related to the decorrelation rate of neighboring particles. Mixing is analyzed by three measures: finite-time Lyapunov exponent, which is more suitable for short-time description because it automatically reduces the mixing information to a simple exponential fit; two-particle dispersion, which is robust measure on all timescales; and fractal correlation dimension, which is a more large-scale measure.

#### 1) Finite-time Lyapunov exponents

In general, large and small values of the local, finite-time Lyapunov exponents indicate regions of enhanced and weakened mixing rates, respectively. Among several techniques for estimating the exponents, we choose the one which, although it requires calculating spatial derivatives of velocity, is very robust (Pierrehumbert and Yang 1993). We define the deformation tensor,

along a Lagrangian trajectory, **X**(*t*). The equation for small displacements, **x̂**(*t*), along this trajectory is

The displacement tensor, *M*_{ij}(*t*), evolves according to

where Σ_{k} denotes *k* summation and *δ*_{ij} is the Kroneker symbol. Because the 2D flow is nondivergent, *M*_{ij}(*t*) has a pair of real eigenvalues with opposite signs. The positive eigenvalue, *λ̃*, is the rate of stretching in the direction of the corresponding eigenvector, and the negative one is the rate of contraction in the perpendicular direction. The finite-time Lyapunov exponent at time *t* = *T* is defined as

and it expresses the net Lagrangian straining as an exponential growth rate for small displacements around **X**(*t*).

To calculate *λ* within each fluid layer, we integrate (4) for 200 × 200 uniformly distributed particles and solve (12), (14), and (15). Figure 13a illustrates a typical structure for *λ*(*T, **x, **y*) at *T* = 2 days for a deployment time corresponding to Fig. 1b, and Fig. 13b shows the corresponding histogram. Large values of *λ* (between 0.005 and 0.06 day^{−1}) are found in the subtropical WBC, in both eastward jets, and around edges of rings that detach from them. These are regions of efficient mixing. Small values of the finite-time Lyapunov exponent (*λ* < 0.01 day^{−1}) characterize regions with weak mixing, such as away from the swift currents where planetary waves dominate, between the eastward jets, and in the cores of rings (e.g., Provenzale 1999), which are also regions with slow two-particle dispersion rates [sections 2c(2)]. The subpolar eastward jet has relatively large *λ,* and it is a strong meridional transport barrier [sections 2a and 2d(1)]. This suggests that most of the mixing occurs on either side of the jet rather than across it (e.g., Samelson 1992). The ensemble-mean, finite-time Lyapunov exponent, 〈*λ*(*T, **x, **y*)〉, is calculated by averaging over many deployments. It is spatially smooth, with maxima located along the main currents and with averaged out information about transient trapping events inside rings (Fig. 14). In the deep ocean and at smaller Re, *λ* is qualitatively similar, and its amplitude decreases both with greater depth (at large Re it is 5 times smaller near the bottom) and with smaller Re.

#### 2) Two-particle dispersion

The two-particle (or relative) dispersion,

is the ensemble average of the variance of the relative distance between two particles, and it provides a description of the mixing process (i.e., distortion and homogenization) of the fluid. On intermediate times *D*^{R}(*t*) is a more robust and complete measure of mixing than *λ* because it accounts for deviations from exponential relative dispersal of particles. In the energy-cascading inertial range of isotropic, stationary, homogeneous 3D turbulence a dimensional analysis indicates that *dD*^{(R)}/*dt* ∼ [*D*^{(R)}]^{2/3}, or equivalently, *D*^{(R)} ∼ *t*^{3} (Obukhov 1941); this relationship is referred to as the Richardson law. In homogeneous 2D turbulence, there is an energy cascade toward large scales and an enstrophy cascade toward small scales. In the enstrophy-cascading inertial range, analogous dimensional arguments show that *dD*^{(R)}/*dt* ∼ *D*^{(R)}; hence for each principal direction:

where *θ*_{ii} is the enstrophy cascade timescale (Lin 1972). This regime is observed in the atmosphere with *θ*_{ii} about 5–6 days (Morel and Larcheveque 1974; Er-El and Peskin 1981).

We calculate *D*^{(R)}_{ii}(*t*) for the same dataset described in section 2b. Here *D*^{(R)}(*t, **x, **y*) distinguishes between the western and eastern regions in the basin, as found in ocean observations (LaCasce and Bower 2000); the former includes the region near the southern boundary, and the latter includes most of the subpolar gyre (Fig. 15). In the western region, *D*^{(R)}_{ii}(*t*) is nearly isotropic. It increases rapidly and exponentially at short times. The growth rate gradually decreases at around 30 days as the particles spread over a substantial area. In the eastern region, *D*^{(R)}_{ii}(*t*) is markedly anisotropic, with larger meridional rather than zonal growth rates, and it shows an exponential regime that persists for several *θ*_{ii}. Here *θ*_{ii}(*x, **y*) is estimated from a fit over 20 days after the deployment. Its magnitude varies widely, from several days in the subtropical WBC to several hundred days near the eastern boundary (Fig. 16). It is nearly isotropic in the western and anisotropic in the eastern region. The small (large) values of *θ*_{ii} indicate rapid (slow) mixing process in the western (eastern) regions, and they are consistent with the degree of filamentation in the particle distributions from localized releases (section 2a, Fig. 4), with the distributions of finite-time Lyapunov exponents [section 2c(1)], and with the homogenization rates [section 2c(3)]. Exponential fit (17) requires both longer time intervals and binning the basin (as in section 2d); therefore the geographical map of *θ*_{ii} (Fig. 16) is smoother than that of 〈*λ*〉 (Fig. 14).

In the western basin, *D*^{(R)}_{ij} has similar behavior at all depths, but in the eastern basin *θ* is two to three times larger in the deep rather than upper ocean. At intermediate Re, *θ*_{ii}(*x, **y*) are similarly distributed, but they are larger by an order of magnitude, and this tendency continues into the small-Re case.

#### 3) Homogenization

For each of the localized releases described in section 2a, we calculate the two-particle correlation function (Pierrehumbert 1991). Given *N* particles, there are *N*(*N* − 1)/2 pairs, and the correlation function, *H*(*t, **r*), is defined as the number of pairs separated by a distance less than *r* at time *t.* If *H*(*r*) ∼ *r*^{d}, then *d* is the fractal correlation dimension that measures the degree of homogenization at scale *r.* In a 2D domain, *d* may vary in the interval 0 ≤ *d* ≤ 2. In an unbounded domain, a well-mixed (i.e., homogenized) ensemble of particles has *d* = 2, an ensemble organized in filaments has *d* ≈ 1, and an infinitesimal ensemble has *d* = 0. In the present problem the domain is bounded, and at large time *H*(*t, **r*) converges to *H*_{∞}(*r*).

For localized particle releases complete mixing occurs over several thousand days (by an order of magnitude faster than at small Re; section 3). The homogenization rate is faster in the subtropical rather than in the subpolar gyre (Fig. 17), and in the western rather than in the eastern basin. In the subpolar gyre, the tracer is substantially concentrated in filaments even after 2000 days (cf. Fig. 4), which is consistent with the slow two-particle dispersion rates induced by the planetary waves [section 2c(2)] and with the long residence time allowed by slow time-mean advection. At all depths, particles that cross the intergyre boundary from the northern gyre to the southern [section 2d(1)] are rapidly homogenized there, as indicated by the steeper slope of *H*(*t, **r*) at large *r* and late time (Fig. 17b). In the eastern basin, the homogenization rates are substantially slower in the deep rather than upper ocean, but in the western basin the change with depth is small.

### d. Large-scale transport

Recently, several theoretical works as well as oceanic float observations focused on the Lagrangian transport across the meandering eastward jet extension of the WBC (section 1b). This transport is thought to be one of the important aspects of global meridional heat flux. The key questions that have been asked so far are whether the eastward jets are transport barriers and whether the cross-jet transport is more inhibited in the upper rather than deep ocean. Our work differs from the previous studies by using the circulation regime that is more realistic in terms of the model, dynamics, and much larger Re than before. In the double-gyre circulation at large Re, the WBCs separate in two distinct jets (i.e., as in the North Pacific), and only the weakly meandering subpolar jet, which is the intergyre boundary, is a strong barrier to the transport. On longer than intermediate times, local net Lagrangian fluxes are not completely defined by local *α*_{ii} because they incorporate transport information from much larger areas and the time-mean currents. In the following sections we look separately at mean Lagrangian fluxes across each jet.

#### 1) Flux across the intergyre boundary

Localized releases of particles (section 2a) show that transport across the subpolar eastward jet is inhibited, in agreement with idealized models of a weakly meandering zonal jet (section 1b). To quantify the flux across the intergyre boundary, we release *N* = 10^{4} particles randomly distributed over the basin at 100 different deployment times and follow them for 500 days. A record is kept of where each particle first crosses the intergyre boundary (marked in Figs. 3–6) and the initial and final positions of particles that originated in one gyre and are found at time *t* in the other gyre. The ensemble-mean intergyre flux in the *i*th layer, averaged over the time *t* since deployment, is

where *V*_{i} = *L*^{2}*H*_{i}/*N* is the fluid volume corresponding to each of *N* 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 is from the northern or southern gyre. The *n* and *s* components of *F*_{i} are equal because of time-mean mass conservation within each layer and gyre. The intergyre boundary is an isoline of time-mean *ψ*_{1}(*x, **y*), thus there can be no net Eulerian flux across it. However, *N*^{(n)}_{i}(*t, **x*) ≠ *N*^{(s)}_{i}(*t, **x*) because locally a net intergyre Lagrangian flux, as defined above, can occur.

The percentage of particles that have crossed the intergyre boundary over time *t* and the values of the time-average flux, *F*_{i}(*t*), are listed in Tables 1 and 2. During 500 days *F*_{i}(*t*) changes smoothly, and *N*^{(n,s)}_{i}(*t, **x*) do not change qualitatively, indicating that their definitions are time dependent but robust. Tables 1 and 2 show that the intergyre flux strongly increases with Re, especially in the deep ocean. In the large-Re solution, the percentage of crossing particles decreases with increasing depth, but the deep-ocean transport remains large due to large layer thicknesses.

In the upper layer, the particles that cross from the south are released either just to the south of the intergyre boundary in the western region between the eastward jets or in the vicinity of the subtropical WBC (Figs. 18a,b). Upper-layer particles crossing in the opposite direction originate in the far northern part of the subpolar gyre (Figs. 18c,d). In the lower layers, the particles that cross come from and stay within a narrow strip near the intergyre boundary (Figs. 18e,f). At intermediate Re the crossing particles are more confined to near the western boundary.

The *N*^{(n,s)}_{i}(*t, **x*) distributions have strong maxima near the western boundary and several smaller maxima along the subpolar jet (Fig. 19). These are the primary intergyre transport pathways. Substantial differences between *N*^{(n)}_{i}(*t, **x*) and *N*^{(s)}_{i}(*t, **x*) occur in the western part of the basin (Fig. 19c); they correspond to net Lagrangian fluxes in several narrow currents in opposite directions. Enhanced, superdiffusive transport rates along lateral boundaries do not show up in *N*_{i}(*x*) because the velocity variance decreases to zero on the walls, but they do appear in normalized *Ñ*_{i} as it is defined below. The smaller *N*_{i} values in the interior are associated with the subdiffusive, single-particle dispersion found there (section 2b). Crossing particles in the subdiffusive region tend to come from very near the intergyre boundary and arrive at very near the intergyre boundary (Fig. 18).

Meridional single-particle dispersion rates along the intergyre boundary are subdiffusive both in the upper and deep ocean (Fig. 8), and the corresponding velocity variances change with depth (Fig. 2)—both factors make it hard to quantify how the flux efficiency changes with depth. We nondimensionalize the fluxes in such a way that their efficiencies may be compared. Local depth-average meridional rms velocity is

and its integral along a contour *C* (*y* = *y*(*x*), *x*_{1} < *x* < *x*_{2}) is

The normalized crossing distributions along the contour are

and our *C* is the intergyre boundary. The ratio of the deep-ocean *Ñ*_{i}(*t, **x*) to the upper-ocean *Ñ*_{1}(*t, **x*) is larger than unity along most of the eastward jet extension of the subpolar WBC (Fig. 20), suggesting that the cross-jet flux is more efficient in the deep rather than upper ocean, although the upper-ocean flux is more intensive (Table 1). That is also true for the transport across the subtropical eastward jet [section 2d(2)]. Enhanced efficiency of the deep-ocean flux may be explained by the ability of fluid particles to preserve their potential vorticity (PV) and by the presence of the sharper PV front in the upper ocean, but this hypothesis has yet to be verified by calculating PV balances following Lagrangian particle trajectories.

#### 2) Flux across the subtropical eastward jet

The eastward jet extension of the subtropical WBC is characterized by intense meandering and large-scale interannual variability associated with changes in the position and intensity of the jet core and its neighboring recirculation currents (Berloff and McWilliams 1999a). This region has the strongest velocity fluctuations (Fig. 1) and the most intensive transport. We calculate the cross-jet flux by the method presented in the previous section.

To define the subtropical jet boundary, we choose two segments of the time-mean, upper-layer streamline in the core of the jet since no streamline in the jet originates and ends on the eastern or western boundaries. Both segments start 60 km away from the western boundary and terminate in the basin interior at either *x* = 1300 or 1900 km. The cross-jet fluxes are listed in Table 3 for *t* = 20, 40, 60, 80, and 100 days using the shorter of the two segments. The flux values across the longer segment are larger by about 20% since most of the flux occurs within ∼600 km in the western sector where *N*^{(n,s)}_{i}(*t, **x*) is larger.

In the upper layer most of the crossing particles are near the time-mean core of the jet, which also includes the WBC (Fig. 21). About 40% of the southward flux is due to particles that arrive from a broad region near the subtropical WBC (Fig. 21c) and move northward through the narrow gap between the western boundary and the beginning of the jet segment; then cross the jet from the north and return to the south (Fig. 21d). (The analogous percentage in the deep ocean is much lower, about 10%.) The compensating flux (Figs. 21a,b) consists of particles that originate in the near-boundary region of the subtropical WBC, are ejected from the WBC by vigorous meanders and eddies (Berloff and McWilliams 1999b), and then cross the jet in the northward direction. As a result of stronger fluctuations, the cross-jet flux (Table 3) is 1.3–1.5 times larger than the whole intergyre flux (Table 1), although both boundaries are locally subdiffusive in meridional direction (section 2b). Away from the western boundary, both eastward jets are qualitatively similar in terms of enhanced cross-jet flux efficiency in the deep rather than upper ocean (Fig. 20), however the subtropical jet extends farther eastward.

Experiments with particles initialized near the Gulf Stream separation point suggest that such particles more efficiently cross the current in the deep rather than upper ocean (Bower and Lozier 1994). A hypothesis is that this indicates efficient deep-ocean transport pathway below the upper-ocean transport barrier coinciding with the jet, however the floats deployed in a small single location represent only a fraction of the total cross-jet transport. We deploy floats near the WBC separation point, in the core of the current (Fig. 22), to examine how the cross-jet flux is inhibited in the upper ocean compared to the deeper one (cf. section 1b). There are synoptic situations when the upper-ocean trajectories remain within the instantaneous eastward jet core (Fig. 22a), but in other cases some particles continue to follow the core while many others escape either in rings or by other fluctuations (Fig. 22b). Velocities generally decrease with depth, therefore deep particles close to the release point lag far behind upper-ocean ones and rapidly disperse (cf. Figs. 3 and 5). For a more synchronous comparison, we release deep particles at the instantaneous (*x, **y*) positions of the upper-ocean particles as they move along the eastward jet (Figs. 22c,d). This ensemble shows a larger degree of dispersal away from the instantaneous jet, compared to the upper-ocean ensemble. Thus, we do see enhanced deep-ocean crossing of the instantaneous jet by particles deployed near the jet core but, as it is demonstrated in this section, that does not indicate that the total cross-jet flux is more intensive there.

#### 3) Recurrence time

The Poincare recurrence time, *τ*_{r}, is the most diverse transport measure among those used in this paper: it combines information about the large-scale transport [alternatively to the methods in sections 2d(1) and 2d(2)], spreading, and mixing aspects. It is defined as the time between a particle leaving from and returning to a subdomain, *A.* The Poincare distribution function is *P*_{r}[*τ*_{r}(*x, **y*)] (Zaslavsky 1999). In a bounded domain with nondivergent 2D flow, *P*_{r} is well defined for any *A,* and the mean recurrence time,

is finite (Kac 1957). Subdomains with short 〈*τ*_{r}〉 are associated with flow structures that trap the fluid (i.e., if a particle escapes *A,* it quickly returns), and long 〈*τ*_{r}〉 indicates dynamical flights.

We calculate *P*_{r}[*τ*_{r}(*x, **y*)] following 10 000 randomly initialized particles. The *A* are defined by dividing the basin into 40 × 40 bins. Figures 23a,b show 〈*P*_{r}〉 averaged over all bins. At large *τ*_{r}, 〈*P*_{r}〉 decays with a slope steeper than −1, indicating that the integral in (22) converges. In each layer 〈*P*_{r}〉 has a maximum at 30–40 days, the timescale of dominant planetary waves (cf. section 2b). The relatively small values of 〈*τ*_{r}(*x, **y*)〉 in the eastern basin and in the deep ocean (Figs. 23c,d) correspond to this type of motion, and the weak local minima in the upper WBCs and central subtropical gyre indicate the presence of fluid-trapping vortices (section 2b). The upper-ocean 〈*P*_{r}〉 has a minimum at around 500 days and a broad secondary maximum in the interannual range, but in the deep ocean it decays monotonically with *τ*_{r}. This secondary maximum is associated with particles that migrate into the opposite gyre and remain there for a long time. The origins of such particles correspond to large 〈*τ*_{r}(*x, **y*)〉 in the WBCs and near-wall regions upstream of the WBCs (Figs. 23c,d; see also section 2d(1) and Figs. 18a–d). In the deep ocean, 〈*τ*_{r}(*x, **y*)〉 is generally smaller because of the larger contribution by planetary waves. In both layers, there is a local maximum in 〈*τ*_{r}(*x, **y*)〉 at the end of the subtropical eastward jet (Figs. 23c,d) that is related to particles injected into the eastern basin and slowly carried away by the gyre advection. Here 〈*τ*_{r}〉 increases as Re decreases, due both to the general weakening of fluctuations and to the increasing timescale of the planetary waves, but the qualitative shape of *P*_{r}[*τ*_{r}(*x, **y*)] is not Re dependent.

## 3. Transport at small Re

At a marginally supercritical value of Re, the steady gyre circulation (Fig. 27a) becomes unstable in a Hopf bifurcation whose fluctuations have appreciable amplitude only in the upper ocean.^{4} Here the Eulerian flow is time periodic (with period *T* = 221 days), but Lagrangian trajectories are generally chaotic (e.g., Ottino 1989). In this small-Re solution the transport rates are, of course, smaller than at large Re because the velocities are smaller, both the time-mean and the fluctuating components. This effect can be factored out, to some degree, by examining evolving particle distributions from localized releases over longer time intervals (Fig. 24).

For a relatively long time after localized particle releases (i.e., for times of the order of *T*; Figs. 24a,c,e), the particle distributions show a significant amount of spatial fine structure, indicative of weak mixing relative to spreading, and in contrast to analogous large-Re distributions (cf. Figs. 3–6). This fine structure is associated with the invariant manifolds connected to a few saddle points of the flow (see below). However, over longer times the fine structure wanes and the particles become rather well mixed throughout the gyre in which they are released (Figs. 24b,d,f). As in the large-Re regime, intergyre fluxes are inhibited at all timescales and spreading occurs more rapidly in the subtropical (Figs. 24a–d) rather than subpolar gyre (Figs. 24e–f). Particles released in the outer parts of the WBCs cross into the other gyre more often than particles released in the inner sites (Figs. 24c–f). This distinction is more pronounced at small rather than large Re because, in the first case, fluctuations in the WBC are much weaker than in the offshore jet. Thus, over intermediate time intervals the evident fine structure for small Re is a manifestation of transient transport barriers related to the invariant manifolds, which are not pronounced at large Re.^{5} But over long time intervals the large-scale transport barriers—between the inner and outer parts of gyres and between the subpolar and subtropical gyres—are qualitatively similar in the different Re regimes. At small Re the intermediate-time single-particle dispersion, *D*^{′}_{ii}(*t*) (not shown), is more often subdiffusive than at larger Re (section 2b); a subdiffusive regime is recovered everywhere in the basin except for a small superdiffusive region around the WBC separation point and the first offshore meander of the eastward jet (Fig. 27). The latter region is responsible for a big portion of the intergyre flux (see below). Homogenization and two-particle dispersion rates, as measured by *H*(*r*) and *θ*_{ii}(*x, **y*), are very small at small Re, but their spatial patterns are qualitatively similar to those calculated for the large-Re solution [sections 2c(2)–(3)].

The simplicity of the flow at small Re permits us to analyze the intergyre flux not only using the statistical methods described in section 2d(1), but also with the geometrical method derived from dynamical system theory (see appendix). The initial and final positions of particles crossing the intergyre boundary (Fig. 25; cf. Figs. 18a–d) indicate that the northward flux comes from a narrow strip of fluid from the outer part of the subtropical gyre, and the southward flux brings material from the outer part of the subpolar gyre and then delivers about two-thirds of it directly into the central part of the subtropical gyre and the rest into the outer part near the intergyre boundary. This conclusion is also achieved with the geometrical analysis. The crossing probability distribution (Fig. 26) shows that, unlike at large Re (Fig. 19a), the intergyre flux is strongly confined near the western boundary and in a narrow region around *x* = 0.1*L* in the offshore jet (cf. Fig. 28). The sequence of net Lagrangian currents with opposite signs at large Re (Fig. 19c) is replaced with only two currents at small Re (Fig. 26). The associated time-average fluxes are *F*_{1}(*T*) = 0.63, *F*_{1}(2*T*) = 0.49, and *F*_{1}(3*T*) = 0.44 Sv (Sv ≡ 10^{6}m ^{3} s^{−1}). These values are systematically larger than the 0.39 Sv obtained by the geometrical method for the chosen primary point on the invariant manifolds. These estimates differ because of the difference between the positions of the Eulerian time-mean and the invariant-manifold inter-gyre boundaries, where the latter oscillates substantially around the former in both space and time. We conclude that both methods yield similar information about the intergyre transport. The statistical method, however, allows for an easier calculation of the intergyre transport because of the complexity in tracking the transient geometry of the invariant manifolds, even for relatively simple flows. On the other hand, the statistical method gives a less detailed characterization of the relevant transient currents.

## 4. Conclusions and discussion

We analyze the Lagrangian transport (i.e., spreading and mixing of advected material or, equivalently, passive tracer) in an idealized double-gyre, wind-driven, quasigeostrophic model of the midlatitude ocean in statistical equilibrium forced by a steady wind. We focus on the large-Re solution (Re ≈ 100), although the Re dependence of the transport properties is explored in two other cases with Re ≈ 12 and 6. The spreading (i.e., propagation) of material from its initial location is described by the evolving ensemble-mean distribution function and by the single-particle dispersion, *D*(*t*). The mixing process (i.e., disappearance of the patchiness in individual realizations) is described by the finite-time Lyapunov exponent, *λ,* the two-particle dispersion, *D*^{(R)}(*t*), and the homogenization rate. Finally, the large-scale transport is explored by analyzing the fluxes across the intergyre boundary, which is the principal transport barrier, and the subtropical eastward jet, and by calculating the recurrence time distributions.

Using the local homogeneity assumption, we geographically partition the basin in terms of the intermediate-time exponents of the power-law fit, *D*′(*t*) ∼ *t*^{α}. In the short-time limit, the spreading is always ballistic (i.e., *D*′ ∼ *t*^{2}), but the diffusive regime (i.e., *D*′ ∼ *t*) that occurs in homogeneous and stationary turbulence in an unbounded domain (Taylor 1921) is often irrelevant because the local homogeneity assumption is violated long before that regime is reached. The *α*(*x, **y, **z*) show widespread anisotropy and extensive regions of both subdiffusion (*α* < 0.8) and superdiffusion (*α* > 1.2) due to coherent structures such as mesoscale vortices, meanders of swift currents, and planetary waves. The nondiffusive behavior is due to the persistence of intermediate-time Lagrangian velocity correlations described by the autocorrelation function, *R*(*τ*), which slowly decays and strongly oscillates. The subdiffusive (superdiffusive) regime is associated with pronounced first negative (second positive) lobe of *R*(*τ*). Regions away from the WBCs with their eastward jet extensions and from the lateral boundaries are typically subdiffusive because the local dynamics is dominated by planetary waves that induce oscillatory trajectories. Another subdiffusive region is in the subtropical WBC and its eastward jet extension where material is trapped inside strong vortices (rings) generated by the current (e.g., Provenzale 1999); the vortex cores are also characterized by slow mixing rates associated with small values of *λ.* Superdiffusion occurs in the currents along the boundaries and in the central part of the subtropical gyre in the zonal direction. In the deep layers and in the solutions at intermediate and small Re, the subdiffusive behavior is more pronounced and only small regions are characterized by superdiffusivity. A geographical partitioning analogous to the one provided by the single-particle dispersion is obtained from the recurrence time, *τ*_{r}(*x, **y, **z*), which is shorter where the subdiffusion occurs.

Thus, the horizontally isotropic eddy diffusion process, which provides the almost universally used in oceanic GCMs subgrid-scale parameterization for material transport by mesoscale eddies, is not confirmed in our eddy-resolving model and needs reconsideration [similar conclusions are made by Figueroa (1994) and Griffa (1996) among other works]. New transport models and the subgrid-scale parameterizations based on them have to account for gradually decaying and oscillating *R*(*τ*) and for associated deviations from the diffusion regime. One approach is to use an hierarchy of stochastic transport models that, on the simplest level, yields the diffusion model but, on more sophisticated levels, takes into account nondiffusive spreading behaviors associated with slowly decaying, oscillating, and asymmetric *R*(*τ*) due to the presence of coherent fluid structures (Berloff and McWilliams 2002). Transport models that, in addition to temporal correlations, account for spatial correlations between neighboring Lagrangian particles, and therefore simulate not only the spreading but also the mixing process, are to be developed as well (e.g., Sawford 2001).

Analyzing the mixing process, we find that *D*^{(R)}(*t*) grows exponentially over a time interval of several weeks in the WBCs and years near the eastern boundary. The exponential growth is in agreement with a dimensional analysis of 2D turbulence in the enstrophy-cascading inertial range (Lin 1972). In terms of the exponential growth timescale, *θ*(*x, **y, **z*); the finite-time Lyapunov exponent, *λ*(*x, **y, **z*); and the recurrence time, *τ*_{r}(*x, **y, **z*), the gyres can be partitioned into well-distinguished western and eastern regions. In the west *θ* is appreciably shorter, *λ* larger, and *τ*_{R} longer, indicating faster mixing and more effective dispersing by strong vortices than by weak planetary waves, consistent with oceanic mixing (LaCasce and Bower 2000). In general, the homogenization of tracer proceeds from large to small scales and it is faster in the upper rather than deep ocean; a fully mixed tracer field is achieved only after thousands of days.

The weakly meandering eastward extension of the subpolar WBC is a long-time barrier to the intergyre transport because the time-mean flow does not cross it, the local meridional single-particle dispersion is subdiffusive, and the local velocity variance is moderate compared to the other, subtropical, eastward jet. Most of the fluid participating in the upper-ocean, intergyre fluxes arrives to the intergyre boundary through the WBCs and from remote parts of the basin. This behavior is also manifested by large values of *τ*_{r}(*x, **y*). The subtropical eastward jet is also a transport barrier because of its subdiffusive behavior but that is compensated for by local, energetic flow fluctuations; therefore, the barrier property is not pronounced. Fluxes across both jets are more intensive in the upper ocean and at larger Re, as a result of stronger fluctuations, but they are more efficient (as it is evident after normalizing by the local velocity variance) in the deep ocean, in agreement with the weaker potential vorticity front there (Bower and Lozier 1994). The cross-jet flux calculations presented are to be extended by finding the fluxes across instantaneous rather than time-mean cores of the jets. That extension is both conceptually and technically challenging because the instantaneous core is quite transient and sometimes incoherent.

At small Re, in addition to smaller overall transport rates because of weaker fluctuations, the fluid in the central parts of gyres is relatively isolated from the gyre peripheries. However, in our solutions with strong gyre asymmetry, the intergyre transport barrier is periodically broken by an eddy that penetrates deeply into the subtropical gyre. The simplicity of the flow in this regime permits analysis of the intergyre transport by a geometrical method of invariant manifolds (e.g., Beigie et al. 1994). The geometrical approach is more detailed but comparable to more common statistical approaches. However, the feasibility of identifying the invariant manifolds at large Re is problematic.^{6}

The greatest limitation of the present theoretical study is its dynamical simplification and geographical idealization in relation to nature. Therefore, the following extensions of the present work are needed. Wind fluctuations on synoptic, seasonal, and climatic scales will change Lagrangian velocity correlations and, therefore, will influence the spreading and mixing rates. Disturbances induced by the bottom topography and irregular boundaries will likely modify the transport, especially at depth and near coasts. Also, it is unknown to what extent the presence of bottom friction influences abyssal dynamics and transport. The transport at even larger Re may be modified by the increased vortex population (Siegel et al. 2001). Instabilities arising in the WBCs are likely exaggerated in the double-gyre model, because commonly used Newtonian friction parameterization does not well represent unresolved coastal dynamics. Also, oceanic gyres have finite Rossby number, hence dynamical behaviors beyond quasigeostrophy, and the transport by 2D nondivergent geostrophic velocities is only part of the 3D transport.

## Acknowledgments

We acknowledge S. Kravtsov and A. Provenzale for many interesting discussions. A. B. thanks UCLA for their hospitality. Funding was provided by NSF Grant OCE 96-33681 and by ONR Grant N00014-98-1-0165.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Geometrical Analysis of Transport at Small Re

The upper-ocean, small-Re flow (Fig. 27a) has several fixed points (i.e., points where **u** = 0), which are either hyperbolic (saddle points) or elliptic (centers) depending on whether linearization of the velocity field around the point yields eigenvalues with nonzero or zero real part, respectively. The generalized Poincare–Hopf index formula for a 2D nondivergent flow in a simply connected domain (Ma and Wang 1998) states that

where Ce, Sa, and Bo are the numbers of centers, interior saddles, and boundary saddles, respectively. In our situation, the values of Ce and Sa vary with time (with Ce = 7 and Sa = 5 the maxima), but the value of Bo is fixed (Bo = 2) with one saddle point continuously moving on each boundary, western and eastern. The transport analysis here is based on finding the template of invariant manifolds for a time-periodic flow (e.g., Beigie et al. 1994; Coulliette and Wiggins 2001). The stable (unstable) invariant manifold of a hyperbolic point is given by the curve containing the full set of Lagrangian particles that approach the point as *t* → +∞ (*t* → −∞). In steady flow the stable and unstable manifolds of the same or different points coincide and form the homoclinic and heteroclinic manifolds, respectively, which separate regions of fluid. In unsteady flow (as in our case), the manifolds split apart and intersect each other infinitely many times. A simply connected patch of fluid bounded by segments of the manifolds is called a lobe. For a pair of hyperbolic points (here the boundary saddles), the intergyre transport is associated with the unstable manifold that emanates from the western point and intersects infinitely many times with the stable manifold of the eastern point. We designate a particular intersection point as the primary one and partition the domain with a line that consists of the unstable manifold segment from the western to the primary point and the stable manifold segment from the primary to the eastern point. This choice of a primary point is, of course, nonunique, and the associated transport across its partition line differs with different choices. The two lobes on each side of the primary intersection point are called the entrainment and detrainment turnstile lobes and because the flow is nondivergent, the lobe areas are proportional to the fluid volumes inside them, and therefore yield the cross-partition transport over an oscillation period. When the time-dependent perturbations are small, the area of the lobes may be estimated analytically (Melnikov 1963), but in general it is done computationally.

The unstable (stable) manifold of the western (eastern) boundary point (Fig. 27b) is calculated by integrating particle trajectories forward (backward) in time. A first Lagrangian particle is kept on the fixed point, a second particle is initialized nearby, and new particles are added as the manifold stretches and folds. When the distance between particles *n* and *n* + 1 becomes larger than small distance, *d*_{max}, the manifold is interpolated cubically, and a new particle is added in between. If the distance between particles *n* − 1 and *n* + 1 becomes smaller than *d*_{min} < *d*_{max}, the *n*th particle is removed. The manifolds (Figs. 27b, 28) are always plotted at the same phase of the time-periodic oscillation, and their geometrical structure shows that fluid participating in the intergyre (i.e., cross-partition) transport arrives from the WBCs and then penetrates either into the center of the subtropical gyre or into the outer part of the subpolar gyre near the intergyre boundary (Fig. 27b). A closer view on the evolution of the turnstyle lobe associated with a primary point near the western boundary (Fig. 28) reveals the following: the area labeled A is occupied by fluid entrained in the subtropical from the subpolar gyre; the fluid in B (two small fragments) and C arrives from and remains in the subtropical and subpolar gyres, respectively; and the fluid in D moves in opposite direction to A. About two-thirds of fluid in A is carried by the strong eddy in the partition line that deeply penetrates inside the subtropical gyre, and about one-third of it, as well as fluid in D and C, is carried by the outer part of the subpolar gyre. This process is associated with substantial stretching and folding of material lines (Fig. 28). The deeply penetrative, intensive transport induced by the eddy does not occur (hence the intergyre transport is much smaller) when the circulation is more symmetrical because of the wind symmetry (Coulliette and Wiggins 2001). The lobe areas are calculated by the Monte Carlo method (Press et al. 1992). We find that the area A, equal to D, is about 2.475 × 10^{4} km^{2}, which, multiplied by *H*_{1} and divided by the period *T,* yields the transport rate of about 0.39 Sv (a small number compared to the mean wind-driven Sverdrup transport).

We find that the lobe analysis, successfully applied here and in other small-Re situations (e.g., Poje and Haller 1999; Coulliette and Wiggins 2001), becomes very complicated and arguably impractical at large Re, for the following reasons. First, the fixed points are very mobile, and tracking them requires computationally expensive root- or extrema-searching algorithms. Second, the lifetime of a transient saddle can be so short that the associated manifolds do not have time to accomplish much transport. Third, in turbulent flows the invariant manifolds experience strong stretching and folding, hence calculating them is a very delicate and computationally expensive task. Fourth, Lagrangian floats limited to the manifolds yield less accurate and efficient estimates of other Lagrangian transport properties, such as one- and two-particle dispersions. Finally, in temporally complex flows many realizations of individual events and their associated manifolds must be analyzed to determine mean Lagrangian fluxes.

## Footnotes

*Corresponding author address:* Pavel S. Berloff, Institute of Geophysics and Planetary Physics, University of California, Los Angeles, CA 90095-1567. Email: pavel@atmos.ucla.edu

^{1}

Substantial asymmetry of the wind forcing is found crucial for achieving generic dynamical behavior of the double-gyre circulation (Berloff and McWilliams 1999a).

^{2}

We advect test particles with frozen instantaneous streamfunction fields and find that the particles closely follow streamlines, even where the velocity gradients are large.

^{3}

The overline applied to localized Lagrangian quantities such as *R*(*τ, ***x**) and *D*(*t, ***x**) involves ensemble averaging over a set of Lagrangian trajectories.

^{4}

Because of the approximate confinement of mean and fluctuating currents to the upper ocean, a three-layer solution is essentially indistinguishable from the analogous, slightly supercritical, 1.5-layer solution analyzed in Berloff and McWilliams (1999a).

^{5}

In almost every 2D nondivergent flow there are instantaneous saddle points and connected manifolds that can be related to nearby particle trajectories over a short-time interval. However, in a temporally complex flow these structures have quite different geometries at different times and therefore do not imprint their pattern on the mean particle distributions.

^{6}

Ide et al. (2001) propose an extension of the geometrical method for time-dependent flows.