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 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.

1 Introduction

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.

2 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

$$\begin{array}{}\text{(1)}& E=\left(\mathit{\nu }/\mathrm{\Omega }{H}^{\mathrm{2}}\right)={\left({\mathit{\delta }}_{\mathrm{E}}/H\right)}^{\mathrm{2}},\end{array}$$

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).

Figure 1 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₀ is the offset between the lid and cylinder rotational centers, as set for the Navier–Stokes simulations.

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 Ψ:

$$\begin{array}{}\text{(2)}& \mathrm{\Psi }=-E^{\mathrm{1}/\mathrm{2}}R\left(r\right)F\left(z\right),\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }F\left(z\right)& {\displaystyle }=A\left[\sin \right(\mathit{\zeta }\left)\sinh \right(\mathit{\zeta })-\cos (\mathit{\zeta }\left)\cosh \right(\mathit{\zeta }\left)\right]\\ \text{(3)}& {\displaystyle }& {\displaystyle }+B\left[\sin \right(\mathit{\zeta }\left)\sinh \right(\mathit{\zeta })+\cos (\mathit{\zeta }\left)\cosh \right(\mathit{\zeta }\left)\right]-D,\end{array}$$

where ζ is a transformed vertical coordinate,

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\zeta }={\displaystyle \frac{z-\mathrm{1}/\mathrm{2}}{E^{\mathrm{1}/\mathrm{2}}}},\end{array}$$

and the constants are defined by

$$\begin{array}{ll}{\displaystyle }A& {\displaystyle }={\displaystyle \frac{-\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{cS}{s^{\mathrm{2}}{C}^{\mathrm{2}}+c^{\mathrm{2}}{S}^{\mathrm{2}}}},B={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{sC}{s^{\mathrm{2}}{C}^{\mathrm{2}}+c^{\mathrm{2}}{S}^{\mathrm{2}}}},\\ {\displaystyle }D& {\displaystyle }=A(sS-cC)+B(sS+cC),\\ {\displaystyle }s& {\displaystyle }=\sin \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{E}^{\mathrm{1}/\mathrm{2}}}}\right),c=\cos \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{E}^{\mathrm{1}/\mathrm{2}}}}\right),S=\sinh \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{E}^{\mathrm{1}/\mathrm{2}}}}\right),\\ \text{(5)}& {\displaystyle }C& {\displaystyle }=\cosh \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{E}^{\mathrm{1}/\mathrm{2}}}}\right).\end{array}$$

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

$$\begin{array}{}\text{(6)}& R\left(r\right)=r(a-r{)}^{\mathrm{2}}/\mathrm{2},\end{array}$$

giving velocities

$$\begin{array}{}\text{(7)}& {\displaystyle }U& {\displaystyle }={\displaystyle \frac{-\partial \mathrm{\Psi }}{\partial z}}=r(a-r{)}^{\mathrm{2}}[A\sin \left(\mathit{\zeta }\right)\cosh \left(\mathit{\zeta }\right)+B\cos \left(\mathit{\zeta }\right)\sinh \left(\mathit{\zeta }\right)],\text{(8)}& {\displaystyle }W& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial r\mathrm{\Psi }}{\partial r}}=-(a-r)(a-\mathrm{2}r)E^{\mathrm{1}/\mathrm{2}}F\left(z\right),\end{array}$$

where U is radial and W is vertical.

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

$$\begin{array}{}\text{(9)}& V(r,z)=r(a-r{)}^{\mathrm{2}}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}+B\sin \left(\mathit{\zeta }\right)\cosh \left(\mathit{\zeta }\right)-A\cos \left(\mathit{\zeta }\right)\sinh \left(\mathit{\zeta }\right)\right].\end{array}$$

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⁻¹ 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⁻¹, 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.

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

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).

Figure 2 (a–c) Background overturning streamfunction for a=1; (a) E=0.125, (b) E=0.02, and (c) E=0.0005. Blue dots are rz-fixed points. (d) Horizontal perturbation streamfunction for γ=2, $x_{\mathrm{0}}=-\mathrm{0.5}$. Note that the center of rotation in the perturbation streamfunction is not at the origin.

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\mathrm{\Psi }}=\mathit{ϵ}{\displaystyle \frac{-\sinh (z/E^{\mathrm{1}/\mathrm{2}})}{\mathrm{2}\sinh (\mathrm{1}/E^{\mathrm{1}/\mathrm{2}})}}(a^{\mathrm{2}}-r^{\mathrm{2}})({\mathit{\gamma }}^{\mathrm{2}}{a}^{\mathrm{2}}-s^{\mathrm{2}}),\\ \text{(10)}& {\displaystyle }& {\displaystyle }s=\sqrt{(x-x_{\mathrm{0}}{)}^{\mathrm{2}}+y^{\mathrm{2}}}.\end{array}$$

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₀. If x₀=0, the disturbance is axially symmetric; if it is nonzero, the disturbance has an azimuthal variation in amplitude ϵx₀. 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

$$\begin{array}{ll}\text{(11)}& {\displaystyle }\stackrel{\mathrm{̃}}{u}& {\displaystyle }=\partial \stackrel{\mathrm{̃}}{\mathrm{\Psi }}/\partial y=\mathrm{4}y\mathit{ϵ}{\displaystyle \frac{\sinh (z/\sqrt{E})}{\sinh (\mathrm{1}/\sqrt{E})}}\left[(a^{\mathrm{2}}-r^{\mathrm{2}})+({\mathit{\gamma }}^{\mathrm{2}}{a}^{\mathrm{2}}-s^{\mathrm{2}})\right],{\displaystyle }\stackrel{\mathrm{̃}}{v}& {\displaystyle }=-\partial \stackrel{\mathrm{̃}}{\mathrm{\Psi }}/\partial x\\ \text{(12)}& {\displaystyle }& {\displaystyle }=-\mathrm{4}y\mathit{ϵ}{\displaystyle \frac{\sinh (z/\sqrt{E})}{\sinh (\mathrm{1}/\sqrt{E})}}\left[(x-x_{\mathrm{0}})(a^{\mathrm{2}}-r^{\mathrm{2}})+x({\mathit{\gamma }}^{\mathrm{2}}{a}^{\mathrm{2}}-s^{\mathrm{2}})\right].\end{array}$$

The corresponding azimuthal and radial velocity perturbations are

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{V}& {\displaystyle }=-\mathrm{2}\mathit{ϵ}{\displaystyle \frac{\sinh (z/\sqrt{E})}{\sinh (\mathrm{1}/\sqrt{E})}}\left[(a^{\mathrm{2}}-r^{\mathrm{2}})+({\mathit{\gamma }}^{\mathrm{2}}{a}^{\mathrm{2}}-s^{\mathrm{2}})\right.\\ \text{(13)}& {\displaystyle }& {\displaystyle }\left.-{\displaystyle \frac{x_{\mathrm{0}}}{r}}\cos \left(\mathit{\theta }\right)(a^{\mathrm{2}}-r^{\mathrm{2}})\right],\text{(14)}& {\displaystyle }\stackrel{\mathrm{̃}}{U}& {\displaystyle }=\mathrm{2}\mathit{ϵ}{x}_{\mathrm{0}}{\displaystyle \frac{\sinh (z/\sqrt{E})}{\sinh (\mathrm{1}/\sqrt{E})}}\sin \left(\mathit{\theta }\right)(a^{\mathrm{2}}-r^{\mathrm{2}}).\end{array}$$

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₀ of the lid's center (labeled in Fig. 1; not to be confused with x₀ 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₀, 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}$.

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₀ 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 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.

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.

3 Lagrangian Batchelor scale

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²time⁻¹), and the average exponential contraction rate λ₃, which quantifies the thinning of a filament due to chaotic advection (and has units of time⁻¹):

$$\begin{array}{}\text{(15)}& \mathit{\delta }=\sqrt{\mathit{\kappa }/\left|{\mathit{\lambda }}_{\mathrm{3}}\right|}.\end{array}$$

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,

$$\begin{array}{}\text{(16)}& \mathrm{\Delta }x=\mathrm{\Delta }{x}_{\mathrm{0}}{e}^{\mathit{\lambda }\mathrm{\Delta }t}.\end{array}$$

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, λ₁≥0, λ₃≤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., λ₃. FTLEs are most commonly computed as

$$\begin{array}{}\text{(17)}& {\mathit{\lambda }}_i=\mathrm{1}/\left|T\right|\ln \sqrt{{\mathit{\sigma }}_i},\end{array}$$

where σ_i values are the eigenvalues of the right Cauchy–Green deformation tensor,

$$\begin{array}{}\text{(18)}& G=[\mathrm{\Delta }{x}_i/\mathrm{\Delta }{x}_{\mathrm{0}}j{]}^T[\mathrm{\Delta }{x}_i/\mathrm{\Delta }{x}_{\mathrm{0}}j].\end{array}$$

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 λ₃.

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 λ₃, is denoted by σ. As shown in the Appendix, the differential equation for σ has a fixed point at

$$\begin{array}{}\text{(19)}& \mathit{\sigma }=\sqrt{\mathit{\kappa }/\left|{\mathit{\lambda }}_{\mathrm{3}}\right|},\end{array}$$

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 λ₃ direction will converge to the Lagrangian Batchelor scale. Mathematically there are also fixed points with negative λ₃ and with negative σ for positive λ₃, but neither corresponds to a real positive tracer distribution. The full solution for σ is

$$\begin{array}{}\text{(20)}& \mathit{\sigma }=\sqrt{\mathit{\kappa }/\left|{\mathit{\lambda }}_{\mathrm{3}}\right|}{\left(({\mathit{\lambda }}_{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}/\mathit{\kappa }-\mathrm{1})e^{\mathrm{2}{\mathit{\lambda }}_{\mathrm{3}}t}+\mathrm{1}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

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

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=\mathrm{1.5}+x$, v=0.5y, $w=-\mathrm{1.5}z$. Color shows z values at $t=\mathrm{0}.$.

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

$$\begin{array}{}\text{(21)}& \mathit{\kappa }=\mathrm{0.0103}{l}^{\mathrm{1.15}},\end{array}$$

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:

$$\begin{array}{}\text{(22)}& \mathit{\delta }=\mathrm{0.0046}|{\mathit{\lambda }}_{\mathrm{3}}{|}^{-\mathrm{1.1765}},\end{array}$$

where λ₃ (in s⁻¹) yields δ (in cm s⁻¹).

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⁻¹; 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, λ₃, 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² s⁻¹ at 1 m to 0.06 m² s⁻¹ at 140 m, which are considerably smaller than diffusivities on the horizontal scale of eddies themselves, about 0.5–8.2 m² s⁻¹ 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.

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² s⁻¹) are between 10⁻⁴ and $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{3}}$ for the three larger Ekman numbers and between $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ and $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{2}}$ for E=0.0005.

4 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

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\mathrm{d}{x}_i}{\mathrm{d}t}}=U_{\mathrm{b}i}\left(\mathit{x}\right)+u_i^{\prime },\end{array}$$

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² s⁻¹ 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₁, weight 1∕6), the halfway point estimated from the current position (v₂, weight 1∕3), the halfway point estimated using v₂ (v₃, weight 1∕3), and the final point estimated using v₃ (v₄, weight 1∕6). Only v₁ and v₄ include stochastic additions, leading to the 1∕3 factor. Together, these give

$$\begin{array}{}\text{(24)}& \mathit{\kappa }={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\Delta }t}{\mathrm{18}}},\end{array}$$

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.

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

Figure 8 Range in ψ for ensembles of trajectories started from a sphere of radius 0.002. Steady perturbation ($\mathit{ϵ}\in \mathit{\{}\mathrm{0.01},\mathrm{0.08}\mathit{\}}$), stochastic perturbations ($\mathit{\kappa }\in \mathit{\{}{\mathrm{10}}^{-\mathrm{5}},{\mathrm{10}}^{-\mathrm{6}},{\mathrm{10}}^{-\mathrm{7}}\mathit{\}}$), or both ($\mathit{\kappa }={\mathrm{10}}^{-\mathrm{7}}$, ϵ=0.01) are added to the background flow. (a, c) Initial sphere in the chaotic sea region, away from fixed points, at $(r,z)=(\mathrm{0.1},\mathrm{0.5})$. (b, d) Initial sphere centered on $(r,z)=(\mathrm{0.4},\mathrm{0.5})$, a resonant region. In (a), the dashed black line is ${\mathrm{10}}^{-\mathrm{5}}\sqrt{t}.$

5 Tracer simulations and Nakamura effective diffusivity

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

$$\begin{array}{}\text{(25)}& {\mathit{\kappa }}_{\mathrm{eff}}\left(C\right)=k{\displaystyle \frac{\mathrm{1}}{(\partial C/\partial V{)}^{\mathrm{2}}}}{\widehat{\left|\mathrm{\nabla }C\right|}}^{\mathrm{2}},\end{array}$$

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² s⁻¹), multiplied by m⁴, 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:

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\partial }{\partial V}}\left({\mathit{\kappa }}_{\mathrm{eff}}{\displaystyle \frac{\partial C}{\partial V}}\right).\end{array}$$

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, χ² (Pattanayak, 2001):

$$\begin{array}{}\text{(27)}& {\mathit{\chi }}^{\mathrm{2}}=\underset{V}{\int }|\mathrm{\nabla }C{|}^{\mathrm{2}}\mathrm{d}V/\underset{V}{\int }|C{|}^{\mathrm{2}}\mathrm{d}V,\end{array}$$

where V here is simple volume. Stirring increases the variance of a tracer, while mixing decreases it. When χ² is increasing, stirring is dominant and the slope of χ²(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₀ from the cylinder axis so that X₀ 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⁻⁴ or 10⁻⁶. 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$.

Figure 9 (a–c) Tracer variance, χ²; (d–f) κ_eff integrated over volume. (a, d) $k={\mathrm{10}}^{-\mathrm{6}}$ and E=0.02, (b, e) $k={\mathrm{10}}^{-\mathrm{4}}$ and E=0.02, and (c, f) $k={\mathrm{10}}^{-\mathrm{4}}$ and E=0.125. Solid blue lines include the deterministic perturbation which induces chaos, $X_{\mathrm{0}}=-\mathrm{0.02}$, green dashed lines are unperturbed, and solid red lines include the deterministic perturbation with $X_{\mathrm{0}}=-\mathrm{0.16}$. Black dashed lines indicate κ_eff integrated over volume in the case of nested circular tori.

The set of simulations performed allows for an examination of the effects of changing E, k, and X₀. 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₀ 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 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 χ² (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 χ². At longer times, the integrated effective diffusivity reaches a nearly constant positive value unlike χ², 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π⁶∕8, which we expect to be the minimum for ∫κ_effdV 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 ∫κ_effdV 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 ∫κ_effdV 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₀ 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) 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 (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.

Figure 10 Results from three Navier–Stokes simulations with E=0.02: (a, d, g, j) $x_{\mathrm{0}}=-\mathrm{0.02}$ and $k={\mathrm{10}}^{-\mathrm{6}}$, (b, e, h, k) x₀=0 and $k={\mathrm{10}}^{-\mathrm{4}}$, and (c, f, i, l) $x_{\mathrm{0}}=-\mathrm{0.02}$ and $k={\mathrm{10}}^{-\mathrm{4}}$. The x₀=0, $k={\mathrm{10}}^{-\mathrm{6}}$ case is not shown but is qualitatively similar to the x₀=0, $k={\mathrm{10}}^{-\mathrm{4}}$ case. (a–c) Dye, t=39. (d–f) κ_eff, t=39. (g–i) Dye, t=299. (j–l) κ_eff, t=299.

Figure 11 Results from Navier–Stokes simulations for E=0.125 and $k={\mathrm{10}}^{-\mathrm{4}}$, with three deterministic perturbation levels: (a, d, g, j) X₀=0, (b, e, h, k) $X_{\mathrm{0}}=-\mathrm{0.02}$, and (c, f, i, l) $X_{\mathrm{0}}=-\mathrm{0.16}$. (a–c) Dye, t=39. (d–f) κ_eff, t=39. (g–i) Dye, t=299. (j–l) κ_eff, t=299.

When the imposed diffusivity k is decreased by 2 orders of magnitude, with X₀ 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₀ 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₀ 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 χ²) and highly elevated enhancement (much larger κ_eff).

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.

6 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 (λ₃), 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, λ₃, 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 λ₃ is calculated. Finally, it is not yet possible to calculate λ₃ 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 λ₃ 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, χ². 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 χ² 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 χ², reach peak values. This is followed by a long period in which χ² 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.

Appendix A: Bifurcation analysis of fixed points of the 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 rz-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 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)=(\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 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.

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 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.

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 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:

$$\begin{array}{}\text{(B1)}& C=c_{\max }\left(t\right)\exp \left({\displaystyle \frac{-x^{\mathrm{2}}{\mathit{\alpha }}^{\mathrm{2}}\left(t\right)}{\mathrm{2}}}+{\displaystyle \frac{-y^{\mathrm{2}}{\mathit{\beta }}^{\mathrm{2}}\left(t\right)}{\mathrm{2}}}+{\displaystyle \frac{-z^{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{2}}\left(t\right)}{\mathrm{2}}}\right),\end{array}$$

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 λ₃ 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

$$\begin{array}{}\text{(B2)}& {\displaystyle }u& {\displaystyle }={\mathit{\lambda }}_{\mathrm{3}}x({\mathit{x}}_{\mathrm{0}},t),\text{(B3)}& {\displaystyle }v& {\displaystyle }={\mathit{\lambda }}_{\mathrm{2}}y({\mathit{x}}_{\mathrm{0}},t),\text{(B4)}& {\displaystyle }w& {\displaystyle }={\mathit{\lambda }}_{\mathrm{1}}z({\mathit{x}}_{\mathrm{0}},t),\text{(B5)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{1}}& {\displaystyle }>{\mathit{\lambda }}_{\mathrm{2}}>{\mathit{\lambda }}_{\mathrm{3}},\text{(B6)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{1}}& {\displaystyle }>\mathrm{0},{\mathit{\lambda }}_{\mathrm{3}}<\mathrm{0},\end{array}$$

with x(x₀,t) indicating the initial position x₀ of the water parcel at t=0. The Lagrangian tracer evolution equation is

$$\begin{array}{}\text{(B7)}& {\displaystyle \frac{\partial C}{\partial t}}+{\mathit{\lambda }}_{\mathrm{3}}x{\displaystyle \frac{\partial C}{\partial x}}+{\mathit{\lambda }}_{\mathrm{2}}y{\displaystyle \frac{\partial C}{\partial y}}+{\mathit{\lambda }}_{\mathrm{1}}z{\displaystyle \frac{\partial C}{\partial z}}=\mathit{\kappa }{\mathrm{\nabla }}^{\mathrm{2}}C,\end{array}$$

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

$$\begin{array}{}\text{(B8)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{c_{\max }}}{\displaystyle \frac{\mathrm{d}{c}_{\max }}{\mathrm{d}t}}=-\mathit{\kappa }\left({\mathit{\alpha }}^{\mathrm{2}}+{\mathit{\beta }}^{\mathrm{2}}+{\mathit{\gamma }}^{\mathrm{2}}\right),\text{(B9)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\alpha }}{\mathrm{d}t}}=-{\mathit{\lambda }}_{\mathrm{3}}\mathit{\alpha }-\mathit{\kappa }{\mathit{\alpha }}^{\mathrm{3}},\text{(B10)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\beta }}{\mathrm{d}t}}=-{\mathit{\lambda }}_{\mathrm{2}}\mathit{\beta }-\mathit{\kappa }{\mathit{\beta }}^{\mathrm{3}},\text{(B11)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\gamma }}{\mathrm{d}t}}=-{\mathit{\lambda }}_{\mathrm{1}}\mathit{\gamma }-\mathit{\kappa }{\mathit{\gamma }}^{\mathrm{3}}.\end{array}$$

The width parameters' equations are nonlinear, but when rewritten in terms like α⁻², they give the following:

$$\begin{array}{}\text{(B12)}& {\displaystyle \frac{\mathrm{d}{\mathit{\alpha }}^{-\mathrm{2}}}{\mathrm{d}t}}& {\displaystyle }=\mathrm{2}{\mathit{\lambda }}_{\mathrm{3}}{\mathit{\alpha }}^{-\mathrm{2}}+\mathrm{2}\mathit{\kappa },\text{(B13)}& {\displaystyle \frac{\mathrm{d}{\mathit{\beta }}^{-\mathrm{2}}}{\mathrm{d}t}}& {\displaystyle }=\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\beta }}^{-\mathrm{2}}+\mathrm{2}\mathit{\kappa },\text{(B14)}& {\displaystyle \frac{\mathrm{d}{\mathit{\gamma }}^{-\mathrm{2}}}{\mathrm{d}t}}& {\displaystyle }=\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\gamma }}^{-\mathrm{2}}+\mathrm{2}\mathit{\kappa },\end{array}$$

which are Bernoulli equations, solvable with integrating factors, giving

$$\begin{array}{}\text{(B15)}& {\displaystyle }\mathit{\alpha }& {\displaystyle }=\sqrt{\left|{\mathit{\lambda }}_{\mathrm{3}}\right|/\mathit{\kappa }}{\left(({\mathit{\lambda }}_{\mathrm{3}}{\mathit{\alpha }}_{\mathrm{0}}^{-\mathrm{2}}/\mathit{\kappa }-\mathrm{1})e^{\mathrm{2}{\mathit{\lambda }}_{\mathrm{3}}t}+\mathrm{1}\right)}^{-\mathrm{1}/\mathrm{2}},\text{(B16)}& {\displaystyle }\mathit{\beta }& {\displaystyle }={\left(({\mathit{\beta }}_{\mathrm{0}}^{-\mathrm{2}}+\mathit{\kappa }/{\mathit{\lambda }}_{\mathrm{2}})e^{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}t}-\mathit{\kappa }/{\mathit{\lambda }}_{\mathrm{2}}\right)}^{-\mathrm{1}/\mathrm{2}},\text{(B17)}& {\displaystyle }\mathit{\gamma }& {\displaystyle }=\sqrt{{\mathit{\lambda }}_{\mathrm{1}}/\mathit{\kappa }}{\left(({\mathit{\lambda }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{0}}^{-\mathrm{2}}/\mathit{\kappa }+\mathrm{1})e^{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}t}-\mathrm{1}\right)}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

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 λ₂ 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

$$\begin{array}{}\text{(B18)}& c_{\max }\left(t\right)=c_{\mathrm{0}}\mathit{\alpha }\left(t\right)\mathit{\beta }\left(t\right)\mathit{\gamma }\left(t\right).\end{array}$$

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

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{c}_{\max }}{\mathrm{d}t}}& {\displaystyle }={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(c_{\mathrm{0}}\mathit{\alpha }\mathit{\beta }\mathit{\gamma }\right)=c_{\mathrm{0}}\left(\mathit{\beta }\mathit{\gamma }{\displaystyle \frac{\mathrm{d}\mathit{\alpha }}{\mathrm{d}t}}+\mathit{\alpha }\mathit{\gamma }{\displaystyle \frac{\mathrm{d}\mathit{\beta }}{\mathrm{d}t}}+\mathit{\alpha }\mathit{\beta }{\displaystyle \frac{\mathrm{d}\mathit{\gamma }}{\mathrm{d}t}}\right)\\ {\displaystyle }& {\displaystyle }=c_{\mathrm{0}}\left(-\mathit{\alpha }\mathit{\beta }\mathit{\gamma }({\mathit{\lambda }}_{\mathrm{3}}+\mathit{\kappa }{\mathit{\alpha }}^{\mathrm{2}})\right.\\ {\displaystyle }& {\displaystyle }\left.-\mathit{\alpha }\mathit{\beta }\mathit{\gamma }({\mathit{\lambda }}_{\mathrm{2}}+\mathit{\kappa }{\mathit{\beta }}^{\mathrm{2}})-\mathit{\alpha }\mathit{\beta }\mathit{\gamma }({\mathit{\lambda }}_{\mathrm{1}}+\mathit{\kappa }{\mathit{\gamma }}^{\mathrm{2}})\right)\\ {\displaystyle }& {\displaystyle }=-c_{\mathrm{0}}\mathit{\alpha }\mathit{\beta }\mathit{\gamma }\left({\mathit{\lambda }}_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}+\mathit{\kappa }[{\mathit{\alpha }}^{\mathrm{2}}+{\mathit{\beta }}^{\mathrm{2}}+{\mathit{\gamma }}^{\mathrm{2}}]\right)\\ {\displaystyle }& {\displaystyle }\Rightarrow {\displaystyle \frac{\mathrm{1}}{c_{\max }}}{\displaystyle \frac{\mathrm{d}{c}_{\max }}{\mathrm{d}t}}=-\mathit{\kappa }\left({\mathit{\alpha }}^{\mathrm{2}}+{\mathit{\beta }}^{\mathrm{2}}+{\mathit{\gamma }}^{\mathrm{2}}\right).\end{array}$$

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.

Appendix C: Long-time limit of effective diffusivity for the axially symmetric rotating cylinder flow

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 ∫κ_effdV, 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=kA², the background diffusivity multiplied by the squared surface area of a torus.

Recall that the volume of a circular torus is

$$\begin{array}{}\text{(C1)}& V_{\mathrm{ct}}=\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}{r}^{\mathrm{2}}R,\end{array}$$

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

$$A_{\mathrm{ct}}=\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}rR.$$

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

$$\begin{array}{ll}{\displaystyle }\int {\mathit{\kappa }}_{\mathrm{eff}}\mathrm{d}V& {\displaystyle }=\iiint k{A}^{\mathrm{2}}\mathrm{d}V\\ {\displaystyle }& {\displaystyle }=k\underset{\mathrm{0}}{\overset{r_{\max }}{\int }}{A}^{\mathrm{3}}\mathrm{d}r\\ {\displaystyle }& {\displaystyle }=k\underset{\mathrm{0}}{\overset{r_{\max }}{\int }}(\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}rR{)}^{\mathrm{3}}\mathrm{d}r\\ {\displaystyle }& {\displaystyle }={\mathrm{4}}^{\mathrm{3}}{\mathit{\pi }}^{\mathrm{6}}{R}^{\mathrm{3}}k\underset{\mathrm{0}}{\overset{r_{\max }}{\int }}{r}^{\mathrm{3}}\mathrm{d}r\\ \text{(C2)}& {\displaystyle }& {\displaystyle }={\mathrm{4}}^{\mathrm{2}}{\mathit{\pi }}^{\mathrm{6}}{R}^{\mathrm{3}}k{r}^{\mathrm{4}}{|}_{\mathrm{0}}^{r_{\max }}=k{\mathit{\pi }}^{\mathrm{6}}/\mathrm{8},\end{array}$$

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.