## Abstract

Oceanic surface submesoscale currents are characterized by anisotropic fronts and filaments with widths from 100 m to a few kilometers; an *O*(1) Rossby number; and large magnitudes of lateral buoyancy and velocity gradients, cyclonic vorticity, and convergence. We derive an asymptotic model of submeoscale frontogenesis—the rate of sharpening of submesoscale gradients—and show that in contrast with “classical” deformation frontogenesis, the near-surface convergent motions, which are associated with the ageostrophic secondary circulation, determine the gradient sharpening rates. Analytical solutions for the inviscid Lagrangian evolution of the gradient fields in the proposed asymptotic regime are provided, and emphasize the importance of ageostrophic motions in governing frontal evolution. These analytical solutions are further used to derive a scaling relation for the vertical buoyancy fluxes that accompany the gradient sharpening process. Realistic numerical simulations and drifter observations in the northern Gulf of Mexico during winter confirm the applicability of the asymptotic model to strong frontogenesis. Careful analysis of the numerical simulations and field measurements demonstrates that a subtle balance between boundary layer turbulence, pressure, and Coriolis effects (e.g., turbulent thermal wind; Gula et al. 2014) leads to the generation of the surface convergent motions that drive frontogenesis in this region. Because the asymptotic model makes no assumptions about the physical mechanisms that initiate the convergent frontogenetic motions, it is generic for submesoscale frontogenesis of *O*(1) Rossby number flows.

## 1. Introduction

Surface oceanic submesoscale currents consist of coherent fronts, filaments, and eddies with spatial scales on the order of 0.1–10 km, and time scales on the order of hours to days, that are dynamically characterized by *O*(1) Rossby and Richardson numbers (Thomas et al. 2008; McWilliams 2016). Due to their significant ageostrophic vertical velocities, they are thought to be instrumental in transferring properties from the surface of the ocean to the interior, thereby affecting a wide range of physical and biogeochemical processes (Thomas et al. 2008; McWilliams 2016; Mahadevan 2016). In addition, they are believed to influence the dispersion of contaminants (Poje et al. 2014) and to provide an important route for kinetic energy dissipation (Capet et al. 2008c; Molemaker et al. 2010; Barkan et al. 2015), thus impacting the general circulation of the oceans.

Submesoscale fronts and filaments are identified as anisotropic flow features with large magnitudes of vorticity, horizontal divergence, buoyancy gradients, and velocity gradients. Their generation is explained by a variety of processes, including mixed layer instabilities (Boccaletti et al. 2007; Fox-Kemper et al. 2008), deformation frontogenesis (Hoskins and Bretherton 1972), and boundary layer turbulence (Gula et al. 2014; McWilliams et al. 2015; Sullivan and McWilliams 2018) all of which rely on the conversion of available potential energy (APE) to kinetic energy (McWilliams 2016). Because of the importance of local advection to their dynamical evolution; their plausible interaction with boundary layer turbulent processes (Bachman and Taylor 2016; McWilliams 2017), surface waves (McWilliams and Fox-Kemper 2013; Hamlington et al. 2014; Suzuki et al. 2016), and internal waves (Whitt and Thomas 2015; Barkan et al. 2017c; Thomas 2017); and their small temporal and spatial scales, they are challenging to study theoretically, numerically, and observationally.

Hoskins (1982) studied the deformation frontogenesis mechanism in an atmospheric context by evaluating the evolution equation for the horizontal buoyancy gradient squared

where $D/Dt$ is the material derivative associated with the three-dimensional velocity field $\u2061(u,\upsilon ,w)$; $j,k=1,2$; $N2$ is the stratification; summation over repeated indices is assumed; commas denote derivatives; and $Vmixb$ and $Hdiffb$ denote viscous forcing and dissipation processes associated with vertical mixing and horizontal diffusion, respectively. Flow regions with positive $Fhorb$ correspond to the sharpening of lateral buoyancy gradients, or frontogenesis (Hoskins 1982; Capet et al. 2008b), where it is traditionally assumed that the horizontal velocity leading to the sharpening is associated with a balanced straining field (Hoskins and Bretherton 1972; Lapeyre and Klein 2006; McWilliams et al. 2009). In quasigeostrophic theory, the horizontal velocity in $Fhorb$ is in geostrophic balance and, consequently, the rate of lateral buoyancy gradient sharpening in the inviscid limit is exponential. In semigeostrophic theory (Hoskins 1975) the horizontal velocity in $Fhorb$ has, in addition, an ageostrophic component that is associated with the induced secondary circulation around the front. This secondary circulation is convergent near the surface and therefore enhances the rate of lateral buoyancy gradient sharpening, leading to a finite time singularity in the inviscid limit (Hoskins and Bretherton 1972; Hoskins 1982). It has been argued that to adequately quantify the rate of frontal sharpening and subsequent vertical fluxes the ageostrophic contribution must be accounted for (Thomas et al. 2008). The $Fvertb$ term, during frontogenesis, is associated with the restratification that accompanies the lateral sharpening process.

In the ocean a variety of mechanisms, other than classical deformation frontogenesis, can produce an ageostrophic secondary circulation in the mixed layer, and strong lateral buoyancy gradients may be driven by a combination of balanced and unbalanced motions (Srinivasan et al. 2017). Often, the strong lateral buoyancy gradients at fronts and filaments are accompanied by strong velocity gradients, large cyclonic vorticity, and large convergence values, as illustrated in Fig. 1. Therefore, to understand the underlying physics of frontogenesis it is instructive to also examine the corresponding evolution equations:

Above, $\zeta =\u2212\epsilon ijui,j$ denotes the vertical vorticity component; $\delta =ui,i$ is the horizontal divergence; $\Phi $ is the geopotential; *f* is the Coriolis frequency; $i,j,k=1,2$; $\Lambda i=\u2061(uz,\upsilon z)$ denotes the vertical shear components; $\epsilon ij$ is the Levi-Civita symbol; commas denote derivatives; and summation over repeated indices is assumed. For all fields in Eqs. (2)–(4) $Vmix$ and $Hdiff$ denote viscous forcing and dissipation processes associated with vertical mixing and horizontal diffusion, respectively. The interpretation of $Fhoru,Fhor\zeta $, and $Fhor\delta $ is equivalent to that of $Fb$, where flow regions with positive $Fhoru$ and $Fhor\zeta $ correspond to the enhancement of lateral velocity gradients and of cyclonic vorticity, respectively. Similarly, flow regions with negative $Fhor\delta $ correspond to the reduction in positive divergence and the enhancement of convergent motions. Throughout this manuscript we will refer to the simultaneous increase of all gradient fields in Eqs. (1)–(4) generically as submesoscale frontogenesis, where it is implied that classical deformation frontogenesis (Hoskins and Bretherton 1972) is only one possible triggering mechanism. Irrespective of the motions driving submesoscale frontogenesis, the resulting flow structures—whether fronts, filaments, or eddies—have a preferred orientation and hence are anisotropic (Fig. 1). The ageostrophic secondary circulation that develops around the flow structures releases APE and acts to restratify the mixed layer, and therefore the $Fvert$ terms will oppose the gradient sharpening induced by the $Fhor$ terms. As long as the $Fhor$ terms remain larger than the corresponding $Fvert$ terms, frontogenesis will persist until, eventually, additional instabilities and dissipative three-dimensional turbulent processes will arrest it and prevent the hypothetical finite-time singularity, a process typically referred to as frontal arrest (McWilliams and Molemaker 2011; Sullivan and McWilliams 2018).

In the current paper we propose an asymptotic model that highlights the role of the horizontal convergence induced by the secondary circulation in governing submesoscale frontogenesis (section 2). Under the proposed model, analytical solutions for the inviscid evolution of all gradient fields in a Lagrangian reference frame can be obtained. The validity of the asymptotic model is verified with respect to submesoscale-resolving numerical simulations in the northern Gulf of Mexico (section 4) and drifter observations (section 5) collected as part of the Lagrangian Submesoscale Experiment (LASER) during January–February 2016 (details of the numerical setup and data analysis are provided in section 3). The asymptotic model is then applied to estimate vertical buoyancy fluxes during frontogenesis (section 6), and a summary and discussion of the results is given in section 7.

## 2. The asymptotic model

Oceanic flows with energetic submesoscale currents have velocity and buoyancy power spectral densities $E\u02dc[uh\u2061(kh)]$ and $E\u02dc[b\u2061(kh)]$ with spectral slopes of $kh\u2212\beta $, where $\beta ~2$ (Capet et al. 2008a). To link these spectral slopes to the observed anisotropic coherent structures (Fig. 1), it is insightful to examine the corresponding gradient spectra:

where $kh$ denotes the horizontal wavenumber vector, $E\u02dc$ denotes the power spectral density (PSD), and when $\beta ~2$ the spectral slopes of $E\u02dc[\u2207huh\u2061(kh)]$ and $E\u02dc[\u2207hb\u2061(kh)]$ are nearly white.

If the white spectral slopes of the gradient fields are assumed to be associated with randomly spaced Dirac delta functions in physical space, then the fields themselves will be associated with Heaviside functions. These one- or two-sided Heaviside functions are a simple functional representation of anisotropic fronts or filaments, respectively.^{1} An important spectral identity relates the PSD of the velocity gradient with the PSDs of the vorticity and the divergence fields, namely,

This spectral identity reveals how, in general, the variances of both vorticity and divergence contribute to the velocity gradient variance.^{2} Indeed, the velocity PSD and the corresponding PSDs of the velocity gradient, vorticity, and divergence in our submesoscale-resolving numerical simulations (Fig. 2) illustrate the importance of the divergence field at submesoscales (i.e., horizontal scales smaller than 10 km), even though the velocity PSD is somewhat steeper than $kh\u22122$ (Fig. 2a). Specifically, the ratio of the vorticity PSD to the divergence PSD is about six at *k*_{h} = 10^{−4} m^{−1} and gradually decreases at higher wavenumbers, suggesting that the magnitudes of *ζ* and *δ* are asymptotically similar at these (submeso) scales.

The phenomenological structure of the submesoscale currents (Fig. 1) and the similar magnitudes of vorticity and divergence (Fig. 2) guide us in developing an asymptotic model for submesoscale frontogenesis. We choose the following scaling, which emphasizes the anisotropic structure of submesoscale currents within a surface mixed layer of depth $hml$ and local Coriolis frequency *f*,

Without loss of generality we set *u* and *υ* to be the cross-front and alongfront velocities, which correspond to the alongfront length scale *L* in the *y* direction and the cross-front length scale *l* in the *x* direction. Variable *T* is the time scale, and $b=\u2212g\rho /\rho 0$ is the buoyancy scale associated with the hydrostatic geopotential $\Phi $. The vertical stratification $N2$ is set to be smaller than $b/hml$ by a parameter $\eta \u226a\u20091$ to highlight the well-mixed structure of the surface mixed layer. For the same reason, the vertical shears $uz,\upsilon z$ are set to be smaller than $RoV/hml$ and $V/hml$, respectively, by a parameter $\xi \u226a1$. Nondimensional parameters are

Here Ro is the Rossby number, *ε* is the anisotropy ratio, and *λ* is the aspect ratio. We further choose the distinguished limit^{3}

The main novelty of the proposed scaling is that for $Ro~O\u2061(1)$ (i.e., submesoscale circulations) the leading-order cross-front and alongfront velocities have similar magnitudes, meaning that the cross-front velocity *u* is purely ageostrophic ( appendix A). Furthermore, $\zeta ~\delta $, and the velocity gradient magnitude squared

has the same form as the spectral identity in Eq. (6). As will be shown in the following sections, this is more appropriate for submesoscale frontogenesis than the semigeostrophic scaling of Hoskins and Bretherton (1972), where the cross-front velocity is much smaller than the alongfront (geostrophic) velocity, and where $\zeta \u226b\delta $ and $|\u2207huh|2~\zeta 2$.

The nondimensional gradient fields subject to the scaling [Eq. (7)] are

where subscripts denote derivatives. Under the assumption $\lambda \u226a1,\epsilon \u226a1,Ro~O\u2061(1)$, which is valid for hydrostatic, anisotropic, submesoscale currents, the leading-order gradient equations [Eqs. (1)–(4)] simplify to ( appendixes A and B)

where the simplification in Eqs. (12a)–(12d) makes the use of index notation unnecessary, $\u2207$ is the gradient operator, the *h* subscripts denote horizontal vector components, and the equations are written in dimensional form. The Laplacian of the pressure $\u2207h2\Phi $ is written as $f\zeta g$ in Eqs. (12b) and (12d), with the *g* and ag subscripts denoting the geostrophic and ageostrophic vorticity components, respectively. Note that *ζ* × Eq. (12c) + *δ* × Eq. (12d) = Eq. (12b), in agreement with Eq. (10).

### a. Frontogenetic tendencies and rates

The most striking simplifications in the asymptotic model, in comparison to the original equation set [Eqs. (1)–(4)], are in the frontogenetic tendency terms $Fhorb,Fhoru,Fhor\zeta ,and\u2009Fhor\delta $, and we therefore begin by examining the physical insight that can be gained from their simplified form. If we define the corresponding advective frontogenetic rates as

and

we obtain the nontrivial relationship

The above equality is a major result of this study that shows that the advective tendency rates of all fields are the same and equal the horizontal convergence $\u2212\delta $. Subsequent sections will test its applicability in models and data.

It is insightful to contrast the simple form of the various frontogenetic tendencies (and rates) above, which are applicable for strong anisotropic submesoscale flows, with the dynamical regime $\lambda \u226a1,Ro~\epsilon \u226a1$, which describes the initiation of frontogenesis. In this dynamical regime ( appendix B) the leading-order advective frontogenetic tendencies $Fhoru$ and $Fhor\zeta $ remain the same as in Eqs. (12b) and (12c), that is,

and

but

and

The nondimensional tensors $Skj$ and $b,kb,j$ are written in Eq. (A4), where the strain-rate tensor $Skj$ contains contributions from both normal and shear strains ( appendix A). This illustrates that at the initial stage of deformation frontogenesis buoyancy gradient sharpening may be driven by both confluent and convergent flow components. At later stages, as Ro increases, the convergent component gradually becomes the dominant one in $Fhorb$ (the convergent normal strains dominate the shear strains), as shown by Eq. (12a), even if the frontogenesis is initiated by a horizontally nondivergent deformation flow with small Ro. On the other hand, the advective rate of velocity gradient variance is still governed by the convergence rate $\u2212\delta $, even in the small Ro limit, as long as the anisotropic assumption is made.^{4}

### b. Analytical solutions in a Lagrangian reference frame

Equations (12a)–(12d) can be solved in the inviscid limit as long as we adopt a Lagrangian reference frame, that is, $D/Dt\u2192\u2202/\u2202\tau $. To gain insight into the nature of the equations, we first investigate solutions while ignoring the Coriolis and pressure terms. Such solutions may apply to situations where the vertical mixing terms largely balance out the Coriolis and pressure terms, a balance previously referred to as turbulent thermal wind (TTW; Gula et al. 2014; McWilliams et al. 2015; Sullivan and McWilliams 2018). In this case

where the *τ* notation is to remind the reader that the solutions are in a Lagrangian reference frame. The 0 subscript denotes the initial values of the different fields at time $\tau 0$, when the strong anisotropic flow regime is reached [i.e., when $\epsilon \u226a1,Ro~O\u2061(1)$]. The solutions above demonstrate that for any initially convergent flow a finite time singularity is expected at $\tau =\tau 0\u2212\delta 0\u22121$. Because $\delta 0~f$, this blowup time is faster than the semigeostrophic blowup time of $log2/\alpha $ (Hoskins and Bretherton 1972; Hoskins 1982), where $\alpha \u226af$ is a horizontally nondivergent strain rate.

Next, we will examine solutions to the full asymptotic inviscid equations, keeping the Coriolis and pressure terms in Eqs. (12b)–(12d). To make progress we have to assume a functional form for the unknown ageostrophic vorticity $\zeta ag$ and for simplicity we set $\zeta ag=\xb1a2f$, where *a* is assumed to be constant. It can be shown that any cyclonic ageostrophic vorticity will prevent a finite time singularity from occurring,^{5} and we therefore pick

We can now readily solve for all gradient fields, which become

where $C=arctan[\u2212\delta 0\u2061(af)\u22121]$, $af=|\zeta ag|f$, and Eq. (18b) is repeated [see Eq. (10)] to remind the reader that knowledge of the divergence and vorticity is sufficient to determine the Lagrangian evolution of $|\u2207huh|2$ . As for the case without rotation and pressure terms, any initially convergent flow will develop a finite time singularity with a modified blow up time $\tau =\tau 0+(\pi /2\u2212C)\u2061(af)\u22121$. In appendix C we generalize the above solutions for a time-varying ageostrophic anticyclonic vorticity and, for illustration, show the solution to

The analytical solutions for the divergence field [Eqs. (16d) and (18d) and appendix C] are shown in Fig. 3 (similar Lagrangian evolutions are found for the other gradient fields; not shown). These solutions illustrate how sensitive submesoscale frontogenesis is to even small values of ageostrophic vorticity (black and reds curve compared with the blue curve), which can substantially speed up the gradient sharpening. In all cases, the Lagrangian evolution and inviscid blowup times are independent of the physical mechanisms that initiate submesoscale frontogenesis. This implies that these analytical solutions [Eqs. (16d) and (18d) and appendix C] are generic with respect to the gradient evolution during submesoscale frontogenesis, as long as the scaling assumptions [Eq. (7)] are valid. In the following sections, we test and validate the asymptotic model [Eqs. (12a)–(12d)] against submesoscale resolving numerical simulations and drifter measurements in the northern Gulf of Mexico.

## 3. Models and measurements

In this section we describe the numerical setup and drifter measurements we used to validate the asymptotic model. The results follow in sections 4 and 5.

### a. Model setup

The numerical simulations used in this study are carried out with the Regional Oceanic Modeling System (ROMS; Shchepetkin and McWilliams 2005) using a one-way nesting procedure (Mason et al. 2010) and focusing on solutions in the northern Gulf of Mexico (GoM) region with an approximately 500-m, nearly isotropic, horizontal resolution. These 500-m solutions are gradually nested down from a 7-km ROMS solution of the entire Atlantic basin (Barkan et al. 2017a). The atmospheric forcing is climatological with a QuikSCAT-based daily product of scatterometer wind stresses (Risien and Chelton 2008), CORE (Large and Yeager 2009) monthly heat-flux atmospheric forcing, and HOAPS (Andersson et al. 2010) monthly freshwater atmospheric forcing. No tidal forcing is used, and a daily river forcing is applied based on daily river volume flux data from the USGS (http://waterdata.usgs.gov/nwis/rt) for the year 2010. This forcing methodology largely eliminates the generation of internal gravity waves and allows us to focus on the dynamics of submesoscale currents. The vertical mixing of momenta and tracers are parameterized using the *K*-profile parameterization (KPP; Large et al. 1994), and the horizontal diffusion of tracer and momenta are associated with the third-order upstream-biased advection scheme ( appendix D). The analysis in Figs. 5 and 6 is carried out based on two years of equilibrated solutions for winter (January–March) and summer (July–September) months, using a twice-per-day output frequency. Figures 7–9, 11, D1, E1, and E2 are based on offline particle advection, using a 15-min output frequency. Additional information about the numerical setup is provided in Barkan et al. (2017a).

### b. Drifter measurements

LASER took place in northern GoM near the site of the Deepwater Horizon oil spill during the winter of 2016, as part of the Consortium for Advanced Research on the Transport of Hydrocarbons in the Environment (CARTHE). Over a thousand surface drifters were deployed during LASER, along with aerial and shipboard observations to help define the dynamics (D’Asaro et al. 2018). The analysis shown in section 5 is based on drifter observations collected between 20 January and 20 February 2016.

The mostly biodegradable drifters consisted of a floater extending 5 cm below the surface, which contained a GPS tracking and satellite communication system, and a drogue extending 60 cm below the surface, which was hanging beneath the floater on a flexible tether (Novelli et al. 2017). On some floats, the drogue separated from the floater during storms. The detection of drogue loss for the drifters is described in Haza et al. (2018), where a combination of algorithms based on position transmission data and comparison of nearby drifter tracks were used. Here, we apply a more restrictive approach to exclude undrogued drifters by exclusively examining changes in data transmission intervals. The motivation for this is that surface velocity gradient computations based on drifter tracks are very sensitive to velocity differences between nearby drifters. Using both drogued and undrogued drifters leads to anomalously large vorticity and divergence values. When the drifters are without a drogue, they have a tendency to flip over in the water, temporarily placing the antenna of the satellite communication system underneath other components. This leads to detectable increases in transmission intervals (Haza et al. 2018). The default position transmission interval is 5 min for the LASER drifters. Therefore, if the transmission interval averaged over 60 position updates exceeds 6 min the drifter is rejected from further analysis. While this approach may remove some drifters that were still drogued, it works well to avoid noise and anomalous values in the computed surface velocity gradient values. After the undrogued drifter data are removed we apply local polynomial regression [locally weighted scatterplot smoothing (LOWESS)] curve fitting (Cleveland and Devlin 1988) to the raw GPS position drifter tracks to reduce GPS position noise. The span used for the curve fitting is set to be 10 points. The curve fitting also removes the occasional *O*(100) m errors in GPS fixes that occurred a limited number of times due to unknown reasons [see Fig. S2 in the supplemental material of D’Asaro et al. (2018)]. Drifter positions are subsequently interpolated to integer 15-min intervals.

Using the cleaned and interpolated drifter tracks described above, at each time, all possible combinations of three drifters are considered. A triplet of drifters is the smallest cluster from which velocity gradients in two directions can be determined, thus maximizing the number of velocity gradient values that can be applied for statistical analysis.^{6} Every drifter triplet forms a triangle from which we can obtain scale and aspect ratio information. The scale is defined as the root-mean-square (RMS) of the distance between the triangle drifters and the aspect ratio is defined as the ratio of the triangle height to triangle base (an equidistant triangle has an aspect ratio of $3/2$). Smaller-scale triangles are susceptible to position errors and low-aspect-ratio triangles are susceptible to undersampling errors (Ohlmann et al. 2017) and we therefore constrain our analysis to triangle scales of 1–3 km and triangle aspect ratios of $3/4\u20133/2$. For these triangle shape constraints, the gradient errors are dominated by undersampling errors, which are approximated to be 3.4 × 10^{−5} s^{−1} ($\u22480.5\u2009f$) for both vorticity and divergence, and 5.3× 10^{−9} s^{−2} for velocity gradient variances. These estimates are based on an assumed velocity RMS error of 4 cm s^{−1} at 2-km scales (Ohlmann et al. 2017) and are consistent with the RMS errors found in the model data ( appendix E).

## 4. Analysis of simulated fronts and filaments

In this section we test the asymptotic model developed in section 2 against submesoscale-resolving numerical simulations in the northern GoM (section 3a). These numerical simulations are carried out in the same region where the LASER field campaign took place (section 3b) and provide dynamical context for the analysis of drifter measurements that follows (section 5). The general circulation of the GoM consists of a meandering Loop Current to the south, which frequently sheds mesoscale warm-core Loop Current eddies, and active river inflow associated with the Mississippi–Atchafalaya River system (Oey et al. 2005). The corresponding submesoscale currents in the northern part of the basin are affected by the mesoscale strain associated with the Loop Current eddies, the Mississippi–Atchafalaya River system, and are seasonally and spatially variable (Barkan et al. 2017a,b; Choi et al. 2017). This intrinsic dynamical variability makes this region suitable to test our theoretical predictions in a rich mesoscale- and submesoscale-turbulent environment.

### a. Eulerian analysis of model solutions

The asymptotic model [Eqs. (12a)–(12d)] and the inviscid solutions [Eqs. (16) and (18) and appendix C] derived in section 2 are most naturally evaluated in a Lagrangian reference frame. Nevertheless, before we carry out a thorough comparison using Lagrangian particles (section 4b), it is instructive to examine a few Eulerian computed statistical measures. We focus here on the advective frontogenetic tendencies and rates, which are most commonly analyzed in the context of frontogenesis (e.g., Hoskins 1982; Capet et al. 2008b). The analysis of the remaining terms is delayed to the next subsection.

First, we would like to demonstrate the advantage of the frontogenetic rates [*T*; Eqs. (13a)–(13d)] over the frontogenetic tendencies [$Fhor$; Eqs. (1)–(4)] in quantifying the rates of sharpening and blurring of lateral buoyancy and velocity gradients, and the spatial resemblance with the divergence field. Snapshots of $Fhorb$, $Tb$, and $|\u2207hb|2$ illustrate the sensitivity of the frontogenetic tendencies to the local gradient magnitudes (Fig. 4). The seemingly strong frontogenetic tendency near the shelf that extends far offshore in the eastern part of the domain (Fig. 4a) reflects the lateral buoyancy gradient signature of the Mississippi–Atchafalaya River system and the Mississippi River jet (Barkan et al. 2017b) (Fig. 4b). Similarly, the seemingly weak frontogenetic tendency west of the Mississippi River jet reflects the relatively weak buoyancy gradients in the interior of a mesoscale eddy. The frontogenetic rate (Fig. 4c) is unbiased with respect to weak/strong local buoyancy gradients and therefore more accurately represents frontogenetic (and frontolytic) regions. The divergence field (Fig. 4d) is predominantly convergent in frontogenetic regions and is strikingly similar to $Tb$, in agreement with Eq. (14), with typical magnitudes of $O\u2061(f)$. A similar visual resemblance is found for $Tu$ (not shown).

The probability density functions (PDFs) of $Tb/f$, $Tu/f$, and $\u2212\delta /f$, normalized by their corresponding RMS values, demonstrate the agreement with Eq. (14) in both shelf and offshore regions, during winter and summer (Fig. 5). The overall positive skewness in all cases demonstrates the generic propensity for buoyancy and velocity frontogenesis in submesoscale resolving simulations. The RMS values of the convergence field are in good agreement with those of the frontogenetic rates; however, the resemblance in PDF shapes is better for frontogenetic regions compared to frontolytic regions. As a result, the skewness values of the convergence field are lower than those of the frontogenetic rates. Because the physics that leads to frontolysis involves *unbalanced* turbulent processes and secondary instabilities, we do not expect the asymptotic model, which is *balanced* and focuses on inviscid dynamics, to adequately describe frontolysis.

To further assess the validity of Eq. (14) we choose a more stringent measure and examine the equality using scatterplots of $Tb$ and $Tu$ versus $\u2212\delta $ (Fig. 6). The corresponding least squares fits and $R2$ coefficient of determination show a near-perfect agreement with the scaling for $Tu$ during both seasons, as expected based on the weak and strong anisotropic scalings [Eqs. (12b), (15a)]. For $Tb$, a better agreement is found in winter compared with summer. This is explained by the reduction in Rossby number compared with winter [$RMS\u2061(Ro)$ values reduce from 0.72 in winter to 0.4 in summer, where $Ro=\zeta /f$], which suggests that the buoyancy frontogenetic tendency rate will be determined more often by the weak anisotropic scaling [Eq. (15c)] rather than by the strong anisotropic scaling [Eq. (12a)]. Furthermore, the agreement with the scaling improves in flow regions with at least marginally strong frontogenesis [i.e., where $Tb/f>0.2$ and the corresponding $RMS\u2061(Ro)=0.99and0.68$ for winter and summer, respectively]. In fact, the slope gets closer to one and the coefficient of determination improves when the analysis is carried out for even stronger frontogenetic regions (i.e., where $Tb/f>1$, not shown).

### b. Lagrangian analysis of model solutions

Figures 5 and 6 suggest that our asymptotic theory should be most applicable to strong frontogenetic regions, which in our numerical simulations are most often found during wintertime. This is expected because submesoscale currents are most energetic during winter, when the mixed layer is deepest (Callies et al. 2015). We therefore focus our Lagrangian analysis on the month of February, during which the LASER experiment took place (section 5). Specifically, we aim to carefully quantify how the various terms in the asymptotic model [Eqs. (12a)–(12d)] agree with the corresponding terms in the full equation set [Eqs. (1)–(4)]. To this end 5625 equally spaced Lagrangian particles are seeded at the surface within the region shown in Fig. 1 (rectangular box in Fig. 4, top). The different terms in Eqs. (1)–(4) are interpolated onto the particle locations as they evolve with the horizontal surface velocities, and tracked for 48 h beginning on 3 February. The same seeding and tracking procedure is repeated six more times until 17 February. The time period and seeding location are chosen to match the time and place of the LASER field experiment (section 3b). To quantify the dynamics during submesoscale frontogenesis an algorithm is developed to identify events during which all gradient fields $|\u2207hb|2,|\u2207huh|2,\zeta ,\u2009and\u2009\delta $ increase simultaneously ( appendix D). The results shown in Figs. 7–9 and 11 are based on the identification of 656 such events.

The different terms in Eqs. (1)–(4) averaged over all of the frontogenetic events are shown in Figs. 7 and 8, where the breakup to two figures is done to facilitate the presentation. The agreement with our asymptotic model for the $Fhor$ terms (solid and dashed black lines in Fig. 7), the $Fpresu$ term (dashed black and solid red lines in Fig. 8c), and the relationship between velocity gradient variance, vorticity variance, and divergence variance [Eq. (10); dashed black and magenta lines in Fig. 9c] is evident during the entire event duration. The same agreement is found when we pick the median instead of the mean or when we bin the results based on magnitude ranges of the different gradient terms (not shown) as long as we average over *O*(10) events ( appendix D). In addition, Fig. 7 shows that the increase in all gradient fields ($D/Dt$ terms; thick gray lines) is largely dominated by the $Fhor$ terms (solid black lines) with little contribution from the other terms. This is consistent with the analytical solutions in Eqs. (16a)–(16d), which ignore Coriolis and pressure effects. Indeed, the evolution of all gradient fields agrees well with these analytical solutions up until the last hour (green and magenta lines in Fig. 9). This dominance of the advective frontogenetic tendencies is explained by the partial cancellation of the remaining terms in the gradient evolution equations (Figs. 8a–c). Specifically, the pressure ($Fpres$; solid red lines), Coriolis ($FCor$; solid black lines), and vertical mixing ($Vmix$; solid cyan lines) terms do not contribute directly to frontogenesis because they are partially balanced (dashed blue lines in Figs. 8a–c). This TTW balance (Gula et al. 2014; McWilliams et al. 2015; Sullivan and McWilliams 2018), describes a situation in which the convergent ageostrophic velocities and secondary circulation are forced by the vertical momentum flux in boundary layer turbulence. In turn, these convergent ageostrophic motions drive frontogenesis of a preexisting front or filament. In the TTW model, the ageostrophic vorticity is expected to be anticyclonic in response to the primarily cyclonic geostrophic vorticity in a front of filament (McWilliams et al. 2015), as is indeed shown in Fig. 9b (dashed and dot–dashed magenta curves). It therefore appears that in our numerical solutions TTW is the main mechanism responsible for the surface convergent motions that drive frontogenesis. This is another major result of this study, which will be compared with the drifter measurements in section 5.

The $Fvert$ terms are frontolytic (negative) in all cases (green curves in Fig. 8) and are largest during the last hour (see section 6 for more discussion). It is during this time that horizontal diffusion (gray curves) and vertical mixing (cyan curves) are most frontolytic as well, implying that frontal arrest is taking place. Unfortunately, in the last hour the Lagrangian error is also largest [dot–dashed gray lines in Fig. 7; Eq. (D1)], and we cannot reliably quantify the dynamics during this time.

## 5. Measured fronts and filaments

In the previous section we demonstrated a good agreement between the asymptotic model [section 2; Eqs. (12a)–(12d)] and the numerical solutions, during winter. We further identified TTW as the main mechanism responsible for the generation of the surface convergent motions that lead to submesoscale frontogenesis. In this section, we test the applicability of the asymptotic regime to the surface drifter measurements collected during LASER.

The triangle method described in section 3b allows us to compute velocity gradients from drifter observations and to validate the asymptotic model with respect to the $Fhor$ terms in Eqs. (2)–(4), and of the inviscid analytical solutions [Eqs. (16) and (18) and appendix C] with respect to the Lagrangian evolution of $\delta ,\zeta ,$ and $|\u2207huh|2$. To this end we apply the same algorithm used to identify frontogenetic events in the numerical solutions ( appendix D) and quantify all the terms that can be computed using the triangle method (Fig. 10). Because surface drifters rapidly collect on convergent frontal lines the triangle scales and aspect ratios are limited, to minimize position and undersampling errors (section 3b), and we cannot track gradient sharpening as far in time as in the model ( appendix E). As a result, the corresponding event duration is shortened from seven to five hours and the magnitudes of all the terms are reduced. To estimate this magnitude reduction bias we deploy Lagrangian particles in the model solutions and compare Lagrangian derived terms and triangle derived terms during frontogenetic events ( appendix E). Because the position and undersampling errors in the model solutions are calculable we can quantify the bias (dashed versus solid lines in Figs. E2a–c and black versus blue lines in Figs. E2d–f) and correct for it in the measured fields (cf. Figs. 10 and E3).

A good agreement is found between the frontogenetic tendencies of divergence and velocity gradient variance $Fhor\delta ,Fhoru$ [Eqs. (4) and (2)], and the corresponding asymptotic prediction [Eqs. (12d) and (12b); solid green and dashed blue lines in Figs. 10a,c]. Furthermore, $|\u2207huh|2~\delta 2+\zeta 2$ (dashed green and solid black lines in Fig. 10f) in agreement with the asymptotic scaling in Eq. (10). Because the pressure and vertical mixing terms cannot be measured directly from the LASER drifters, we are unable to quantify how well the TTW balance holds, as was done in the model solutions (dashed blue lines in Fig. 8). Nevertheless, there is a good match between the analytical solutions that ignore pressure and Coriolis effects [Eqs. (16b)–(16d)] and the Lagrangian evolution of $\zeta ,\delta ,$ and $|\u2207huh|2$ (dashed blue and solid black lines in Figs. 10d–f), in agreement with the model results (Figs. 9a–c). This strongly suggests that some degree of cancellation between the vertical mixing, Coriolis, and pressure terms must have been at play during LASER, and that TTW is a plausible mechanism for the generation of the surface convergent motions that led to frontogenesis in the field measurements as well.

## 6. Vertical buoyancy fluxes during submesoscale frontogenesis

An important aspect of buoyancy gradient sharpening during frontogenesis is that it is accompanied by positive vertical buoyancy fluxes $w\u2032b\u2032$ that tend to restratify the mixed layer, convert APE to kinetic energy, and inject buoyancy into the ocean interior. It is therefore plausible to relate buoyancy gradient sharpening during frontogenesis with APE evolution. Mechanistically, the dynamical connection is between $w\u2032b\u2032$ and the $Fvertb$ term [Eq. (1)]. To see this we evaluate the leading-order vertical tendency rate

subject to the chosen scaling and coordinate orientation [Eqs. (7), (A3)], where $sb=\u2212bx/bz$ is the isopycnal slope. Equation (20) highlights that a classical ageostrophic secondary circulation around a stably stratified submesoscale feature, that is, downward on the dense side and upward on the light side, must have $w\u2032b\u2032>0$ and $Tvertb~Fvertb<0$ (in agreement with Fig. 8d, green curve). In this section the asymptotic model developed in section 2 is used to estimate the vertical buoyancy fluxes that release APE during submesoscale frontogenesis.

The amount of APE stored in an idealized front of width *l*, with a constant lateral buoyancy gradient $|\u2207hb|$, that extends through a mixed layer of depth $hml$ and constant stratification $N2$, can be computed with respect to an adiabatically sorted reference state (i.e., Winters et al. 1995) to be exactly $\u2061(1/4)|\u2207hb|\u2009hml\u2009l$. If viscous effects are ignored, all of the APE will be converted to kinetic energy during frontogenesis and the mixed layer will restratify completely. Assuming that the buoyancy gradient sharpening evolution follows the solution in Eqs. (16), as suggested by Fig. 9, the mixed layer depth and frontal width will decrease in size like $hml=hml0\u2061(1+\delta 0\tau ),l=l0\u2061(1+\delta 0\tau )$ and

where $l0$ and $hml0$ are the initial frontal width and mixed layer depth, respectively, and $w\u2032b\u2032>0$ for any initial convergence magnitude ($\delta 0<0$). For realistic fronts, the mixed layer does not completely restratify during frontogenesis and, therefore, $hml$ varies much slower than *l* and $|\u2207hb|$ ($Fvertb\u226aFhorb$, Figs. 7d, 8d). Nevertheless, this idealized physical model suggests that a bulk estimate for the vertical buoyancy fluxes of a front or filament undergoing submesoscale frontogenesis should scale like

To quantify the relationship between buoyancy gradient sharpening and APE and to test the scaling [Eq. (22)], we repeat the seven particle seeding experiments described in section 4a at a depth of 10 m, where the vertical velocity does not vanish. Applying the same methodology as in the preceding sections, frontogenetic events are identified and the Eulerian-computed terms in Eq. (1) are interpolated onto the particle positions ( appendix D). In addition to the horizontal and vertical advective tendencies $Fhorb$ and $Fvertb$, we also keep track of the vertical buoyancy fluxes $w\u2032b\u2032$, where primes denote perturbation from a time mean taken over the 48-h particle release durations. Finally, we sort the events based on the averaged magnitudes of $|\u2207hb|2$ during the last hour (hour 7).

The Lagrangian evolution of $w\u2032b\u2032$, $Fvertb$, and the scaling [Eq. (22)] are shown in Fig. 11 for the hundred strongest (solid lines) and hundred weakest (dashed lines) fronts. The buoyancy fluxes and vertical advective tendencies are indeed larger in magnitude for stronger fronts (solid versus dashed gray lines), with $w\u2032b\u2032>0$ and $Fvertb<0$. This demonstrates that the frontolytic vertical tendency term, which is associated with the vertical component of the secondary circulation [Eq. (20)], is closely tied with the release of APE. Most importantly, the buoyancy flux scaling estimate [Eq. (22)] that relies on our asymptotic model (section 2) agrees well with the observed $w\u2032b\u2032$ (Fig. 11).

## 7. Summary and discussion

Submesoscale-current turbulence consists of flow structures with spatial scales of *O*(0.1–10) km, temporal scales on the order of hours to days, in a dynamical regime characterized by *O*(1) Rossby and Richardson numbers (Thomas et al. 2008; McWilliams 2016). Phenomenologically, the flow exhibits anisotropic patterns of lines or streaks with large magnitudes of buoyancy and velocity gradients, cyclonic vorticity, and convergence (Fig. 1). The associated velocity power spectral density (PSD) is shallow ($~kh\u22122$) and can be explained by sharp velocity gradients that are concentrated in the fronts and filaments. In this scale range, the velocity gradient PSD (Fig. 2) has roughly equal contributions from the divergence and vorticity PSDs that comprise it [$\delta ~\zeta $; Eq. (6)].

These dynamical, phenomenological, and spectral characteristics motivate us to develop an asymptotic theory that focuses on the formation of the fronts and filaments, that is, the Lagrangian evolution of $|\u2207hb|2,|\u2207hu|2,\zeta ,\delta $, which we generally refer to as submesoscale frontogenesis (section 2). A major result of our asymptotic theory is that convergent horizontal advection dominates the frontogenetic tendencies ($Fhor$ terms) that drive frontogenesis [Eq. (14)]. Mechanistically, as fronts and filaments sharpen and a secondary circulation ensues, a strong ageostrophic cross-frontal flow develops with a similar magnitude to the alongfrontal flow. This ageostrophic cross-frontal flow is convergent near the surface and dominates the dynamics of frontal sharpening. This is different than classical deformation-driven frontogenesis (Hoskins and Bretherton 1972), where the cross-front velocity is assumed to be weaker than the alongfront velocity at all times ($\delta \u226a\zeta $) and where horizontally nondivergent confluent motions govern gradient sharpening. We argue that although the “classical” assumptions may be correct for the early stage of deformation-driven frontogenesis, in later stages the cross-front velocity strengthens and the convergent motions dominate the confluent motions. Analytical solutions for the inviscid Lagrangian evolution of the fronts and filaments [Eqs. (16) and (18) and appendix C] further highlight the importance of horizontal convergence and of ageostrophic vorticity in determining the flow evolution during frontal sharpening and, specifically, the hypothetical finite time singularity. The asymptotic model can further be used to determine the positive vertical buoyancy fluxes that accompany submesoscale frontogenesis [Eq. (22), Fig. 11], and that are presumably important for mixed layer restratification and for buoyancy (and tracer) injection into the ocean interior.

Submesoscale resolving numerical simulations and drifter measurements collected during the LASER field campaign confirm the applicability of the asymptotic regime to strong frontogenesis in the northern GoM, during winter (Figs. 7, 10a–c). Furthermore, through careful investigation of all viscid and inviscid terms in the evolution equations of the modeled gradient fields [Eqs. (1)–(4)], we demonstrate that the dominant mechanism responsible for the generation of the ageostrophic convergent motions is boundary layer turbulence (TTW; Gula et al. 2014; McWilliams et al. 2015; Sullivan and McWilliams 2018). In the TTW model the vertical mixing term associated with boundary layer turbulence, which can be substantial, is largely balanced by the pressure and Coriolis terms (Fig. 8), and generates ageostrophic motions that drive frontogenesis. As a result the corresponding Lagrangian evolution of the frontogenetic fields is dominated by the frontogenetic tendencies [$Fhor$ terms in Eqs. (1)–(4)], with negligible contribution from the other terms ($Fpres$ and $FCor$), and is in good agreement with the analytical solutions that ignore Coriolis and pressure effects [Eqs. (16a)–(12d), Fig. 9]. Although we cannot explicitly compute the pressure and vertical mixing terms from the data, the fact that the same analytical solutions match the Lagrangian evolution of the drifter-based gradient fields (Figs. 10d–f) strongly suggests that an approximate TTW balance was at play during LASER as well.

Our asymptotic model and solutions make no assumptions about the underlying physical processes that produce the ageostrophic convergent motions that lead to gradient sharpening (e.g., deformation, mixed layer instabilities, boundary layer turbulence), and it can therefore be regarded as a generic model for submesoscale frontogenesis. If TTW is an important frontogenetic mechanism globally, as it is in the Gulf of Mexico, then it implies that boundary layer turbulence plays a significant role in determining the transport of materials from the mixed layer to the ocean interior.

The clearest limitation of our model is ignoring the onset of frontal instabilities that are most likely the cause of frontal arrest (McWilliams and Molemaker 2011; Sullivan and McWilliams 2018). We also note that our theory ignores the effects of near-inertial motions and internal gravity waves, which might overlap in scales (Callies and Ferrari 2013) and can alter the dynamics at fronts (Whitt and Thomas 2015; Barkan et al. 2017c; Thomas 2017). In addition, the theory does not consider the dynamical effects of surface gravity waves, which has been previously shown to affect frontal evolution (McWilliams and Fox-Kemper 2013; Hamlington et al. 2014; Suzuki et al. 2016; McWilliams 2018). Nevertheless, we argue that it is instructive to first develop a theory for the underlying “balanced,” submesoscale-current dynamics, and leave the incorporation of frontal instabilities and wave interactions to future work.

## Acknowledgments

We thank three anonymous referees for comments that greatly improved the presentation of this work. This work was made possible by a grant from The Gulf of Mexico Research Initiative through the CARTHE consortium. Data are publicly available through the Gulf of Mexico Research Initiative Information and Data Cooperative (GRIIDC) at https://data.gulfresearchinitiative.org (doi: 10.7266/N7PK0DK2, doi: 10.7266/N7F18X4S, doi: 10.7266/N75H7DQ5, doi: 10.7266/N7JS9NVS, and doi: 10.7266/N79885FW). RB is further supported by ONR N000141812697 and JCM is further supported by ONR N000141410626. The Extreme Science and Engineering Discovery Environment (XSEDE) provide support for computing.

### APPENDIX A

#### Derivation of the Scaled Frontogenetic Tendency Equations

We begin by applying the scaling in Eq. (7) and the distinguished limit in Eq. (9) to the inviscid Boussinesq equations of motion, resulting in the following nondimensional form

where all fields are nondimensional; $\u2061(u,\upsilon ,w)$ are the velocities in the $\u2061(x^,y^,z^)$ directions; *b* is the buoyancy; $\Phi $ is the geopotential; subscripts denote derivatives; and the nondimensional numbers $Ro,\epsilon ,$ and *λ* are described in Eq. (8).

The buoyancy frontogenetic tendency equation Eq. (1) is derived by taking the horizontal gradient $\u2207h$ of Eq. (A1d) and then multiplying it by $\u2207hb$, where the subscript *h* denotes horizontal vector components. The index notation is chosen to concisely write the terms on the right hand side. Because the velocity gradient tensor $um,j$ can be decomposed into the symmetric part $Smj$ and antisymmetric part $\Omega mj$, namely,

the buoyancy frontogenetic equation can be simplified to

where $j,k=1,2$, summation over repeated indices is assumed, and commas denote derivatives. Only $Skj$, typically referred to as the (horizontal) strain rate tensor, contributes to $Fhorb$, whereas the contribution by $\Omega kj$, typically referred to as the (horizontal) rotation tensor, cancels out. Because $Fhorb$ in Eq. (A3) is written as a tensor dot product, it is in coordinate invariant form. Furthermore, Eq. (A3) illustrates that the components of the horizontal strain-rate tensor determine the rate of frontogenesis of frontolysis. The nondimensional $Skj$ and $b,kb,j$ subject to the scaling in Eq. (7) are

where the diagonal components of $Skj$ are often attributed to normal strains and the off-diagonal components are often attributed to shear strains. Note that because the horizontal velocity is divergent the normal strains are not $\xb1(ux\u2212\upsilon y)$. The corresponding nondimensional $Fhorb$ and $Fvertb$ are

and

where $N2=bz$.

The velocity frontogenetic tendency equation in Eq. (2) is derived by taking the horizontal gradient of Eqs. (A1a) and (A1b), adding them up, and multiplying by $ui,j$. Applying the decomposition in Eq. (A2), the velocity frontogenetic tendency equation can be simplified to

where, as for $Fhorb$, only the strain rate tensor contributes to $Fhoru$, and $\Lambda i=(uz,\upsilon z)$. The components of the nondimensional $Skj$ are written in Eq. (A4) and

The corresponding nondimensional $Fhoru$, $Fvertu$, and $Fpresu$ become

and

The vertical vorticity and horizontal divergence equations [Eqs. (3) and (4)] are derived by taking the curl and divergence of the momentum equations [Eqs. (A1a) and (A1b)], respectively. The corresponding nondimensional $Fhor\zeta ,Fvert\zeta ,Fhor\delta ,Fvert\delta $, and $pressure\delta $ are

### APPENDIX B

#### The Asymptotic Model in Coordinate Invariant Form

In the asymptotic limit $\epsilon \u226a1,\lambda \u226a1,Ro~O\u2061(1)$ the various frontogenetic tendency terms ($F$) and the pressure term in the velocity gradient evolution equation ($Fpresu$) take the form

where H.O.T denotes higher-order terms; det denotes the determinant; $fug=k^\xd7\Phi $; the Jacobian is defined as $J\u2061(a,b)=axby\u2212bxay$; and $Q\u02dcj=\u2212[J\u2061(\upsilon ,b),J\u2061(u,b)]$. The terms to the right of the arrows in Eqs. (B1) are in coordinate invariant form. To see that the $ux$ terms in Eqs. (B1) generalize to the horizontal divergence $\delta =Roux+\epsilon \upsilon y$, we decompose the nondimensional horizontal velocity field into its divergent and rotational components (i.e., Helmholtz decomposition), namely,

where *ϕ* is a nondimensional velocity potential; *ψ* is a nondimensional streamfunction; and $\varphi ~Ro\psi $ under our scaling [Eq. (7)]. Combining Eqs. (B1) and (B2) yields

which demonstrates that the leading-order advective frontogenetic rates [Eqs. (13a)–(13d)] are governed by the divergent part of the horizontal velocity field only, and that $ux$ must generalize to the coordinate invariant *δ* (the coordinate invariant strain rate is $SkjSkj=ux2+1/2\upsilon x2\u2260ux$).

It is useful (section 2a) to contrast Eqs. (B1) with the dynamical regime $\lambda \u226a1,Ro~\epsilon \u226a1$. In this dynamical regime

### APPENDIX C

#### Solution to Eqs. (12a)–(12d) for Time-Varying Ageostrophic Vorticity

It is possible to generalize the solution in Eqs. (18) and allow for a time variable ageostrophic vorticity by replacing Eq. (17) with

where we have made the transformation to a Lagrangian reference frame. Equation (C1) is a Riccatti equation and by setting $\delta =w\tau /w$, it becomes

If we assume that the increase in ageostrophic vorticity is algebraic in nature (which is sensible in order for the total cyclonic vorticity to blow up) such that $g\u2061(\tau )=a2\tau n$ then Eq. (C2) can be solved using power series solutions, by setting $w=\u2211k=0\u221eAk\tau k$. After some algebra we get

with the recursion relation

and with the initial conditions $A2=w\u2061(0)/2,A3=w\u2032\u2061(0)/6$. The solution for *δ* based on the recursion relation [Eq. (C5)] is plotted in Fig. 3 for different *a* values (red curves). The solution for the other fields can then be determined by separation of variables because these power series are integrable.

### APPENDIX D

#### Closing the Balance and the Selection Algorithm for Frontogenetic Events

##### a. Closing the balance

The different terms in the momentum and tracer equations are computed online in a way that is exactly consistent with ROMS’s discretization, and are saved every 15 min.

The required derivative operations are applied to the balances computed in step 1 to arrive at Eqs. (1)–(4). We carefully compute the various advection terms because ROMS is written in volume flux form. The $Vmix$ terms are computed based on the KPP (Large et al. 1994) mixing profiles and values. The $Hdiff$ terms in ROMS are fourth-order hyperviscous/diffusive and are associated with the advection schemes for tracers and momenta. They are computed as the difference between a third-order upstream-biased and a fourth-order central advection schemes.

- Two-dimensional Lagrangian trajectories are computed offline using the surface velocities output at the same fifteen minute intervals. The different terms computed in steps 1 and 2 are then interpolated onto the particle positions. Note that the different Lagrangian tendency terms consist of the interpolated three-dimensional material derivatives. The Lagrangian errors shown in Fig. 7 are computed as the difference between the Lagrangian and Eulerian $Fhor$ terms ($FhorLag\u2212FhorEul$), where the $FhorEul$ terms are computed directly from the model and are interpolated onto the particle locations. The $FhorLag$ terms are computed as$FhorLag\delta =\u2202\u2202\tau \u2009\delta \u2212\u2207h\u22c5DhadvDtuh,FhorLag\zeta =\u2202\u2202\tau \u2009\zeta \u2212\u2207h\xd7DhadvDtuh,$(D1)

and similarly for $|\u2207huh|2$ and $|\u2207hb|2$. The first term on the right-hand side in Eq. (D1) corresponds to the Eulerian-computed fields that are interpolated and *τ*-differentiated following the Lagrangian trajectories. The second term on the right-hand-side corresponds to the divergence/curl of the Eulerian-computed horizontal material derivative (subscript *h*) in advective form (subscript adv) that is interpolated onto the particle locations. Because any particle advection code is subject to the discretization errors of the modeled Eulerian velocities, a Lagrangian advection error is always expected. As shown in Fig. 7, these errors are largest at the final stage of frontogenesis when the gradient fields are largest and, hence, have spatial scales closest to the grid resolution.

##### b. Selection algorithm

The algorithm used to select for frontogenetic events is based on the gradient fields $\delta ,\zeta ,|\u2207hb|2$, and $|\u2207huh|2$ that are interpolated onto the particle locations following the procedure described above. On every particle $pi$ we identify the maxima $maxi$ and associated time $\tau maxi$ for each one of the gradient fields (for the divergence field it was maxima of convergence). As it turns out, all fields have a large number of common maxima and events of simultaneous increase in all fields are easily found. As long as the gradient field values prior to every maxima $\tau maxi$ are gradually increasing for at least 7 h (i.e., allowing for 28 data points), we define the corresponding Lagrangian time series as a frontogenetic event. We allow for the possibility of having local maxima during an event as long as their associated magnitudes are smaller than 40% of $maxi$. We further require that the divergence values during every frontogenetic event are always negative (i.e., convergent). Because we are interested in relatively strong frontogenetic events, we discard events with maximal convergent values smaller than $1.5f$ and require every gradient field magnitude to at least double during an event. We experimented with events shorter and longer than 7 h (between 4 and 9 h) and found that our results did not vary much. We choose 7 h for the model solutions because it allows the largest number of events to be identified. The same algorithm is applied to the LASER drifter data using the triangle-based gradient computation (section 3b). Because of the position error and higher signal to noise ratio, we only make gradient estimates at scales greater or equal to 1 km. As a result, the highest $\delta ,\zeta ,$ and $|\u2207huh|2$ values are not measured, and we find that 5 h allows for the largest number of events to be identified in the data.

Figure D1 illustrates the events selected in the model solutions using the above algorithm for seven particles (gray curves) that collected onto the same submesoscale feature. Although the variability among the particles is considerable, their mean evolution agrees well with the theoretical prediction [Eq. (16d)]. A similar agreement was found for other submesoscale structures as long as we averaged over *O*(10) particles (not shown).

### APPENDIX E

#### Error Biases in the Triangle-Based Gradient Computation Method

To test the ability of the triangle method (section 3b) to quantify dominant balances during frontogenetic events we apply it to one of our particle release experiments in the model. Because the modeled particles have no position error and the subgrid-scale energy is small, it allows us to extend the triangle aspect ratios to $3/40$ (instead of $3/4$ in the observations, section 3b), and we carefully examine the undersampling error and any other biases of the method with respect to a ground “truth.” Figure E1 demonstrates how the triangles deform (aspect ratio decreases) as the front sharpens and particles collect onto convergence lines. It further indicates what can be expected from the aspect-ratio limit used in the observations (Fig. E1, dashed black line). Figure E2 demonstrates the biases associated with the triangle method. In all measured fields the triangle method underestimates the gradient magnitudes (Figs. E2d–f), with the biases larger for convergence and velocity gradient variance than for vorticity. As expected, the biases decrease when we allow for lower aspect ratios than the ones applied to the drifter measurements (Figs. E2d–f, thin vs thick blue lines), although they do not disappear completely even for the lowest aspect ratios used. Because the particles aggregate in nearly straight lines the biases are expected to vanish only in the limit that the aspect ratio goes to zero. Having quantified the biases using same scale and aspect ratio limits as in the drifter measurements (Figs. E2a–c, dashed vs solid lines; Figs. E2d–f, thick blue vs thick black lines) we correct the drifter measurements using the ratios between the triangle derived quantities to the drifter derived quantities (Fig. 10). For example, the corrected divergence field in the LASER measurements (solid black line in Fig. 10d) $\delta \u2061(\tau )corrected=\delta \u2061(\tau )\xd7\delta \u2061(\tau )particle/\delta \u2061(\tau )triangle$, where $\delta \u2061(\tau )$ is the uncorrected divergence field (solid black line in Fig. E3d), and $\delta \u2061(\tau )particle$ and $\delta \u2061(\tau )triangle$ are the model-based divergence fields computed on particles and using the triangle methods, respectively (thick black and blue lines in Fig. E2d). All the plotted terms in Fig. 10 are corrected in the same manner.

Finally, we compute the undersampling errors of the triangle method in the model based on the RMS differences between the triangle-based gradient computation and the “true” Eulerian values. These come out to be $0.5f$ for both vorticity and divergence, where we focus on frontogenetic regions. These model-based error estimates are in agreement with the observationally based estimates discussed in section 3b and are the basis for the shaded error bars in Figs. 10a–c and E3a–c.

## REFERENCES

*Nat. Commun.*,

**6**, 6862, https://doi.org/10.1038/ncomms7862.

*Circulation in the Gulf of Mexico: Observations and Models*,

*Geophys. Monogr.*, Vol. 161, Amer. Geophys. Union, 31–56.

*Ocean Modeling in an Eddying Regime*,

*Geophys. Monogr.*, Vol. 177, Amer. Geophys. Union, 17–38.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Note that the Garrett–Munk (Garrett and Munk 1972) internal wave continuum spectral slope is also close to $kh\u22122$ at submesoscales, demonstrating that knowledge of the spectral slope is generally insufficient to determine the underlying physical structures that lead to it.

^{2}

For quasigeostrophic flow, which is horizontally nondivergent, $E\u02dc[\u2207huh\u2061(kh)]=E\u02dc\u2061[\zeta \u2061(kh)]=|kh|2E\u02dc\u2061[u\u2061(kh)]~kh\u22121$, in agreement with the enstrophy cascade in geostrophic turbulence (Salmon 1980).

^{4}

In this scaling regime $|\u2207huh|2~\zeta 2$, and so $Tu=T\zeta =\u2212\delta $ by construction.

^{5}

When $\zeta ag=+a2\u2009f$, the solution to Eq. (17) does not have a finite time singularity.

^{6}

Increasing the number of drifters in a cluster will reduce the error in the computed gradients, but only at a rate proportional to the square root of the drifter number (Ohlmann et al. 2017).