Trajectory encounter volume as a diagnostic of mixing potential in fluid flows

Fluid parcels can exchange water properties when coming into contact with each other, leading to mixing. The trajectory encounter mass and a related simplified quantity, the encounter volume, are introduced as a measure of the mixing potential of a flow. The encounter volume quantifies the volume of fluid that passes close to a reference trajectory over a finite time interval. Regions characterized by a low encounter volume, such as the cores of coherent eddies, have a low mixing potential, whereas turbulent or chaotic regions characterized by a large encounter volume have a high mixing potential. The encounter volume diagnostic is used to characterize the mixing potential in three flows of increasing complexity: the Duffing oscillator, the Bickley jet and the altimetry-based velocity in the Gulf Stream extension region. An additional example is presented in which the encounter volume is combined with the u approach of Pratt et al. (2016) to characterize the mixing potential for a specific tracer distribution in the Bickley jet flow. Analytical relationships are derived that connect the encounter volume to the shear and strain rates for linear shear and linear strain flows, respectively. It is shown that in both flows the encounter volume is proportional to time.

1 Encounter volume

Main idea
Mixing is an irreversible exchange of properties between different water masses.This process is important for maintaining the oceanic large-scale stratification and general circulation, and it plays a key role in the redistribution of biogeo-chemical tracers throughout the world oceans.Mixing occurs between different water masses when they come into direct contact with each other.Thus, the mixing potential of the flow, i.e., the opportunity for mixing to occur, is generally enhanced in regions where water parcels meet or encounter many other water parcels and are thus exposed to a large amount of fluid passing by them as the flow evolves.This would be the case, for example, for a parcel within a chaotic zone, which is a region of the flow that is in a state of chaotic advection.There, the separation between initially nearby water parcels grows exponentially in time and, in the infinite time limit, each water parcel encounters all the other water parcels within the same zone and comes into contact with the entire volume of the chaotic zone.Similarly, high encounter volumes will exist in turbulent regions.In contrast, the mixing potential and encounter volume is expected to be smaller in regions where water parcels do not experience many encounters with other water parcels and remain close to their initial neighbors as the flow evolves.This would be the case, for example, for a water parcel that is located inside a coherent eddy.If the eddy is in a state of solid body rotation, the water parcel would forever stay close to its initial neighbors and will not have any new encounters at all.If some amount of azimuthal shear is present within the eddy, then for a water parcel located at a radius r from the eddy center, the encounters will be limited to those water parcels located within a circular strip centered at the same r.
Of course, the presence of a mixing potential does not guarantee that the mixing of a tracer will occur; it is also essential that the tracer distribution is nonuniform so that irreversible property exchange can take place between the different water parcels during their encounters.This exchange hap-pens through diffusion and therefore relies on a concentration difference between the two parcels.Thus, the intensity of the mixing would depend on both the tracer distribution and the flow, whereas the mixing potential is the property of only the flow field alone.In this work, we introduce the concept of an encounter mass, M, and the encounter volume, V , which serves as a simplified representation of M in incompressible flows, as an objective measure of the encounters between different fluid elements in order to quantify the mixing potential of a fluid flow.There are many existing trajectory-based measures of fluid stirring; ours has the virtue of having a straightforward physical interpretation that is easy to implement and immediately applicable to ocean float and drifter data.Our method does not require sophisticated bookkeeping as in braid theory (Allshouse and Thiffeault, 2012) or finite-time entropy (Froyland and Padberg-Gehle, 2012).

Definition and numerical implementation
For a given reference trajectory, x(x 0 , t 0 ; T ), the encounter mass, M(x 0 , t 0 ; T ), is defined as the total mass of fluid that passes within a radius R of a reference trajectory over a finite time interval t 0 < t < t 0 +T .One might imagine a sphere that has a radius R and that is centered at and moves with the reference trajectory.The encounter mass then consists of the mass of the fluid that is initially located within the sphere along with the mass of all the fluid that passes through the sphere over the time interval t 0 < t < t 0 + T .Note that it is generally not possible to compute the latter by simply integrating the mass flux into the sphere over t 0 < t < t 0 + T since some fluid may leave and then reenter the sphere and would be counted more than once; Lagrangian information is therefore required to keep track of the history of each fluid parcel trajectory entering the sphere.
To this end, we subdivide the entire fluid at t = t 0 into small compact fluid elements with masses δM i = ρ i δV i , where ρ i is the density of a fluid element and δV i is its volume.We wish to follow the motion of each fluid element over time interval t 0 < t < t 0 + T , and we assume that the elements remain compact over such time so that the motion of each fluid element can be well represented by one trajectory.If the fluid elements stretch and deform too much, we can evoke the continuum hypothesis and make δM sufficiently small that such compactness is assured.In the limit of infinitesimal fluid elements, δM i → dM i , we can associate a unique trajectory with each infinitesimal fluid element.The encounter mass is then For an incompressible flow, the density and volume of each fluid element, ρ i and δV i , remain constant following a trajectory, although different fluid elements are still allowed to have different densities such as, for example, in stratified three-dimensional geophysical flows.If the flow is unstratified, the densities of all fluid elements are equal, ρ i = ρ, and the encounter mass becomes is the encounter volume defined as the total volume of fluid that passes within a radius R of the reference trajectory over a finite time interval t 0 < t < t 0 +T .When all volume elements are equal, dV i = dV , the encounter volume can be further simplified to where the encounter number, N (x 0 , t 0 ; T ), is the number of trajectories that come within a radius R of the reference trajectory over a time interval t 0 < t < t 0 + T .We will refer to t 0 as the starting time, T as the trajectory integration time and x 0 as the trajectory initial position; i.e., x (x 0 , t 0 ; T = 0) = x 0 .For practical applications with geophysical flows, the limit in the definition of the encounter volume can be dropped, and one can simply use the approximation with the dense grid of initial positions x 0 .Mathematically, the encounter number can be written as where the indicator function I is 1 if true and 0 if false, and K is the total number of Lagrangian particles released.The encounter volume depends on the starting time, the integration time, the encounter radius and the number of trajectories (i.e., the grid spacing); all of these parameter dependences will be discussed below.Once the encounter volume is estimated, regions of space with large or small V would then be associated with enhanced or inhibited mixing potential.For the remainder of this paper, we will focus on incompressible fluid flows and will be concerned with the encounter volume rather than the encounter mass.We define V (x 0 , t 0 ; T ) and N (x 0 , t 0 ; T ) in Eq. (1a, b) based on the number of encounters with different trajectories, not the total number of encounter events; even if a trajectory first comes close to the reference trajectory and then moves away and reapproaches it again later, it is only counted once.In a flow field with no sources or sinks of tracer variance, where the variance is therefore decaying, it is reasonable to expect that the most property exchange between two parcels will often occur during their first encounter; this is thus the motive for counting only the first encounter.Note that this assumption may not hold if the parcels reacquire different properties after their first encounter due to encountering and exchanging properties with other parcels.In this case, or in the case when tracer variance is continuously introduced, it may be more reasonable to count the total number of encounters.
For a numerical implementation of the trajectory encounter volume-based mixing characterization, one would need to start at a chosen time t 0 with a grid of initial positions spanning the flow domain, and then evolve the trajectories under the flow field over the time interval T .This time interval should be chosen based on the physical properties of the flow and with specific scientific questions in mind.For example, if the research focus is on ocean submesoscale dynamics, the timescale T would be on the order of hours to days, whereas the corresponding timescale for mesoscale dynamics would be on the order of weeks to months.
V (x 0 , t 0 ; T ) is a Lagrangian quantity that characterizes the mixing potential of a flow over a time interval from t 0 to t 0 + T .As the flow field evolves in time, its mixing characteristics can change, and V (x 0 , t 0 ; T ) will reflect this change.For example, if a coherent eddy with weak mixing potential embedded in a chaotic zone with enhanced mixing potential was present in the flow from time t 1 to time t 2 but dispersed and disappeared afterwards, then V (x 0 , t 0 ; T ) is expected to be small at those locations x 0 that correspond to the interior of an eddy for t 0 ≥ t 1 and t 0 + T ≤ t 2 .For t 0 > t 2 , when the eddy is no longer present, V (x 0 , t 0 ; T ) would increase.Dependences on T and t 0 are similarly expected to be present within a chaotic zone.
In the infinite time limit, T → ∞, when all parcels within a chaotic zone (or turbulent region) of a finite extent encounter all other parcels within the same chaotic zone, the encounter volume V (x 0 , t 0 ; T → ∞) approaches a constant equal to the volume (or the area in two dimensions) of the chaotic zone.For a two-dimensional, incompressible flow, the encounter rates over finite T are locally the largest near a hyperbolic trajectory and along the segments of its associated stable manifolds.The stable manifolds serve as pathways that bring water parcels from remote regions into the vicinity of the hyperbolic trajectory, where parcels stay for extended periods of time and where many encounters occur.Note that the unstable manifolds, on the other hand, will rapidly remove a particle from a hyperbolic region, thus limiting its exposure to the high-encounter region near the hyperbolic trajectory.For this reason, the unstable manifolds are not revealed by the encounter volume calculation performed in forward time and require a backward-time calculation instead.This exclusive link between the forward-and backward-time calculations of the trajectories and the stable and unstable manifolds, respectively, is not specific to the encounter volume diagnostic.It is rather typical for many finitetime methods from the dynamical systems theory, including finite-time Lyapunov exponents (FTLEs), which in forward time approximate segments of the stable manifold as max-imizing ridges (Haller, 2002;Shadden et al., 2005;Lekien and Ross, 2010).
Since the locations of the hyperbolic trajectories and manifolds generally evolve in time, V (x 0 , t 0 ; T ) is also expected to vary with t 0 .As the trajectory integration time T increases, the water parcels initially located further from the hyperbolic trajectory will have the opportunity to come into its vicinity along the stable manifold.Such parcels, as they approach the hyperbolic trajectory, are expected to have more encounters than their neighbors that are initially located off the manifold and thus bypass the vicinity of the hyperbolic trajectory where many encounters occur.Thus, V (x 0 , t 0 ; T ) reveals longer segments of stable manifolds for longer integration time T , as will be illustrated numerically in the next section.In the long integration time limit when each manifold, either stable or unstable, densely fills the entire chaotic zone forming a dense homoclininc or heteroclinic tangle, the whole tangle will be characterized by high encounter volumes in both forward and backward time.Again, this is similar to how the maximizing ridges of the forward-time FTLEs elongate and sharpen with increasing integration time.
The radius R, which defines how close to a reference trajectory another trajectory should come in order to be counted as an encounter, is an important parameter for the calculation of the encounter volume V .Generally, R should be small compared to the spatial scale of the smallest features of interest.Specifically, for the V field to delineate a flow feature like an eddy, the trajectories within the eddy interior should not encounter those on its exterior.The boundary region near the eddy perimeter where such encounters can occur has the width 2R.So, if that width is comparable to or larger than the eddy size, then the eddy would become completely smeared out and not be resolved.From a practical viewpoint, however, using a very small R would require very dense grids of trajectories to be computed.Otherwise, zero or a very small number of trajectory encounters will occur in the entire flow domain.The numerical examples in the next section suggest that choosing R to be a fraction, up to about half of the radius of the smallest features of interest, works best.
Finally, the approximation V ≈ N δV breaks down for sparse grids of initial positions with the insufficient number of Lagrangian particles when N is small and δV is large.It also works poorly when applied to two-dimensional divergent flows due to δV changing following trajectories.The numerical simulations in the next section suggest that a grid spacing of ≤ R/2 is sufficient and that the method can also be applied to characterize the mixing potential in slightly divergent two-dimensional flows.
Once the timescale T is identified, the grid of initial positions is chosen, the trajectories are computed, the radius R is defined and the number of encounters, N (x 0 , t 0 ; t), is counted for each trajectory, then the encounter volume can be estimated as V ≈ N δV and plotted as a function of the trajectory initial position x 0 .The resulting V field delineates the flow regions with different mixing properties as subdomains having different values of V .

Examples
We proceed to test the performance of the encounter volume technique in quantifying the mixing potential of several geophysically relevant sample flows of increasing complexity, starting from a simple analytically prescribed periodically perturbed double-gyre Duffing oscillator system, followed by a dynamically consistent solution of the PV conservation equation on a beta plane known as the Bickley jet and finishing with an observationally based geostrophic velocity field in the North Atlantic derived from the sea surface height altimetry.

Duffing oscillator
The Duffing oscillator flow with its figure-eight geometry has become a standard test case for emerging techniques related to the dynamical systems theory.This flow consists of two gyres with the same sign of rotation (clockwise in our case) with elliptic centers that oscillate in time around their mean positions.A hyperbolic point is located at the origin between the two gyres, and a pair of stable and unstable manifolds emanate from it, forming a figure eight in the absence of the time-dependent perturbation or forming a classic homoclinic tangle in the presence of the perturbation.The velocity field is two-dimensional, incompressible and given by u = y and v = (x − ax 3 )(1 + cos(ωt + φ)) with a = 1, ω = 3π/2, φ = π/4 and = 0.1.With these parameters, the Poincaré section (Fig. 1b) shows the presence of two main regular elliptical regions with a O(1) radius corresponding to the interiors of the gyres, which are embedded in a chaotic zone shaped like a figure eight, within which a number of island chains with smaller regular islands are present.The winding time for most of the trajectories in the system is on the order of 5T pert with T pert = 2π ω , except for the trajectories near the hyperbolic point, for which the winding time is much longer (Fig. 1a).
The encounter volume was computed for a range of trajectory integration times from T = T pert (which is significantly shorter than the trajectory winding time) to T = 50T pert (significantly longer than the trajectory winding time) and for a range of encounter radii from R = 0.01 R eddy (significantly smaller than the eddy core radius) to R = 1 ≈ R eddy (comparable to the eddy core radius).The results in Fig. 2 suggest that the encounter volume method works best for integration times longer than the trajectory winding time and with an encounter radius about one-third to one-half of the gyre radius (the right three panels of the middle row).For a very small encounter radius (top row), V is noisy because the trajectories simply do not encounter many neighbors.Thus, delineating the domain into regions with different mixing po- tentials, as in the top right panel, requires a long integration time.For T = 50T pert , good agreement with the Poincaré section is observed, and the use of a small encounter radius allows for a precise identification of smaller regular island chains, such as the chains of four islands located just outside of the perimeter of both the left and right eddy cores.Note that the noise in the V field can be suppressed by using a denser initial grid of trajectories, but at the cost of a more expensive computation.For very short integration times (left column) when the trajectory segments are very short, the encounter volume does not capture the difference between the regular and chaotic regions.This is not surprising as velocity shear is probably a dominating factor over such small times.As the integration time increases, the difference in the encounter volume becomes more pronounced between the trajectories that remain within the eddy cores and the trajectories that are free to move around the chaotic zone.Over a timescale of approximately one winding period (or about five periods of the perturbation; second column), the two regular eddy cores (the blue regions with the small V ) and a segment of the stable manifold (the red curve emanating from the origin with the largest V ) become clearly visible for R = 0.2 and R = 1.The revealed manifold segment becomes longer, narrower and more tangled, eventually filling up the whole chaotic zone.At the same time, the shape of the core region becomes more exact and approaches the "true" core in the Poincaré section as the integration time increases to 50 periods of the perturbation.The agreement with the Poincaré section is excellent in the right middle panel, although the smaller island chains are not as visible as in the top right panel because of the use of a larger encounter radius that is comparable to their size (see Fig. 3).Finally, for the large encounter radius that is comparable to the size of the eddy (bottom row), the boundary region near the perimeter of an eddy, within which the trajectories on the inside of the eddy can encounter the trajectories passing by on the outside, is as wide as the eddy itself, essentially wiping out all small scales from the V field.All of these trends are in agreement with the theoretical expectations described in Sect. 1.
In order to more clearly highlight the link between high values of V and stable (rather than unstable) manifolds, we have computed both the stable and unstable manifolds for the Duffing oscillator flow using a direct method, where we grew manifolds from a small segment starting at the hyperbolic trajectory.For the Duffing oscillator, this computation is straightforward since the the hyperbolic trajectory stays at the origin at all times.Both the stable and unstable directly computed manifolds were then superimposed on a forwardtime encounter volume plot in Fig. 4. The comparison shows that, as anticipated, the encounter volume diagnostic clearly highlights the stable manifolds as maximizing ridges of V computed in forward time.
With a variety of dynamical systems techniques available, it is important to understand the advantages and limitations of the different methods.We compared the encounter vol- ume to two well-established and commonly used methods, the Poincaré section (Fig. 3) and the FTLEs (Fig. 5).Since the Poincaré section requires stroboscopic sampling of the trajectories in time, it can only be applied to time-periodic flows and requires that the trajectories are computed over a long integration time, typically thousands of periods of the perturbation.On the other hand, it generally requires only a few parcels to be released at some key locations rather than releasing a dense grid of initial positions to map out the entire phase space.The encounter volume and FTLEs, on the other hand, are not limited to time-periodic flows and also work with significantly shorter segments of trajectories (the longest integration time in our simulations in Fig. 2 is only 50 periods of the perturbation).They are also better suited for identifying the manifolds than the Poincaré sectioning as they do not require any a priori knowledge about the location of the hyperbolic trajectory.On the other hand, they require many more parcels to be released in order to map out the phase space.When applied to the same set of trajectories (the same initial positions and integration times), the FTLEs and the encounter volume methods produced similar results (Fig. 5), with V being arguably better suited for (1) identifying the coherent core regions of the eddies where the FTLEs have spiraling patterns that complicate the analysis and (2) producing more continuous segments of manifolds at intermediate integration times when the FTLE-based ridges become discontinuous near the turning points of a manifold.The advantage of FTLEs, on the other hand, is that they have fewer parameters (T and grid spacing), whereas V also depends on R, and they are less expensive computationally.The more expensive computational cost of V compared to FTLEs is due to two reasons: first, the FTLEs only depend on the initial and final positions of the trajectories, whereas V depends on the entire trajectory history; second, the FTLEs depend on the relative distance between a trajectory and its closest neighbors, whereas V keeps tracks of the encounters with all trajectories, not just the neighboring trajectories.Thus, the cost of evaluating the FTLE for each particle is indepen-dent of the total number of particles released, and the cost of evaluating V for each particle increases in proportion to the number of particles (since one needs to keep track of the encounters with all particles).The calculation of V is still feasible for realistic geophysical flows, as illustrated below.Also note that, depending on the physical question being studied, the information about the entire trajectory and not only about the final and initial position might, in fact, be advantageous.
Related to the issue of computational cost is the question of a sufficient grid size.We have carried out numerical simulations (Fig. 6) to investigate the dependence of the encounter volume on the grid size and to come up with a rule-of-thumb recommendation regarding the appropriate grid spacing.Our simulations suggest that the encounter volume values (approximated by V ≈ NdV ) are relatively insensitive to the variations in the grid spacing between one-tenth and one-half of the encounter radius (with the encounter radius being a fraction of the size of the feature of interest, as suggested by Fig. 2) and that the major effect of a coarser grid is the degraded resolution of the resulting V map rather than incorrect V values.

Bickley jet
The meandering Bickley jet flow is an idealized, but linearly dynamically consistent, model for the eastward zonal jet in the Earth's stratosphere (del-Castillo-Negrete and Morrison, 1993;Rypina et al., 2007Rypina et al., , 2011)).This flow consists of a steady eastward zonal jet on which two eastward propagating Rossby-like waves are superimposed.All flow parameters used here are identical to those used in our previous 2007 and 2011 papers.In the reference frame moving at a speed of one of the waves, the flow consists of a steady background velocity subject to a time-periodic perturbation.The background looks like a meandering jet, with three recirculation gyres to the north and south of the jet core.Between the recirculation gyres, there are three hyperbolic points with the associated stable and unstable manifolds.Under the influence of the time-periodic perturbation imposed by the second wave, heteroclinic tangles are formed by the manifolds emanating from different hyperbolic regions between the recirculations, and a chaotic zone emerges on either side of the jet.The manifolds, however, cannot penetrate through the jet core, which remains regular and acts as a transport barrier separating the northern and southern chaotic zones.All of these features are clearly visible in the Poincaré section shown in Fig. 7 (top).The bottom subplot shows the V field computed using the encounter radius R = 5 × 10 5 , which is about half of the recirculation region radius, and using a trajectory integration time on the order of a few winding times within the recirculations.As expected, the encounter volume identified six recirculation regions and the jet core as zones with a small mixing potential (blue).Six blue recirculation regions are embedded in two distinct chaotic zones with an enhanced mixing potential (yellow and red) on either side  of the jet.The mixing potential is the largest (red) along the segments of stable manifolds emanating from the hyperbolic trajectories between the recirculations.

Altimetry-based velocity in the meandering Gulf Stream region
Past its separation point from the coast at Cape Hatteras, the strong and narrow Gulf Stream current turns offshore, where it loses its coherence, broadens and weakens and then starts to meander.Some of the meanders then grow and eventually detach from the current, forming strong mesoscale eddies known as the Gulf Stream rings.On 11 July 1997, a number of such Gulf Stream rings of various strengths and sizes at different stages of their lifetimes were clearly present both north and south of the Gulf Stream extension current (Fig. 8).
The flow in the Gulf Stream extension region, with a nonsteady meandering jet, the Gulf Stream rings and the recirculations to the north and south of the jet core, has a lot in common, at least qualitatively, with the Bickley jet example.Unlike the idealized model, however, the real Gulf Stream rings have finite lifetimes, and the jet is not periodic in the zonal direction.Nevertheless, many of the qualitative features of the Bickley jet's V field hold in this example.Specifically, the trajectories inside the coherent eddy cores have smaller encounter volumes than the eddy peripheries, and the jet centerline has a smaller encounter volume than the flanks.
The velocity field that we used was downloaded from the AVISO website (http://www.aviso.altimetry.fr/en/data/products/sea-surface-height-products/global.html) and corresponds to their gridded product with 0.25 • spatial resolution and a temporal step of 1 day.This velocity is based on the altimetric sea surface height measurements made from satellites.The heights were converted into velocities using a geostrophic approximation.For the encounter volume estimation, the trajectories were seeded on a regular grid with dx = dy ∼ = 0.06 • on 11 July 1997 and were integrated forward in time for 90 days using a fifth-order, variable-step Runge-Kutta integration scheme with a bilinear interpolation between grid points in space and time.The encounter radius was chosen to be 0.3 • , which is about one-third of the radius of a typical 200 m wide Gulf Stream ring.for the Bickley jet flow.For the V calculation, the trajectories were released on a regular grid spanning the entire domain with a grid spacing of about 10 5 in both the x and y directions.
The encounter volume was estimated for three different integration times, T = 30, 60 and 90 days (Fig. 8).The V field clearly indicates that a number of Gulf Stream rings were present on both sides of the meandering jet.Among those, the two strongest ones can be seen at 54 • W, 36 • N and 52 • W, 41 • N with a low-V (blue) core and a high-V (red) periphery.As the integration time increases from 30 to 90 days, the Gulf Stream rings generally start to leak fluid.Their cores start to lose coherence and the encounter volume within the eddy cores starts to increase as more and more trajectories escape into the eddy surroundings over time.After a 90 day integration time, only a few Gulf Stream rings still possess coherent cores, whereas others become leaky throughout.Even for the two strongest rings, the coherent Lagrangian cores (the bluish regions with V ≈ 0) shrink down in size and, importantly, become significantly smaller than what the Eulerian velocity field would suggest.The core of the northern eddy also becomes slightly shifted to the east from the corresponding Eulerian stagnation point and becomes deformed into a non-convex, sickle-like shape.
The overall leakiness of the Gulf Stream rings and the small extent of their coherent Lagrangian core regions suggests that the coherent transport by the Gulf Stream rings (and maybe by mesoscale eddies in general) over time intervals of a few months or longer may be significantly smaller than what is generally anticipated from Eulerian diagnostics based on closed streamlines or Okubo-Weiss criteria.Inter-estingly, the prominent red rings (the large V values) around the eddy cores in Fig. 8 indicate that the significant contribution to transport by Lagrangian eddies may be due to the peripheries with a high mixing potential rather than the coherent cores themselves.
To visualize the Lagrangian evolution of the core regions and to illustrate the eddy leakiness, we extracted the trajectories from the core of the northern eddy in Fig. 8a (i.e., the trajectories with V < 6000 km 2 from the 30-day V field) and plotted their subsequent positions after 30, 60 and 90 days.The results in Fig. 9 confirm that the eddy core stays completely coherent over 30 days (i.e., all trajectories stay together), but starts to deteriorate at 60 days with only a small fraction of the initial patch still staying together and the rest of the patch dispersing and forming long and narrow filaments.
The jet region, although noisy, seems to suggest higher V near the flanks and smaller V near the centerline.The center region is not as well defined as in the Bickley jet example, possibly because the Gulf Stream inhibits but does not fully prevent the meridional transport in this region, and because our encounter radius might have been too large to reveal the central region if the true center region was narrower than 2R (0.6 • ).Finally, the V field suggests that the mixing potential of the flow is not symmetric with respect to the jet centerline and is higher on the northern side.It would be interesting to see if this is a general property of the flow in this region or if this phenomenon is specific to the time interval chosen.This investigation is left for a future study.

Encounter volume for some simple flow regimes
By analogy with molecular diffusion, the eddy diffusivity, K, is often used to characterize the eddy-induced downgradient tracer transfer in realistic geophysical fluid flows (LaCasce, 2008;Vallis, 2006;Rypina et al., 2012;Kamenkovich et al., 2015).Because of the simplicity of this approach, the majority of existing non-eddy-resolving oceanic numerical models are diffusion based, despite the somewhat questionable assumptions underlying this approach.An analytical connection between the encounter volume and the diffusivity would thus be useful for the parameterizations in numerical models.
Although we have not been able to find an analytical expression connecting V and K, we outline below some steps in that direction that help frame the problem.Let us start by considering a simple, diffusive, random walk particle motion in two dimensions, where particles take steps of a fixed length L in random directions at time intervals t.For such a process, the single particle dispersion is which characterizes the mean square displacement from the particle's initial position (x 0 , y 0 ) and grows in proportion to  the number of steps, n: with the proportionality coefficient, K = L 2 / t, denoting the diffusivity.The angular brackets denote the ensemble average.We are interested in finding an analytical expression for the encounter number, i.e., the number of particles that pass within radius R from a reference particle over time T as a function of K and T .
It is convenient to move to a reference frame that is tied to a reference particle, which would then always stay at the origin, while the other particles would be involved in a random walk motion.The problem of finding the encounter number is then reduced to counting the number of particles that come within the radius R from the origin over time T in the moving frame.The properties of the random walk process in the moving reference frame are different from those in the stationary frame.Specifically, the direction of each step in the moving reference frame still remains random (since it is a sum of two random variables, each uniformly distributed within an interval [0; 2π ]), but the step size is no longer fixed.Instead, the step size can be written as This change in the step size between the stationary and moving frames leads to a doubling of the diffusivity in the moving reference frame.To show this, we write the dispersion in the moving frame as where x = x − x 0 is the deviation from the initial position in the stationary frame (similarly for y, x ref and y ref ).We have used < x >=< y >= 0 to obtain the last equality.When averaged over many reference trajectories, < x 2 ref + y 2 ref >= D. This is because in the stationary reference frame, the reference particle is doing a random walk just like all the other particles, so that < D m >= 2D; equivalently, < K m >= 2K.
We thus seek an expression for the number of particles involved in a random walk process with a diffusivity of 2K that come within an encounter radius R from the origin during their first n steps (n plays the role of discretized integration time).This quantity is related to the first-passage time density, which characterizes the probability that a particle has first reached an absorbing boundary (often referred to as a cliff in statistics) at time t, and its integral quantity, called the survival probability, which characterizes the probability that a particle has not come into contact with an absorbing boundary over time t (i.e., it survived after time t without falling off a cliff).So far, however, we have not been able to complete the derivation, and we leave this development for a future investigation.
The numerical Monte Carlo simulations of a random walk process suggest that the dependence of the encounter number (and the encounter volume) on the integration time T is neither a linear nor a square root function.The power-law least square fit of the form V ∼ T α returns α values between 0.64 and 0.78 for a wide variety of R and K, each spanning an interval of values equaling an order of magnitude.Similarly, the power-law least square fits V ∼ K β and V ∼ R γ yield β ∼ = 0.664 and γ ∼ = 0.69.
The ballistic spreading that is dominated by a local velocity shear is another commonly encountered spreading regime.The separation between particles grows in proportion to time.Ballistic spreading can often be observed in nonsteady realistic oceanic flows at timescales that are much shorter than the onset of diffusive spreading (which develops after a trajectory samples multiple different eddies or other flow features).To derive a connection between the encounter volume and the velocity shear, we consider a trajectory that is advected by a flow field with a constant meridional velocity shear, γ , of the zonal velocity.In a reference frame moving with a reference trajectory, the velocity profile is u(y) = γ y, where u denotes the x component of velocity and the encounter volume becomes (2) This suggests a linear growth with time for a ballistic regime.Note that expression (2) quantifies the encounter volume as a volume of fluid that is initially located outside of the encounter sphere and that passes through the sphere over time T .To include the volume of fluid that is initially located within the encounter sphere (or within the encounter circle in this two-dimenstional case), one needs to add π R 2 to expression (2).The contribution of this term becomes negligibly small as T → ∞.Expression (2) has been tested numer- ically and shows good agreement with the numerically estimated encounter volume for a linear shear flow (Fig. 10b).
The steady linear saddle flow with a constant strain rate α and the following velocities is another commonly considered example that is often used to approximate the vicinity of a hyperbolic trajectory in more complicated, nonsteady and nonlinear situations.A unique property of this flow is that the velocity profile is unchanged in any reference frame moving with a trajectory.This can be shown by applying the coordinate transformation x = x − x tr (t); ŷ = y − y tr (t), where (x; y) are coordinates in a stationary frame, ( x; ŷ) are coordinates in a moving frame and (x tr (t); y tr (t)) is the trajectory.The velocity in a moving frame is then where the last equality holds because dx tr dt = αx tr ; dy tr dt = −αy tr .Thus, without a loss of generality, we can consider a flow in a reference frame, which is moving with a reference trajectory that is located at the origin.The encounter volume that comes within a radius R of the origin over the time interval T can be written as where dx and dy denote the grid spacing between neighboring trajectories, and the flux in the trajectories entering the circle is given by Nonlin.Processes Geophys., 24, 189-202, 2017 www.nonlin-processes-geophys.net/24/189/2017/ Again, as in our treatment of the linear shear flow, expression (5) does not include the volume of fluid that is initially located within the encounter sphere (or the encounter circle in this two-dimensional case), but only the volume that was initially located outside that passes through the sphere over time T .The contribution of that fixed volume (π R 2 ) becomes negligibly small as T → ∞.Here, u ⊥ is the inwardlooking normal component of the velocity at a circle of radius R, and ds is an infinitesimal segment of the circle arc.From symmetry, the flux is the same in each of the four quadrants, so we will consider the first quadrant only.From geometry (Fig. 11), and ds = rdβ, leading to and Adding the other three quadrants then gives The numerical simulations of the encounter volume in a linear strain flow show excellent agreement with theoretical expression (9) (Fig. 10a).
The linear growth of the encounter volume with time in the linear shear and linear strain flows could be anticipated by noting that both flows are steady in a reference frame moving with a reference trajectory, and all particles only encounter the origin once and never come back.Thus, the flux through the encounter circle is constant in time, and the encounter volume, which is a time integral of flux, is proportional to time.The random walk flow seems to be different because the particles can encounter the reference trajectory more than once, leading to a nonsteady flux of first encounters and a nonlinear time dependence of the encounter volume.
4 Mixing potential for a specified tracer: the u * approach The above examples are centered on the mixing potential of a flow field, but there may be value in computing the encounter volume for swarms of trajectories of biological organisms, drifting sensors and other non-Lagrangian trajectories.For example, if one is interested in the actual transport of scalar properties such as heat, salt or vorticity, then it may be useful to calculate V using a velocity field that is directly linked to the vector flux of the scalar of interest.This approach has been used in connection with heat transport in advective and diffusive flows (Bejan, 1995;Costa, 2006;Mahmud and Fraser, 2007;Mukhopadhyay et al., 2002, andSpeetjens, 2012) and more recently with the transport of more general scalars in forced and dissipative (and possibly turbulent) flows (Pratt et al., 2016).The central idea is to a define velocity field u * based on the (known) flux F of a scalar with concentration C. Here, the bold italic quantities denote vectors.The concentration is assumed to obey a conservation equation of the form where S contains the sources and sinks of C. The velocity u * is defined as the velocity of a hypothetical flow in which the flux of C is purely advective: F = Cu * .Pratt et al. (2016) show that, in the absence of sources or sinks of C, the total amount of C contained within any material boundary advected by this hypothetical flow is conserved: d dt V CdV = 0. Thus, u * is linked to scalar property fluxes, while u is limited to fluid volume (or area) fluxes.
If indeed F is due entirely to advection by the actual fluid velocity field u, then u * = u.More generally, F will contain contributions from eddy fluxes, molecular or subgrid diffusion, and even forcing and dissipation terms that can be expressed as the divergence of a flux.In addition, F may be augmented by the addition of any non-divergent vector without altering Eq. ( 10).As shown by Speetjens (2012), this lack of uniqueness can be dealt with by defining a physically relevant reference scalar distribution and then focusing on the flux of the scalar anomaly, which is an approach that we adapt below.Thus, by estimating the encounter volume V for trajectories of the u * field, one is quantifying the rate at which different "parcels" of tracer anomaly are brought into contact with each other.An example is presented next.We applied the encounter volume diagnostic to quantify the mixing potential for a specific tracer in the Bickley jet flow.Our goal is to describe an example where the u * field for a given tracer is significantly different from the flow velocity u and where the corresponding encounter volume field for a given tracer, V * , is significantly different from the water particle trajectory-based encounter volume V .
Consider the Bickley jet flow with the same parameters as in Sect.2.2 and assume that one is interested in a tracer that, at initial time t 0 , has a uniform value c 0 south of the jet and a constant meridional gradient north of the jet; i.e., C 0 = c 0 + 0.5y(sign(y − 5 × 10 5 ) + 1) with c 0 = 1.Ignoring the diffusive terms, the tracer evolution is governed by the advection equation ∂C ∂t = −∇ • (uC), where u is the Bickley jet flow velocity.Since the jet core acts as a transport barrier separating the northern and southern chaotic zones, this tracer will rapidly filament and develop high property gradients north of the jet, but will remain uniform south of the jet.So, despite the fact that the mixing potential of the Bickley jet flow is exactly the same on both sides of the jet (Fig. 7b), stirring will not lead to mixing for this particular tracer distribution south of the jet where the tracer gradient is zero, thus leading to zero mixing potential for this particular tracer.We seek to capture this effect by applying a mixing diagnostic based on the encounter volume to the corresponding u * field for this tracer.
In the spirit of Speetjens (2012), we regard c 0 as the reference concentration, here constant, and define F to be the flux of a tracer anomaly: ) is zero south of the jet where C = c 0 and approximately equal to u north of the jet where C c 0 , leading to the u * -based encounter number V * = 0 south of the jet and V * ≈ V north of the jet.This behavior was further numerically validated in Fig. 12, where we first numerically simulated the evolution of this tracer in the Bickley jet flow and then estimated u * , counted N * and estimated V ∼ = NdV for the trajectories advected by the u * field.The result confirms that the mixing potential for this tracer is zero south of the jet, V * = 0. North of the jet, however, V * is very close to V from Fig. 7b.Thus, by combining the u * approach with the encounter volume idea, we were able to correctly capture the mixing potential for a specific tracer.

Summary and discussion
When water parcels come into direct contact with each other, they can exchange water properties, leading to mixing.The trajectory encounter volume, V , quantifies the volume of fluid that passes close to a reference trajectory over a time interval t 0 < t < t 0 + T .Thus, the encounter volume is proportional to, and can be used as a measure of, the mixing potential of a flow.For incompressible flows densely seeded with particles, the encounter volume can be approximated by V ∼ = N δV , where N is the encounter number, i.e., the number of trajectories that come come within radius R from the reference trajectory over time t 0 < t < t 0 + T , and δV is a small volume element.
The encounter volume diagnostic was tested in three flows with increasing complexity: the Duffing oscillator, the Bickley jet and the altimetry-based velocity in the Gulf Stream extension region.In all cases, V was smaller within the cores of the coherent eddies and jets, where the mixing potential was low, and V was larger in the chaotic zones near the peripheries of the eddies and at the flanks of the meandering jets, where the mixing potential of the flow was high.
Similar to finite-time Lyapunov exponents (FTLEs) that are commonly used to delineate regions with qualitatively different motion (Haller, 2002;Shadden et al., 2005;Lekien and Ross, 2010), V depends on the trajectory starting time, t 0 , which allows us to track the evolution of oceanic features by repeating the calculation at different t 0 , and on the trajectory integration time, T , revealing different structures that impact the mixing potential of the flow from time t 0 to time t 0 + T .Specifically, longer segments of stable or unstable manifolds emanating from the hyperbolic regions are revealed for longer T in forward or backward time.In the long T limit, when both the stable and unstable manifolds densely fill the entire chaotic zone, V approaches a constant equaling the volume of the chaotic zone.
V also depends on the encounter radius R, which defines how close two trajectories need to be in order to be counted as an encounter.The analytic arguments and numerical simulations both suggest that R on the order of a fraction (onethird) of the radius of the smallest feature of interest should work well in most cases.
Finally, while V was initially introduced in the continuous limit of infinitely many and infinitely small fluid elements (i.e., an infinitely dense grid of initial positions), its approximation V ∼ = N δV depends on the initial spacing between the neighboring trajectories.The numerical simulations suggest that this approximation works well for a grid spacing as large as R/2 (with the appropriately chosen R as discussed above) and that the major effect of increasing the grid spacing is the degraded resolution of the resulting V map rather than incorrect V values.
As with FTLEs, complexity measures (Rypina et al., 2011), Lagrangian descriptors (Mendoza et al., 2014) and other techniques from the dynamical systems theory (Beron-Vera et al., 2013;Budisic and Mezic, 2012;Froyland et al., 2007;Haller et al., 2016), V can be computed for forwardand backward-time trajectories, with the backward computation revealing the unstable manifolds.Our encounter number could plausibly be related, in a limiting case, to the mixing geometry of Karrash and Keller (2017).
For a ballistic spreading regime dominated by the velocity shear γ and for the linear saddle flow with a constant strain α, V was shown to be proportional to γ T and αT , respectively.The linear growth of the encounter number with integration time for the linear shear and the linear strain flows is a consequence of the steady flux of first encounters through the encounter circle.
An analytical connection between the encounter volume and a widely used measure of mixing, the diffusivity K, would be a desirable result for parameterizing the effects of eddies in numerical models.Some initial developments towards deriving such a formula were outlined for a diffusive random walk process.It was numerically shown that the dependence of V on time is nonlinear, but the numerical simulations were too inconclusive to make further inferences.
The mixing potential is the property of the flow field and characterizes the intensity of stirring, whereas the actual tracer mixing depends on both the flow and the tracer.For example, no tracer mixing will occur if the tracer gradient is zero, even if the mixing potential of the flow is high.To address this, we have proposed combining the encounter number diagnostic with the u * approach of Pratt et al. (2016) for characterizing the mixing potential for a specific tracer C. The u * depends on, and includes information about, the tracer fluxes.In the absence of sources and sinks of C, the amount of tracer is conserved within any Lagrangian volume advected by u * , so the encounter volume V * computed for the trajectories advected by u * can be used to quantify the mixing potential for a specific tracer.An example was presented where V * for a specified tracer distribution in the Bickley jet flow was significantly different from V in a large part of the domain.
The encounter volume is a frame-independent quantity because it is based on the relative distances between water parcel trajectories rather than on the properties of isolated trajectories.The encounter volume values do not change under the orthogonal transformations of the coordinates, i.e., under the rotations and translations of a reference frame.This is a desirable property because the ability of a flow to mix tracers should not depend on the reference frame.
The encounter volume and, more generally, the encounter mass ideas presented in this paper are not restricted to two dimensions and can be used to quantify the mixing potential in three-dimensional flows.This framework also does not require incompressibility and can work with unstructured irregular grids.The investigation of the performance of the method in quantifying the mixing potential of a flow in more complicated cases is left for a future study.

Figure 1 .
Figure 1.The trajectory segments for the different integration times (a) and the Poincaré section (b) for the Duffing oscillator.

Figure 2 .
Figure2.The encounter volume for the Duffing oscillator for the various integration times from T = 0.1T pert (on the left) to T = 50T pert (on the right) and for the various encounter radii from R = 0.01 (on the top) to R = 1 (on the bottom).The trajectories were released on a regular grid spanning the entire domain with a grid spacing of 0.013 in both the x and y directions.

Figure 3 .
Figure3.The Poincaré section (the black dots; same as in Fig.1b) superimposed onto the encounter volume (in color; same as the top and middle right panels in Fig.2).Only select trajectories from the Poincaré section are shown.

Figure 4 .
Figure4.The encounter volume (in color; the same as the second row and the second column subplot of Fig.2) and the stable (black) and unstable (white) manifolds for the Duffing oscillator flow computed using the direct method.

Figure 5 .
Figure 5.A comparison between the FTLEs (a) and the encounter volume (bottom; same as the middle row of Fig. 2) for the Duffing oscillator flow for various integration times from T = 0.1T pert = 0.13 (on the left) to T = 50T pert = 66.67 (on the right).The same set of trajectories deployed on a dense initial grid with a 0.02 grid spacing is used in all simulations.(b) R = 0.2.

Figure 6 .
Figure 6.The encounter volume, V , for the Duffing oscillator flow for various grids of the initial positions from a dense grid spacing of 0.02 (a), to an intermediate grid spacing of 0.04 (b) to a coarse grid spacing of 0.1 (c).The encounter radius, R = 0.2, and integration time, T = 6.67, are the same in all three simulations.

Figure 7 .
Figure 7.The Poincaré section (a) and the encounter volume V (b)for the Bickley jet flow.For the V calculation, the trajectories were released on a regular grid spanning the entire domain with a grid spacing of about 10 5 in both the x and y directions.

Figure 8 .
Figure 8.The encounter volume for the AVISO velocities in the Gulf Stream extension region for the trajectories released on 11 July 1997 and integrated over 30 days (a), 60 days (b) and 90 days (c).The trajectories were released on a regular grid spanning the domain from 65 to 35 • W and 30 to 50 • N with a grid spacing of about 0.06 • in both longitude and latitude.Additional simulations were performed to ensure that the release domain was sufficiently large and that a further increase of the release domain does not lead to changes in the encounter volume for the trajectories starting in the subdomain shown.

Figure 9 .
Figure 9.The positions of the trajectories that were initially located within the eddy core on 11 July 1997 (blue patch) after 30 days (green), 60 days (red) and 90 days (yellow) of integration.The background shows the kinetic energy of the flow as a snapshot on 11 July 1997.
where dx and dy correspond to the displacements of a particle in the x and y directions at an instance in time, and the subscripts m and ref denote the moving reference frame and the reference trajectory, respectively.Denoting the angle in which the step is taken by ϕ, the displacements are dx = L cos ϕ, dy = L sin ϕ, dx ref = L cos ϕ ref and dy ref = L sin ϕ ref , leading to L m = 2L sin α, where α = ϕ−ϕ ref 2 .Since both ϕ and ϕ ref are random variables uniformly distributed between 0 and 2π , α is a random variable with a flat pdf distribution ∈ [0; π ].

Figure 10 .
Figure 10.A comparison between the numerically computed encounter volume (blue) and the analytical predictions (Eqs.8 and 9; red) for the linear strain (a) and linear shear flows (b).For the linear shear flow, α = 0.1, r = 5 and dx = dy = r/25; for the linear strain flow, γ = 0.1, r = 5 and dx = dy = r/25.The other parameter choices show good agreement as well.

Figure 11 .
Figure 11.A schematic diagram for estimating the encounter number for a linear saddle.

Figure 12 .
Figure 12.The u * -based encounter volume, V * , for a tracer with an initial distribution south of the jet and a constant meridional gradient north of the jet.