**Research article**
05 Apr 2019

**Research article** | 05 Apr 2019

# Competition between chaotic advection and diffusion: stirring and mixing in a 3-D eddy model

Genevieve Jay Brett Larry Pratt Irina Rypina and Peng Wang

^{1},

^{2},

^{2},

^{3}

**Genevieve Jay Brett et al.**Genevieve Jay Brett Larry Pratt Irina Rypina and Peng Wang

^{1},

^{2},

^{2},

^{3}

^{1}IPRC, University of Hawai`i at Mānoa, Honolulu, HI, USA^{2}Woods Hole Oceanographic Institution, Department of Physical Oceanography, Woods Hole, MA, USA^{3}Department of Atmospheric and Oceanic Sciences, University of California, Los Angeles, CA, USA

^{1}IPRC, University of Hawai`i at Mānoa, Honolulu, HI, USA^{2}Woods Hole Oceanographic Institution, Department of Physical Oceanography, Woods Hole, MA, USA^{3}Department of Atmospheric and Oceanic Sciences, University of California, Los Angeles, CA, USA

**Correspondence**: Genevieve Jay Brett (brett33@hawaii.edu)

**Correspondence**: Genevieve Jay Brett (brett33@hawaii.edu)

Received: 28 Nov 2018 – Discussion started: 04 Dec 2018 – Revised: 13 Mar 2019 – Accepted: 16 Mar 2019 – Published: 05 Apr 2019

The importance of chaotic advection relative to turbulent diffusion is investigated in an idealized model of a 3-D swirling and overturning ocean eddy. Various measures of stirring and mixing are examined in order to determine when and where chaotic advection is relevant. Turbulent diffusion is alternatively represented by (1) an explicit, observation-based, scale-dependent diffusivity, (2) stochastic noise, added to a deterministic velocity field, or (3) explicit and implicit diffusion in a spectral numerical model of the Navier–Stokes equations. Lagrangian chaos in our model occurs only within distinct regions of the eddy, including a large chaotic “sea” that fills much of the volume near the perimeter and central axis of the eddy and much smaller “resonant” bands. The size and distribution of these regions depend on factors such as the degree of axial asymmetry of the eddy and the Ekman number. The relative importance of chaotic advection and turbulent diffusion within the chaotic regions is quantified using three measures: the Lagrangian Batchelor scale, the rate of dispersal of closely spaced fluid parcels, and the Nakamura effective diffusivity. The role of chaotic advection in the stirring of a passive tracer is generally found to be most important within the larger chaotic seas, at intermediate times, with small diffusivities, and for eddies with strong asymmetry. In contrast, in thin chaotic regions, turbulent diffusion at oceanographically relevant rates is at least as important as chaotic advection. Future work should address anisotropic and spatially varying representations of turbulent diffusion for more realistic models.

Chaotic advection (Aref, 1984; Shepherd et al., 2000) is a process by which rapid stirring, as manifested by the stretching and folding of material, is produced within a smooth and well-organized Eulerian velocity field. The enhancement of stirring can be attributed to chaotic fluid parcel trajectories and their rapid separation from nearby trajectories. There are many examples, ranging from simple models of purely laminar flow (e.g., Rom-Kedar et al., 1990; Samelson, 1992; Pierrehumbert, 1994; Malhotra et al., 1998; Poje and Haller, 1999; Coulliette and Wiggins, 2001, and other work reviewed in the texts of Ottino, 1990; Samelson and Wiggins, 2006) to modeled or observed, oceanographically or atmospherically relevant flows (e.g. Rogerson et al., 1999; Miller et al., 2002; Deese et al., 2002; Olascoaga and Haller, 2012; Sayol et al., 2013; Rypina et al., 2007, 2009, 2011a, 2012). In most cases the flow fields are two-dimensional and time-dependent, and when observed, they often occur at the sea surface or within the stratosphere (Polvani et al., 1995; Ngan and Shepherd, 1997). Three-dimensional examples also exist (e.g., Fountain et al., 2000; Rypina et al., 2015; Solomon and Mezić, 2003; Yuan et al., 2004; Branicki and Kirwan Jr., 2010, and Pratt et al., 2014, hereafter P2014) and often involve numerically modeled velocity fields, due to the limitations of observational methods.

A feature that is intriguing and quite common in these studies is that Lagrangian chaos is confined to certain sub-regions of the flow field, separated from each other by bands of material curves or surfaces that contain no chaotic Lagrangian motion. The chaotic regions are rapidly stirred as a result of the signature rapid separation of nearby trajectories, but the non-chaotic bands act as barriers that confine the stirring. In textbook examples, including area-preserving maps of time-periodic 2-D or steady 3-D velocity fields, the chaotic and non-chaotic regions form a fractal geometry, with bounded chaotic regions imbedded in larger chaotic seas, themselves bounded and imbedded in even larger chaotic regions (Chirikov, 1971, 1979; Casati and Ford, 1979; Gromeka, 1881; Dombre et al., 1986). In finite-time systems or systems with arbitrary time dependence, the distinction between chaotic and regular trajectories is difficult to define. A great deal of recent work in the field has resulted in the development of methods for identifying material barriers based on the notion of Lagrangian coherence. These methods include, for instance, finding sets of trajectories that experience the fastest separation rates from their close neighbors, identifying contours that undergo minimal stretching, locating sets of trajectories that remain compact in some sense and/or share a common property, or identifying trajectories that encounter the largest number of other trajectories (see Haller, 2002; Shadden et al., 2005; Froyland et al., 2007, 2012; Rypina and Pratt, 2017; Rypina et al., 2018, 2011b; Hadjighasem et al., 2017; Haller and Beron-Vera, 2012, 2013, as well as the review by Haller, 2015, and references contained therein). Applications of these methods often result in the identification of material contours and surfaces that act as barriers over finite time, thus allowing for partitioning between strongly and weakly stirred regions of the flow field.

Completely impenetrable material barriers only exist because of the deterministic nature of the trajectories. Even a low level of background turbulence at small scales, if represented as a diffusive process, would cause the barriers to become permeable or fuzzy over sufficiently long periods of time and perhaps nonexistent in any practical sense if the timescale of interest is long enough. The relevance of chaotic advection for the stirring of material within geophysical flows would appear to rest on several criteria. The first is that the flow field must contain persistent, long-lived (on the timescale of interest) features such as gyres, eddies, and jets, which by themselves generate regions of elevated stirring as well as separating barriers. Secondly, the stirring within these regions should be at least as important as that due to smaller scale, intermittent features (i.e., small-scale turbulence). Third, the barriers that exist in the absence of small-scale turbulence should retain meaning as suppressors of exchange between the rapidly stirred regions in the presence of the small-scale turbulence. For the flow considered in this paper the first aspect has been investigated and shown to be true (P2014; Rypina et al., 2015); this work concentrates on investigating the second and third aspects.

The terms “important” and “relevant” are somewhat subjective, and a particular aspect, such as the existence of barriers, that is of interest for one scientific question may not be so for another. We examine several measures of stirring and mixing in a particular case of a three-dimensional flow field: an idealized representation of an isolated eddy with horizontal swirl and vertical overturning. This idealized eddy is most likely to be similar to a submesoscale eddy within a surface mixed layer of the ocean, although the velocities of such eddies have not been well observed. The effects of stirring and mixing at these smaller scales, where vertical velocities become important, is increasingly under study (e.g. Mahadevan, 2016). Generally, increased resolution improves ocean model behavior (Griffies et al., 2015), so at lower resolutions, an ongoing challenge is parameterizing sub-grid-scale processes (e.g. Hallberg, 2013).

Our three-dimensional flow contains Ekman layers at the top and bottom of a cylindrical domain, and their thickness relative to the full depth is measured by an Ekman number. The Lagrangian structure of the steady as well as time-periodic, deterministic versions of this flow has previously been explored (P2014; Fountain et al., 2000; Rypina et al., 2015). This deterministic flow field can be approximated by an analytically described velocity field (Sect. 2), favorable for the efficient calculation of large numbers of trajectories. In this paper, we will add a stochastic disturbance representing small-scale turbulent diffusion to the deterministic flow. In addition, some of our calculations are done using velocity fields from a direct numerical integration of the Navier–Stokes equations (used in Sect. 5).

In order to examine the relevance and importance of stirring and mixing due to large-scale Lagrangian chaos compared to that due to small-scale turbulent diffusion, we use several distinct measures applied to our isolated eddy model. The first measure is a Lagrangian version of the Batchelor scale (Sect. 3), a measure of the smallest tracer filament width that can be produced by chaotic advection before small-scale turbulent diffusion arrests the progression to smaller scales. The second measure (Sect. 4) involves the dispersion of ensembles of initially closely spaced trajectories. The final measure (Sect. 5) is a bulk or “effective” diffusivity (Nakamura, 1996) that indicates the rate of irreversible mixing between volumes with different tracer concentrations. The analyses in Sects. 3–4 are based on a “kinematic” analytical model with and without stochastic perturbation; the analysis in Sect. 5 is based on a “dynamical” numerical solution of the Navier–Stokes equations.

We will consider the steady flow of a homogeneous and incompressible fluid in
a rotating cylinder of depth *H*, driven at the top by the stress due to a
differentially rotating lid. The resulting circulation has Ekman layers at
the top and bottom, and thus a central parameter is the Ekman number

where *ν* is the kinematic viscosity, Ω is the angular rate of
rotation of the cylinder, and *δ*_{E} is the thickness of the Ekman
layers. Much oceanographic literature has been devoted to the case in which
the differential lid rotation *δ*Ω is small $(\mathit{\delta}\mathrm{\Omega}/\mathrm{\Omega})\ll \mathrm{1})$, and the Ekman layers are relatively thin, *E*≪1. In
this case, a linear, asymptotic solution is available
(Greenspan, 1968, and Appendix A of P2014). According to this
solution (with *δ*Ω>0), fluid is drawn up into the top Ekman
layer from an inviscid and vertically rigid interior region that rotates at
half the angular velocity of the lid. The fluid is carried radially outward
and then downward within thin, viscous side-wall layers. When it reaches the
bottom, the fluid flows radially inward in a bottom Ekman layer and is expelled
upward into the interior region. Fluid trajectories thus spiral upwards in
the interior, outwards in the top Ekman layer, downwards near the side walls,
and inward in the bottom Ekman layer; Fig. 1 is a diagram of this
flow (see also Fig. 1 of P2014).

Although the setup described above and its linear asymptotic treatment have provided a foundation for a wide variety of models with geophysical and industrial applications (e.g., Lopez and Marques, 2010), it is not the most convenient for Lagrangian studies. One difficulty is that all fluid trajectories pass through small corner regions at the top and bottom of the cylinder. These regions are not resolved by the asymptotic solution and can be difficult to resolve numerically, particularly when the velocity field is to be used to accurately calculate trajectories that are cycling through the cylinder numerous times. For this reason it is advantageous to modify the forcing at the upper surface to conform to a stress that still acts in the azimuthal direction and is zero at the cylinder axis but approaches zero at the cylinder boundary as well. P2014 used one such forcing distribution to create a flow in which the downwelling occurs over a broad outer region of the inviscid interior, no longer confined to the thin, viscous side-wall layers. We will use the same velocities (obtained from a numerical model) for the tracer release experiments discussed later in this work.

Since numerical solutions are required to get a complete, dynamically consistent velocity field for the rotating cylinder, Lagrangian calculations requiring long integration times can become cumbersome, making it difficult to explore the variations in the governing parameters. As a compromise, past investigators have developed phenomenological models in which an incompressible Eulerian velocity field containing the qualitative features of the dynamically consistent fields is specified analytically and fluid trajectories are computed from it. Many of the calculations described below are based on such a model, hereafter referred to as the kinematic model. This new model is an improvement of the phenomenological model used by P2014 and Rypina et al. (2015) in terms of its more realistic portrayal of Ekman layers and inclusion of the Ekman number as a parameter.

The kinematic model specifies an analytically prescribed background velocity
field that is steady, incompressible, and has no azimuthal structure. Under
these conditions, all trajectories are regular, or non-chaotic. When
perturbed through the addition of an analytically prescribed
symmetry-breaking disturbance, one with azimuthal structure, Lagrangian
chaos arises in portions of the three-dimensional flow field. To see the
qualitative behavior of the flow, examine Fig. 1. The velocity
field is specified in nondimensional cylindrical coordinates $(r,\mathit{\theta},z)$, with $(\mathrm{1}\ge z\ge \mathrm{0})$ and (*r*≤*a*), where *a* is the
width-to-height ratio of the domain. The background flow has
$\partial /\partial \mathit{\theta}=\mathrm{0}$ and can be expressed as the sum of an azimuthal
velocity *V*(*r*,*z*) and an overturning circulation with radial and vertical
velocity components *U*(*r*,*z*) and *W*(*r*,*z*). The latter are specified by the
streamfunction Ψ:

where *F*(*z*) is the vertical portion of the streamfunction, and *R*(*r*) is the
radial portion of the streamfunction. The streamfunction relates to the
velocities by the negative *z* derivative of Ψ being the radial velocity
and the radial derivative being the vertical velocity. The vertical portion
of the streamfunction is

where *ζ* is a transformed vertical coordinate,

and the constants are defined by

In the limit of infinite cylinder radius, *a*→∞, the radial
portion of the streamfunction, $R\left(r\right)={r}^{\mathrm{2}}/s$, yields a dynamically consistent
solution for flow between two differentially rotating, horizontal plates.
Fluid flows radially inward within the bottom Ekman layer and is expelled
upward and eventually into the top Ekman layer, where it moves radially
outward. When *a* is finite the velocity needs to vanish at the cylinder
walls, and this can be accomplished by choosing *R* as

giving velocities

where *U* is radial and *W* is vertical.

The axisymmetric azimuthal velocity *V*, satisfying the incompressibility condition in 3-D, is defined as

This velocity leads to typical nondimensional trajectory rotation times of
20–200 for all Ekman numbers examined; the central orbit at
$(r,z)=(\mathrm{0.5},\mathrm{0.5})$ has a period of 16*π*≈50. At the maximum azimuthal
velocity, which occurs at $r=aH/\mathrm{3}$, the period is about 20. Model horizontal
velocities are typically between 0.01 and 0.1 in magnitude, which are
reasonable ocean velocities in meters per second. This choice of the velocity
scale being 1 m s^{−1} gives rotation times of several hours,
assuming the eddy radius is equal to its height (*a*=1). Using the same
scaling for vertical velocities, whose nondimensional values are *E*^{1∕2}
smaller, gives overturning times of 7 h to 2 months; although eddies with
this structure have not been carefully observed, vertical velocities near
submesoscale fronts reach 30 m day^{−1}, which is in line with
these rates. These and all other relationships between nondimensional model
values and their dimensional equivalents are listed in Table 1.
For all parameter values, there is upwelling in the center
(*r*=0) and weaker downwelling near the sides of the cylinder (strongest at
*r*=0.75*a*). There is horizontal convergence near the bottom and divergence
near the top; for *E* near 1, these are true for the full bottom and top
halves of the system.

As the Ekman number varies, the overturning streamfunction changes
qualitatively (Fig. 2). For $E>\mathrm{1}/\mathrm{60}$ the overturning circulation
is rounded and has a single internal fixed point corresponding to the
horizontal, circular trajectory described above as the central orbit (Fig. 2a, b).
For $E<\mathrm{1}/\mathrm{60}$ additional fixed points in the overturning
circulation arise at *r*=0.5 (Fig. 2c). These fixed points in
Fig. 2c are again circular periodic trajectories in 3-D, and the increasing number
arises through pitchfork bifurcations as *E* decreases (see Appendix A for
more details). The additional circular trajectories are associated with
smaller overturning cells imbedded in the larger cell (detailed example in
Appendix A, Fig. A2). The overturning streamfunction also
exhibits more vertical rigidity as *E* decreases, analogous to deeper oceanic
columns, in accordance with the Taylor–Proudman theorem
(Greenspan, 1968).

## 2.1 Symmetry-breaking perturbation

In the kinematic, axially symmetric, and analytically prescribed background flow described above all trajectories move along toroidal surfaces and are thus non-chaotic. In order to use this system to study the interplay of chaotic advection and turbulent diffusion, we must perturb the system to break the axial symmetry, which will introduce chaotic trajectories. The applied perturbation, approximating the flow produced by a lid rotating off-center, is a horizontal flow that decays in strength with depth and is described by the following streamfunction:

This general form allows for an *r*- and *z*-dependent adjustment to the
strength of the azimuthal velocity, with amplitude *ϵ*, and a symmetry-breaking component governed by the offset parameter *x*_{0}. If *x*_{0}=0, the
disturbance is axially symmetric; if it is nonzero, the disturbance has an
azimuthal variation in amplitude *ϵ**x*_{0}. The parameter *γ* can
be used to make adjustments in the radial structure of the disturbance. This
streamfunction is for velocities in the *x* and *y* directions, unlike *r*-
and *z*-dependent background overturning streamfunction; the velocities from
the two are added together. The perturbation velocities in *x* and *y* are

The corresponding azimuthal and radial velocity perturbations are

The perturbation streamfunction's
overall strength decays with depth and goes to 0 at the bottom (*z*=0). For
the rest of the work, we will use *a*=1 and *γ*=2 (Fig. 2d).
We note that the total, i.e., background plus perturbation,
azimuthal velocity can be zero at some locations in the domain for certain
choices of *ϵ*, but with *ϵ*<0.05 these locations are all very
close to the boundaries of the cylinder.

## 2.2 Comparison to dynamic model

In this section we compare our kinematic model to the Navier–Stokes (NS)
simulation of a rotating cylinder flow by P2014. We will use the kinematic
model for the analyses in Sect. 3.1 and 3.2 and the NS simulation for the
analysis in Sect. 3.3. We are interested in comparing the qualitative
features of the two model flows under steady symmetry-breaking perturbation.
It is important to note that the parameters of the two systems are slightly
different. The parameters that arise in the NS simulation are the Ekman
number, *E*, the aspect ratio, *α*, the displacement *X*_{0} of the lid's
center (labeled in Fig. 1; not to be confused with *x*_{0} in the
kinematic model), and the Rossby number, $\mathit{\text{Ro}}=\mathit{\delta}\mathrm{\Omega}/\mathrm{\Omega}$. The
kinematic model parameters are the Ekman number, *E*, the aspect ratio, *a*,
the perturbation offset parameter, *x*_{0}, and the strength of the
perturbation, *ϵ*. For matching the kinematic model to the NS
simulation, we set $\mathit{\alpha}=a=\mathrm{1}$ and examine four Ekman numbers used in P2014,
$E\in \mathit{\{}\mathrm{0.25},\mathrm{0.125},\mathrm{0.02},\mathrm{0.0005}\mathit{\}}$. The displacement and strength of the
kinematic perturbation are adjusted to match the behavior for a given Rossby
number and displacement of the lid in the dynamic simulation. The chosen
values are maintained throughout the rest of the work unless otherwise noted.
We do this rather than attempting a mathematical equivalence because the
kinematic perturbation has a different form than that describing a physical
lid rotating off-center. Our model mimics a flow with a small Rossby number,
so we compare our results to those from P2014's *Ro*=0.2, with lid
displacement ${X}_{\mathrm{0}}=-\mathrm{0.02}$.

Figures 3–4 show some examples of Poincaré
maps from the NS simulation (panels a, b, reproduced from P2014)
with maps from the kinematic model (panels c, d). It is important for our
purposes to achieve qualitative agreement in terms of the depth of the Ekman
layers, the vertical rigidity of the interior regions, and the overall layout
of regular, chaotic, and resonant regions. For the choice of the parameters
described above, there is a good match of these qualitative features. Each
case is marked by the presence of a substantial chaotic region that extends
from the radial center around the top and bottom boundaries and to our
largest radii near the perimeter of the cylinder. We henceforth refer to this
region as the “chaotic sea”. Also, in all cases there are many more
points near the surface than near the bottom; this is due to the higher
azimuthal velocities near the surface, and is seen in both the dynamic and
kinematic model. In *E*=0.25, both Poincaré sections show a series of nested
closed curves centered around $(r,z)=(\mathrm{0.5},\mathrm{0.5})$ corresponding to
quasiperiodic trajectories on nested tori. Between these are some thin
resonant layers with high numbers of small islands. For *E*=0.125, the main
feature is a series of larger islands between a set of nested tori and the
chaotic sea. For *E*=0.02, there is one large island with a number of
resonant layers surrounding it, including small islands. For *E*=0.0005, the
vertical structure of both models is more rigid, the kinematic model more so
than the NS simulation. Altogether, the kinematic model reproduces the
general features of the NS simulations, though there are often differences in
details, such as the number and widths of islands.

In this section, we examine the relative importance of chaotic advection and
turbulent diffusion for tracer distribution using a Lagrangian Batchelor
scale. The Batchelor scale, *δ*, is the length scale at which advection
and diffusion balance in their respective thinning and widening of a patch of
tracer. Chaotic advection thins tracer patches through averaged exponential
contraction in the contracting direction or directions, decreasing the relevant
length scale towards small scales where turbulent diffusion is dominant. In
this section, we represent turbulent diffusion as a scale-dependent
diffusivity. This diffusion widens tracer patches by moving tracer down its
gradient, spreading it out from its maximum. Below *δ*, diffusion
dominates tracer behavior, while above *δ* advection dominates. If
*δ* is larger than the structures in the flow induced by chaos, then
diffusion will overcome advection and wipe out these structures. The
structures of interest, induced by the deterministic, symmetry-breaking
perturbation (see Figs. 3–4), are the bands of
chaos, called resonant layers, surrounding regular island chains (see blue
diamond in Fig. 3d), and the chaotic sea region (outside the
red oval in Fig. 3d) located near the cylinder perimeter and
central axis, which are identified by visual inspection of Poincaré
sections. When we compare *δ* to these structures, we define their
widths as the difference between distances from the central orbit,
$(r,z)=(\mathrm{0.5},\mathrm{0.5})$, to the outermost and innermost part of the structure, measured
in Poincaré sections like Figs. 3–4.

In principle, the width of a tracer filament should approach the Batchelor scale regardless of initial conditions. If we consider an initial patch of tracer that is far from the Batchelor scale, advection and diffusion will not balance. If the patch is larger than the Batchelor scale, chaotic advection constricts the patch in the direction of fastest contraction so that it approaches the Batchelor scale. If the patch of tracer is smaller than the Batchelor scale, diffusion widens the patch to approach the Batchelor scale. When the width of a filament is at the Batchelor scale, the width will be steady in time but the concentration will continue falling.

Traditional formulations of the Batchelor scale use the Eulerian quantity – strain rate – to quantify advection and to find the scale at which advective and diffusive effects balance. Several rigorous derivations of a Lagrangian Batchelor scale have been presented (e.g. Thiffeault, 2004; Fereday and Haynes, 2004; Son, 1999), and a few papers have used less-rigorous scaling arguments to estimate the importance of chaotic advection (Rypina et al., 2010; Ledwell et al., 1993, 1998). Below we present a simple explanation for the Lagrangian Batchelor scale to gain intuition about this quantity, followed by a rigorous derivation of a Lagrangian Batchelor scale for a Gaussian tracer in a 3-D linear strain flow. The latter extends the work of Flierl and Woods (2015) from 2-D to 3-D.

The first formulation of the Lagrangian Batchelor scale uses dimensional
arguments to construct a quantity that has units of length from the
diffusivity *κ*, which quantifies the intensity of diffusion (and has
units of length^{2}time^{−1}), and the average exponential contraction rate
*λ*_{3}, which quantifies the thinning of a filament due to chaotic
advection (and has units of time^{−1}):

In a flow field with uniform steady strain, one could simply use the Eulerian
strain rate as the filament thinning rate. However, in flows with
a non-constant strain rate, the tracer will feel different strain as it is
advected by the flow, so a Lagrangian quantity such as the finite-time
Lyapunov exponent (FTLE) would be more appropriate. The FTLE quantifies the
average exponential separation rate between a trajectory and its close
neighbors over a finite-time interval Δ*t*,

Since separation rates between trajectories are generally different in
different directions, in 3-D flows there are three FTLEs that can be ordered
${\mathit{\lambda}}_{\mathrm{1}}\ge {\mathit{\lambda}}_{\mathrm{2}}\ge {\mathit{\lambda}}_{\mathrm{3}}$ and can be thought of as the
stretching and contraction rates of the three major axes of an infinitesimal
spherical blob of fluid as it deforms into an ellipsoid under the influence
of the flow field (see Fig. 5). For incompressible flows,
*λ*_{1}≥0, *λ*_{3}≤0 and ${\mathit{\lambda}}_{\mathrm{1}}+{\mathit{\lambda}}_{\mathrm{2}}+{\mathit{\lambda}}_{\mathrm{3}}=\mathrm{0}$. For the Batchelor scale in Eq. (16), the appropriate FTLE is that for
the most contracting direction, i.e., *λ*_{3}. FTLEs are most commonly
computed as

where *σ*_{i} values are the eigenvalues of the right Cauchy–Green deformation
tensor,

Here Δ*x*_{i} and Δ*x*_{i0} are the final and initial displacements
in the *i*th direction between initially nearby particles that are advected
by the flow over time interval Δ*t*. *G* can be calculated directly
from dense grids of simulated Lagrangian trajectories. We use the latter
method in our calculations to estimate *λ*_{3}.

As an alternative motivation of the Lagrangian Batchelor scale, we show
analytically that the width of a Gaussian tracer distribution asymptotically
approaches the Batchelor scale in a simple flow field. This derivation is an
extension to three dimensions of a two-dimensional calculation by
Flierl and Woods (2015). The main steps of the derivation are described
below, with more details in Appendix B. First, we assume that in the
Lagrangian frame, the velocity field is a steady linear strain with rates
*λ*_{i} in each direction such that the sum of the *λ* is zero,
giving an incompressible flow. Second, we assume that the tracer
concentration *C* initially has a Gaussian distribution in each direction,
and we look for a solution to the tracer evolution equation where it remains
Gaussian. In this case we can use the standard deviation of the Gaussian
distribution to measure the width of the filament in each direction. The
width in the most contracting direction, which is shrinking with rate
*λ*_{3}, is denoted by *σ*.
As shown in the Appendix, the differential equation for *σ* has a fixed point at

meaning that the width of the Gaussian patch in the fastest contracting
direction has a fixed point at the Batchelor scale, as expected from the
physical arguments about the balance between advection and diffusion. This
fixed point is attracting, meaning that for any initial width, the width in
the *λ*_{3} direction will converge to the Lagrangian Batchelor scale.
Mathematically there are also fixed points with negative *λ*_{3} and with
negative *σ* for positive *λ*_{3}, but neither corresponds to a real
positive tracer distribution. The full solution for *σ* is

More details and the full solution for *C* are in Appendix B.

## Results of Batchelor scale analysis

In order to calculate the Lagrangian Batchelor scale, *δ*, we use the
oceanic diffusivity estimates from Okubo (1971). In the ocean,
diffusivity is scale-dependent, increasing with size, as described by Okubo.
He used observations of horizontal dye diffusion at various scales ranging
between about 20 m and 100 km to find the empirical relationship

where *l* is the horizontal length scale of the dye patch (in cm) and *κ*
is in centimeters squared per second. Consistent with the lack of
density stratification in our model, we assume an isotropic
three-dimensional diffusivity. This assumption is supportable in the upper
ocean mixed layer and is consistent with our assumption of shallow eddies.

The variable nature of Okubo's *κ* makes determination of the
Batchelor scale a bit more subtle. In the case of a spatially variable
*κ*, the thinning of an initially large tracer patch will occur as
before, but as the filaments decrease in width, the corresponding *κ*
decreases as well. Following Rypina et al. (2010), we hypothesize that
equilibration will occur if during this process, the tracer scale *L*
approaches
$\left(\mathit{\kappa}\right(L)/|{\mathit{\lambda}}_{\mathrm{3}}|{)}^{\mathrm{1}/\mathrm{2}}=(\mathrm{0.0103}{L}^{\mathrm{1.15}}/\left|{\mathit{\lambda}}_{\mathrm{3}}\right|{)}^{\mathrm{1}/\mathrm{2}}$. Solving
for *L* yields the Batchelor scale:

where *λ*_{3} (in s^{−1}) yields *δ* (in cm s^{−1}).

To relate our dimensionless kinematic model FTLEs to Okubo's diffusivities,
we need to set dimensional time and diffusivity scaling factors. We
previously discussed the winding times and associated velocity scaling of
1 m s^{−1}; our desired scaling factors can be computed with
this velocity scaling and a length scale. The main parameter of the background
model is the Ekman number, the square of the ratio of Ekman layer thickness
to eddy depth. Due to the unstratified nature of our flow, we focus on two
intermediate Ekman numbers: *E*=0.125 and *E*=0.02. Assuming an Ekman depth
of about 40 m, which is within the range of open-ocean observations
(see Lenn and Chereskin, 2009, and references therein), our shallower eddy is
about 110 m deep, whereas *E*=0.02 would correspond to an eddy depth of
about 280 m. Depending on region and season, it is possible for either of
these to be within the surface mixed layer of the ocean, which can reach
500 m in subpolar regions in the winter but may decrease to a few meters in
the summer. Since the aspect ratio of the width to depth of our eddy is 1,
the corresponding eddy radius is also between roughly 100 and 300 m. Using the
product of the dimensional depth of the eddy and the chosen velocity scale,
Okubo's diffusivities can be made nondimensional. In contrast, the FTLEs could
be made dimensional using the time step in seconds. These scalings are
explicitly given in Table 1; we will discuss the results in
nondimensional terms.

The calculated *δ* values are shown in Fig. 6 next to the
widths of chaotic regions; both widths are made dimensional using the eddy
depths. The range of *δ* values is due to the spatial variation in the
most contracting FTLE, *λ*_{3}, in the region (see Figs. 3–4 for most stretching FTLEs, which are of
the same magnitude). FTLEs were estimated over an integration time of 400;
the range of FTLE magnitudes does not noticeably change when the integration
time is decreased by half. The widths of the chaotic sea and smaller resonant
regions were estimated from inspection of Poincaré sections. The Batchelor
scale is generally about 0.01–0.08, which is similar to the resonant layer
widths and smaller than the chaotic sea widths. The dimensional diffusivities
at these scales range from $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m^{2} s^{−1} at 1 m to 0.06 m^{2} s^{−1} at 140 m, which are considerably
smaller than diffusivities on the horizontal scale of eddies themselves,
about 0.5–8.2 m^{2} s^{−1} for 1–10 km. The
Batchelor scale results imply that chaotic advection is expected to influence
tracer distribution throughout the system but dominates only in the wider
chaotic sea region.

In this section, we quantify the relative effects of turbulent diffusion and chaotic advection using the dispersion (or spread) of sets of initially nearby trajectories in the kinematic model. We consider chaotic advection to be dominant compared to diffusion when the ensemble spread is greater for the deterministic perturbation that induces chaos than for the stochastic perturbation that simulates turbulent diffusion. Ensembles of 100 to 300 trajectories that begin inside a small sphere have been examined for their behavior under various perturbations. Other initial conditions, on a torus or axial circle, give similar results (not shown). The spread of trajectories is measured in terms of Ψ values; the streamfunction of the background flow is given by Eq. (2). Examining the spread in Ψ is convenient because it leads to zero spread for particles following the background flow. However, it is important to note that this interpretation limits the directions of chaotic stretching that are considered – it is possible for the fastest-spreading direction to be along the background streamlines, which would not be visible in the coordinates chosen.

To simulate turbulent diffusion, we add a stochastic velocity perturbation to
the background model flow. The stochastic perturbation is a random flight
model created by adding small pseudorandom values with a Gaussian
distribution to the velocity at fixed intervals of time Δ*t*. The
equation governing a fluid particle trajectory is then

where *i* is a direction index, *U*_{bi} is the background velocity, and
${u}_{i}^{\prime}$ are the stochastic additions. These velocity additions are
uncorrelated and lead to a Gaussian random walk behavior
(Zambianchi and Griffa, 1994). Using the described stochastic perturbation,
although it is quite simple, with *U*_{bi}=0 or a constant, the variance of a
set of trajectories grows linearly in time, while the standard deviation
grows linearly with the square root of time, as expected for diffusion. The
diffusivity, *κ*, is computed from the 1-D relationship for a Gaussian
random walk: $\mathit{\kappa}={s}^{\mathrm{2}}/\mathrm{2}\mathrm{\Delta}t$, where *s* is the standard deviation of
step size in the random walk. To choose the level of diffusivity for the
stochastic perturbation, we consider the turbulent diffusivities near the
Batchelor scale as computed in the previous section. The Okubo diffusivities
at the Batchelor scale are in the range $\mathit{\kappa}\in [{\mathrm{10}}^{-\mathrm{4}},{\mathrm{10}}^{-\mathrm{2}}]$ m^{2} s^{−1} across the four Ekman
numbers examined, which is nondimensionally $\mathit{\kappa}\in \mathit{\{}{\mathrm{10}}^{-\mathrm{6}},\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}\mathit{\}}$. As our primary example, we will discuss the level of
diffusivity $\mathit{\kappa}={\mathrm{10}}^{-\mathrm{6}}$. This diffusivity requires a certain step size
*s* for the stochastic perturbation, which relates to the distribution of
*u*^{′} by $s=\mathit{\sigma}\mathrm{\Delta}t/\mathrm{3},$ with *σ* being the standard deviation of *u*^{′}
and Δ*t* being the numerical time step (0.01) and having the factor of 3 due to the
details of a fourth-order Runga–Kutta integration. The next position, using
this method, is estimated using the weighted sum of estimates of the velocity
at the current position (*v*_{1}, weight 1∕6), the halfway point estimated from
the current position (*v*_{2}, weight 1∕3), the halfway point estimated using
*v*_{2} (*v*_{3}, weight 1∕3), and the final point estimated using *v*_{3} (*v*_{4},
weight 1∕6). Only *v*_{1} and *v*_{4} include stochastic additions, leading to
the 1∕3 factor. Together, these give

so *σ*=0.042. We will also discuss a smaller stochastic
perturbation, $\mathit{\kappa}={\mathrm{10}}^{-\mathrm{7}}$ and *σ*=0.013, and a larger one,
$\mathit{\kappa}={\mathrm{10}}^{-\mathrm{5}}$ and *σ*=0.13. The stochastic perturbation with
$\mathit{\kappa}={\mathrm{10}}^{-\mathrm{6}}$ has kinetic energy (integrated over the cylinder) about the
same as the background flow: $\int ({\mathit{u}}^{\prime}{)}^{\mathrm{2}}\approx \int ({\mathit{U}}_{\mathrm{b}}^{\mathrm{2}})\approx \mathrm{0.63}$. The perturbation with $\mathit{\kappa}={\mathrm{10}}^{-\mathrm{7}}$ has kinetic energy
about the same as the deterministic perturbation with *ϵ*=0.01 and ${x}_{\mathrm{0}}=-\mathrm{0.5}$, such that $\int ({\mathit{u}}_{\mathrm{d}}^{\prime}{)}^{\mathrm{2}}\approx \int ({\mathit{u}}_{\mathrm{s}}^{\prime}{)}^{\mathrm{2}}\approx \mathrm{0.075}$, where ${\mathit{u}}_{\mathrm{d}}^{\prime}$ is the deterministic perturbation velocity
and ${\mathit{u}}_{\mathrm{s}}^{\prime}$ is stochastic.

We begin with an example for *E*=0.125, showing the spread of trajectories
(measured in terms of the background streamfunction Ψ) in the presence
of either the deterministic or the stochastic perturbation. Trajectories are
started on a small sphere located entirely in the chaotic sea region centered
on $(r,z)=(\mathrm{0.1},\mathrm{0.5})$ (see Fig. 3 for the Poincaré section).
For the deterministic perturbation at early times, trajectories oscillate
through the background streamfunction because the perturbation velocities
form an azimuthal wave (Fig. 7a). The frequency of this
oscillation depends on the exact location of the trajectory, so with time,
trajectories move out of phase due to the cumulative effect of their slightly
different oscillatory frequencies. It takes a few cycles of overturning to
develop noticeable spreading, but then the spread grows quickly.

For the stochastic perturbation (Fig. 7b), trajectories are uncorrelated as they spread across the background streamfunction. There are no oscillations in time because the perturbation acts separately on each trajectory at each time step, leading to continuous and monotonic spreading of the ensemble. This spreading is similar to diffusion, but the increase in the range of trajectories does not depend on the gradients of concentration – Fick's law does not apply. If both perturbations are included (Fig. 7c), trajectory ensembles maintain some of their oscillatory behavior but spread out in a more continuous fashion due to the stochastic perturbation. In this example, and over timescales considered, we conclude that the stochastic perturbation dominates at early times but that chaotic spreading takes over at times larger than about 1000. Over an even longer time period, turbulent diffusive spreading is expected to overtake chaotic spreading.

We next compare the spreading of trajectory ensembles in Ψ with a variety of perturbations for the same initial conditions as in Fig. 7 using the range over time (Fig. 8); results are similar when the variance in Ψ is used for comparison (not shown). Chaotic advection dominates when the spread in Ψ for an ensemble under deterministic perturbation is larger than the spread under stochastic perturbation. The spread from the deterministic perturbation is very fast, appearing to be qualitatively exponential, for a period of time, as expected for a region with high FTLEs, which indicates exponential growth on average but is limited to the width of the chaotic region in which the ensemble begins (e.g., red curve in Fig. 8a). In contrast, the stochastic perturbation will spread with the square root of time until it reaches the cylinder boundaries (e.g., dark blue curve in Fig. 8a). Therefore, the time when the deterministic perturbation has greater spread will be limited to between when fast chaotic-advection-induced separation starts in the deterministic perturbation, which requires sufficient interaction with hyperbolic regions, and when the stochastic perturbation spreads the ensemble to the width of the chaotic region.

In the chaotic sea region (Fig. 8a, c), ensembles with
stochastic perturbations all have their ranges in Ψ grow in a manner
similar to the square root of time, and the spreading is faster for larger
*κ*. The ensembles with deterministic, chaos-inducing perturbations
experience an initial delay before they begin quickly growing. Once rapid
growth sets in, they spread to the width of the chaotic region between the times
500 and 3000. Larger deterministic perturbations lead to earlier and faster
spreading as well as wider chaotic regions. For the weaker deterministic
perturbation *ϵ*=0.01, there are some time intervals over which
chaotic spreading in the chaotic sea dominates stochastic spreading. These
instances occur more readily in the case of the shallower eddy (*E*=0.125;
Fig. 8a) and less so for the deeper eddy (*E*=0.02; Fig. 8c).
However, larger deterministic perturbations (e.g.,
*ϵ*=0.08) produce chaos that is dominant over longer times, an extreme
example being the pink curve in Fig. 8a.

We can also consider the timescales over which diffusive and advective
processes with similar kinetic energy (red and light blue curves in Fig. 8)
dominate over each other. The ensembles released in the
chaotic sea show that over the first few hundred time steps, turbulent
diffusion dominates the spread (Fig. 8a and c at *t*<1000), as
chaotic advection does not yet show significant growth. After that we see a
period of fast growth due to chaotic advection, which quickly overtakes the
slower diffusive spreading. This rapid growth stops when the advective spread
reaches the width of the chaotic region, and the diffusive spreading, which
is not limited by the chaotic region width, is then able to catch up and
exceed chaotic advection. Of course, these processes will be acting at the
same time, not separately; the green curves in Fig. 8 are
examples when small perturbations of both types are present. In this case,
spreading of the ensemble begins immediately, as in simulations with only
stochastic perturbations, but then has a time period of pronounced growth and
some oscillations, as seen in simulations with only steady perturbations.

We also examined the behavior of trajectories beginning at $(r,z)=(\mathrm{0.4},\mathrm{0.5})$,
a small distance from the central fixed orbit, within the region containing
resonant layers (Figs. 3–4). In these cases,
the same behavior as in the chaotic sea region occurs for the spreading under
stochastic perturbations (Fig. 8b, d). The spreading under
deterministic perturbations is much slower than in the outer chaotic sea
region for *ϵ*=0.01 (red curves in Fig. 8b, d), and
diffusion dominates at all times for all values of *κ* shown. With
*ϵ*=0.08, the chaotic region is larger, and growth due to the
deterministic perturbation is generally more rapid than that due to
diffusion, at least within the time window when chaotic advection begins and
until saturation occurs (pink curves in Fig. 8b, d).

From the spreading of ensembles of trajectories, we see that the wider chaotic regions are where chaotic advection dominates over turbulent diffusion (at least over some time intervals), as expected from our scaling arguments. However, those scalings did not include considerations of time including considerations of when fast chaotic-advection-induced stretching begins, as FTLEs are time averages; the delay in chaotic stretching decreases the period of time when chaotic advection is important. This time period begins when fast advective stretching is first apparent and ends when turbulent diffusion has spread across the region under consideration. From these ensembles, we would expect a set of passive 3-D drifters or an injected tracer beginning in a blob to spread out diffusively, be stretched and folded throughout the chaotic sea, producing strong filamentation, then gradually diffuse across the barriers of the chaotic sea and into the remainder of the eddy. During the later stage, tracer variance due to the formation of filaments by chaotic advection would be gradually eroded by turbulent diffusion. This sequence of events will be apparent in tracer simulations shown in the next section.

In this section we analyze the effects of the symmetry-breaking,
chaos-inducing deterministic velocity perturbation on the stirring and
mixing of a diffusive tracer in a dynamically consistent numerical model of a
rotating cylinder flow. Dye experiments are often used in both the ocean and
the laboratory to understand the stirring and mixing in a fluid (examples
include Fountain et al., 2000; Ledwell et al., 1993, 1998).
The distributions of passive tracers like dye are created by the advective
and diffusive patterns without the feedback onto the flow that would occur
with temperature or salinity, allowing for insight into those processes. For
our simulations we turn away from the kinematic model and take advantage of
the existing numerical model that solves Navier–Stokes equations
corresponding to the rotating cylinder flow accompanied by integration of the
advection–diffusion equation with diffusivity *k* for a passive tracer, both
described in P2014. As discussed earlier, these simulations have the
advantage of being dynamically consistent at the cost of being
computationally expensive, whereas economy of the kinematic model allows us
to explore a wider range of parameters.

Our main quantification tool is Nakamura's effective diffusivity: a background diffusivity scaled by a representation of the stretching of dye concentration contours by advection. Two-dimensional and quasi-three-dimensional analyses of effective diffusivity have been applied to the atmosphere and ocean (Nakamura, 1996; Nakamura and Ma, 1997; Haynes and Shuckburgh, 2000; Abernathey et al., 2010). For our fully three-dimensional system with constant density, the effective diffusivity can be written as

where *C* is tracer concentration, *V* is volume, and $\widehat{f}$ indicates an
average of function *f* over the area of a concentration surface. The imposed
small-scale diffusivity *k* is constant and is thus more closely related to
the *κ* used in Sect. 4 for the stochastic perturbation than the
scale-dependent Okubo *κ* in Sect. 3. (It is not clear how one would
incorporate a scale-dependent diffusivity into Nakamura's formulation.)
The volume *V* is a one-to-one mapping of tracer concentration and volume
such that *V*(*C*) is the volume occupied by fluid with concentrations greater
than *C*. The derivation leading to the above definition for *κ*_{eff}
can be found in Shuckburgh and Haynes (2003), who perform the algebra in
2-D but note that the 3-D development is identical.
Equation (16)
describes an effective diffusivity that is amplified from the small-scale
diffusivity by a factor of the degree of contortion of the concentration
contour. The units of the effective diffusivity are those of *k* (typically
m^{2} s^{−1}), multiplied by
m^{4}, or volume squared divided by length squared, which is
the same as surface area squared. Larger effective diffusivity leads to
larger diffusive fluxes of tracer. This amplification can be understood as
being caused by advective stretching and folding of tracer contours which
increases the area of surfaces of constant *C*, thereby amplifying gradients
of *C* and speeding up diffusive fluxes. This amplification factor is
precisely the surface area squared in the rare situation where $\left|\mathrm{\nabla}C\right|$
is constant on a *C* surface.

Both advection and diffusion redistribute tracer concentration and influence effective diffusivity. The effective diffusivity allows the effects of advection to be included in a diffusive term:

As advection stretches and folds the initial tracer, creating filaments, the
surface area of a contour and gradients of the tracer increase, leading to
larger *κ*_{eff}. Then, as diffusion smooths the tracer field, wiping
away the filaments, gradients decrease and contours become smoother, with a
lower surface-area-to-volume ratio. We compare the effective diffusivity with
a deterministic perturbation to that without; any increase is due to
increased stirring, which gives a quantitative measure of how important that
stirring is for the distribution of tracer in each region of the flow.

As a secondary quantification tool, we use the volume-integrated tracer variance function, *χ*^{2} (Pattanayak, 2001):

where *V* here is simple volume. Stirring increases the variance of a tracer,
while mixing decreases it. When *χ*^{2} is increasing, stirring is dominant
and the slope of *χ*^{2}(*t*) quantifies the stirring rate. The tracer
variance function was used to relate the Ekman number, perturbation strength, and
stirring rate for the rotating cylinder in P2014; the authors found that
stirring increased with larger perturbations and was nonmonotonic with *E*,
peaking near *E*=0.01.

The numerical simulations are run using the solver NEK5000 for several
diffusivities and strengths of the symmetry-breaking deterministic
perturbation. This model solves the incompressible Navier–Stokes equations
using a spectral element method (see https://nek5000.mcs.anl.gov, last access: 3 September 2018, P2014,
Fischer, 1997). The domain has an identical radius and height,
matching the aspect ratio assumed in our kinematic model. The
symmetry-breaking perturbation is created by moving the central axis of the
imposed surface lid stress a fraction of the radius *X*_{0} from the cylinder
axis so that *X*_{0} becomes the primary parameter determining the
perturbation strength. The ${X}_{\mathrm{0}}=-\mathrm{0.02}$ case is what was used to compare
Poincaré sections with the kinematic model, so qualitative features match
the *ϵ*=0.01 cases. The ${X}_{\mathrm{0}}=-\mathrm{0.16}$ case is a significantly larger
perturbation, similar to the *ϵ*=0.08 case in the previous section.
The nondimensional imposed tracer diffusivity, *k*, is 10^{−4} or
10^{−6}. Using Okubo's scaling, the lower diffusivity is appropriate for
scales near 1 m, while the larger is appropriate for scales near 50 m. After
the simulated velocity field is spun up, the tracer concentration, *C*, is
initialized with a constant vertical gradient, $C=\mathrm{1}-z$.

The set of simulations performed allows for an examination of the effects of
changing *E*, *k*, and *X*_{0}. They are *E*=0.125 and $k={\mathrm{10}}^{-\mathrm{4}}$,
${X}_{\mathrm{0}}\in \mathit{\{}\mathrm{0},-\mathrm{0.02},-\mathrm{0.16}\mathit{\}}$ and *E*=0.02, and $k\in \mathit{\{}{\mathrm{10}}^{-\mathrm{4}},{\mathrm{10}}^{-\mathrm{6}}\mathit{\}}$ and
${X}_{\mathrm{0}}\in \mathit{\{}\mathrm{0},-\mathrm{0.02}\mathit{\}}$ for a total of seven simulations. Each simulation is run
for a time of 300 after the tracer is initialized. The evolution in time of
the tracer variance function and Nakamura effective diffusivity integrated
over the volume of the cylinder are described first; we then discuss the
evolution of the dye and finally the spatial characteristics of the Nakamura
effective diffusivity.

The tracer variance function over time (Fig. 9a–c) initially grows nearly linearly as stirring creates filaments and large gradients. The function then has a single maximum that occurs at the time when diffusive mixing starts to overcome stirring so that the variance of the tracer begins to decrease. The maximum occurs earlier when either the imposed diffusivity or the strength of the deterministic perturbation increases. Increasing the diffusivity makes the maximum occur earlier by increasing the strength of the mixing (Fig. 9a, b). Increasing the deterministic perturbation also makes the maximum occur earlier, as faster stirring creates larger gradients, in turn increasing diffusive fluxes (Fig. 9c, red curve).

The maximum of the tracer variance function increases with decreased
diffusivity, as more filamentation can occur before diffusion wipes the
filaments out. This change in maximum is most evident in the difference
between $k={\mathrm{10}}^{-\mathrm{4}}$ and $k={\mathrm{10}}^{-\mathrm{6}}$ for *E*=0.02, where the decrease in
diffusivity increases the maximum of the tracer variance function by an order
of magnitude (Fig. 9a, b). Changes in the maximum as the size
of *X*_{0} is increased from 0 to 0.02 are small and negative, because the
slightly earlier time of the maximum combined with similar stirring rates leads
to a slightly smaller maximum with the perturbation. In the case of
*E*=0.125 and ${X}_{\mathrm{0}}=-\mathrm{0.16}$, the maximum is larger than with either *X*_{0}=0 or
${X}_{\mathrm{0}}=-\mathrm{0.02}$ due to faster stirring and a different spatial pattern of the
dye, which will be discussed later.

The effective diffusivity, *κ*_{eff}, integrated over the total volume
shows an overall progression similar to the tracer variance function, which
indicates the dominance of the gradient term over both the $\partial C/\partial V$ term in *κ*_{eff} and the $|C{|}^{\mathrm{2}}$ term in *χ*^{2} (Fig. 9d–f).
The initial slope and details of the maximum can be
understood as relating to perturbation and diffusivity strengths in the same
manner as for *χ*^{2}. At longer times, the integrated effective diffusivity
reaches a nearly constant positive value unlike *χ*^{2}, which approaches
zero. This constant value can be estimated by using the surface area
representation of *κ*_{eff}. At long times, here meaning after many
overturns but before diffusion removes all gradients, the shape of tracer
surfaces is distorted nested tori (look ahead to Fig. 10h). If
the *C* surfaces were nested circular tori, $\left|\mathrm{\nabla}C\right|$ would be constant
along the surfaces, and then ${\mathit{\kappa}}_{\mathrm{eff}}=k{A}^{\mathrm{2}},$ where *A* is the surface
area of a given toroidal tracer contour. The volume integral of the squared
surface area of circular tori nested around $(r,z)=(\mathrm{0.5},\mathrm{0.5})$ multiplied by
the background diffusivity is *k**π*^{6}∕8, which we expect to be the minimum
for ∫*κ*_{eff}d*V* in this system while gradients are nonzero (see
Appendix C for details). This value is shown as black dashed lines in Fig. 9d–f
and is just below the lowest ∫*κ*_{eff}d*V* value
seen. The higher values for *κ*_{eff} with steady perturbations at long
times correspond to persistent asymmetries in the tracer field which result
in larger constant concentration surface areas. The extreme case is
*E*=0.125, ${X}_{\mathrm{0}}=-\mathrm{0.16}$, and $k={\mathrm{10}}^{-\mathrm{4}}$, which has the most asymmetric dye
contours (Fig. 11i); here, the long-time value of ∫*κ*_{eff}d*V* is about twice as large as for circular tori.

Further insight can be gained by perusal of vertical sections of *C* and
*κ*_{eff} (Figs. 10 and 11). A caveat is that
*κ*_{eff} is a nonlocal property, so plots show the values
*κ*_{eff}(*C*) mapped onto the locations on the sections with
corresponding dye concentrations, *C*, while they are calculated using the
distribution of *C* over the whole volume at that time. These mappings are
noisier than the sections of *C* because the numerically computed
*κ*_{eff}(*C*) is nonmonotonic and can have large changes with small
changes in *C*. Nevertheless, these plots can yield some insights into the
time histories shown in Fig. 9. Figure 10 is restricted to
cases with *E*=0.02, while Fig. 11 is restricted to *E*=0.125. The
two are laid out differently, with the former designed to emphasize the
effects of varying *k* and the latter designed to explore variations in the
strength *X*_{0} of the perturbation. Both figures contain snapshots from an
early time (*t*=39) in the simulation, before diffusion has had a chance to
arrest growth in the tracer variance function, and at a late time (*t*=299)
when *κ*_{eff} has reached a quasi-steady value. In all cases, the *C*
sections become smoother and their range decreases between the snapshots, due
to continued mixing. The high *κ*_{eff} values are enhanced over much of
the sections' area at the early time and localized to mostly the chaotic sea
region at the late time.

The early development (*t*=39) of the tracer field, *C*, and of
*κ*_{eff} can be seen in Fig. 10a–f. With no
disturbance present (*X*_{0}=0) and $k={\mathrm{10}}^{-\mathrm{4}}$ (Fig. 10b), the
initially horizontal lines of constant *C* have been advected by the axially
symmetric overturning circulation such that contours of constant *C* are
roughly aligned with the overturning streamfunction. For an initial broad
gradient in any direction, we expect the same realignment after the first few
overturnings as the tracer is passively advected by the background velocity
field. We believe, then, that the tracer distribution that exists at later
times is somewhat independent of initial distribution. The corresponding
*κ*_{eff} at *t*=39 (Fig. 10e) exhibits high values at the
edges of filaments created by the straining motion of the symmetric
background flow, despite the fact that no trajectories are chaotic. When a
disturbance is added (${X}_{\mathrm{0}}=-\mathrm{0.02}$, Figs. 10c, f) the axial
symmetry is broken and the peak values of *κ*_{eff} are reduced. The
latter is somewhat surprising, since we have already seen (Fig. 9b)
that the volume-integrated values of *κ*_{eff} are nearly the same for
the disturbed and undisturbed case. The situation is made clearer if one
notes that moderate values of *κ*_{eff} (light blue in Fig. 10f)
are more widely distributed in the disturbed case. A similar result can be
seen by comparing the case *X*_{0}=0 (Fig. 11a, d) to ${X}_{\mathrm{0}}=-\mathrm{0.02}$
(Figs. 11b, e), all for *E*=0.125. Again, the unperturbed
(symmetric case) has larger peak values while the perturbed case has more
locations with moderate values of *κ*_{eff}, resulting in a similar
volume-integrated value of *κ*_{eff} (Fig. 9f). It is possible
that slight increases in stirring in the perturbed cases have caused more
mixing than in the unperturbed cases, even over the short interval before
these snapshots, leading to a lower range of *C* and smaller average
gradients in the perturbed cases. However, the volume-integrated measures
(Fig. 9) do not show any clear indications of that process occurring.

When the imposed diffusivity *k* is decreased by 2 orders of magnitude,
with *X*_{0} fixed at −0.02, the results are remarkably different. To begin
with, a comparison of Fig. 9d with Fig. 9e shows that
*κ*_{eff} is generally larger at any particular time when *k* takes the
smaller value. As Fig. 10a and c show, the tracer field
contains much finer filaments when $k={\mathrm{10}}^{-\mathrm{6}}$, consistent with the reduction
of the Batchelor scale. The distribution of *κ*_{eff} is broader and
with larger peak values for this lower numerical diffusivity (compare Fig. 10d and f).
The higher *κ*_{eff} indicates that despite
the decrease in *k*, the effects of stirring on the contours, as measured by
*κ*_{eff}∕*k*, have more than compensated, resulting in a higher rate of
irreversible property exchange. Thus the combined effect of smaller
diffusivity and finer filaments (i.e., stronger tracer gradients) leads to
more rapid mixing across tracer contours.

The results that have just been described occur early (*t*=39) in the
evolution of the tracer field, at a time when fluid parcels have overturned
just a few times and the perturbation amplitude *X*_{0} has been small. For
this weakly perturbed flow, Lagrangian chaos requires many overturns to be
significant, so we now turn our attention to the results for *t*=299 (Figs. 10g–l and 11g–l). Here a comparison between the
unperturbed and perturbed cases (contrast Fig. 10h and k with
Fig. 10i and l and also Fig. 11g and j with Fig. 11h and k)
reveal only modest differences in the spatial distribution and magnitude of
*C* and *κ*_{eff}. As in the early snapshots, there is a tendency for
the unperturbed flows to have higher peak values of *κ*_{eff}, while the
perturbed flows produce moderate values over a larger area. Decreasing the
value of *k* again has the effect of creating a more fine structure (Fig. 10g)
and of increasing the peak values of *κ*_{eff} by an
order of magnitude (Fig. 10j).

So far, the consequences of the symmetry-breaking disturbance are modest.
However, dramatic differences occur when *X*_{0} is increased from −0.02 to
−0.16 for *E*=0.125. The tracer distribution is markedly distorted at early
times (compare Fig. 11b with Fig. 11c), and strong tracer
gradients remain present even at *t*=299, at a time when the gradients in the
unperturbed and weakly perturbed cases have been strongly eroded (compare
Fig. 11g and h with Fig. 11i). The peak values of
*κ*_{eff} at *t*=299 (Fig. 11l) remain comparable to those
of the weakly perturbed case (Fig. 11k) but occupy a much larger
volume, making the volume-integrated *κ*_{eff} much larger, in agreement
with Fig. 9f.

For a different perspective, we examine the mean *κ*_{eff} in subdomains
of the system corresponding to a regular island and a region of the chaotic
resonant layer of roughly the same size. The cross sections of the cylinder
along the *x* and *y* axes are broken into different regions using the
matching Poincaré sections of the perturbed flow (Fig. 12).
Demarcation of these subdomains was most straightforward for the case
*E*=0.02, due to its large island and extended resonant region. While we used
Poincaré sections as guidance for defining regular and chaotic regions,
other methods (for example, Haller et al., 2018) could be used instead
for the more precise delineation of the phase space. The mean *κ*_{eff}
in the chosen subdomains gives a clear result in the *E*=0.02 and $k={\mathrm{10}}^{-\mathrm{4}}$
case (Fig. 12c), where at long times, when the overall gradients
have smoothed out, the resonant regions have about twice the effective
diffusivity as the islands. The islands' *κ*_{eff} at that time
approximately matches the value from the same region in the unperturbed
simulation, indicating that chaos has not affected this area. In the
*E*=0.02 and $k={\mathrm{10}}^{-\mathrm{6}}$ case (Fig. 12d), the mean *κ*_{eff} is
typically higher in the resonant region than in the island, but the
differences are less pronounced. It is notable that at *t*>130,
*κ*_{eff} is larger in the island than in the same unperturbed region,
perhaps because islands are not completely regular and contain smaller
chaotic resonant regions within them.

Overall, these dye experiments show that chaotic advection enhances Nakamura
effective diffusivity within the chaotic sea at some times in all cases
examined. The amount of enhancement is controlled by both the size of the
perturbation and the imposed diffusivity. A larger perturbation leads to
greater enhancement (higher *κ*_{eff}). A smaller diffusivity leads to
more filamentation (higher *χ*^{2}) and highly elevated enhancement (much
larger *κ*_{eff}).

The main goal of this work is to establish whether the stirring due to chaotic advection in an idealized model of an upper ocean eddy remains relevant in the presence of levels of background turbulent diffusion that are consistent with observations. The answer is that chaotic advection can indeed be relevant, and in some cases dominant, within certain regions of the flow field and over certain time intervals. The region most likely to feel the effects of chaotic advection is the extensive chaotic sea that exists in all simulations, and this is especially pronounced when the eddy is shallow. Chaotic stirring in the smaller and more isolated resonant regions is less likely to be important. This conclusion comes with many caveats related to idealizations (e.g., homogeneous turbulent diffusion) and uncertain parameter values (e.g., background diffusivity and strength of perturbation).

A second focus of the work has been to explore different bases for comparison
of the effects of chaotic advection and homogeneous turbulent diffusion. To
this end we have identified three metrics for comparison and are now in a
position to discuss their advantages and disadvantages. The first metric is
the Lagrangian Batchelor scale (Sect. 3), an estimate of the equilibrium
width of a passive tracer filament. Equilibrium is achieved when transverse
compression due to advection, as quantified by the negative Lyapunov exponent
with the largest magnitude (*λ*_{3}), is balanced by the diffusive
spreading of the tracer. Below the Batchelor scale, diffusion is stronger
than advection. When this width is smaller than that of the chaotic
regions, advection dominates; when it is larger, diffusion dominates. We fixed the turbulent diffusivity
using Okubo's empirical formula and calculated the Batchelor scale *δ*
using the rate of chaotic filament stretching, *λ*_{3}, computed
numerically as the largest negative finite-time Lyapunov exponent for the
kinematic model. The resulting Batchelor scale varies from O(1 m) for *E*=0.25
to O(100 m) for *E*=0.0005. These values of *δ* are smaller than the
spatial extent of the chaotic sea over all *E* values considered (0.25,
0.125, 0.02, and 0.0005) but of similar magnitude to the widths of the
resonant regions.

Interpretation of the Lagrangian Batchelor scale analysis would appear to be
straightforward, but it does not comprehend the fact that chaotic advection
may only be dominant over a finite-time interval, which is averaged in the
FTLEs. Even when the level of background turbulent diffusion is weak, it will
eventually spread beyond the region of Lagrangian chaos. There is also a
level of uncertainty due to the choice of integration time over which
*λ*_{3} is calculated. Finally, it is not yet possible to calculate
*λ*_{3} from ocean data with contemporary float or drifter technology.
Vertical velocities are typically very weak, and Lagrangian drifters that are
able to follow water parcels in 3-D are expensive and have only been deployed
in small numbers (D'Asaro et al., 1996; D'Asaro, 2015).

As a second basis for comparison, we computed the dispersion over time of initially small clusters of trajectories (Sect. 4) as they spread across isosurfaces of the background streamfunction. Background turbulent diffusion is simulated as a Lagrangian random walk based on spatially uniform diffusivity. We consider the dispersion characteristics that arise when this representation of turbulent diffusion is added to a background flow with no chaotic advection and compare it to flows that are undergoing chaotic advection but lack turbulent diffusion. Since the chaotic regions occupy sub-volumes of the entire eddy, spread of trajectories or tracers due to turbulent diffusion will eventually surpass that due to chaotic advection: chaos alone cannot distribute parcels across Lagrangian boundaries. However, it remains meaningful to compare the rate of spreading of parcels at earlier times. One immediate observation is that the character of ensemble spreading is qualitatively different for advective as opposed to diffusive perturbations. For the former, the spreading rate is significantly enhanced at some key times when trajectories pass near strong hyperbolic regions. In the latter case, the spread grows similarly to the square root of time at all times.

When the eddy is moderately shallow (*E*=0.125), there are many instances in
which chaotic advection in the chaotic sea dominates turbulent diffusion,
even at the higher ranges of turbulent diffusivity. When the perturbation
strength is moderately large (*ϵ*=0.08, ${x}_{\mathrm{0}}=-\mathrm{0.02}$), chaotic advection
produces more rapid spreading than diffusion for two of three diffusivities
considered (pink curve in Fig. 8a). Even when the
perturbation strength is small (*ϵ*=0.01), spread due to chaotic
advection in the chaotic sea (red curve in Fig. 8a) is of
a comparable order to turbulent diffusion at the lowest *k* values considered
(light blue curve in Fig. 8a). These results are in agreement
with the Batchelor scale analysis.

When the eddy is deeper (*E*=0.02) spreading due to turbulent diffusion in the
chaotic sea and resonant regions generally dominates over spreading due to
chaotic advection. This holds even when the perturbation strength is
moderately large (*ϵ*=0.08). These results are not in strict agreement
with the Batchelor scale analysis (Fig. 6) result that the
dimension of the chaotic sea is greater or equal to that of the Lagrangian
Batchelor scale for deeper eddies. To reconcile these inconsistencies, note
that as *E* gets small, a greater percentage of the eddy volume becomes
occupied by an inviscid, vertically rigid interior. For a very small *E*,
parcels experience relatively low levels of strain while rising or descending
through the region. When a fluid parcel nears the top or bottom boundary,
however, it become vertically squashed and horizontally stretched, suggesting
that the main contribution to *λ*_{3} comes from close encounters with
these boundaries. A Batchelor scale that is based only on a single parameter
measuring the time-averaged contraction over several overturning cycles may
be too simplistic when a parcel divides its time between kinematically
distinct regions.

This method of comparison based on parcel spreading has several advantages over the Batchelor scale. First, it offers a direct measure of fluid stirring. Also, it reveals information about the time history of dispersion that is hidden in the Lagrangian Batchelor scale analysis. Disadvantages include the fact that the analysis, as presented, does not account for scale-dependent diffusivity. Also, like the Batchelor scale analysis, it requires the tracking of fluid parcels in 3-D, something that is currently difficult in the ocean.

The third method for comparison (Sect. 5) differs
from the first two in that it is based on metrics of irreversible property
exchange (mixing). These metrics consist of the Nakamura effective
diffusivity, *κ*_{eff}, and a volume-integrated tracer variance
function, *χ*^{2}. We consider a flow with a given background turbulent
diffusivity, *k*, and calculate how much the irreversible property exchange
is amplified as a result of chaotic stirring. The volume-integrated
*κ*_{eff} and *χ*^{2} both depend on time and show rapid initial
growth, a result of filamentation of an initially smooth tracer distribution.
Growth is arrested when diffusion begins to dominate due to the enhanced
gradients produced by the filamentation process, at which time both measures,
*κ*_{eff} and *χ*^{2}, reach peak values. This is followed by a long
period in which *χ*^{2} slowly diminishes to zero and the volume integral of
*κ*_{eff} reaches a nearly constant value. In most cases, chaotic
advection leads to more rapid initial growth, a lower peak value for both
measures, and a larger long-term, near-equilibrium value of *κ*_{eff}.
In weakly perturbed cases, the differences in initial growth and peak value
of *κ*_{eff} are minor, usually on the order of 10 % or 20 %, while
differences in the longer term, near-equilibrium value of *κ*_{eff} are
more significant. For strongly perturbed cases the initial growth is an order
of magnitude larger and the amplification in the long-term value of
*κ*_{eff} is larger by a factor of 2 than in the unperturbed case.

The spatial structure of *κ*_{eff} also yields interesting information,
though one must be aware of the caveat that the local value is due to
nonlocal processes. The chaotic sea region generally has enhanced values
compared to the interior and its resonant regions. Under weak perturbation,
maximum values of *κ*_{eff} were smaller than in the unperturbed case,
but the spatial extent of the intermediate values was larger, leading to the
enhanced volume-integrated values discussed above. Larger changes in
*κ*_{eff} are evident for lower *k* due to the occurrence of more
numerous small-scale filaments. With a larger perturbation, chaotic
advection dramatically changes the effective diffusivity, but there are also
stronger barriers present, evident from isolated areas with different tracer
concentrations. We conclude that the spatial structures of chaotic and regular
regions can play an important role in how a tracer is distributed.

The use of effective diffusivity as a metric has several advantages. First of all, it provides a direct measure of irreversible property exchange between regions with different dye concentration. Its time history leads to insights about the evolution of mixing and, in particular, the time periods when chaotic advection is most relevant. Also, it can be measured, at least in principle, by performing an ocean dye release and measuring the dye concentration along sections that cut through the dye plume at different depths or angles, all in an attempt to recreate a concentration map in 3-D. Of the three methods proposed herein, it would appear to be the one most testable by ocean observations. The main disadvantage of effective diffusivity is that it requires the background diffusivity to be constant, which is strictly true only if the diffusivity is interpreted as the molecular diffusivity.

In this work, we examined the relative strengths of advection and diffusion for the redistribution of a passive tracer in a rotating cylinder flow as an analogue for an overturning submesoscale eddy. Since a major challenge of this work was developing ways of thinking about the competition between chaotic advection and turbulent diffusion, the numerical experiments described in this paper have been necessarily idealized. Although the focus of this current paper is on the behavior of a steady 3-D eddy flow subject to a turbulent diffusion, similar results are expected to hold for 3-D eddy flows with time-periodic and time-quasi-periodic behavior. Exploration with models that are more realistic for the ocean presents a number of challenges, including the development of more anisotropic and spatially varying representations of turbulence to account for differences between the ocean surface mixed layer and the stratified fluid underneath. In addition, finite eddy lifetimes must be confronted, as a separation of timescales between feature lifetimes and the periods of trajectories within them is needed for these analyses.

Matlab codes for performing trajectory integrations and for making most figures are available at Github (Brett, 2018, https://doi.org/10.5281/zenodo.1560663). Data (code outputs) for making the figures are available through Zenodo (Brett and Wang, 2018, https://doi.org/10.5281/zenodo.1560204). Outputs of NS simulations using NEK5000 are included, but the source code is available instead at the following link: https://nek5000.mcs.anl.gov, last access: 3 September 2018.

Here we provide detail about the fixed points, and their bifurcations, of the background velocity field in the kinematic model of the rotating cylinder. Then we present the bifurcation diagram and an example of the flow with many fixed points in the overturning streamfunction.

The overturning streamfunction is described by Eqs. (2)–(7), with radial and
vertical velocities (Eqs. 8–9) and azimuthal velocity (Eq. 10). All three velocity
components are zero at *z*=0 and *r*=*a*. The azimuthal velocity *V* only has
one other zero at *r*=0. However, additional points with zero
vertical and radial velocity exist, which correspond to circular periodic orbits in
the horizontal plane and which we refer to as *r**z*-fixed points.

All *r**z*-fixed points in the interior occur at *r*=0.5*a* because this is the
only place where *W*=0. Finding the *r**z*-fixed points is thus equivalent to
finding points in *z* where $U(r=\mathrm{0.5}a,z)=\mathrm{0}$. One such point exists for all
*E* at *z*=0.5. Additional *r**z*-fixed points appear through pitchfork
bifurcations, where new pairs split from *z*=0.5 and move apart in *z* as *E*
decreases from 1 (Fig. A1).

It is possible to classify the *r**z*-fixed points as elliptic or hyperbolic
according to their behavior in the *r*–*z* plane: the overturning
streamfunction is a local maximum in both *z* and *r* at elliptic points and
a saddle, i.e., a minimum in *r* but a maximum in *z*, at hyperbolic points.
At *E*=1, the only stationary point is at $(r,z)=(\mathrm{0.5},\mathrm{0.5})$, and it is
elliptic. As *E* decreases to about 1∕62, the first bifurcation creates two
elliptic points above and below the now-hyperbolic central point at
$(r,z)=(\mathrm{0.5},\mathrm{0.5})$. As *E* decreases, the newly created points move away
vertically from the central point, until the next bifurcation creates two new
hyperbolic points, and the central fixed point becomes elliptic again. This
process continues; the number of fixed points increases as *E* decreases
through a repeated pitchfork bifurcations of the $(r,z)=(\mathrm{0.5},\mathrm{0.5})$ fixed
point. As these bifurcations occur, their effects remain within a region
bounded by trajectories between the first pair of hyperbolic points, meaning
that their effects are quite local. The spreading of the first pair of
hyperbolic points, and not the total increase in *r**z*-fixed points, causes the
increasing vertical homogeneity of the flow with decreasing *E*, which appears
qualitatively similar to Taylor columns. An example with nine *r**z*-fixed
points is shown in Fig. A2 for *E*=0.00125; the central point
is now elliptic. Trajectories in the vertical plane are level curves of the
streamfunction; these show the elliptic and hyperbolic nature of the
*r**z*-fixed points, where trajectories near an elliptic point remain nearby
but trajectories near a hyperbolic point may travel a long distance before
returning or may move toward another hyperbolic point.

In this Appendix, we present the derivation of the evolution of a three-dimensional tracer in a steady linear strain flow. This result was used in the main text to show that the thinnest width of the Gaussian tracer distribution will asymptotically approach the Lagrangian Batchelor scale. We start with the definitions of the velocity field, the tracer evolution equation, and the form of the solution. Then we derive the full time-dependent solution for the tracer distribution.

We are solving for the evolution of tracer concentration, *C*, with a solution in the form of a Gaussian
function:

where
*c*_{max} is the maximum concentration and *α*, *β*, and *γ* are
the reciprocal of the standard deviations in each direction. In the
Lagrangian frame of reference that is moving with the center of mass of the
tracer, these four parameters are dependent on time but not space. The
smallest width of the distribution is $\mathit{\sigma}=\mathrm{1}/\mathit{\alpha}$, and in the main
text, we have used the fact that it has a stable fixed point
$\mathit{\sigma}=\sqrt{\mathit{\kappa}/\left|{\mathit{\lambda}}_{\mathrm{3}}\right|}$, where *λ*_{3} is the contraction rate
of the velocity field. We are now going to formally prove it.

The velocities are defined in the Lagrangian frame by

with ** x**(

*x*_{0},

*t*) indicating the initial position

*x*_{0}of the water parcel at

*t*=0. The Lagrangian tracer evolution equation is

where *κ* is the diffusivity.

The form of *C* and the tracer evolution equation allow us to find
differential equations for each of our four parameters, which are

The width parameters' equations are nonlinear, but when rewritten in terms like *α*^{−2}, they give the following:

which are Bernoulli equations, solvable with integrating factors, giving

where subscript 0 indicates the value at *t*=0. The differences in these
equations are due to the different signs of each *λ*, with the ambiguity
of the sign of *λ*_{2} preventing its factoring.

The *c*_{max} equation depends on the width parameters and is not simple to
solve directly. However, a careful inspection shows that *c*_{max}∕(*α**β**γ*) is conserved, so we can write

For anyone in doubt, we plug in this solution to check it:

The full solution for the tracer concentration *C* then has been fully solved by
Eq. (B1),
with $\mathit{\alpha},\mathit{\beta},\mathit{\gamma}$, and *c*_{max} given by Eqs. (B15)–(B18).

For a three-dimensional Gaussian tracer advected by a linear strain field in the presence of constant diffusivity, in the Lagrangian frame the width of the tracer distribution will increase in the stretching direction or directions forever but reach a fixed value in the contracting direction or directions.

For the axially symmetric rotating cylinder flow at long times, the dye
contours resemble nested tori, although with cross sections that are
somewhat between a circle and a square. Here, we derive the expected limit of
∫*κ*_{eff}d*V*, assuming that the dye iso-contours at late times are
nested tori with a circular cross section and that the gradient of the dye
concentration is constant along each torus. In this case the effective
diffusivity on each torus is *κ*_{eff}=*k**A*^{2}, the background diffusivity
multiplied by the squared surface area of a torus.

Recall that the volume of a circular torus is

where *r* is the radius of the circular cross section and *R* is the
distance from the center of mass of the torus to the center of the
cross section. The surface area is

Noting that ${A}_{\mathrm{ct}}=\mathrm{d}{V}_{\mathrm{ct}}/\mathrm{d}r$, we can calculate the volume-integrated effective diffusivity as

using *R*=0.5 and *r*_{max}=0.5. This circular-torus-based result gives a
lower bound, because there is still volume outside the largest torus that
fits in the cylinder, and the final cross sections are somewhat square, thus
having a larger surface area per volume.

GJB carried out the main investigation and formal analyses of this work and wrote the initial draft. LP and IR acquired funding, formulated the kinematic model, and supervised the project. PW performed the NS simulations and contributed to their design and analysis. LP, IR, and PW revised the written work, with LP contributing significant portions of the final text.

The authors declare that they have no conflict of interest.

This work was supported by DOD (MURI) grant no. N0004110087 as well as the US National Science Foundation grants OCE-1154641 and OCE-1558806 and the National Aeronautics and Space Administration (NASA) grant no. NNX14AH29G. Genevieve Jay Brett received additional support from the Woods Hole Oceanographic Institution Academic Programs Office. We also thank Marianna Linz, Dan Collura, and four anonymous reviewers for their constructive input on how to present this material.

This paper was edited by Vicente Perez-Munuzuri and reviewed by four anonymous referees.

Abernathey, R., Marshall, J., Mazloff, M., and Shuckburgh, E.: Enhancement of mesoscale eddy stirring at steering levels in the Southern Ocean, J. Phys. Oceanogr., 40, 170–184, 2010. a

Aref, H.: Stirring by chaotic advection, J. Fluid Mech., 143, 1–21, 1984. a

Branicki, M. and Kirwan Jr., A.: Stirring: the Eckart paradigm revisited, Int. J. Eng. Sci., 48, 1027–1042, 2010. a

Brett, G.: 3D rotating cylinder eddy codes, https://doi.org/10.5281/zenodo.1560663, 2018. a

Brett, G. and Wang, P.: 3D rotating cylinder eddy data, https://doi.org/10.5281/zenodo.1560204, 2018. a

Casati, G. and Ford, J.: Stochastic behavior in classical and quantum hamiltonian systems, in: Stochastic Behavior in Classical and Quantum Hamiltonian Systems, 93, Springer-Verlag, Berlin Heidelberg, 1979. a

Chirikov, B. V.: Research concerning the theory of non-linear resonance and stochasticity, Tech. rep., CM-P00100691, original report for the Nuclear Physics Institute of the Siberian section of the USSR Academy of Sciences, 1969, location Novosibirsk, translation by: Sanders, A. T., CERN Translation 71-40, 1971, Geneva, 1971. a

Chirikov, B. V.: A universal instability of many-dimensional oscillator systems, Phys. Rep., 52, 263–379, 1979. a

Coulliette, C. and Wiggins, S.: Intergyre transport in a wind-driven, quasigeostrophic double gyre: An application of lobe dynamics, Nonlin. Processes Geophys., 8, 69–94, https://doi.org/10.5194/npg-8-69-2001, 2001. a

D'Asaro, E.: Surface wave measurements from subsurface floats, J. Atmos. Ocean. Tech., 32, 816–827, 2015. a

D'Asaro, E. A., Farmer, D. M., Osse, J. T., and Dairiki, G. T.: A Lagrangian float, J. Atmos. Ocean. Tech., 13, 1230–1246, 1996. a

Deese, H. E., Pratt, L. J., and Helfrich, K. R.: A laboratory model of exchange and mixing between western boundary layers and subbasin recirculation gyres, J. Phys. Oceanogr., 32, 1870–1889, 2002. a

Dombre, T., Frisch, U., Greene, J. M., Hénon, M., Mehr, A., and Soward, A. M.: Chaotic streamlines in the ABC flows, J. Fluid Mech., 167, 353–391, 1986. a

Fereday, D. and Haynes, P.: Scalar decay in two-dimensional chaotic advection and Batchelor-regime turbulence, Phys. Fluids, 16, 4359–4370, 2004. a

Fischer, P. F.: An overlapping Schwarz method for spectral element solution of the incompressible Navier–Stokes equations, J. Comput. Phys., 133, 84–101, 1997. a

Flierl, G. R. and Woods, N. W.: Copepod aggregations: influences of physics and collective behavior, J. Stat. Phys., 158, 665–698, 2015. a, b

Fountain, G., Khakhar, D., Mezić, I., and Ottino, J.: Chaotic mixing in a bounded three-dimensional flow, J. Fluid Mech., 417, 265–301, 2000. a, b, c

Froyland, G., Padberg, K., England, M. H., and Treguier, A. M.: Detection of coherent oceanic structures via transfer operators, Phys. Rev. Lett., 98, 224503, https://doi.org/10.1103/PhysRevLett.98.224503, 2007. a

Froyland, G., Horenkamp, C., Rossi, V., Santitissadeekorn, N., and Gupta, A. S.: Three-dimensional characterization and tracking of an Agulhas Ring, Ocean Model., 52, 69–75, 2012. a

Greenspan, H. P.: The theory of rotating fluids, CUP Archive, Cambridge, UK, 1968. a, b

Griffies, S. M., Winton, M., Anderson, W. G., Benson, R., Delworth, T. L., Dufour, C. O., Dunne, J. P., Goddard, P., Morrison, A. K., Rosati, A., Wittenberg, A. T., Yin, J., and Zhang, R.: Impacts on ocean heat from transient mesoscale eddies in a hierarchy of climate models, J. Climate, 28, 952–977, 2015. a

Gromeka, I.: Some cases of incompressible fluids motion, Scientific notes of the Kazan University, Kazan, Russia, pp. 76–148, 1881. a

Hadjighasem, A., Farazmand, M., Blazevski, D., Froyland, G., and Haller, G.: A critical comparison of Lagrangian methods for coherent structure detection, Chaos: An Interdisciplinary Journal of Nonlinear Science, 27, 053104, https://doi.org/10.1063/1.4982720, 2017. a

Hallberg, R.: Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects, Ocean Model., 72, 92–103, 2013. a

Haller, G.: Lagrangian coherent structures from approximate velocity data, Phys. Fluids, 14, 1851–1861, 2002. a

Haller, G.: Lagrangian coherent structures, Annu. Rev. Fluid Mech., 47, 137–162, 2015. a

Haller, G. and Beron-Vera, F. J.: Geodesic theory of transport barriers in two-dimensional flows, Physica D, 241, 1680–1702, 2012. a

Haller, G. and Beron-Vera, F. J.: Coherent Lagrangian vortices: The black holes of turbulence, J. Fluid Mech., 731, R4, https://doi.org/10.1017/jfm.2013.391, 2013. a

Haller, G., Karrasch, D., and Kogelbauer, F.: Material barriers to diffusive and stochastic transport, P. Natl. Acad. Sci. USA, 115, 9074–9079, 2018. a

Haynes, P. and Shuckburgh, E.: Effective diffusivity as a diagnostic of atmospheric transport: 2. Troposphere and lower stratosphere, J. Geophys. Res.-Atmos., 105, 22795–22810, 2000. a

Ledwell, J. R., Watson, A. J., and Law, C. S.: Evidence for slow mixing across the pycnocline from an open-ocean tracer-release experiment, Nature, 364, 701–703, 1993. a, b

Ledwell, J. R., Watson, A. J., and Law, C. S.: Mixing of a tracer in the pycnocline, J. Geophys. Res.-Oceans, 103, 21499–21529, 1998. a, b

Lenn, Y.-D. and Chereskin, T. K.: Observations of Ekman currents in the Southern Ocean, J. Phys. Oceanogr., 39, 768–779, 2009. a

Lopez, J. M. and Marques, F.: Sidewall boundary layer instabilities in a rapidly rotating cylinder driven by a differentially corotating lid, Phys. Fluids, 22, 114109, https://doi.org/10.1063/1.3517292, 2010. a

Mahadevan, A.: The impact of submesoscale physics on primary productivity of plankton, Annu. Rev. Mar. Sci., 8, 161–184, 2016. a

Malhotra, N., Mezić, I., and Wiggins, S.: Patchiness: A new diagnostic for Lagrangian trajectory analysis in time-dependent fluid flows, Int. J. Bifurcat. Chaos, 8, 1053–1093, 1998. a

Miller, P. D., Pratt, L. J., Helfrich, K. R., and Jones, C. K.: Chaotic transport of mass and potential vorticity for an island recirculation, J. Phys. Oceanogr., 32, 80–102, 2002. a

Nakamura, N.: Two-dimensional mixing, edge formation, and permeability diagnosed in an area coordinate, J. Atmos. Sci., 53, 1524–1537, 1996. a, b

Nakamura, N. and Ma, J.: Modified Lagrangian-mean diagnostics of the stratospheric polar vortices: 2. Nitrous oxide and seasonal barrier migration in the cryogenic limb array etalon spectrometer and SKYHI general circulation model, J. Geophys. Res.-Atmos., 102, 25721–25735, 1997. a

Ngan, K. and Shepherd, T. G.: Chaotic mixing and transport in Rossby-wave critical layers, J. Fluid Mech., 334, 315–351, 1997. a

Okubo, A.: Oceanic diffusion diagrams, in: Deep Sea Research and Oceanographic Abstracts, Elsevier, Amsterdam, the Netherlands, 18, 789–802, 1971. a

Olascoaga, M. J. and Haller, G.: Forecasting sudden changes in environmental pollution patterns, P. Natl. Acad. Sci. USA, 13, 4738–4743, https://doi.org/10.1073/pnas.1118574109, 2012. a

Ottino, J.: Mixing, chaotic advection, and turbulence, Annu. Rev. Fluid Mech., 22, 207–254, 1990. a

Pattanayak, A. K.: Characterizing the metastable balance between chaos and diffusion, Physica D, 148, 1–19, 2001. a

Pierrehumbert, R.: Tracer microstructure in the large-eddy dominated regime, Chaos Soliton. Fract., 4, 1091–1110, 1994. a

Poje, A. and Haller, G.: Geometry of cross-stream mixing in a double-gyre ocean model, J. Phys. Oceanogr., 29, 1649–1665, 1999. a

Polvani, L. M., Waugh, D., and Plumb, R. A.: On the subtropical edge of the stratospheric surf zone, J. Atmos. Sci., 52, 1288–1309, 1995. a

Pratt, L. J., Rypina, I. I., Özgökmen, T., Wang, P., Childs, H., and Bebieva, Y.: Chaotic advection in a steady, three-dimensional, Ekman-driven eddy, J. Fluid Mech., 738, 143–183, 2014. a, b

Rogerson, A. M., Miller, P. D., Pratt, L. J., and Jones, C. K. R. T.: Lagrangian motion and fluid exchange in a barotropic meandering jet, J. Phys. Oceanogr., 29, 2635–2655, 1999. a

Rom-Kedar, V., Leonard, A., and Wiggins, S.: An analytical study of transport, mixing and chaos in an unsteady vortical flow, J. Fluid Mech., 214, 347–394, 1990. a

Rypina, I., Brown, M., Beron-Vera, F., Kocak, H., Olascoaga, M., and Udovydchenkov, I.: On the Lagrangian dynamics of atmospheric zonal jets and the permeability of the stratospheric polar vortex, J. Atmos. Sci., 64, 3595–3610, 2007. a

Rypina, I., Pratt, L., Wang, P., Özgökmen, T., and Mezic, I.: Resonance phenomena in a time-dependent, three-dimensional model of an idealized eddy, Chaos: An Interdisciplinary Journal of Nonlinear Science, 25, 087401, https://doi.org/10.1063/1.4916086, 2015. a, b, c, d

Rypina, I. I. and Pratt, L. J.: Trajectory encounter volume as a diagnostic of mixing potential in fluid flows, Nonlin. Processes Geophys., 24, 189–202, https://doi.org/10.5194/npg-24-189-2017, 2017. a

Rypina, I. I., Brown, M. G., and Koçak, H.: Transport in an idealized three-gyre system with application to the Adriatic Sea, J. Phys. Oceanogr., 39, 675–690, 2009. a

Rypina, I. I., Pratt, L. J., Pullen, J., Levin, J., and Gordon, A. L.: Chaotic advection in an archipelago, J. Phys. Oceanogr., 40, 1988–2006, 2010. a, b

Rypina, I. I., Pratt, L. J., and Lozier, M. S.: Near-surface transport pathways in the North Atlantic Ocean: Looking for throughput from the subtropical to the subpolar gyre, J. Phys. Oceanogr., 41, 911–925, 2011a. a

Rypina, I. I., Scott, S. E., Pratt, L. J., and Brown, M. G.: Investigating the connection between complexity of isolated trajectories and Lagrangian coherent structures, Nonlin. Processes Geophys., 18, 977–987, https://doi.org/10.5194/npg-18-977-2011, 2011b. a

Rypina, I. I., Kamenkovich, I., Berloff, P., and Pratt, L. J.: Eddy-induced particle dispersion in the near-surface North Atlantic, J. Phys. Oceanogr., 42, 2206–2228, 2012. a

Rypina, I. I., Llewellyn Smith, S. G., and Pratt, L. J.: Connection between encounter volume and diffusivity in geophysical flows, Nonlin. Processes Geophys., 25, 267–278, https://doi.org/10.5194/npg-25-267-2018, 2018. a

Samelson, R.: Fluid exchange across a meandering jet, J. Phys. Oceanogr., 22, 431–444, 1992. a

Samelson, R. M. and Wiggins, S.: Lagrangian transport in geophysical jets and waves: The dynamical systems approach, vol. 31, Springer Science & Business Media, Berlin, Germany, 2006. a

Sayol, J.-M., Orfila, A., Simarro, G., López, C., Renault, L., Galán, A., and Conti, D.: Sea surface transport in the Western Mediterranean Sea: A Lagrangian perspective, J. Geophys. Res.-Oceans, 118, 6371–6384, 2013. a

Shadden, S. C., Lekien, F., and Marsden, J. E.: Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows, Physica D, 212, 271–304, 2005. a

Shepherd, T. G., Koshyk, J. N., and Ngan, K.: On the nature of large-scale mixing in the stratosphere and mesosphere, J. Geophys. Res.-Atmos., 105, 12433–12446, 2000. a

Shuckburgh, E. and Haynes, P.: Diagnosing transport and mixing using a tracer-based coordinate system, Phys. Fluids, 15, 3342–3357, 2003. a

Solomon, T. and Mezić, I.: Uniform resonant chaotic mixing in fluid flows, Nature, 425, 376–380, 2003. a

Son, D.: Turbulent decay of a passive scalar in the Batchelor limit: Exact results from a quantum-mechanical approach, Phys. Rev. E, 59, R3811, https://doi.org/10.1103/PhysRevE.59.R3811, 1999. a

Thiffeault, J.-L.: Stretching and curvature of material lines in chaotic flows, Physica D, 198, 169–181, 2004. a

Yuan, G. C., Pratt, L., and Jones, C.: Cross-jet Lagrangian transport and mixing in a 2 1∕2-layer model, J. Phys. Oceanogr., 34, 1991–2005, 2004. a

Zambianchi, E. and Griffa, A.: Effects of finite scales of turbulence on dispersion estimates, J. Marine Res., 52, 129–148, 1994. a

- Abstract
- Introduction
- Models
- Lagrangian Batchelor scale
- Particle dispersion
- Tracer simulations and Nakamura effective diffusivity
- Conclusions
- Code and data availability
- Appendix A: Bifurcation analysis of fixed points of the background streamfunction
- Appendix B: Gaussian tracer in linear strain
- Appendix C: Long-time limit of effective diffusivity for the axially symmetric rotating cylinder flow
- Author contributions
- Competing interests
- Acknowledgements
- Review statement
- References

- Abstract
- Introduction
- Models
- Lagrangian Batchelor scale
- Particle dispersion
- Tracer simulations and Nakamura effective diffusivity
- Conclusions
- Code and data availability
- Appendix A: Bifurcation analysis of fixed points of the background streamfunction
- Appendix B: Gaussian tracer in linear strain
- Appendix C: Long-time limit of effective diffusivity for the axially symmetric rotating cylinder flow
- Author contributions
- Competing interests
- Acknowledgements
- Review statement
- References