Lagrangian study of dispersion and transport by submesoscale currents at an upper-ocean front

Density fronts, ubiquitous in the upper ocean, are an important source of submesoscale dynamics. The dynamics typically occur at length scales of 0.1 - 10 $\,\rm{km}$ and time scale of O(1) day and are characterized by Rossby number, $Ro\equiv U/fL=O(1)$, where $f$ is the Coriolis parameter, and $U$ and $L$ are characteristic velocity and length scales, respectively (Thomas et al., 2008; McWilliams, 2016). The submesoscale dynamics plays a significant role in the restratification of the upper ocean and the vertical transport of tracers such as buoyancy, salinity and carbon from the surface ocean to the interior (Boccaletti et al., 2007; Thomas et al., 2008; Fox-Kemper et al., 2008; Omand et al., 2015). These processes affect the the upper-ocean structure and impact the interactions between the ocean and the atmosphere, thereby influencing the Earth’s climate. The submesoscale dynamics also play a significant role in ocean’s biochemical cycle by aiding phytoplankton growth through supply of nutrients from the upper thermocline into the surface layer (Mahadevan, 2016).

Many of the upper-ocean processes driven by the submesoscale dynamics are possible because of their ability to develop large vertical velocity (Mahadevan and Tandon, 2006), presumably with spatial and temporal coherence. This is in contrast to the small-scale turbulent motions, which are relevant for the local mixing, or the balanced mesoscale motions in which the vertical velocity is orders of magnitude smaller. The lateral transport is believed to be dominated by the mesoscale currents and eddies, but the role of submesoscales can be significant as they can provide interconnections between the mesoscale transport barriers and enhance horizontal spread (Haza et al., 2016). The submesoscales are also important for predicting the dispersion of buoyant pollutants such as oil (D’Asaro et al., 2018). An understanding of the organization of vertical velocity and transport pathways is therefore crucial for understanding the submesoscale upper-ocean transport and dispersion processes.

Because of their size and relatively fast dynamics, the submesoscale motions have been difficult to investigate using conventional observational methods: ships surveys and satellite remote sensing. However, recent observations employing innovative techniques have uncovered some interesting features of the submesoscale dynamics. By measuring horizontal velocity synchronously along two-parallel tracks, Shcherbina et al. (2013) were able to calculate the velocity gradient tensor at $O(1)\,\rm{km}$ in the North Atlantic Mode Water region where there is an active submesoscale. Their observations were consistent with dynamics associated with a predominance of filaments of $O(f)$ cyclonic vorticity in a soup of relatively weak anti-cyclonic vorticity. The filament structures with cyclonic vorticity are known to develop through frontogenesis that can occur due to straining of the front by a large scale confluent flow. A front can also undergo frontogenesis through non-linear evolution of baroclinic instability (BI) (Hoskins and Bretherton, 1972; Hoskins, 1982). The initial stages of the frontogenetic development of BI at an atmospheric front has been studied in detail by Mudrick (1974). Recent studies have shown that the interaction of a cold filament in thermal wind balance with boundary layer turbulence can drive secondary circulations in the lateral-vertical plane that is frontogenetic and restratifies the filament within a few hours (McWilliams et al., 2015; Sullivan and McWilliams, 2018). The ageostrophic circulation in the case of especially strong fronts can lead to nonlinear bores (Pham and Sarkar, 2018). Filament structures with cyclonic vorticity were also observed in the northern Gulf of Mexico in an observational campaign utilizing a large number of satellite-tracked surface drifters (D’Asaro et al., 2018). The structures were smaller than $1\,\rm{km}$ in width, separated dense water mass from the light water mass, and were found to be convergent, attracting surface drifters into a line which then wrapped into a cyclonic eddy. The convergence of water mass implies downwelling, and the measured vertical velocity was as large as $1-2\,\rm{cm\,s^{-1}}$. In comparison, the typical vertical velocity at a mesoscale front is $O(0.01)\,\rm{cm\,s^{-1}}$ (Rudnick, 1996).

The evolution of BI in upper-ocean density fronts is an important mechanism for generating submesoscale currents. The problem has been studied extensively using large-scale ocean models (Capet et al., 2008) and turbulence resolving models (Skyllingstad and Samelson, 2012; Hamlington et al., 2014; Stamper and Taylor, 2017; Verma et al., 2019). Simulating a density front that is initially in thermal wind balance, Verma et al. (2019) (hereafter VPS19) find that the evolution of BI generates long, thin vortex filaments with cyclonic vorticity and downwelling vertical velocity that roll into coherent submesoscale eddies. These submesoscale filaments and the large vertical velocity inside them are similar to the submesoscale filament-like features observed during the the surface drifter measurements of D’Asaro et al. (2018). VPS19 showed that the coherent structures, i.e., vortex filaments and eddies, provide a 3D organization to the secondary circulation whose velocity field suggests that water is transported laterally and vertically across the front. Although there are organized 3D structures, the actual paths followed by the fluid parcels over time are not apparent from the instantaneous velocity field as the dynamics is transient. Furthermore, the spatial pattern of the velocity field changes when the coherent structures are transported by the mean down-front jet. A Lagrangian framework is better suited for a study of material transport by the submesoscale, and is the subject of this paper. A related problem is about the time scale of subduction and restratification of the front. The vertical velocity observed in the filaments can be so large as to produce vertical displacement of $O(1)\,\rm{km}$ in a day if sustained in magnitude and direction. However, the restratification is likely to progresses on the time scale of baroclinic instability which is $O(2\pi/f)$ (Stone, 1966). Here, we show evidence of a slow restratification at near-inertial time scale emerging from relatively fast motions in the filament structures.

Lagrangian drifters and floats have been widely used in the ocean for understanding local flow properties and dynamics (e.g. see review article of LaCasce (2008)). Single-particle metrics are used for calculating the mean flow and eddy kinetic energy (Richardson, 1983; Fratantoni, 2001; Jakobsen et al., 2003) and the eddy diffusivities (Zhurbas and Oh, 2003) in different parts of the ocean and have been utilized for investigating the local transport of tracers (Davis, 1985). The metric of single-particle dispersion is also a convenient tool for predicting the spread of a particle from the point of release by a velocity field which has coherent unsteady currents and turbulence.

Particle-pair statistics are often used to probe the scale dependence of dynamics. In a recent analysis of the trajectories of surface drifters deployed during the Grand Lagrangian Dynamics (GLAD) experiment in the Gulf of Mexico, Poje et al. (2014) found that the second-order structure function of the velocity field showed power-law behavior from 200 m to 100 km, including the submesoscale range, suggesting the dominance of local advective dynamics. Balwada et al. (2016) applied a Helmholtz decomposition of the second-order structure function computed from the GLAD drifters into divergent and rotational components finding that the divergent component dominated at scales below $5\,\rm{km}$, and also computed the third-order structure function. From their analysis, they inferred forward 3-D energy cascade below $5\,\rm{km}$, 2-D enstrophy cascade between 5 to 40 km (the deformation radius), and an inverse energy cascade between 40 - 100 km. Beron-Vera and LaCasce (2016) examined pair-separation statistics in the submesoscale range, using synthetic drifter trajectories from data-assimilated NCOM simulations conducted with 1 km horizontal resolution. They found that the pair separation grew exponentially in accord with non-local dynamics. They further attributed the discrepancy of their result with the results from GLAD observational drifter trajectories to the strong inertial oscillations experienced by the GLAD drifters and their limited number of independent samples with possibly low statistical significance. In the ocean, internal gravity waves are likely to further complicate the statistical measure arising due to submesoscale dynamics.

Two-particle statistics also reveals the spread of the particles about the center of mass (COM) of a particle cloud (Batchelor, 1952). Multiparticle studies have also been used, mainly to measure flow properties such as relative vorticity and horizontal divergence (Molinari and Kirwan Jr, 1975). Multiparticle statistics using groups of four particles (tetrads) can be used to investigate the changes in the shape of the particle clusters, which result from straining by the large-scale flows and dispersion by the finescale motions, as shown by Pumir et al. (2000).

In this study, we employ the model front of VPS19 to investigate the transport and dispersion characteristics of the submesoscale currents, including finescale turbulence, during the evolution of BI when the dynamics is dominated by the coherent vortex filaments and eddies with localized finescales. VPS19 does not contain strong inertial motions, internal gravity waves or surface forcing, thus providing an ideal setup for studying the dispersion characteristics of submesoscale currents in isolation. Owing to the high resolution (2 m in all three directions) of the simulation, we are able to capture 3D turbulence generated during the evolution of BI. The domain size of 4 km captures a wide range of the submesoscale but not the mesoscale. The study is performed in the Lagrangian framework by releasing a large number of tracer particles, which move with the local fluid velocity.

The paper is structured as follows. In section 2, the numerical method and the setup of the model front of this work is described. In section 3, the generation of coherent structures such as vortex filaments and eddies by the baroclinic instability and their organization in 3D are summarized. The details about the tracer-particle simulation are given in section 4. In section 5, the 3D organization of the particles, the typical features of the transport pathways, and the correlation between the lateral and vertical motions of the particles circulating at the front are examined. In section 6, the motion of particle clouds, each cloud containing fluid of similar density, is studied by monitoring the vertical motion of their centers of mass and the dispersion of the particles about the cloud centers. Additionally, the transport of mean properties such as temperature, subgrid viscosity and kinetic energy associated of with constituent particles are studied. Section 7 focuses on the dispersion characteristics of the submesoscale currents including the localized three-dimensional finescale associated with the currents. In this section, single- and pair-particle dispersion statistics are investigated, and shape changes following groups of four particles (tetrads) are reported. Finally in section 8, conclusions are drawn based on the results, with a brief discussion about their implications.

The model used here is the same as the one employed by VPS19. Nevertheless, we describe the model for completeness. The model consists of an upper-ocean front in thermal wind balance with a surface jet. The width of the front is $L=1.2\,\rm{km}$ and is confined in a surface layer of depth $H=50\,\rm{m}$ situated above a strongly stratified thermocline. The density variation is due to the changes in temperature, and both the quantities are assumed to be related by a linear equation of state, ${\rho}/{\rho_{0}}=-\alpha T$, where $\alpha=2\times 10^{-4}\,\rm{K}^{-1}$ is the coefficient of thermal expansion, $\rho$ is the density deviation from a reference density $\rho_{0}=1028\,\rm{kg\,m^{-3}}$, and $T$ is the temperature deviation from a reference temperature. The temperature profile is given by

Here, the front is assumed to be aligned with the $x$ direction (along-front), and the temperature variation is in the $y$ direction (cross-front) and the $z$ direction (vertical), which also coincides with the axis of rotation; $M_{0}^{2}$ is the value of $M^{2}=-(g/\rho_{0})\partial\rho/\partial y$ at the center $y=0$, where $M^{2}$ is defined analogous to the square of buoyancy frequency associated with the vertical density gradient, $N^{2}=-(g/\rho_{0})\partial\rho/\partial z$; $N^{2}_{M}$ and $N^{2}_{T}$ are the squared buoyancy frequencies in the surface layer and the thermocline, respectively; $\delta_{H}=5\,\rm{m}$ is a thin region between the surface layer and the thermocline where the temperature profile joins smoothly from its value in the surface layer to that in the thermocline; $\rm{g}=9.81\,\rm{m}\,\rm{s}^{-2}$ is the gravitational acceleration. The values of the temperature profile parameters used in the simulation are $M^{2}_{0}=1.5\times 10^{-7}\,\rm{s^{-2}}$, $N^{2}_{M}=3.0\times 10^{-7}\,\rm{s^{-2}}$ and $N^{2}_{T}=10^{-5}\,\rm{s^{-2}}$.

The surface jet, $U(y,z)$, is constructed from the density field by integrating the thermal wind relation, $\partial U/\partial z=-M^{2}/f$, where $f=1.4\times 10^{-4}\,\rm{s}^{-1}$ is the Coriolis parameter. Additionally, a broadband velocity noise with amplitude of $10^{-4}\,\rm{m}\,\rm{s}^{-1}$ is superimposed to the frontal jet for instigating the instabilities.

The contours of initial velocity, temperature and potential vorticity over a $y$-$z$ plane are shown in Fig. 1. The potential vorticity is defined as $\Pi=(\boldsymbol{\omega}+f\mathbf{k})\cdot\boldsymbol{\nabla}b$, where $\boldsymbol{\omega}$ is the relative vorticity, $\mathbf{k}$ is a unit vector in the vertical direction and $b=\alpha Tg$ is the buoyancy. As shown in Fig. 1, the potential vorticity at the front is initially negative and the setup is unstable to symmetric perturbations.

The evolution of the model front is studied by numerical means, utilizing the large eddy simulation (LES) approach and solving the non-hydrostatic Navier-Stokes equations under the Boussinesq approximation. Along-front velocity $u_{1}$, cross-front velocity $u_{2}$, vertical velocity $u_{3}$, temperature $T$ and dynamic pressure $p$ are advanced in time $t$ as follows:

where $i,\,j,\,k=1,\,2,\,3$, and a repeated index implies summation; $\nu$ is the molecular viscosity and $\kappa$ is the molecular diffusivity; $\tau^{sgs}_{ij}=-\nu^{sgs}(\partial u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i})$ is the modeled LES subgrid stress tensor and $q^{sgs}_{j}=-\kappa^{sgs}(\partial T/\partial x_{j})$ is the modeled LES subgrid heat flux, with $\nu^{sgs}$ and $\kappa^{sgs}$ representing the subgrid viscosity and diffusivity, respectively. Parameters $\nu$ and $\kappa$ are related by the Prandtl number $Pr=\nu/\kappa$; the value of molecular viscosity used is $\nu=10^{-6}\,\rm{m^{2}s^{-1}}$, and the Prandtl number, $Pr=7$. The subgrid scale $\nu^{sgs}$ and $\kappa^{sgs}$ are related by a turbulent Prandtl number taken to be $Pr^{sgs}=1$. An alternate notation for the velocity components is also used wherein the along-front, cross-front, and vertical velocity components are expressed as $u$, $v$, and $w$, respectively.

When Eq. 2 is scaled by the velocity scale $U_{0}=M_{0}^{2}H/f$, the maximum geostrophic jet velocity at the surface, and the buoyancy scale $N^{2}_{M}H$, the non-dimensional parameters are as follows: the Ekman number, $Ek=\nu/fH^{2}$, the non-dimensional lateral buoyancy gradient, $M_{0}^{2}/f^{2}$, and the Richardson number, $Ri=N^{2}_{M}f^{2}/M_{0}^{4}$. In the present study, $Ri=0.26$ and $Ek=2.86\times 10^{-6}$. The ratio $M_{0}^{2}/f^{2}=7.65$ is comparable to the values used in the studies of Skyllingstad and Samelson (2012) and Hamlington et al. (2014). Also, note that the Rossby number, $Ro=U_{0}/fL$, based on the initial horizontal shear is 0.32 and the Reynolds number, $Re=U_{0}H/\nu$, is $2.67\times 10^{6}$.

The subgrid fluxes are parametrized following Ducros et al. (1996). First, the subgrid viscosity, $\nu^{sgs}$, is calculated, and then the subgrid diffusivity of temperature, $\kappa^{sgs}$, is predicted knowing the turbulent Prandtl number, $Pr^{sgs}$. The subgrid viscosity, $\nu^{sgs}$, is computed dynamically at every grid point $(i,j,k)$ using a local velocity structure function $F$:

where $C_{K}=0.5$ is the Kolmogorov constant, $\Delta=||\Delta x_{i}||$ is the magnitude of the filter grid spacing, and

For calculating $F(x,\Delta x_{i},t)$, the velocity field $\tilde{\mathbf{u}}_{i,j,k}$ is obtained by passing the LES velocity through a discrete high-pass Laplacian filter. The parametrization is efficient in predicting $\nu^{sgs}$, and the values are substantial only at grid points with large velocity fluctuation. Once $\nu^{sgs}$ is known, the subgrid diffusivity $\kappa^{sgs}$ is calculated assuming $Pr^{sgs}=1$. Note that the dynamic Ducros model used here has been employed in several previous studies, including the oceanic examples of turbulent baroclinic eddies (Skyllingstad and Samelson, 2012) and the formation of gravity currents from strong fronts (Pham and Sarkar, 2018).

The computational domain is also the same. It is a rectangular box bounded by $0\leq x\leq 4098\,\rm{m}$, $-3073\,\rm{m}\leq y\leq 3073\,\rm{m}$ and $-130\,\rm{m}\leq z\leq 0$. The domain is discretized in two different ways during during the solution. A uniform grid with $2050\times 3074\times 66$ points in employed initially, providing a grid resolution of 2 m in each direction. Later, during the evolution of baroclinic instability, the solution is obtained using a grid which is exactly the same in the horizontal, but has $98$ grid points in the vertical, with uniform stretching such that the grid spacing changes from $0.5\,\rm{m}$ near the surface to $1.5\,\rm{m}$ near the bottom of the surface layer. The finer grid resolution in the vertical is needed near the surface to resolve the surface intensified turbulence in the vortex filaments that develops during the nonlinear evolution of BI. As in VPS19, the domain size is chosen to accommodate the growth of the most unstable baroclinic mode (Stone, 1966) whose wavelength, $L_{b}$, and the time scale, $\tau_{b}$, are:

With the parameters used in the present study, the chosen domain is large enough to accommodate at least two wavelengths of the most unstable baroclinic mode.

For obtaining the numerical solution, one also needs to specify the boundary conditions appropriately. We consider the domain to be periodic in the along-front direction. Free-slip on the velocity and no-flux on the temperature are used as the boundary conditions at the surface ($z=0$) and the lateral boundaries. The bottom boundary conditions are free slip for the velocity and a constant heat flux for the temperature corresponding to the vertical gradient in the thermocline. Sponge layers are employed at the lateral and bottom boundaries to prevent reflection of spurious waves. The sponge layers at the lateral boundaries have a thickness of $64\,\rm{m}$ and the sponge layer at the bottom boundary is $20\,\rm{m}$ thick. The governing equations (Eq. 2) are advanced in time using a mixed third-order Runge-Kutta (for advective fluxes) and Crank-Nicolson (for diffusive fluxes). Second-order finite difference discretization is used to compute spatial derivatives. The dynamic pressure is obtained by solving the Poisson equation with a multi-grid iterative method.

The evolution of the front is discussed in detail in VPS19. Here, we describe the evolution briefly to motivate the Lagrangian studies performed in the remainder of this paper. The front evolves through symmetric and baroclinic instabilities. Initially, the front is unstable to symmetric perturbations as the potential vorticity is negative (Hoskins, 1974). The symmetric instability (SI) grows and forms convection cells nearly aligned with the isopycnals. However, SI does not persist for long. The vertical shear in the convection cells become unstable to secondary Kelvin-Helmholtz instabilities (Taylor and Ferrari, 2009), which break down into turbulence through tertiary instabilities (Arobone and Sarkar, 2015). This leads to restratification of the front, making it stable to SI.

The evolution of the front, to a large extent, is controlled by baroclinic instability (BI), which becomes dominant after SI subsides. The non-linear growth of BI spawns submesoscale coherent structures such as vortex filaments and eddies. Understanding these coherent structures provides deeper insights into the dynamics at the front, specifically the vertical and lateral transport, as well as dispersion of tracer particles. Here, we briefly explain how the submesoscale coherent eddies and filaments evolve and how their spatial organization changes in time. The evolution is readily noticeable from the submesoscale flow component obtained by applying a 2D, low-pass Lanczos filter in $x$- and $y$-coordinate directions, with a cutoff wavenumber, $k_{c}=0.04\,\rm{rad\,m^{-1}}$, and wavelength, $\lambda_{c}\equiv 2\pi/k_{c}=157\,\rm{m}$. The application of the filter separates the coherent submesoscale from the small scales.

The time evolution of the coherent structures is illustrated in Fig. 2, which shows the submesoscale vertical vorticity at $10\,\rm{m}$ (Figs. 2a-c) and $30\,\rm{m}$ (Figs. 2d-f) depth at different times, $t=57.2\,\rm{h},75\,\rm{h},$ and $84.9\,\rm{h}$. The filaments of cyclonic vorticity connected to the heavy edge of the front and wrapping into coherent eddies in the central region can be observed in Fig. 2a plotted at $t=57.2\,\rm{h}$. At this time, the vorticity filaments have just begun to roll up, and eddies are small in size, slightly larger than the width of the filaments, $O(100)\,\rm{m}$. The eddies are vertically coherent and can be identified at $30\,\rm{m}$ depth in Fig. 2d. As the instabilities evolve, the vorticity filaments grow in length and the eddies become larger in diameter (Figs. 2b,c). Moreover, the structures are advected by the mean jet velocity, which is in the negative-$x$ direction in the present model. In Fig. 2a, three small eddies connected to the heavy-edge vortex filaments can be observed. However, Fig. 2b, plotted at $t=75\,\rm{h}$, shows two relatively larger eddies, whereas the one situated between them does not grow. This represents a merger of the two cyclonic eddies on the left side of the panel as well as amalgamation of surrounding cyclonic vorticity into the growing vortices. The vortex merger is nearly complete at $t=84.9\,\rm{h}$ where the imprint of the eddy, which was initially between the other two, is weak at both $10\,\rm{m}$ and $30\,\rm{m}$ depth (Figs. 2c,d). At depth, there are vorticity filaments attached to the light edge of the front that wrap around the eddies. The light-edge filaments can be observed in Figs. 2e,f plotted at $30\,\rm{m}$ depth, and they wrap around the eddies at the side opposite to where the heavy-edge filaments join with the eddies.

The three-dimensional organization of the coherent structures at a late time ($t=84.9\,\rm{h}$) is visualized in Fig. 3, where the iso-surfaces of Q are plotted. In this figure Q is calculated using the submesoscale velocity fields, i.e., $\tilde{Q}=(\overline{\Omega}_{ij}^{2}-\overline{S}_{ij}^{2})/2$, and the iso-surfaces are plotted at $\tilde{Q}/f^{2}=-0.4$ and $0.4$; here, $\overline{\Omega}_{ij}=(\partial\overline{u}_{i}/\partial x_{j}-\partial\overline{u}_{j}/\partial x_{i})/2$ and $\overline{S}_{ij}=(\partial\overline{u}_{i}/\partial x_{j}+\partial\overline{u}_{j}/\partial x_{i})/2$ are the second-order rotation and strain-rate tensors, respectively. Thus, regions with positive $\tilde{Q}$ represent rotation-dominated flow at the front, whereas those with negative $\tilde{Q}$ represent strain-dominated flow. The vertical coherence of the submesoscale coherent eddies is evident from the figure, as these structures dominated by cyclonic vorticity appear columnar. The vortex filaments connected to the heavy edge (on the negative-$y$ side) of the front can also be observed. The filaments are shallow at the end connected to the heavy edge; however, they grow deeper as one moves along their length towards the core of the eddies. Surrounding the eddies, the filaments appear to span the entire depth of the front. The filaments have regions which are dominated by strain as well as those where rotation dominates. In the neighborhood of the eddies, there are strain- and rotation-dominated vertical layers which are arranged alternately. We also note that the vortex filaments connected to the light-edge of the front are not obvious in the Q-visualization, possibly due to weaker cyclonic vorticity ($\sim f$) at depth.

In this paper, the main focus is to understand dispersion and transport by the submesoscale currents generated by BI. Such a focus necessitates studies of an essentially Lagrangian nature which require the knowledge of how material points, i.e., the tracer particles, move under the influence of submesoscale currents. To this end, tracer particles were introduced in the flow at $t=57.2\,\rm{h}$. At this time, the vortex filaments have formed at the front and have begun to wrap into eddies (Figs. 2a,d). The particles are placed at the nodes of a regular lattice over a rectangular subdomain that occupies the entire domain length in the along-front ($x$ direction), $-1.6\,\rm{km}\leq y\leq 1.6\,\rm{km}$ in the lateral ($y$ direction) and $-70\,\rm{m}\leq z\leq-2\,\rm{m}$ in the vertical (z direction) with resolution $16\,\rm{m}\times 16\,\rm{m}\times 2\,\rm{m}$. For multi-particle analysis, additional particles are released at $10\,\rm{m}$ and $30\,\rm{m}$ depth. Before introducing the particles, the simulation run with a uniform vertical grid of resolution $2\,\rm{m}$ is interpolated to the higher vertical resolution grid (0.5 m near the surface to 1.5 m near the bottom of the upper-ocean front) at $t\approx 56\,\rm{h}$. Tracer particles are seeded in the domain at $57.2\,\rm{h}$ at which point the simulation on the grid with high vertical resolution has progressed by about an hour.

The tracer particles are passive and, by definition, move with the local fluid velocity. Thus, the position of a tracer particle $\mathbf{x}_{p}=(x_{p},y_{p},z_{p})$ is expressed as

where $\mathbf{u}_{f}(\mathbf{x}_{p},t)$ is the fluid velocity at the particle’s position. The trajectories of the particles are computed by integrating Eq. 6. The time integration is performed following a third-order Runge-Kutta (RK3) scheme, and the particle velocity $\mathbf{u}_{f}(\mathbf{x}_{p},t)$ is obtained by the fourth-order Lagrange interpolation of a cell-centered velocity field. The Navier-Stokes solver stores the velocity components at the edge centers, and the cell-centered velocity is obtained by the linear interpolation. A CFL value smaller than one is ensured for the particle advection to get a stable numerical trajectory.

The trajectories followed by individual tracer particles moving in a time-varying velocity field at the front are complex and can differ considerably even for particles released close to each other. Nevertheless, the particles are strongly influenced by the coherent structures, which leads to an overall organization in their motions as elaborated below. Anticipating differences in the transport of particles in the frontal zone and at the edges, it will be convenient at times to distinguish among particle groups according to their cross-front ($y$) locations as follows: (i) central-region particles released in $-500\,\rm{m}\leq y\leq 500\,\rm{m}$, (ii) heavy-edge particles released in $y<-500\,\rm{m}$ and (iii) light-edge particles released in $y>500\,\rm{m}$.

The influence of coherent vortex filaments and eddies on the transport and organization of tracer particles at the front is illustrated in Figs. 4 and  5, using particles released at 10 m and 40 m depth. Figure 4 shows temperature, lateral velocity, and vertical velocity of the particles at $t=84.9\,\rm{h}$, which is $\sim 28\,\rm{h}$ after they were released. An examination of the particles reveals that they are organized into two lobes (LB1 and LB2) within the front, each associated with a coherent eddy. As the eddies move along the front, so do the lobes. We also notice that the lobes are stratified (Figs. 4a,b), with warm particles constituting the upper portions of the lobes, facing the lighter side of the front, and the cold particles constituting the undersides. Recall that the along-front velocity is vertically sheared, and the near-surface particles in the lobes have larger magnitude of along-front velocity ($u_{p}$) compared to those near the bottom, whose velocity magnitude is nearly zero. The along-front motions of the lobes therefore leads to a complex circulation of the constituent particles. The particles appear to circulate clockwise in the inclined lobes when viewed from the top, having negatively correlated lateral and vertical motions. This circulation has been illustrated by plotting lateral ($v_{p}$) and vertical ($w_{p}$) velocities of the particles: Figs. 4c,e for the particles released at $10\,\rm{m}$ depth and Figs. 4d,f for the particles released at $40\,\rm{m}$ depth. Negatively correlated $v_{p}$ and $w_{p}$ of the lobe particles can be observed in these figures. For example, particles at the front side LBF (marked by circles with dots in Fig. 4c) of the lobe LB2 have positive $v_{p}$. At the same time, the overall vertical motion in LBF is downwelling, as $w_{p}$ is mostly negative (Figs. 4e,f). These particles moving along the lobe surface in the positive-$y$ direction move downwards. Note that the back-to-front direction is down-front, i.e. oriented along the jet velocity, which is the negative-$x$ direction in the present configuration. Eventually, the lateral velocities become negative, when the particles reach to the back side LBB (marked by circles with crosses in Fig. 4c) of the lobe LB2. The overall vertical velocity becomes positive in LBB, as seen in Figs. 4e, and the particles climb up the lobe. In general, the particles circulate within the same lobe, especially when the interaction between the neighboring eddies is weak.

It is worth noting that the correlation between $v_{p}$ and $w_{p}$ is such that the associated circulation in the lateral-vertical plane follows the isopycnal slope. In the present case, the lateral density gradient points in the negative-$x$ direction and, therefore, $v_{p}$ and $w_{p}$ have negative correlation.

The coherent filaments transfer edge particles to the lobes. There are two light-edge filaments (LEF1 and LEF2) that can be identified in Fig. 4b. The filaments have mostly positive $w_{p}$ (Fig. 4f) and lift the warm-edge particles towards the surface and transfer them to the upper layers of the stratified lobes. Although not visible in the plots shown in Fig. 4, there are coherent filaments connected to the heavy edge of the front, as well. The role of the filament structures on the transport is further examined in Fig. 5 by plotting the heavy-edge particles released at $10\,\rm{m}$ depth (Fig. 5a) and the light-edge particles released at $40\,\rm{m}$ depth (Fig. 5b) at $t=84.9\,\rm{h}$, after a flight time of $\sim 28\,\rm{h}$. In Fig. 5a, the effect of two coherent filaments connected to the heavy edge of the front is apparent as downwelling of the heavy-edge particles (HEF1 and HEF2) to the underside of the lobes. Once within the lobes, the particles undergo motions that are characteristic of the lobe as described previously in this section. The light-edge filaments LEF1 and LEF2, also marked in Fig. 4b, can be observed in Fig. 5b. Comparing Figs. 5a,b shows that the light-edge filaments are located between the heavy-edge filaments. Additionally, we find that the particles moving through the light-edge filaments loop back to the bottom of the front after they are lifted upwards. The explanation is that the light-edge filaments transfer the particles to the lobe, upon which they subduct through the downwelling limb of the lobe. Interestingly, some particles in LEF2 detach from the main branch near the surface and spread laterally under the influence of the near-surface circulation. The particles detaching from the main branch near the surface are enclosed within the rectangular box shown in Fig. 5b.

The transport processes mediated through filaments and eddies and the circulation of fluid particles organized into lobes, as discussed in Figs. 4 and Fig. 5, are further illustrated by the schematic shown in Fig. 6. The schematic depicts the overall motions through the heavy- and light-edge filaments, which transfer fluid particles from the edges to the lobe structures in the central region of the front; the circulations within the lobes are also shown.

The flux of edge particles into the central region causes the local particles to adjust, leading to the restratification of the front. In order to understand the subsequent adjustment of the front, it is key to characterize how the particles that were released in the central region redistribute spatially as time progresses. Figure 7 shows the time evolution of the particles released in the central region, with $-500\,\rm{m}<y<500\,\rm{m}$ and $z>-50\,\rm{m}$. For identifying their distribution, the particles are sampled in the cells of a rectangular grid that has lateral resolution $\Delta_{s}y=16\,\rm{m}$ and vertical resolution $\Delta_{s}z=2\,\rm{m}$. The distributions are plotted at three different times: $t=69\,\rm{h},84.9\,\rm{h}$, and $99\,\rm{h}$, about $12\,\rm{h},28\,\rm{h}$ and $40\,\rm{h}$ after their release. As the front evolves through BI, the particles contained in the region $-500\,\rm{m}<y<500\,\rm{m}$ and $z>-50\,\rm{m}$ spread laterally while remaining confined in the central region, indicating slumping of the front. The number of particles in the cells within the particle cloud away from the edges does not change considerably with time, consistent with the incompressibility of the flow. Thus, the edge particles that are brought into the front are primarily organized in regions above or below the central-region particles. The isotherms corresponding to the mean temperature of the sampled particles are also plotted. They reveal that the stratification is maintained during the lateral spread of the particles. Further, the stratification becomes stronger with increasing time.

The aforementioned features of the Lagrangian transport can be identified in individual particle trajectories. In Fig. 8, the trajectories $x_{p}(t)$, $y_{p}(t)$ and $z_{p}(t)$ of a few selected particles are shown; the trajectories correspond to the particles released at different lateral ($y$-direction) locations with fixed $z=-30\,\rm{m}$ and $x=1490\,\rm{m}$. Depending on their initial $y$-coordinates (marked in the middle panel of Fig. 8), the particles can be distinguished as the heavy-edge particles (P1 and P2), the central-region particles (P3-P5), and the light-edge particles (P6 and P7). Typically, the heavy-edge particles downwell and the light-edge particles upwell through filaments as they are transported to the central region of the front, where they circulate with the local particles organized into lobes. The lobes are associated with coherent eddies and move in the along-front direction.

The vertical trajectories of the particles are depicted in Fig. 8c. As expected, the figure shows that the heavy-edge particles (P1 and P2) downwell, while the light-edge particles (P6 and P7) upwell. The trajectory of P2 also shows upwelling after $t\approx 100\,\rm{h}$, which is a result of the particle’s motion within the lobe. The negative correlation between the lateral and vertical motions of the particles is also evident from some of the lateral and vertical trajectories (Figs. 8b and c). Such correlation often occurs for particles moving in the lobes or through the coherent filaments. For example, the $y$ and $z$ trajectories of the central-region particle P4 reveals that the particle moves vertically downward during the time when the lateral motion is in positive $y$ direction, whereas it moves vertically upward when the lateral motion is in negative $y$ direction. Moreover, the correlated lateral and vertical motions exhibit oscillations with a time period of about $25\,\rm{h}$, which is twice the inertial time period ($T=12.5\,\rm{h}$). Similarly, the upwelling/downwelling edge particles (e.g. P1 and P7) exhibit negatively correlated lateral and vertical motions. The central-region particle P3 remains trapped inside an eddy and shows oscillations in its $y$ coordinates at near-inertial time scale while maintaining a nearly constant height in the vertical. The decoupled lateral and vertical motions are also evident for the light-edge particle P6, after it upwells to the surface ($t\gtrapprox 90\,\rm{h}$). The vertical trajectory of P6 also shows a fast time-scale event with remarkably rapid transport in the vertical. This event starts at $t\approx 97\,\rm{h}$ when P6 downwells by approximately $20\,\rm{m}$ (from A to B) over a period of about an hour and then upwells back to the surface (from B to C) in the next two and half hours. The downwelling occurs when the particle gets attracted to a heavy-edge filament, which mainly transports cold fluid downwards. However, P6 eventually upwells when it finds itself in a denser background.

The trajectories of the edge particles reveal that their vertical transport (under the influence of the front) commences at different times. In general, the motion of an edge particle farther away from the central region of the front is delayed compared to the one that is closer. In Fig. 8c, the approximate time when the heavy- and the light-edge particles start moving vertically are marked with solid squares in their trajectories. The starting time is determined when the magnitude of vertical displacement exceeds $2\,\rm{m}$ for a heavy-edge particle and $4\,\rm{m}$ for a light-edge particle; a higher threshold for the vertical displacement is used for the light-edge particles as their vertical trajectories have relatively large amplitude oscillations superposed to their initial positions. Indeed, P1 starts moving vertically after P2 and P6 after P5, since P2 and P5 are closer to the central region of the front than P1 and P6. This suggests that the vortex filaments primarily transport the edge particles adjacent to the slumping front, and those outside are transported after the width of the front increases slowly in the lateral.

The along-front particle trajectories are shown in Fig. 8a. The figure shows that the displacements are generally in the negative-$x$ direction, same as the mean along-front velocity at the front. It is worth noting that the $x$ trajectories of particles P5 and P6 cross the boundary of the computational domain at $x=0$. The trajectories are continued into the negative $x$ region using the streamwise periodicity of the domain. It can also be noticed that the upwelling particles (P6 and P7) have larger negative $x$ displacements compared to the downwelling particles (P1 and P2), as the upwelling particles tend to spend more time near the surface where the along-front velocity is larger. Overall, the displacement of the particles in the negative $x$ direction increases with time, indicating negative along-front velocities; however, a few particles can acquire positive along-front velocities, especially when they are near the bottom of the front (e.g. P2 during $t\approx 96-105\,\rm{h}$ and P4 during $t\approx 64-72\,\rm{h}$)

Oscillations with small amplitudes can be observed in the vertical trajectories (e.g. P4, P5, and P6), indicating the influence of the finescales. Further, the effect of the finescales on the along-front and cross-front trajectories is weak, as the trajectories appear smooth (i.e. the deviations are small) and are dominated by the submesoscale velocity components. Next, we quantify the correlation between the lateral and vertical motions of the particles.

The particle trajectories discussed in this section have demonstrated that particles circulating in the lobes or moving through vortex filaments exhibit a trend towards negative correlation of lateral- and vertical-velocity components, a consequence of the secondary circulation induced by the coherent structures. To statistically quantify this relationship between vertical and lateral motions in the central region of the front, a correlator variable $r_{yz}$ is introduced for each particle trajectory and probability density functions (PDFs) are computed over the ensemble of particles released at specific depths. For each particle trajectory $r_{yz}$ is defined as

Here, $\Delta y_{p}^{n}=y_{p}(t_{n})-y_{p}(t_{n-1})$ and $\Delta z_{p}^{n}=z_{p}(t_{n})-z_{p}(t_{n-1})$ for a particle, and the time superscript $n$ varies to cover the entire simulation time from $t_{0}$ to $t_{N}$. The above definition of the correlator $r_{yz}$ can also be interpreted in terms of the particle velocity weighted with the advection time step, i.e., $\Delta y_{p}^{n}\approx v_{p}^{n}\Delta t_{n}$ and $\Delta z_{p}^{n}\approx w_{p}^{n}\Delta t_{n}$, with $\Delta t_{n}=t_{n}-t_{n-1}$. The weighting with $\Delta t_{n}$ used with the velocity components $v_{p}^{n}$ and $w_{p}^{n}$ accounts for the variable time step in the simulation.

The correlator takes a value of $r_{yz}\in[-1,1]$, a magnitude close to unity represent perfectly correlated vertical and lateral motions and magnitudes close to zero correspond to uncorrelated motions. The probability density function (PDF) of the correlator $r_{yz}$ is computed for groups of particles released in the central region at different depths – $10\,\rm{m}$, $20\,\rm{m}$, $30\,\rm{m}$ and $40\,\rm{m}$ – and shown in Fig. 9 are their PDFs. The particles quickly organize into coherent lobes after being released and and continue moving in these structures thereafter. The figure clearly shows negative $r_{yz}$ for the majority of the particles. The PDF of $r_{yz}$ for the particles released at $10\,\rm{m}$ have a broad peak. The peak sharpens and shifts towards -1 as the depth of the release increases, indicating stronger correlation between the lateral and vertical motions. For the particles released at $40\,\rm{m}$ depth the median value of $r_{yz}$ is about -0.6. There are also many particles with small and moderate values of $r_{yz}$, reflecting some unpredictability in the multiscale, chaotic trajectories executed by the particles as they circulate within the front.

The correlator $r_{yz}$ as defined in Eq. 7 is skewed towards the largest magnitudes of $\Delta y_{p}^{n}\Delta z_{p}^{n}$. Alternatively, we can define a correlator $\tilde{r}_{yz}$ that gives equal weight to the lateral and vertical displacements at each time step, i.e.,

where $|\Delta y^{n}_{p}|$ and $|\Delta z^{n}_{p}|$ are the absolute values of $\Delta y^{n}_{p}$ and $\Delta z^{n}_{p}$, respectively. The modified correlator $\tilde{r}_{yz}\in[-1,1]$ and contains similar information as $r_{yz}$. The PDFs of $\tilde{r}_{yz}$ for the particles released in the central region at different depths are qualitatively similar to those shown in Fig. 9 and are not shown. This suggests that the negative correlation between the lateral and vertical motions is not episodic, dominated by large displacements over a few time steps; instead, it is a typical feature of the particle motion within the front at all times.

Thus, baroclinic instability at the front leads to complex Lagrangian dynamics. The paths followed by each tracer particle vary from one another locally, as well as globally over different regions of the front. Nevertheless, an overall underlying structure can be constructed based on the collective motion of the particles. The tracer particles in the central region organize into lobes, each associated with a coherent submesoscale eddy. The particles in a lobe are stratified with colder particles located at the underside and warmer particles at the top. As the lobes move with eddies, particles in the lobes circulate clockwise when viewed from above, and the sense of rotation is opposite to the cyclonic coherent eddies. On average, the lateral and vertical motions are negatively correlated. Typically, a particle circulating within a lobe moves downwards when the lateral velocity is positive, but moves upwards when it becomes negative. The coherent filaments from the light/heavy side of the front connect the edges with the lobes in the central region and transfer warm/cold edge particles to the central region. The newly deposited particles subsequently undergo the typical circulation in the lobes.

In the previous section, we have shown that the vertical motions of particles at the front exhibit oscillatory components at two widely separated time scales: fast oscillations due to the small scales in the vortex filaments and slow oscillations at a near-inertial time scale due to the circulation in the lobes.Therefore, large $w_{p}$ magnitudes, such as those encountered in the vortex filaments, do not necessarily lead to a net vertical transport, responsible for restratifying the front; the long-time displacements must be examined. Following a single particle over a long time is insufficient since the behavior can differ considerably from one particle to another. In this section, we investigate the collective motion of particle clouds and inquire about the relevant time scale for subduction and restratification at the front. The clouds are created such that the constituent fluid particles have similar densities and, therefore, have similar buoyancy control on the dynamics.

In this section, single- and two-particle dispersion statistics, as well as multiparticle statistics using a group of four particles (tetrads) are studied. All the results included in this section consider only those particles released in the central region.

We investigate dispersion and transport by submesoscale turbulent currents generated by the evolution of baroclinic instability at an upper-ocean front. The study employs a LES model and is performed in the Lagrangian framework by releasing a large number of tracer particles that move with the local fluid velocity. The presence of coherent structures such as vortex filaments and eddies is a typical feature of submesoscale dynamics. From the Lagrangian analysis, we find that these structures provide the primary pathways of three-dimensional transport by submesoscale currents and provide a quantitative assessment of the transport.

The paths followed by individual particles are found to be complex and can differ considerably as time progresses, even for a pair released close to each other. Nevertheless, the motions are strongly influenced by the coherent structures, namely the vortex filaments and eddies. Particles inside the filaments experience rapid motions with displacements of $O(10)\,\rm{m}$ over an hour, as well as slower motions at a near-inertial time scale while moving under the influence of the coherent structures. It is possible to identify some common features that dictate the overall transport. The central-region particles cluster into inclined lobes, each associated with a coherent eddy. The lateral and vertical velocity of these particles reveals a clockwise circulation when viewed from above, which is opposite to the circulation induced by the coherent cyclonic eddies. The vortex filaments connect the heavy and light edges of the front with the central region and play a critical role in vertical and lateral transport and the restratification of the front. The process can be described as follows. The coherent filaments draw the edge particles into the central region and transfer them to the lobes. The lobes are stratified, and the heavy-edge particles downwell to the undersides of the lobes, whereas the light-edge particle upwell to the top. The flux of new edge particles into the central region from the edges causes the central-region particles to adjust, which leads the front to restratify.

We find that the lateral ($v_{p}$) and vertical ($w_{p}$) velocity of the particles moving through the filaments and circulating in the lobes have a near-inertial time scale of $O(2\pi/f)$ and have a correlation which is consistent with the lateral stratification. For the present case where the cross-front density gradient is negative, the correlation is also negative. The correlation between $v_{p}$ and $w_{p}$ is quantified by defining a correlator for each particle using its lateral and vertical displacements over small time intervals and computing the accumulated value over the entire flight time of the particle. The median value of the correlator determined for the central region particles is $\sim-0.5$; the value depends on the depth where the particles are released with the correlation being somewhat larger for particles released near the bottom. It is important to note that the large magnitudes of vertical velocity in the filaments or lobes do not independently lead to a net restratification of the front. The restratification process is also dependent on the transport of particles from the edges to the central region with an appropriate correlation between $v_{p}$ and $w_{p}$.

Further analysis by following the centers of mass of the particle clouds $z_{com}$ released at $10\,\rm{m}$ and $40\,\rm{m}$ depth shows that subduction/upwelling through the vortex filaments occurs over 1-2 inertial time periods, and a slow adjustment follows after the particles are accommodated in the lobes and begin circulating. The near-inertial time scale is consistent with the time scale of the growth of baroclinic instability, which drives the restratification of the front. During the subduction/upwelling through the filaments, the particles disperse, mostly along the sloping isopycnals, and the the root-mean-square vertical displacement ($z_{rms}$) of the constituent particles with respect to the COM saturates to $\sim 15\,\rm{m}$. Near-inertial oscillations in $z_{com}$ and $z_{rms}$ of the particle clouds are observed to result from the circulation of the particles in the lobes.

Fluid and flow properties associated with material points are also transported by the submesoscale currents. The mean subgrid viscosity of the particle clouds released at $10\,\rm{m}$ and $40\,\rm{m}$ depth reveals large magnitudes initially, reflecting the motions of the particles through the vortex filaments. Moreover, the clouds released near the surface experience about 2-3 times larger values of subgrid viscosity compared to those released near the bottom, as the finescale activity near the surface is stronger due to surface-intensified frontogenesis. The mean values typically decrease for both groups of particle clouds as time progresses. The possible exceptions are the upwelling particle clouds, which show elevated mean subgrid viscosity at late times when they reach the near-surface region. The average subgrid viscosity at late time is about $50-100$ times larger than the molecular viscosity for both $10\,\rm{m}$-depth and $40\,\rm{m}$-depth clouds. Because of the turbulent exchange with the surroundings, the fluid properties associated with the particles change. We find that the change in the average temperature of the clouds over a flight time of about 45 h is about $4-6\%$ of the cross-front temperature difference. Typically, the downwelling particles on average tend to become warmer while the upwelling particles tend to become colder. The mean KE of the clouds change with $z_{com}$, which reflects a more energetic flow near the surface than at depth.

The process of restratification in the present model front under BI is considerably different than the restratification process depicted in Spall (1995), which assumes sliding of a fluid parcel from the heavy side of the front across to the light side. In contrast, the process described here is three-dimensional and involves continuous stirring of the central-region fluid by submesoscale coherent eddies and the injection of edge particles into the central region by the coherent vortex filaments. Thus, after being subducted, a heavy-edge fluid parcel continues to move under the influence of the eddies at the front.

We also find that vertical distribution of the particles released at a depth remains confined within $\sim 50\,\rm{m}$ depth from the surface over time (see Fig. 12), which is also the initial depth of the front. As shown in VPS19, the vertical velocity in the thermocline is non-zero, but the particles do not subduct below the surface layer. This behavior is consistent with the fact that coherent structures control the vertical transport of the particles at the front. Since these structures are contained within the front, so are the particles.

The near-inertial oscillations observed here can be contrasted with the inertial oscillations associated with the geostrophic adjustment of a front with an initially unbalanced horizontal density gradient, which was analyzed by Tandon and Garrett (1994). In the present simulation, the near-inertial oscillations result primarily from the anticyclonic circulation of the particles within the lobes. Furthermore, the dynamics here are driven by BI. In contrast, the model investigated by Tandon and Garrett (1994) exhibits inertial oscillations when the unbalanced front tries to slump by releasing potential energy, but the Coriolis force acts on the developed velocity field and provides a restoring tendency towards the original configuration. Tandon and Garrett (1994) find that the oscillatory adjustment continues at inertial time scale and, due to the lack of dissipation in their simplified model, the system does not return to a stable stationary state.

The dispersion characteristics of the submesoscale turbulent flow are also studied here. In both single- and pair-particle dispersion, the vertical component is restricted by the mixed layer depth and its value saturates to $O(10)\,\rm{m}$ at long times. The along-front component of single-particle dispersion shows super-diffusive behavior at late times, and the mean-square displacement increases as $t^{1.8}$; this behavior can be related to the mean jet in the negative $x$ direction. In the ocean, long-time super-diffusive behavior has been observed in coastal regions with mean currents. The particle-pair dispersion in $x$ and $y$ directions show $t^{3}$ behavior during the intermediate times, and the root-mean-square displacement is $O(100)\,\rm{m}$, which is comparable to the lateral width of the vortex filaments. This may indicate a Kolmogorovian inertial range, but the role of horizontal shear on relative dispersion cannot be ruled out. The long-time behavior of particle pairs is consistent with single-particle dispersion. The multiparticle analysis reveals strong filamentogenesis in the vortex filaments, as the tetrads moving through these structures deform into thin, needle-like shapes. Probability density functions of shape metrics $I_{1}$ and $I_{2}$ indicate that there is a strong propensity to form needle-shaped structures, more so than in homogeneous turbulence that is either isotropic or stratified. The filamentogenesis is associated with the strong strain field within the coherent filaments.

We are pleased to acknowledge the support of ONR grant N00014-18-1-2137. We thank Hieu Pham for helpful discussion and his comments on a draft version of this manuscript.