Clustering of buoyant tracer in quasi-geostrophic coherent structures

Abstract

We have investigated the dynamics of floating tracer in an idealised turbulent quasi-geostrophic ocean by advecting Lagrangian particles in a high-resolution velocity field enhanced by the potential flow associated with vortex stretching. At first order in the Rossby number expansion, this component of the ageostrophic circulation can be derived through a diagnostic equation in terms of the geostrophic velocities. Borrowing methods from the theory of Lagrangian coherent structures, we identify coherent material loops around strong vortex cores using the Lagrangian averaged vorticity deviation (LAVD). Building on studies of clustering in kinematic, stochastic velocity fields, we utilise methods from statistical topography to show that the coherent vortices dominate the distribution of extreme values of the concentration field. We find that the presence of clusters and voids in a coherent vortex depends on more than just the sense of rotation, but also on the full evolution of the vorticity over its lifecycle. We identify the mechanism behind the cluster formation that respects the symmetries of the quasi-geostrophic equations but can be expected to hold robustly in more complicated regimes, due to the simple physical description. The association of cluster formation with vortex stretching implies that LAVD is a particularly relevant metric for floating tracer dynamics. The detection of intense clustering also has implications for reaction rates between ocean-borne flotsam, meaning that our results are relevant to understanding the general risk of floating microplastics and marine biological populations.

1. Introduction

Spatial heterogeneity in the distribution of floating material is ubiquitous in observations. Whether the subject matter is plankton populations, microplastic waste or floating drifters, the organisation of these materials into dense clusters and sparse voids is self-evident (Jordi, Basterretxea & Anglès 2009; Andrady 2011; Maximenko, Hafner & Niiler 2012; Olascoaga et al. 2013; Pratt, True & Crimaldi 2017; Little et al. 2018). Understanding the processes behind cluster formation is relevant to a number of open questions, such as microplastic fate modelling (Cózar 2014), and the risk of pollutants to the marine ecosystem (Meacham & Berloff 2024). Structures in the distribution of these materials are often qualitatively linked to flow structures in the ocean currents (Lehahn et al. 2011), suggesting an important geophysical component to clustering.

The dynamics of buoyant material can be investigated at several levels of complexity. Beginning with fully three-dimensional models of inertial particles, such as in Reartes & Mininni (2023), there is an observed tendency for buoyant particles in stratified flows to settle at the level of their neutral buoyancy. This leads to a situation where vertical motions of the particles are inhibited, but the material is still passively carried by the flow in the horizontal plane. For significantly positively buoyant particles, such as some microplastics (Andrady 2011) and plankton populations (Gross, Zeuthen & Yonge 1948), the floating material is then trapped at the surface. At this point, the flotsam is carried by the horizontal surface currents exclusively. This mechanism provides the basis for clustering, since even incompressible ocean currents can have a divergent horizontal component, which leads to a weakly compressible flow of floating tracer. On the global scale, it has been argued that Ekman currents provide the dominant surface convergence, and that these determine the locations of the so called ‘convergence zones’ of plastic pollution (Onink et al. 2019). However, at smaller scales, especially at the lower boundary of the mesoscale, we still observe large gradients of chlorophyll and other tracers (Jordi et al. 2009). At these scales, there is a small but significant ageostrophic component to the surface flow, which is naturally divergent, since it is associated with non-zero vertical velocity below the surface.

Geostrophic flows are manifestly non-divergent in the horizontal plane, so cannot be a direct cause of cluster formation. At the quasi-geostrophic level, which includes all terms up to order one in the Rossby number, there is a necessary divergent ageostrophic flow that facilitates the time evolution of the geostrophic component (Gent & McWilliams 1983). This ageostrophic divergence can be defined through a diagnostic equation in terms of the geostrophic flow, then used to reconstruct the potential component of the flow. We examine the dynamics of cluster formation in the geostrophic flow enhanced with this potential component.

We utilise a Lagrangian approach to model clustering. Advantages of this approach include preferential resolution of clustered regions, due to the convergence of Lagrangian paths. It also avoids the inclusion of hyper-diffusivities that are necessary for numerical stability in grid-based models. These are known to damp clustering and destroy the spatial structure of clusters (Meacham & Berloff 2023). The Lagrangian approach also allows for easier inter-comparability with observational studies. Lagrangian coherent structure metrics, such as the finite time Lyapunov exponent, dilation, and Lagrangian averaged vorticity deviation (LAVD), are increasingly utilised in an observational setting (Olascoaga et al. 2013; Huntley et al. 2015; Abernathey & Haller 2018; Rypina et al. 2022). The LAVD is particularly useful for identifying rotationally coherent vortices, which are expected to have distinguished transport properties (Haller et al. 2016), while the dilation is related to the surface divergence of the velocity, so is particularly relevant to the dynamics of floating tracers.

Building on a series of studies that have investigated clustering in kinematic stochastic fields (Koshel et al. 2019; Meacham & Berloff 2023, 2024) and deterministic mesoscale fields with a kinematic, stochastic divergent sub-mesoscale (Stepanov et al. 2020a,b), we move to a fully dynamical regime. Insights from the kinematic studies will be useful in our investigation. First, metrics from statistical topography, such as cluster mass and area (Klyatskin 2003), have proven incredibly effective in characterising spatiotemporal features of clustering (Koshel et al. 2019). Second, the sensitivity of the clustering process has been investigated (Meacham & Berloff 2023), showing that poor representation of temporal statistics can lead to spurious clustering. This motivated the use of offline velocity fields with very fine time resolution, verified with frequency spectrum analysis. However, the dynamical fields have a character very different to that of their kinematic counterparts, because they are dominated by coherent structures. This turns out to be of great importance for the clustering process.

Combining the clustering metrics with Lagrangian coherent structure techniques reveals rich dynamics for buoyant particles that are strongly dependent on the temporal coherence as well as instantaneous flow features. When considering leading-order contributions to the ageostrophic divergent flow, it becomes evident that the association between clusters and coherent structures is not just a consequence of the long-time stability but is also directly related to quasi-geostrophic dynamics.

2. Methods

2.1. Two-layer quasi-geostrophic ocean model

We study the behaviour of floating tracers in an idealised ocean model where we can generate high-resolution velocity and divergence data. The ocean is a doubly periodic square basin with side length 3600 km. It is a mid-latitude beta plane with local Rossby radius of deformation 25 km. In the vertical, there are two layers. The upper layer is shallow and less dense, with thickness 1 km, while the deeper layer has thickness 3 km. The ocean is bounded by a rigid lid and a flat bottom boundary. We opt for a quasi-geostrophic model to set the dynamics of the ocean currents. Quasi-geostrophy is most commonly derived from an expansion of the fluid variables in increasing orders of the Rossby number (Vallis 2017a). This analysis can be applied in many scenarios, but here it is considered as a limit of the multi-layer shallow-water equations.

At first order in the Rossby number, the equations for the zeroth-order contributions (the geostrophic currents) become a closed set. The geostrophic velocity can be expressed in terms of streamfunctions $\psi _1$ and $\psi _2$ in the upper and lower layers, respectively. The associated potential vorticity equations are

(2.1)$$\begin{gather} \frac{\partial \zeta_1}{\partial t} + J(\psi_1,\zeta_1)+\beta\,\frac{\partial \psi_1}{\partial x} = \nu\,\nabla^4\psi_1, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \zeta_2}{\partial t} + J(\psi_2,\zeta_2) +\beta\,\frac{\partial \psi_2}{\partial x} = \nu\,\nabla^4\psi_2-\gamma\,\nabla^2\psi_2, \end{gather}$$

where $\zeta _2,\zeta _1$ are the perturbation vorticities, given by

(2.3)$$\begin{gather} \zeta_1 = \nabla^2\psi_1 + S_1 (\psi_2-\psi_1), \end{gather}$$

(2.4)$$\begin{gather}\zeta_2 = \nabla^2\psi_2 + S_2 (\psi_1-\psi_2), \end{gather}$$

with the so-called ‘stratification parameters’ $S_1$ and $S_2$, which depend on the density contrast and height difference between layers. The advection terms are expressed using the Jacobian, which is defined as

(2.5)\begin{equation} J(\psi_i,\zeta_j) = \frac{\partial \psi_i}{\partial x}\,\frac{\partial \zeta_j}{\partial y} - \frac{\partial \psi_i}{\partial y}\,\frac{\partial \zeta_j}{\partial x}. \end{equation}

On the left-hand sides of (2.1) and (2.2), we have the usual advective terms. In addition to this, on the right-hand sides, we have an ‘eddy’ viscosity in both layers, and a bottom friction in the lower layer. The viscosity will act to stabilise our numerical solutions, but will otherwise be negligible. On the other hand, bottom friction will lead to a decay in the total energy in the system, unless it is resisted by some forcing mechanism. To maintain a steady turbulent velocity field, we follow the method of Karabasov, Berloff & Goloviznin (2009). A fixed background shear is imposed by adding a small, but constant westward current to the upper-layer velocity. In mathematical notation, this is equivalent to applying the transform

(2.6)\begin{equation} \psi_1\rightarrow\psi_1-Uy, \end{equation}

where $y$ is the Cartesian coordinate in the north–south direction, and $U$ is some negative constant velocity. A chaotic turbulent flow is then generated because the eddy field draws energy from the background shear through baroclinic instability (Berloff, Kamenkovich & Pedlosky 2009). By injecting energy at the largest scale, a sustained energy cascade develops, which leads to the formation of a mix of cyclonic and anticyclonic vortices. The upper-layer streamfunction is initialised as a small Gaussian bump, which provides the simulation with some kinetic energy, allowing baroclinic instability to develop.

The system of (2.1)–(2.4) was integrated using the CABARET numerical scheme as described in Karabasov et al. (2009). This numerical scheme provides many advantages for high Reynolds number (turbulent) modelling, due to its conservation properties and low dispersion. First, there was a 2200 day spin-up period, until a steady state for the various domain averaged kinetic energies was reached. This was followed by 260 days of high temporal resolution simulation, which was output at half-day intervals. Some $4096\times 4096$ grid points were used for the spatial discretisation, such that the deformation radius was well resolved. The parameters used for the numerical simulation are summarised in table 1.

Table 1. Parameters relevant to the two-layer quasi-geostrophic system.

From the quasi-geostrophic streamfunctions, we can calculate the geostrophic velocity. For example, in the upper layer it is given by

(2.7)\begin{equation} (u_g,v_g)=\left(-\frac{\partial \psi_1}{\partial y},\frac{\partial \psi_1}{\partial x}\right). \end{equation}

In the derivation of the quasi-geostrophic system, higher-order corrections to the velocity are specifically neglected. Therefore, the total velocity is not available at this level of approximation. This is problematic from the perspective of a simulation of floating tracer clustering, since a non-zero potential velocity is required. However, the divergence of the ageostrophic velocity at first order in the Rossby number is essential to balance the vertical velocities of the moving layer interface, and can be calculated from the geostrophic velocities as a result. To demonstrate this, consider the displacement of the interface between the layers, which can be expressed in terms of the geostrophic streamfunctions:

(2.8)\begin{equation} \eta=\frac{f_0}{g}\,(\psi_2-\psi_1). \end{equation}

As the interface moves deeper, the depth of the upper layer increases. To conserve mass of fluid in the column, there must be an associated convergent flow. This convergence can be provided by the ageostrophic flow only due to the non-divergence of geostrophic velocities. To order one in the Rossby number expansion, the ageostrophic divergence is given by (Vallis 2017a)

(2.9)\begin{equation} \frac{\partial u_a}{\partial x}+\frac{\partial v_a}{\partial y}=\frac{S_1}{f_0} \left(\frac{\partial }{\partial t}+u_g\,\frac{\partial}{\partial x } +v_g\, \frac{\partial}{\partial y}\right) (\psi_2-\psi_1). \end{equation}

If we further assume that the ageostrophic velocities are purely potential, then this equation diagnoses the velocity field upon solving the equation

(2.10)\begin{equation} \frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}=\frac{S_1}{f_0} \left\{\frac{\partial }{\partial t}(\psi_2-\psi_1)+ J(\psi_1,\psi_2)\right\} \end{equation}

subject to doubly periodic boundary conditions. The potential $\phi$ allows us to recover the ageostrophic velocity:

(2.11)\begin{equation} (u_a,v_a)=\left(\frac{\partial \phi}{\partial x},\frac{\partial \phi}{\partial y}\right) \end{equation}

Figure 1(a) shows the energy spectra of the geostrophic ($u_g$), ageostrophic ($u_a$) and total ($u_g+u_a$) velocities in the upper layer, as a function of the magnitude of the wavevector, k. It can be seen clearly that the geostrophic flow is significantly more energetic than the ageostrophic flow, as would be expected from the assumptions of the quasi-geostrophic equations. The energy spectrum of the geostrophic component demonstrates the classic scaling of two-dimensional turbulence from the Kolmogorov theory (Kolmogorov 1941), which is also shown in figure 1(a), denoted by $\varepsilon$. This consists of an inverse cascade of energy towards large scales, associated with a $-\frac {5}{3}$ spectrum, as well as the direct cascade of enstrophy in the inertial range, with a power-law spectrum of order $-3$ (Vallis 2017b). There is a peak in the spectrum at the Rossby radius of deformation $R_d$, which is the scale of the largest vortices. At the lower end of the inertial range, we see the effects of dissipation and a steep loss of spectral power. The ageostrophic velocity contrasts with this, as it exhibits a very flat spectrum. In spite of this, the total velocity has a spectrum that is almost indistinguishable from the geostrophic flow. However, it also lies below the geostrophic spectrum, implying that the ageostrophic flow acts against the geostrophic currents on average.

Figure 1. Various spectra from the quasi-geostrophic simulation. (a) Temporally averaged energy spectra for geostrophic ($u_g$), ageostrophic ($u_a$) and total ($u_g+u_a$) velocities. (b) The scale-dependent Rossby number, which is the square root of the ratio of the ageostrophic energy spectrum to the geostrophic energy spectrum. The vertical dashed line in (a,b) denotes the wavelength associated with the deformation radius ($R_d$). (c) The frequency power spectrum of the geostrophic velocity and the ageostrophic divergence. (The vertical axis on the left-hand side is for the power spectral density (PSD) of $u_g$, and the vertical axis on the right-hand side is for the divergence $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}_a$.) The vertical dashed line shows the per day frequency. (d–f) Zoomed-in plots of the upper-layer streamfunction $\psi _1$, displacement of the layer interface $\eta$, and the ageostrophic divergence $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}_a$, respectively. The domain is chosen to contain a coherent cyclone (upper half) and coherent anticyclone (lower half).

By taking the square root of the ratios of the energy spectra of the ageostrophic to geostrophic components, we can define a scale-dependent Rossby number $Ro_k$, which is shown in figure 1(b). At the deformation radius and above (for $k<{1}/{R_d}$), this number is vanishingly small. As we move from the ‘mesoscale’ of the system towards the ‘sub-mesoscale’, the Rossby number approaches unity, implying that there is a strong potential component at smaller scales. Averaging over all scales leads to Rossby number $0.21$, which represents a significantly ageostrophic system, contrasting with the global Rossby number of approximately $0.015$.

The energy spectra of turbulent flows as a function of wavenumber have been studied extensively in a range of settings (Kolmogorov 1941; Rhines 1979; Tulloch & Smith 2006; Callies et al. 2015; Cosentino, Simon & Morales-Juberías 2019). Temporal spectra are often neglected, and there is generally less interest in quantifying the typical time scales of variability. Previous studies of buoyant and inertial particle clustering have shown that resolving these time scales is essential to accurate modelling. This is because the concentration of buoyant tracer depends on the whole time history of divergence along Lagrangian paths, so it is sensitive to related errors. In particular, overestimating time correlations by using poor temporal interpolation has been shown to lead to spurious clustering (Meacham & Berloff 2023).

To ensure that we have resolved all relevant time scales, frequency power spectral densities were calculated for the geostrophic velocity and ageostrophic divergence, and these are shown in figure 1(c). For the half-day sampled high-resolution velocity data, the high-frequency power is several orders of magnitude smaller than the slowest frequencies. This implies that all time variability has been resolved satisfactorily at this sampling rate, so the data are sufficient for the clustering study. Notably, the spectrum of the ageostrophic divergence is skewed towards lower frequencies compared to the geostrophic velocity. This is opposite to the trend observed in the spatial spectra, but with a simple explanation. In the quasi-geostrophic model, the large-scale structures, such as the geostrophic vortices, are the ones that excite a significant ageostrophic response. This response is at smaller scales than for the generating structures, but inherits the slow time variability. Hence the mesoscale forcing produces compressible sub-mesoscale velocities. This can be seen qualitatively from snapshots of the various fields of the quasi-geostrophic simulation. Figures 1(d–f) show the streamfunction $\psi _1$, the layer interface displacement $\eta$, and the ageostrophic divergence $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}_a$, respectively, around a pair of coherent vortices. The anticyclone in the lower half has a strong core of rotational geostrophic velocity (which can be seen in figure 1d), and coincides with a deepening of the layer interface (shown in figure 1e). The response in the ageostrophic divergence field is a small region of divergence that lies predominantly on the boundary of the vortex. The cyclone in the upper half is the same but reversed, with a shallow upper layer and a strong convergent region. Compared to the first two fields, the regions of non-zero divergence are clearly separated from the non-divergent background. This demonstrates how the ageostrophic response is tied to the large scales of the geostrophic field.

2.2. Lagrangian particle model

We perform two separate Lagrangian simulations in our quasi-geostrophic velocity field. One simulation represents the dynamics of buoyant tracers, which follow both the geostrophic and leading-order ageostrophic surface flow. The other simulation involves only the geostrophic flow, and is used to detect the coherent vortices in the simulation.

For the purposes of these simulations, we take the upper-layer quasi-geostrophic velocities to also represent the surface velocity. This is consistent with the discretisation of the vertical structure into distinct layers in the derivation of a layered quasi-geostrophic system. Limited vertical variation within each layer of velocity and other characteristics is implicitly assumed for this model reduction to be valid. In the real ocean, limited vertical variation in mesoscale properties is expected within the mixed-layer depth (Kara, Rochford & Hurlburt 2003; Vallis 2017c). We require the surface velocity to model floating tracer, since it is assumed that it resides almost exactly at the surface due to the effect of positive buoyancy.

2.2.1. Buoyant tracers

We consider buoyant particles that are trapped at the top surface of the upper layer, but which follow horizontal flows passively. This represents a reduction of the advection problem to two dimensions, where the buoyant particles follow the surface Lagrangian paths. This is a commonly used approximation for the modelling of such material (Maximenko et al. 2012; Koshel et al. 2019; Stepanov et al. 2020a,b; Vic et al. 2022; Maalouly et al. 2023; Meacham & Berloff 2023, 2024). Since the geostrophic flow is horizontally non-divergent, it cannot lead to clustering of buoyant tracer. We construct a velocity field that contains the leading-order contributions to both the solenoidal and potential flows. The path of a particle in this flow satisfies

(2.12)$$\begin{gather} \frac{{\rm d}}{{\rm d} t}\boldsymbol{x}(t;\boldsymbol{x}_0)= \boldsymbol{u}_g(\boldsymbol{x}(t;\boldsymbol{x}_0),t)+ \boldsymbol{u}_a(\boldsymbol{x}(t;\boldsymbol{x}_0),t), \end{gather}$$

(2.13)$$\begin{gather}\boldsymbol{x}(0;\boldsymbol{x}_0)=\boldsymbol{x}_0, \end{gather}$$

where $u_g$ and $u_a$ are specifically the upper-layer geostrophic and ageostrophic velocities.

The construction of the upper-layer surface velocities $u_g$ and $u_a$ is described in § 2.1. Additionally, we consider the concentration of buoyant tracer in this flow:

(2.14)\begin{equation} \frac{\partial C}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}((\boldsymbol{u}_g+\boldsymbol{u}_a) C)=0, \end{equation}

where $\boldsymbol {\nabla }=({\partial }/{\partial x},{\partial }/{\partial y})$ is the horizontal gradient operator. This is simply the conservation of tracer mass, following the particle paths. It can be expressed as a material derivative along paths:

(2.15)\begin{equation} \frac{{\rm d}}{{\rm d} t}C(\boldsymbol{x}(t;\boldsymbol{x}_0),t) ={-} C(\boldsymbol{x}(t;\boldsymbol{x}_0),t)\,\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_a(\boldsymbol{x}(t;\boldsymbol{x}_0),t). \end{equation}

Solving this equation in conjunction with (2.12)–(2.13) allows us to reconstruct the concentration field of buoyant tracer from an ensemble of Lagrangian particles.

We initialise a regular grid of these buoyant tracer particles in the upper layer of the quasi-geostrophic flow. Each particle begins with unit concentration. We then integrate (2.12) and (2.15) in time for each particle using the pseudo-symplectic Runge–Kutta algorithm.

2.2.2. Virtual particles and coherent structures

Motivated by observations, we know that the spatial distribution of floating material is strongly linked to large-scale flow structures such as geostrophic vortices and jets (Brach et al. 2018). Geostrophic velocities are some of the most readily available observations in oceanography, thanks to satellite altimetry (Park 2004; Sánchez-Reales et al. 2012; Doglioni et al. 2023). As a result, it makes sense to try to understand the observed distributions in terms of the geostrophic velocity alone. Observations of sub-mesoscale currents are always improving, and this is a very active area of research. Measurements of large surface divergences exist, but these measurements are often local, provided by drifter estimates or mooring arrays (Callies, Barkan & Garabato 2020; Zhang et al. 2021) that are naturally limited in scope compared to the global coverage of mesoscale current data.

In the quasi-geostrophic setting that is our focus, the divergent velocity is simply a perturbation to the stronger non-divergent flow. This can be verified from the kinetic energies contained in each component, as shown in figure 1. While the surface divergence is a necessary condition for clustering, we can still expect that the distribution and spatial extent of clusters will largely be set by structures in the non-divergent flow. First, this is a consequence of the relative strength of the two components. Second, it also relates to the fact that the divergent ageostrophic flow in quasi-geostrophy does not undergo separate dynamics, but is instead entirely determined by the geostrophic currents.

Temporal dynamics of the velocity field have been shown to be crucial in determining clustering, even in stochastic kinematic models with limited coherence (Meacham & Berloff 2023). This is because the concentration of a particle depends on the entire time history of divergence along the Lagrangian path, not just the instantaneous value. In the dynamical setting considered here, there are long-lived coherent structures that will have a large imprint on temporally averaged quantities, even if the instantaneous strength of divergence is small.

To identify coherent vortices in the geostrophic flow, we follow the LAVD method of coherent structure detection (Haller et al. 2016). This is an objective (frame-independent) method that identifies closed material curves that are minimally deformed over time.

We solve the following equations, for an ensemble of initial conditions $\boldsymbol {x}_0$:

(2.16)$$\begin{gather} \frac{{\rm d}}{{\rm d} t}\boldsymbol{y}(t;\boldsymbol{x}_0)=\boldsymbol{u}_g(\kern0.09em \boldsymbol{y}(t;\boldsymbol{x}_0),t), \end{gather}$$

(2.17)$$\begin{gather}\boldsymbol{y}(0;\boldsymbol{x}_0)=\boldsymbol{x}_0, \end{gather}$$

(2.18)$$\begin{gather}LAVD(t,\boldsymbol{x}_0)=\int_0^t {\rm d}\tau [\omega(\kern0.09em \boldsymbol{y}(\tau;\boldsymbol{x}_0),\tau)-\bar{\omega}(\tau)], \end{gather}$$

where $\omega ={\partial v_g}/{\partial x}-{\partial u_g}/{\partial y}$ is the relative vorticity of the geostrophic velocity, and $\bar {\omega }(t)$ is the domain averaged vorticity at time $t$. This is the LAVD as defined in Haller et al. (2016); however, we can simplify calculations in our case because the domain averaged vorticity in this scenario is always zero. Once we have constructed the LAVD field as a function of the initial conditions, we look for stationary points in this field. The Lagrangian path $\boldsymbol {y}$ that is initialised at a stationary point is then the vortex centre. Finally, we look for the outermost contour of the LAVD field around each vortex centre that is minimally deformed over the full length of the quasi-geostrophic simulation (260 days in this case). We first find the outermost closed contour of the initial LAVD field. After advecting this contour as a material curve in the geostrophic velocity for the full 260 days, we find the convex hull of the final curve. If the area of the convex hull is significantly larger than the initial area, then this implies that filamentation has occurred, and the contour is not coherent. We find the outermost closed contour for which the change in area is below a tolerance parameter that can be chosen, and identify this contour as the vortex boundary. The method of choosing the boundary is relatively unimportant because the edge of each vortex is clearly associated with a rapid decrease in LAVD. As a result, most methods will identify similar boundaries. Figure 2 shows the LAVD field computed with a $1000\times 1000$ grid of virtual particles, as well as the vortex boundaries of $15$ coherent vortices that were observed in our model. Figure 6 below shows the evolution of a large cyclonic vortex, with the concentration field of particles inside the vortex boundary.

Figure 2. (a) The LAVD over the full 6-month model run. Red curves show the cyclonic vortex boundaries, and black curves show anticyclonic boundaries on the first day. (b) Trajectories of the 15 coherent vortices detected using the LAVD method. Red trajectories correspond to cyclones, and black to anticyclones. Vortices begin at circular markers and finish at crosses. Each vortex has been assigned an index in order of their LAVD magnitude at the vortex core. (c) Deviation of the vortex core depth from the stationary depth of the upper layer, averaged over cyclones (red curve) and anticyclones (black curve).

2.2.3. Lagrangian numerical schemes

In order to integrate the Lagrangian paths and concentration, a six-step pseudo-symplectic Runge–Kutta algorithm was implemented, as constructed in Aubry & Chartier (1998). The six-step method is third-order accurate in time, but preserves the area of material volumes with sixth-order accuracy. This is a desirable property, since it ensures that there is limited numerical diffusion of particles by the non-divergent flow. It is especially relevant to the virtual particles used to construct the LAVD field, since we expect the area of vortices to be conserved.

To evaluate gridded velocity and divergence data at particle positions, we use linear interpolation in space and time. The $4096\times 4096$ resolution of the velocity field means that higher-order interpolations are unnecessary, although cubic interpolation has also been implemented, and it was found that there is a negligible effect on the results of the clustering simulation. Similarly, the temporal sampling time (0.5 days) is much shorter than the temporal variability of the flow, as can be seen in figure 1(c). The spectral power at the sampling frequency is five orders of magnitude smaller than the peak.

For the clustering simulation, $16\times 10^6$ particles were initialised in a square grid, whereas $10^6$ virtual particles were used in the LAVD calculation. A large number of particles is needed once clustering can occur in order to resolve the voids for long periods of time.

2.3. Clustering measures

Methods from statistical topography have been applied to the clustering of floating tracers, mostly in the context of tracers in stochastic kinematic velocity fields (Klyatskin 2003; Koshel et al. 2019; Stepanov et al. 2020a,b). Two important quantities used to characterise the rate of clustering are the cluster mass and area, which are defined as

(2.19)$$\begin{gather} M(t,\bar{C}) = \left\{\int{\rm d}^2{\boldsymbol{x}}_0\,C(\boldsymbol{x}_0,0)\right\}^{{-}1} \left\{\int{\rm d}^2{\boldsymbol{x}_0}\,\varTheta(C(\boldsymbol{x}(t;\boldsymbol{x}_0),t)-\bar{C})\right\}, \end{gather}$$

(2.20)$$\begin{gather}S(t,\bar{C}) = \left\{\int {\rm d}^2\boldsymbol{x}_0\,C(\boldsymbol{x}_0,0)\right\}^{{-}1}\left\{\int {\rm d}^2{\boldsymbol{x}_0}\, \frac{\varTheta(C(\boldsymbol{x}(t;\boldsymbol{x}_0),t)-\bar{C})}{C(\boldsymbol{x}(t;\boldsymbol{x}_0),t)}\right\}, \end{gather}$$

where $\varTheta$ is the Heaviside step function, and $\boldsymbol {x}(t;\boldsymbol {x}_0)$ are the buoyant tracer particles as defined in § 2.2.1. The cluster mass $M(t,\bar {C})$ measures the quantity of tracer that has clustered above concentration $\bar {C}$, whereas the cluster area measures the area bounded by $\bar {C}$ isolines of the concentration field. Expressed as an integral over the initial conditions, these quantities are readily calculated using the Lagrangian data from our simulations. Observations of stochastic potential flow show exponential clustering, and this is also expected from the analytical asymptotic theory. In an exponential clustering process, $M(t,\bar {C}) \rightarrow 1$ and $S(t,\bar {C}) \rightarrow 0$ as $t\rightarrow \infty$ for any $\bar {C}$, leading to the creation of point-like singular clusters. The presence of a solenoidal field can slow the process, but it is still generally exponential in character (Koshel et al. 2019).

To isolate the contribution of the coherent vortices to the clustering process, we can restrict the cluster mass and area to the vortex cores,

(2.21)\begin{align} M_{\mathcal{V}}(t,\bar{C}) &= \left\{\int{\rm d}^2{\boldsymbol{x}}_0\, C(\boldsymbol{x}_0,0)\right\}^{{-}1} \nonumber\\ &\quad \times \left\{\iint{\rm d}^2{\boldsymbol{x}_0}\,{\rm d}^2{\boldsymbol{y}_0}\, \varTheta(C(\boldsymbol{x}(t;\boldsymbol{x}_0),t)-\bar{C})\,\varTheta (\boldsymbol{x}(t;\boldsymbol{x}_0)-\boldsymbol{y}(t;\boldsymbol{y}_0))\right\}, \end{align}

and similarly for the cluster area,

(2.22)\begin{align} S_{\mathcal{V}}(t,\bar{C}) &= \left\{\int{\rm d}^2{\boldsymbol{x}}_0\,C(\boldsymbol{x}_0,0)\right\}^{{-}1}\nonumber\\ &\quad \times \left\{\iint{\rm d}^2{\boldsymbol{x}_0}\,{\rm d}^2{\boldsymbol{y}_0}\, \frac{\varTheta(C(\boldsymbol{x}(t;\boldsymbol{x}_0),t)-\bar{C})}{C(\boldsymbol{x}(t;\boldsymbol{x}_0),t)}\, \varTheta(\boldsymbol{x}(t;\boldsymbol{x}_0)-\boldsymbol{y}(t;\boldsymbol{y}_0))\right\}, \end{align}

where the integral over $\boldsymbol {y}_0$ is over all initial positions inside the vortex boundaries, and $\boldsymbol {y}(t;\boldsymbol {y}_0)$ is a Lagrangian path following just the geostrophic velocity as in § 2.2.2.

To investigate clustering in individual coherent vortices, we can look at the total mass inside the vortex boundary. In terms of Lagrangian particles, this is

(2.23)\begin{equation} \mathcal{M}_{i}(t) = \int_{\mathcal{D}_0}\int_{A_i}{\rm d}^2\boldsymbol{x_0}\,{\rm d}^2\boldsymbol{y}_0\, \varTheta(\boldsymbol{x}(t;\boldsymbol{x}_0)-\boldsymbol{y}(t;\boldsymbol{y}_0)), \end{equation}

where $A_i$ is all initial conditions $\boldsymbol {y}_0$ inside the $i$th vortex boundary at time $t=0$, and $\mathcal {D}_0$ is all initial conditions within the doubly periodic domain $\boldsymbol {x}_0$.

3. Results

After the initial spin-up period, the geostrophic virtual particles were introduced to the developed turbulence in a uniform grid of $1000\times 1000$ particles. The LAVD field was calculated over 260 days of Lagrangian integration in the high temporal resolution velocity field. Using the method described in the previous section, 15 coherent vortices were detected that existed at the beginning of the high-resolution dataset. Of these, there were four anticyclones and eleven cyclones. Figure 2(a) shows the initial boundaries of the vortices superimposed over the LAVD field. Each vortex is assigned an index that corresponds to the size of the LAVD signal at the vortex core (i.e. vortex 0 is the largest vortex in terms of LAVD magnitude, and vortex 14 has the smallest magnitude). While additional peaks in the LAVD field can be seen, they were not associated with a coherent material curve that survived until the end of the simulation. We are restricting our analysis to the most coherent vortices that survive over the whole period of simulation.

The vortex core point is always inside the material vortex boundary by definition. We can use the layer depth at the core as a indicator of the strength of the vortex. Figure 2(c) shows the vortex core depth over the course of the simulation, averaged separately over anticyclonic and cyclonic vortices. From this, we can see that the initially strong vortices decay quickly, in spite of their material coherence. The majority of the vorticity is shed in the first 80 days of the simulation.

In addition to the geostrophic particles, we introduce the buoyant tracer particles to the upper layer in a uniform grid of $4000\times 4000$ particles, each with unit concentration. We then run the Lagrangian scheme for 260 days. We observe dense clusters on cyclonic vortices, whilst the anticyclones tend to clear particles from their orbit, leading to the formation of sparse voids. Figure 3 shows the concentration field after the first 60 days of the floating tracer Lagrangian simulation. Clearly, we can see the formation of spiral clusters centred on point-like cores, which are very reminiscent of vortex structures. It is this distribution of material that motivates an investigation of the relationship to the coherent vortices. We will show through the use of clustering metrics and the results of the LAVD analysis that the densest clusters are directly associated with the most coherent vortices. The largest cyclone (vortex 0) contains almost 50 times its initial mass at the end of the simulation; meanwhile, the largest anticyclone (vortex 1) expels all 2500 particles inside its boundary within the first 60 days of the simulation. This can be seen in figures 4(a–d), which contain all the vortex masses for each individual coherent vortex. Figures 4(a–c) show the cyclonic vortices, whereas figure 4(d) shows the anticyclonic vortices.

Figure 3. The concentration field over the whole domain after 60 days of the floating tracer Lagrangian simulation.

Figure 4. Mass curves for buoyant tracer inside vortex boundary normalised by the initial mass in the vortex, for (a–c) cyclones, and (d) anticyclones.

Figure 5(a) shows the ratio of vortex cluster mass $M_{\mathcal {V}}$ to total cluster mass $M$, as defined in § 2.3. The total area of the coherent vortices is approximately $5\,\%$ of the domain area. However, for larger concentration reference values $\bar {C}$, the vortices dominate the distribution. For the first $100$ days of the simulation, between $80\,\%$ and $100\,\%$ of particles with concentration over $\bar {C}=100.0$ lie in the detected vortices. This share decays over the course of the simulation, but slowly. Even after the full 260 days, more than $75\,\%$ of the areas of highest concentration lie within the material boundaries. This is in spite of the decay of the vortices and the emergence of new vortices that are not considered in our methodology. The overall trend, observable in figure 4, is that cyclones form clusters, whereas anticyclones form voids. This is the same trend that has been observed in other studies (Vic et al. 2022; Maalouly et al. 2023). Examples of this behaviour are shown in the supplementary movies (https://doi.org/10.1017/jfm.2024.1229), which show the evolution of the concentration field following a handful of coherent vortices. However, it cannot be the whole explanation, since there are notable exceptions. Vortex 6 is a cyclone, yet it is initially expelling tracer, and ends the simulation as a weak cluster. Figure 6 shows six snapshots of the concentration field inside vortex 6. The first two snapshots show the creation of a weak void, before material begins to spiral in to the centre, with a chirality that matches the cyclonic rotation. Vortex 11 is a cyclone that undergoes multiple phases of expulsion and clusterisation. Figure 7 demonstrates this. Vortex 11 clearly forms the ubiquitous spiral cluster in the first four snapshots, before rapidly voiding in the fifth (at day 200). It then forms a new cluster in the final two snapshots, in a way that is consistent with the vortex mass calculation (figure 7). For the anticyclones, the majority (three out of four) end up as voids. However, vortex 7 ends up as a weak cluster, and vortex 10 initially clusters up to approximately day 50, before it expels all its particles in the remaining time. Figure 8 shows the evolution of the buoyant tracer concentration inside vortex 7. Material initially accumulates at the boundary, as tracer is pushed out from the centre, before spiralling back into the vortex in a spiral with the opposite chirality to the cyclonic clusters. This cluster is also distinct from the cyclonic counterparts because it retains the local minimum at the centre of the cluster, as opposed to a local maximum. Figure 9 shows vortex 10, which initially attracts a small amount of material to its boundary, before rapidly expelling all particles. What these particular vortices demonstrate is that the association between vortex sign and cluster formation is more complicated than we may initially expect.

Figure 5. Cluster masses and areas as defined in § 2.3. (a) Ratio of vortex cluster mass to total cluster mass as a function of time $t$ and reference concentration $\bar {C}$. (b) Total cluster mass. (c) Same as (a), but for the cluster areas. (d) Total cluster area.

Figure 6. Snapshots of the concentration field around and within vortex 6 (a cyclone) at 10, 50, 100, 150, 200 and 250 days, moving clockwise from (a) at the top left to ( f).

Figure 7. Same as figure 6, but for vortex 11 (a cyclone).

Figure 8. Same as figure 6, but for vortex 7 (an anticyclone).

Figure 9. Same as figure 6, but for vortex 10 (an anticyclone). Some plotting artefacts can be seen within the vortex core in (e, f), because the concentration field is being reconstructed from its values known only at the Lagrangian particle locations, and there are no particles within the vortex. This means that there is no resolution of the anticyclone vortex core.

To understand the relation between vortex evolution and cluster/void formation, we can examine the dynamics behind the potential ageostrophic circulation. Consider a coherent vortex that occupies an area $A(t)$ at time $t$, with boundary $\partial A(t)$. We can define a flux function of buoyant material into the vortex by integrating the ageostrophic velocity around the boundary:

(3.1)\begin{equation} F(t) = \oint_{\partial A(t)} \boldsymbol{u}_a \boldsymbol{\cdot} \boldsymbol{n}\,{\rm d} s, \end{equation}

where the boundary has been parametrised by $s$, and $\boldsymbol {n}$ is the normal to the boundary. From here, we can see that $F>0$ corresponds to floating tracer leaving the boundary (voiding), whereas $F<0$ corresponds to an accumulation of material (clustering). Applying the two-dimensional divergence theorem to this integral, and substituting for the diagnostic equation for the ageostrophic divergence (2.9), leads to

(3.2)\begin{equation} F(t) = \int_{A(t)} \frac{S_1}{f_0}\,\frac{{\rm D}}{{\rm D} t}(\psi_2-\psi_1)\,{\rm d} A, \end{equation}

where ${{\rm D}}/{{\rm D} t}$ is the material derivative with respect to the geostrophic velocity. Additionally, using the potential vorticity equation in the upper layer (ignoring the viscous term) will allow us to eliminate any dependence on the lower-layer velocities, i.e.

(3.3)\begin{equation} F(t)={-}\int_{A(t)}\left[\frac{{\rm D}\omega_1}{{\rm D} t}+\beta\,\frac{\partial \psi_1}{\partial x}\right]{\rm d} A, \end{equation}

which then further simplifies using the fact that $\partial A(t)$ is a material curve with respect to the geostrophic velocity:

(3.4)\begin{equation} F(t)={-}\frac{{\rm d}}{{\rm d} t} \int_{A(t)} (\omega_1+f_0+\beta y)\,{\rm d} A. \end{equation}

As a result of this, it becomes clear that the clustering behaviour of a vortex is largely determined by the rate of change of absolute vorticity in its interior. In particular, strong cyclonic vortices will expel floating tracer as they spin up, but will accumulate material as they spin down. Due to having the opposite sign for vorticity, anticyclonic vortices will accumulate material as they spin up, and expel as they spin down. This clearly explains the observed behaviour, since we can see from figure 2 that all the coherent vortices are weakening on average, so cyclones are in a clustering state, and anticyclones are in a repulsive state. However, this is different to a statement that cyclones are clusters and anticyclones are voids. Instead, the observations are as much a consequence of initialising a uniform concentration field in an established velocity field where the state of the vortices has already been set. By considering the time-reversal symmetry of quasi-geostrophy, we might expect that cyclones and anticyclones will spend a comparable amount of time as clusters and voids over a large enough number of realisations, since they map to each other under the action of the symmetry. The different behaviours of the two structures are simply due to the fact that a time-reversed cluster is a void, and vice versa. If we include more complex geophysical considerations, such as next-order corrections to quasi-geostrophy, an asymmetry between cyclones and anticyclones will develop (Maalouly et al. 2023). While these dynamics may have significant consequences for the clustering process, our results show that this asymmetry is not required to observe the expected behaviour.

The above consideration may also help to explain why LAVD detection methods so effectively identified the regions of most intense clustering, since areas with high LAVD values have experienced large changes in vorticity, associated with intense clustering.

4. Discussion

Our results imply that coherent vortices defined through Lagrangian averaged vorticity deviation (LAVD) will be efficient transporters of floating material for reasons beyond coherence itself. These structures do more than simply move existing inhomogeneities in the distribution of material for long times. They also have the potential to create inhomogeneities by attracting material in one location before shedding it in another. In fact, our results show an incredibly powerful ability of these vortices to ‘vacuum’ up surrounding floating tracer, with interior concentrations often hundreds of times greater than their surroundings.

We have observed a novel relationship between the direction of spin of a coherent vortex and its clustering behaviour. In particular, both cyclones and anticyclones have an equal propensity for forming clusters and voids of material. A single vortex can also change its behaviour depending on its evolution and whether it is growing or decaying. Specifically, anticyclones cluster as they grow and expel as they decay, with the opposite relation for cyclones. We argue that this is a natural consequence of the symmetry of quasi-geostrophy, but it still has implications for more complex models, where the tendency is to look for a direct, rather than dynamic, relationship between vortex sign and clustering. Looking towards more realistic scenarios, this will likely have consequences for the distribution of material in the open ocean. For example, enhanced flow of floating tracer in the western boundary currents could be expected. Boundary currents, such as the Agulhas current, contain strong, energetic and coherent geostrophic vortices (Beron-Vera et al. 2013; Casanova-Masjoan et al. 2017). These are prime candidates for the kind of cluster formation that we have observed here, especially since they form along coasts, which are naturally areas of high concentration of both marine populations and ocean pollutants (Lehahn et al. 2011). Our results would suggest that anticyclonic vortices form at the region of retroflection of the current, just south of the African coastline. Here, they will pick up material in their vicinity. Then, as they move eastwards into the open ocean of the southern Atlantic, they may rapidly expel this material. This would constitute a significantly stronger transport pathway from the plankton- and pollutant-rich coastlines to the ocean interior than with passive transport.

The relevance of LAVD to floating tracers also allows for reinterpretation of existing observations, in terms of the potential for clustering and buoyant tracer transport. For instance, in Rypina et al. (2022), both the LAVD and dilation were calculated for an approximately $11\,{\rm km}\times 11\,{\rm km}$ array of floating drifters. Dilation is defined as the Lagrangian averaged divergence. If we separate variables and integrate (2.15) in time, then we see that the concentration of floating tracer is simply the exponential of the negative of the dilation:

(4.1)\begin{equation} C(\boldsymbol{x}(t;\boldsymbol{x}_0),t)=C_0\exp\left(-\int_0^t{\rm d}\tau\, \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}(\tau;\boldsymbol{x}_0),\tau)\right). \end{equation}

The LAVD and dilation, whilst not entirely coincident in this experiment, were found to identify similar coherent regions. Although there will be an impact of many additional processes not included in our idealised model, this would be consistent with our own observations in the numerical setting.

If the characteristics of clustering observed here remain robust as we move to more realistic modelling, then the theory of how and where these clusters should form could aid in observation and removal efforts of plastic waste in the open ocean. The order of magnitude variations in concentration observed here imply that there is room for more efficient searching strategies for such waste, which could be defined in terms of ocean current data, and perhaps even Lagrangian coherent structure detection will be of use here. Methods from statistical topography, particularly cluster mass and area statistics, have allowed us to exactly quantify how much of the clustering is due to the presence of the coherent structures. That only $5\,\%$ of the domain can contain the majority of regions of concentration over 10 times the initial concentration (and $100\,\%$ of the regions over 100 times the initial concentration) is a remarkable result that indicates the importance of these features in determining the distribution of floating tracers. That many vortices ended the simulation with 10–100 times the mass of tracer with which they were initialised is surprising for a low Rossby number, near geostrophic system. However, it is clear that even the mesoscale ocean weather systems can significantly contribute to the formation of clusters in ocean flotsam.

Our results can motivate the definition of new identification metrics for coherent flow structures that specifically relate to their ability to accumulate floating tracer. For instance, we have used LAVD, which is an integral of vorticity over time, to identify coherent structures. However, (3.4) shows that clustering of floating tracer is more closely related to the change in vorticity along a Lagrangian path, so considering this quantity might lead to a more successful metric. Additionally, we may consider an LAVD metric that includes the contribution of the $\beta$ terms to the absolute vorticity.

We can also consider other flow regimes of the quasi-geostrophic model, such as the jet-dominated (Berloff et al. 2009) or oceanic gyre (Karabasov et al. 2009) solutions. These contain coherent jets and westward boundary currents respectively, and may exhibit unique clustering of floating tracer. Baroclinic vortex pulsar solutions can also provide a more detailed model of vortex–wave interactions, which leads to novel evolution of non-isolated coherent vortices (Berloff & Sutyrin 2024). Given the strong dependence of clustering on the evolution of coherent vortices seen in this work, we may be able to gain more insight into vortex clustering in such a solution. Exploring each of these regimes individually can aid in building a more comprehensive picture of clustering in the marine environment, as each of these quasi-geostrophic coherent structures has analogues in real ocean currents.

Having built up from kinematic, stochastic models of turbulence (Koshel et al. 2019; Meacham & Berloff 2023) to a fully dynamical model of ocean currents, we have observed similar behaviour in terms of the global cluster mass and area statistics. This implies that there is a similar exponential clustering process. However, the existence of the temporal coherence changes the spatial character of the clusters, and leads to more intense clustering at the extreme end of the concentration distribution. Given that these results from the kinematic studies have proved to be robust, we might expect that other properties will also carry over to the dynamic model. In Meacham & Berloff (2024), it was shown that coextensive clustering between different floating species can lead to enhanced reaction rates between them. Even for a weak clustering process, it can alter the global equilibria of a system, motivating further study in this direction.