Competition between Chaotic Advection and Diffusion : Stirring and Mixing in a 3 D Eddy Model

The importance of chaotic advection relative to turbulent diffusion is investigated in an idealized model of a 3D 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) 5 stochastic noise, added to a deterministic velocity field, or 3) explicit and implicit diffusion in a spectral numerical model of 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 depends on factors such as the degree of axial asymmetry of the eddy 10 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. 15 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.

A feature that is intriguing and quite common in these studies is that Lagrangian chaos is confined to certain subregions 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 tra-Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
jectories, but the non-chaotic bands act as barriers that confine the stirring.In textbook examples, including areapreserving 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(Chirikov, , 1979;;Casati and Ford, 1979;Gromeka, 1881;Dombre et al., 1986).In finitetime 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., 2007Froyland et al., , 2012;;Rypina and Pratt, 2017;Rypina et al., 2018Rypina et al., , 2011b;;Hadjighasem et al., 2017;Haller andBeron-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 smallscale 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 timeperiodic, 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.

Models
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 (δ / ) 1), and the Ekman layers are relatively thin, E 1.In this case, a linear, asymptotic solution is available (Greenspan, 1968, andAppendix 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, θ, z), with (1 ≥ z ≥ 0) and (r ≤ a), where a is the width-to-height ratio of the domain.The background flow has ∂/∂θ = 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, www.nonlin-processes-geophys.net/26/37/2019/ Nonlin.Processes Geophys., 26, 37-60, 2019 and the constants are defined by (5) In the limit of infinite cylinder radius, a → ∞, the radial portion of the streamfunction, R(r) = r 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) = (0.5, 0.5) has a period of 16π ≈ 50.At the maximum azimuthal velocity, which occurs at r = aH /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.75a).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 > 1/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 < 1/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).

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 rand 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 rand z-dependent background overturning streamfunction; the velocities from the two are added together.The perturbation velocities in x and y are    = −4y 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 cer-tain choices of , but with < 0.05 these locations are all very close to the boundaries of the cylinder.

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 symmetrybreaking 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, Ro = δ / .The kinematic model parameters are the Ekman number, E, the aspect ratio, a, the perturbation offset parameter, x 0 , and the strength of www.nonlin-processes-geophys.net/26/37/2019/ Nonlin.Processes Geophys., 26, 37-60, 2019 the perturbation, .For matching the kinematic model to the NS simulation, we set α = a = 1 and examine four Ekman numbers used in P2014, E ∈ {0.25, 0.125, 0.02, 0.0005}.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 0 = −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) = (0.5, 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.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, x = x 0 e λ t . (16) Since separation rates between trajectories are generally different in different directions, in 3-D flows there are three FTLEs that can be ordered λ 1 ≥ λ 2 ≥ λ 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 λ 1 + λ 2 + λ 3 = 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 ith 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 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.
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 2 × 10 −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.2m 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.

Particle dispersion
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 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, www.nonlin-processes-geophys.net/26/37/2019/ Nonlin.Processes Geophys., 26, 37-60, 2019 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: κ = s 2 /2 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 κ ∈ [10 −4 , 10 −2 ] m 2 s −1 across the four Ekman numbers examined, which is nondimensionally κ ∈ {10 −6 , 3 × 10 −5 }.As our primary example, we will discuss the level of diffusivity κ = 10 −6 .This diffusivity requires a certain step size s for the stochastic perturbation, which relates to the distribution of u by s = σ t/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, κ = 10 −7 and σ = 0.013, and a larger one, κ = 10 −5 and σ = 0.13.The stochastic perturbation with κ = 10 −6 has kinetic energy (integrated over the cylinder) about the same as the background flow: (u ) 2 ≈ (U 2 b ) ≈ 0.63.The perturbation with κ = 10 −7 has kinetic energy about the same as the deterministic perturbation with = 0.01 and x 0 = −0.5, such that (u d ) 2 ≈ (u s ) 2 ≈ 0.075, where u d is the deterministic perturbation velocity and u s 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) = (0.1, 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 per-turbation 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 chaoticadvection-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 en-sembles 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) = (0.4,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.

Tracer simulations and Nakamura effective diffusivity
In this section we analyze the effects of the symmetrybreaking, 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., 1993Ledwell et al., , 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 advectiondiffusion 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 efwww.nonlin-processes-geophys.net/26/37/2019/ Nonlin.Processes Geophys., 26, 37-60, 2019 fective diffusivity can be written as where C is tracer concentration, V is volume, and 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 ( 25) 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 sur-face area squared in the rare situation where |∇C| 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 volumeintegrated 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 symmetrybreaking 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 0 = −0.02case is what was used to compare Poincaré sections with the kinematic model, so qualitative features match the = 0.01 cases.The X 0 = −0.16case 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 = 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 = 10 −4 , X 0 ∈ {0, −0.02, −0.16} and E = 0.02, and k ∈ {10 −4 , 10 −6 } and X 0 ∈ {0, −0.02} 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 = 10 −4 and k = 10 −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 0 = −0.16, the maximum is larger than with either X 0 = 0 or X 0 = −0.02due 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 ∂C/∂V term in κ eff and the |C| 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, |∇C| would be constant along the surfaces, and then κ eff = kA 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) = (0.5, 0.5) multiplied by the background diffusivity is kπ 6 /8, which we expect to be the minimum for κ eff dV 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 dV 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 0 = −0.16,and k = 10 −4 , which has the most asymmetric dye contours (Fig. 11i); here, the long-time value of κ eff dV 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 be- come 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 = 10 −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 0 = −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 0 = −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 mod-erate 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 volumeintegrated 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 = 10 −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 un-  www.nonlin-processes-geophys.net/26/37/2019/ Nonlin.Processes Geophys., 26, 37-60, 2019 perturbed 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 = 10 −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 = 10 −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 ).

Conclusions
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 subvolumes 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 0 = −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 ac-count 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, nearequilibrium 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.
Code and data availability.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.background streamfunction 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 rz-fixed points.
All rz-fixed points in the interior occur at r = 0.5a because this is the only place where W = 0. Finding the rzfixed points is thus equivalent to finding points in z where U (r = 0.5a, z) = 0.One such point exists for all E at z = 0.5.Additional rz-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 rz-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) = (0.5, 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) = (0.5, 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) = (0.5, 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 rz-fixed points, causes the increasing vertical homogeneity of the flow with decreasing E, which appears qualitatively similar to Taylor columns.An example with nine rz-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 rz-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.

Appendix B: Gaussian tracer in linear strain
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 timedependent 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 σ = 1/α, and in the main text, we have used the fact that it has a stable fixed point σ = √ κ/|λ 3 |, where λ 3 is the contraction rate of the velocity field.We are now going to formally prove it.
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 c max (t) = c 0 α(t)β(t)γ (t). (B18) 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 α, β, γ , 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.

Figure 1 .
Figure1.Sketch of the qualitative velocity field (Eqs.7-9).Ekman layers at the top and bottom are where flow has a larger radial component.is the rotation rate of the system.X 0 is the offset between the lid and cylinder rotational centers, as set for the Navier-Stokes simulations.

Figure 3 .
Figure 3. Structures in the kinematic model and dynamical simulation for Ekman numbers of 0.25 (a, c) and 0.125 (b, d).(a, b) Poincaré maps from Pratt et al. (Fig. 10 in 2014), resulting from a dynamically consistent numerical simulation.(c, d) Poincaré maps (black) and largest FTLEs (color) resulting from our non-dynamically consistent kinematic analytic model, with = 0.01 and x 0 either −0.5 (c) or −0.9 (d); in color are maximum FTLEs calculated for the kinematic model with integration time 400.In (d), red oval approximately separates the resonant and regular layers (inside) from the chaotic sea region (outside), with the blue line segment showing the width of the chaotic sea.The blue diamond shows the width of an island, which is also the width of the resonant layer.

Figure 4 .
Figure 4. Structures in the kinematic model (c, d) and dynamical simulation (a, b) for Ekman numbers of 0.02 (a, c) and 0.0005 (b, d), with same format as Fig. 3.

Figure 5 .
Figure 5.An initial sphere in a linear strain field evolving into an ellipsoid during a time of 1. Ellipsoid axes marked by bars, with figure axes ticks showing their endpoint values.Velocity field u = 1.5 + x, v = 0.5y, w = −1.5z.Color shows z values at t = 0..

Figure 6 .
Figure 6.Layer widths in blue, and Lagrangian Batchelor scale δ (Eq.22) in the same region in yellow.(a) Chaotic resonant region between islands; (b) the chaotic sea region.The diffusivities at the Batchelor scale (in m 2 s −1 ) are between 10 −4 and 6 × 10 −3 for the three larger Ekman numbers and between 1 × 10 −2 and 6 × 10 −2 for E = 0.0005.

Figure 7 .
Figure 7. Grey lines are individual trajectories in ψ starting from a sphere of radius 0.002 at (r, z) = (0.1, 0.5) with E = 0.125.Solid black curves are the mean; black dashed-dotted lines are ±1 standard deviations from the mean.

Figure 12 .
Figure 12.E = 0.02 Poincaré sections in the (a) x-z and (b) y-z planes in black.Polygons show the island (blue) and resonant (red) regions used for analysis (c) and (d), which show mean κ eff over time in these regions under both applied background diffusivities.

Figure
Figure A1.z positions of rz-fixed points.Black indicates elliptic points, blue hyperbolic points, and grey the neutrally stable points at the top and bottom.New fixed point pairs separate symmetrically from z = 0.5 as E decreases.At each bifurcation, the central fixed point changes stability.

Figure A2 .
Figure A2.Trajectories in the vertical plane for E = 0.00125 and a = 1.There are nine rz-fixed points along r = 0.5, marked with red stars.Note the closed curves between the outermost hyperbolic points which surround the interior five rz-fixed points; these limit the effects of those points to the local area.
κ is the diffusivity.The form of C and the tracer evolution equation allow us to find differential equations for each of our four parameters, 1 γ − κγ 3 .(B11)Thewidth parameters' equations are nonlinear, but when rewritten in terms like α −2 , they give the following:dα −2 dt = 2λ 3 α −2

Table 1 .
Nondimensional variables, their scale factors, and their dimensional equivalents.