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

 
DDt12(b,jb,j)=uk,jb,kb,jFhorbN2w,jb,jFvertb+Vmixb+Hdiffb,
(1)

where D/Dt is the material derivative associated with the three-dimensional velocity field (u,υ,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:

 
DDt12(ui,jui,j)=uk,jui,kui,jFhoruΛiw,jui,jFvertuui,jΦ,ijFpresu+Vmixu+Hdiffu,
(2)
 
DDtζ=δζFhorζfδFCorζεijΛjw,iFvertζ+Vmixζ+Hdiffζ,and
(3)
 
DDtδ=uk,iui,kFhorδ+fζFCorδΛiw,iFvertδΦ,iiFpresδ+Vmixδ+Hdiffδ.
(4)

Above, ζ=εijui,j denotes the vertical vorticity component; δ=ui,i is the horizontal divergence; Φ is the geopotential; f is the Coriolis frequency; i,j,k=1,2; Λi=(uz,υz) denotes the vertical shear components; ε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ζ, and Fhorδ is equivalent to that of Fb, where flow regions with positive Fhoru and Fhorζ correspond to the enhancement of lateral velocity gradients and of cyclonic vorticity, respectively. Similarly, flow regions with negative Fhorδ 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).

Fig. 1.

Representative wintertime surface snapshots of (a) |hb|2 (s−4), (b)|huh|2 (s−2), (c) ζ/f, and (d) δ/f from the numerical simulations (section 3a) in the box region shown in Fig. 4a. The overly saturated structure in the southwest corner of (a) is associated with the Mississippi River jet (Barkan et al. 2017b). It has similar spatial patterns to the remaining fields.

Fig. 1.

Representative wintertime surface snapshots of (a) |hb|2 (s−4), (b)|huh|2 (s−2), (c) ζ/f, and (d) δ/f from the numerical simulations (section 3a) in the box region shown in Fig. 4a. The overly saturated structure in the southwest corner of (a) is associated with the Mississippi River jet (Barkan et al. 2017b). It has similar spatial patterns to the remaining fields.

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˜[uh(kh)] and E˜[b(kh)] with spectral slopes of khβ, where β~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:

 
E˜[huh(kh)]=kh2E˜[u(kh)]~kh2β,and
(5a)
 
E˜[hb(kh)]=kh2E˜[b(kh)]~kh2β,
(5b)

where kh denotes the horizontal wavenumber vector, E˜ denotes the power spectral density (PSD), and when β~2 the spectral slopes of E˜[huh(kh)] and E˜[hb(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,

 
E˜[huh(kh)]=E˜[ζ(kh)]+E˜[δ(kh)].
(6)

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 kh2 (Fig. 2a). Specifically, the ratio of the vorticity PSD to the divergence PSD is about six at kh = 10−4 m−1 and gradually decreases at higher wavenumbers, suggesting that the magnitudes of ζ and δ are asymptotically similar at these (submeso) scales.

Fig. 2.

The model-based PSD of (a) the horizontal velocity and (b) the horizontal velocity gradient huh (dashed black), the vorticity (solid gray), and the divergence (solid black) computed in the box region shown in Fig. 4a. All PSDs are time-averaged during February and are displayed as a function of the horizontal wavenumber kh. The horizontal wavelength λ is also marked.

Fig. 2.

The model-based PSD of (a) the horizontal velocity and (b) the horizontal velocity gradient huh (dashed black), the vorticity (solid gray), and the divergence (solid black) computed in the box region shown in Fig. 4a. All PSDs are time-averaged during February and are displayed as a function of the horizontal wavenumber kh. The horizontal wavelength λ is also marked.

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,

 
x~l,y~L,z~hml,υ~V,u~RoV,w~RoVhmll,T~lRoV,Φ~fVl,b~fVlhml,Λ1=uz~ξRoV/hml,Λ2=υz~ξV/hml,N2~ηb/hml.
(7)

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=gρ/ρ0 is the buoyancy scale associated with the hydrostatic geopotential Φ. The vertical stratification N2 is set to be smaller than b/hml by a parameter η1 to highlight the well-mixed structure of the surface mixed layer. For the same reason, the vertical shears uz,υz are set to be smaller than RoV/hml and V/hml, respectively, by a parameter ξ1. Nondimensional parameters are

 
Ro=Vfl,ε=lL,λ=hmll.
(8)

Here Ro is the Rossby number, ε is the anisotropy ratio, and λ is the aspect ratio. We further choose the distinguished limit3

 
ξ~η~ε1.
(9)

The main novelty of the proposed scaling is that for Ro~O(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, ζ~δ, and the velocity gradient magnitude squared

 
|huh|2~ζ2+δ2
(10)

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 ζδ and |huh|2~ζ2.

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

 
b,jb,j=bx2+ε2by2,ui,jui,j=υx2+Ro2ux2+ε2υy2+Ro2ε2uy2,ζ=υx+εRouy,δ=Roux+ευy,
(11)

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

 
DDt12|hb|2=δ|hb|2Fhorb+Vmixb+Hdiffb,
(12a)
 
DDt12|huh|2=δ|huh|2FhorufδζgFpresu+Vmixu+Hdiffu,
(12b)
 
DDtζ=δζFhorζfδFCorζ+Vmixζ+Hdiffζ,and
(12c)
 
DDtδ=δ2Fhorδ+fζagFCorδ+Fpresδ+Vmixδ+Hdiffδ,
(12d)

where the simplification in Eqs. (12a)(12d) makes the use of index notation unnecessary, is the gradient operator, the h subscripts denote horizontal vector components, and the equations are written in dimensional form. The Laplacian of the pressure h2Φ is written as fζ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ζ,andFhorδ, 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

 
Tb=Fhorb|hb|2,
(13a)
 
Tu=Fhoru|huh|2,
(13b)
 
Tζ=Fhorζζ,
(13c)

and

 
Tδ=Fhorδδ,
(13d)

we obtain the nontrivial relationship

 
Tb=Tu=Tζ=Tδ=δ.
(14)

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 δ. 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 λ1,Ro~ε1, which describes the initiation of frontogenesis. In this dynamical regime ( appendix B) the leading-order advective frontogenetic tendencies Fhoru and Fhorζ remain the same as in Eqs. (12b) and (12c), that is,

 
Fhoru=δ|huh|2
(15a)

and

 
Fhorζ=δζ,
(15b)

but

 
Fhorb=Skjb,kb,j,
(15c)

and

 
Fhorδ=uk,iui,k.
(15d)

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 δ, 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/τ. 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

 
|hb(τ)|2=|hb0|2[1+δ0(ττ0)]2,
(16a)
 
|huh(τ)|2=|huh0|2[1+δ0(ττ0)]2,
(16b)
 
ζ(τ)=ζ01+δ0(ττ0),and
(16c)
 
δ(τ)=δ01+δ0(ττ0),
(16d)

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 τ0, when the strong anisotropic flow regime is reached [i.e., when ε1,Ro~O(1)]. The solutions above demonstrate that for any initially convergent flow a finite time singularity is expected at τ=τ0δ01. Because δ0~f, this blowup time is faster than the semigeostrophic blowup time of log2/α (Hoskins and Bretherton 1972; Hoskins 1982), where αf 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 ζag and for simplicity we set ζag=±a2f, 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

 
DDtδ=δ2Fhorδa2f2FCorδ+Fpresδ.
(17)

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

 
|hb(τ)|2=|hb0|2cos(C)cos[af(ττ0)+C],
(18a)
 
|huh(τ)|2=ζ(τ)2+δ(τ)2,
(18b)
 
ζ(τ)=ζ0cos(C)2f sin[af(ττ0)+2C2]sin[af(ττ0)2]cos[af(ττ0)+C],and
(18c)
 
δ(τ)=aftan[af(ττ0)+C],
(18d)

where C=arctan[δ0(af)1], af=|ζ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 |huh|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 τ=τ0+(π/2C)(af)1. In  appendix C we generalize the above solutions for a time-varying ageostrophic anticyclonic vorticity and, for illustration, show the solution to

 
τδ=δ2Fhorδa2f2τ2FCorδ+Fpresδ.
(19)

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.

Fig. 3.

Analytical solutions for the divergence field evolution without Coriolis and pressure effects [solution to Eq. (17) with FCorδ+Fpresδ=0, blue], with constant anticyclonic ageostrophic vorticity [solution to Eq. (17), black], and with increasing anticyclonic ageostrophic vorticity [solution to Eq. (19), red]. The initial divergence value in all cases is δ0=0.55f. For the solutions with Coriolis and pressure effects (black and red) FCorδ+Fpresδ=fζag, with the ageostrophic vorticity ζag=a2f and a2=0.1 (dot–dashed), a2=0.2 (solid), and a2=0.4 (dashed).

Fig. 3.

Analytical solutions for the divergence field evolution without Coriolis and pressure effects [solution to Eq. (17) with FCorδ+Fpresδ=0, blue], with constant anticyclonic ageostrophic vorticity [solution to Eq. (17), black], and with increasing anticyclonic ageostrophic vorticity [solution to Eq. (19), red]. The initial divergence value in all cases is δ0=0.55f. For the solutions with Coriolis and pressure effects (black and red) FCorδ+Fpresδ=fζag, with the ageostrophic vorticity ζag=a2f and a2=0.1 (dot–dashed), a2=0.2 (solid), and a2=0.4 (dashed).

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 79, 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/43/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 (0.5f) 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 |hb|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(f). A similar visual resemblance is found for Tu (not shown).

Fig. 4.

Snapshots of (a) surface buoyancy frontogenetic tendency Fhorb [Eq. (1)], (b) surface lateral buoyancy gradient magnitude squared, (c) surface buoyancy frontogenetic rate Tb [Eq. (13a)] normalized by the Coriolis frequency f, and (d) normalized surface divergence field δ/f, from the numerical solutions during winter. All fields are plotted in frontogenetic regions only (i.e., where Fhorb,Tb>0), and in (a) and (b) with a logarithmic color bar. The probability distribution functions of frontogenetic and frontolytic regions are displayed in Fig. 5. The rectangular box in (a) indicates the location where Lagrangian particles are seeded.

Fig. 4.

Snapshots of (a) surface buoyancy frontogenetic tendency Fhorb [Eq. (1)], (b) surface lateral buoyancy gradient magnitude squared, (c) surface buoyancy frontogenetic rate Tb [Eq. (13a)] normalized by the Coriolis frequency f, and (d) normalized surface divergence field δ/f, from the numerical solutions during winter. All fields are plotted in frontogenetic regions only (i.e., where Fhorb,Tb>0), and in (a) and (b) with a logarithmic color bar. The probability distribution functions of frontogenetic and frontolytic regions are displayed in Fig. 5. The rectangular box in (a) indicates the location where Lagrangian particles are seeded.

The probability density functions (PDFs) of Tb/f, Tu/f, and δ/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.

Fig. 5.

PDFs of modeled surface Tb/f,Tu/f, and δ/f normalized by their corresponding RMS values in water depth shallower than 150 m (shelf) and deeper than 500 m (offshore) for winter and summer. RMS and skewness values for Tb/f,Tu/f, and δ/f are shown in the top-right corner of each panel. The maximal standard errors for the reported RMS and skewness values are at most 10% and 0.3% of the signal, respectively, where the number of degrees of freedom is based on an assumed decorrelation time scale of 3 days and an assumed decorrelation length scale of 10 km.

Fig. 5.

PDFs of modeled surface Tb/f,Tu/f, and δ/f normalized by their corresponding RMS values in water depth shallower than 150 m (shelf) and deeper than 500 m (offshore) for winter and summer. RMS and skewness values for Tb/f,Tu/f, and δ/f are shown in the top-right corner of each panel. The maximal standard errors for the reported RMS and skewness values are at most 10% and 0.3% of the signal, respectively, where the number of degrees of freedom is based on an assumed decorrelation time scale of 3 days and an assumed decorrelation length scale of 10 km.

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 δ (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(Ro) values reduce from 0.72 in winter to 0.4 in summer, where Ro=ζ/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(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).

Fig. 6.

Scatterplots from model solutions of (a) Tb/f vs δ/f in winter, (b) Tu/f vs δ/f in winter, (c) Tb/f vs δ/f in summer, and (d) Tu/f vs δ/f in summer. Diagonal lines indicate the least squares fits over the entire range of frontogenetic rates (solid dark gray) and over the range Tb/f,Tu/f>0.2 (dashed light gray). Legends indicate the corresponding least squares fit slopes and, in parentheses, the R2 coefficient of determination. The p value in all cases is less than 0.005.

Fig. 6.

Scatterplots from model solutions of (a) Tb/f vs δ/f in winter, (b) Tu/f vs δ/f in winter, (c) Tb/f vs δ/f in summer, and (d) Tu/f vs δ/f in summer. Diagonal lines indicate the least squares fits over the entire range of frontogenetic rates (solid dark gray) and over the range Tb/f,Tu/f>0.2 (dashed light gray). Legends indicate the corresponding least squares fit slopes and, in parentheses, the R2 coefficient of determination. The p value in all cases is less than 0.005.

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 |hb|2,|huh|2,ζ,andδ increase simultaneously ( appendix D). The results shown in Figs. 79 and 11 are based on the identification of 656 such events.

Fig. 7.

The various terms in the evolution Eqs. (1)(4) for (a) divergence, (b) vorticity, (c) velocity gradient variance, and (d) buoyancy gradient variance, computed from the model solutions. The asymptotic predictions for the Fhor terms are plotted with dashed black lines. The residual of all the terms (dashed gray) in Eqs. (1)(4) is denoted by “residual” (the remaining terms are shown in Fig. 8) and the Lagrangian error (dotted gray) denotes the error associated with Lagrangian particle advection model ( appendix D). The fields in (a) and (b) are normalized by f2 and are in s−3 and in s−5 in (c) and (d), respectively. All quantities are averaged over the selected 656 events ( appendix D). For information about the variability between events and release experiments refer to Figs. 9 and D1. The similarity in (a)–(d) between the Fhor terms (solid black), the corresponding asymptotic predictions [Eqs. (12a)(12d); dashed black), and the D/Dt terms (thick solid gray) demonstrates that the asymptotic model (section 2) adequately describes the advective tendency terms and that these advective tendencies are the main frontogenesis driver.

Fig. 7.

The various terms in the evolution Eqs. (1)(4) for (a) divergence, (b) vorticity, (c) velocity gradient variance, and (d) buoyancy gradient variance, computed from the model solutions. The asymptotic predictions for the Fhor terms are plotted with dashed black lines. The residual of all the terms (dashed gray) in Eqs. (1)(4) is denoted by “residual” (the remaining terms are shown in Fig. 8) and the Lagrangian error (dotted gray) denotes the error associated with Lagrangian particle advection model ( appendix D). The fields in (a) and (b) are normalized by f2 and are in s−3 and in s−5 in (c) and (d), respectively. All quantities are averaged over the selected 656 events ( appendix D). For information about the variability between events and release experiments refer to Figs. 9 and D1. The similarity in (a)–(d) between the Fhor terms (solid black), the corresponding asymptotic predictions [Eqs. (12a)(12d); dashed black), and the D/Dt terms (thick solid gray) demonstrates that the asymptotic model (section 2) adequately describes the advective tendency terms and that these advective tendencies are the main frontogenesis driver.

Fig. 8.

As in Fig. 7, but for the remaining terms in Eqs. (1)(4). TTW denotes the turbulent thermal wind balance for (a) divergence, FCorδ+Fpresδ+Vmixδ [Eq. (4)]; (b) vorticity, FCorζ+Vmixζ [Eq. (3)]; and (c) velocity gradient variance, Fpresu+Vmixu [Eq. (2)]. The dashed black line in (c) denotes the asymptotic prediction for the Fpresu term in Eq. (12b). The orange line in (d) is associated with the temperature relaxation to climatological SST applied in ROMS and is required to close the buoyancy gradient variance balance. The fields in (a) and (b) are normalized by f2 and are in s−3 and s−5 in (c) and (d), respectively. For information about the variability between events and release experiments refer to Figs. 9 and D1. The dominance of the Fhor terms in driving frontogenesis (Fig. 7) is a result of the partial cancelation between the Coriolis, pressure, and vertical mixing terms (TTW balance) in (a)–(c).

Fig. 8.

As in Fig. 7, but for the remaining terms in Eqs. (1)(4). TTW denotes the turbulent thermal wind balance for (a) divergence, FCorδ+Fpresδ+Vmixδ [Eq. (4)]; (b) vorticity, FCorζ+Vmixζ [Eq. (3)]; and (c) velocity gradient variance, Fpresu+Vmixu [Eq. (2)]. The dashed black line in (c) denotes the asymptotic prediction for the Fpresu term in Eq. (12b). The orange line in (d) is associated with the temperature relaxation to climatological SST applied in ROMS and is required to close the buoyancy gradient variance balance. The fields in (a) and (b) are normalized by f2 and are in s−3 and s−5 in (c) and (d), respectively. For information about the variability between events and release experiments refer to Figs. 9 and D1. The dominance of the Fhor terms in driving frontogenesis (Fig. 7) is a result of the partial cancelation between the Coriolis, pressure, and vertical mixing terms (TTW balance) in (a)–(c).

Fig. 9.

The Lagrangian evolution of (a) normalized divergence, (b) normalized vorticity, (c) velocity gradient variance (s−2), and (d) buoyancy gradient variance (s−4) computed from the model solutions. The mean values over the selected 656 events ( appendix D) and the corresponding analytical solutions [Eqs. (16a)(16d)] are shown in magenta and green, respectively. The normalized geostrophic (dashed magenta) and ageostrophic (dot–dashed magenta) vorticities are also plotted in (b). The dashed black line in (c) indicates the asymptotic prediction for the velocity gradient variance in terms of the divergence and vorticity variances [Eq. (10)]. Solid gray lines indicate the mean values for each of the seven particle release experiments, and the shaded blue envelopes indicate a representative event spread (one standard deviation) for each experiment. Similar variabilities between and during release experiments were also found for the various terms in Figs. 7, 8, and 11 (not shown).

Fig. 9.

The Lagrangian evolution of (a) normalized divergence, (b) normalized vorticity, (c) velocity gradient variance (s−2), and (d) buoyancy gradient variance (s−4) computed from the model solutions. The mean values over the selected 656 events ( appendix D) and the corresponding analytical solutions [Eqs. (16a)(16d)] are shown in magenta and green, respectively. The normalized geostrophic (dashed magenta) and ageostrophic (dot–dashed magenta) vorticities are also plotted in (b). The dashed black line in (c) indicates the asymptotic prediction for the velocity gradient variance in terms of the divergence and vorticity variances [Eq. (10)]. Solid gray lines indicate the mean values for each of the seven particle release experiments, and the shaded blue envelopes indicate a representative event spread (one standard deviation) for each experiment. Similar variabilities between and during release experiments were also found for the various terms in Figs. 7, 8, and 11 (not shown).

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 δ,ζ, and |huh|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).

Fig. 10.

(a)–(c) As in Figs. 7 and 8, but for the terms that can be computed from the LASER drifter observations. (d)–(f) As in Fig. 9, but computed based on the LASER drifter observations. The fields in (a) and (b) are normalized by f2 and in (d) and (e) by f. The fields in (c) and (f) are in s−3 and s−2, respectively. All quantities are averaged over the selected 705 events ( appendix D). The shaded envelopes in (a)–(c) denote the estimated sampling error based on vorticity/divergence errors of 0.5f and velocity gradient variance error of 5.3×109s2 (section 3b and  appendix E), dτ of 15 min, and assuming 10% of the total selected events were independent. The shaded envelopes in (d)–(f) denotes the spread (one standard deviation) of all selected events. A bias correction was applied to all fields ( appendix E). Note that the D/Dt terms do not include the vertical advection because the LASER drifters were designed to stay near the surface. The similarity in (a) and (c) between the Fhor terms (solid green) and the corresponding asymptotic predictions [Eqs. (12b), (12d); dashed blue] demonstrates that the asymptotic model (section 2) adequately describes the advective tendency terms in the LASER measurements. The good agreement in (d)–(f) between the Lagrangian evolution of δ,ζ, and |huh|2 (solid black) and the analytical solutions that ignore pressure and Coriolis effects [Eqs. (16b)(16d); dashed blue] suggest that TTW is a plausible dominant mechanism for the initiation of submesoscale frontogenesis during LASER.

Fig. 10.

(a)–(c) As in Figs. 7 and 8, but for the terms that can be computed from the LASER drifter observations. (d)–(f) As in Fig. 9, but computed based on the LASER drifter observations. The fields in (a) and (b) are normalized by f2 and in (d) and (e) by f. The fields in (c) and (f) are in s−3 and s−2, respectively. All quantities are averaged over the selected 705 events ( appendix D). The shaded envelopes in (a)–(c) denote the estimated sampling error based on vorticity/divergence errors of 0.5f and velocity gradient variance error of 5.3×109s2 (section 3b and  appendix E), dτ of 15 min, and assuming 10% of the total selected events were independent. The shaded envelopes in (d)–(f) denotes the spread (one standard deviation) of all selected events. A bias correction was applied to all fields ( appendix E). Note that the D/Dt terms do not include the vertical advection because the LASER drifters were designed to stay near the surface. The similarity in (a) and (c) between the Fhor terms (solid green) and the corresponding asymptotic predictions [Eqs. (12b), (12d); dashed blue] demonstrates that the asymptotic model (section 2) adequately describes the advective tendency terms in the LASER measurements. The good agreement in (d)–(f) between the Lagrangian evolution of δ,ζ, and |huh|2 (solid black) and the analytical solutions that ignore pressure and Coriolis effects [Eqs. (16b)(16d); dashed blue] suggest that TTW is a plausible dominant mechanism for the initiation of submesoscale frontogenesis during LASER.

A good agreement is found between the frontogenetic tendencies of divergence and velocity gradient variance Fhorδ,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, |huh|2~δ2+ζ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 ζ,δ, and |huh|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 wb 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 wb and the Fvertb term [Eq. (1)]. To see this we evaluate the leading-order vertical tendency rate

 
Tvertb=Fvertb/|hb|2=εwxsb,
(20)

subject to the chosen scaling and coordinate orientation [Eqs. (7), (A3)], where sb=bx/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 wb>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 |hb|, 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 (1/4)|hb|hmll. 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(1+δ0τ),l=l0(1+δ0τ) and

 
τAPE=14τ(|hb|hmll)=14hml0l0δ0|hb0|wb,
(21)

where l0 and hml0 are the initial frontal width and mixed layer depth, respectively, and wb>0 for any initial convergence magnitude (δ0<0). For realistic fronts, the mixed layer does not completely restratify during frontogenesis and, therefore, hml varies much slower than l and |hb| (FvertbFhorb, 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

 
wb~hml0l(τ)δ(τ)|hb(τ)|.
(22)

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 wb, 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 |hb|2 during the last hour (hour 7).

The Lagrangian evolution of wb, 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 wb>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 wb (Fig. 11).

Fig. 11.

The model-based Lagrangian evolution of the vertical buoyancy fluxes wb (W kg−1), the vertical advection [Fvertb in Eq. (1); s−5], and the scaling [Eq. (22); W kg−1]. Quantities are averaged over the hundred frontogenetic events ( appendix D) with the largest and weakest final |hb| magnitudes (strong and weak, respectively). The δ(τ) and |hb(τ)| are interpolated directly on the particles, and l(τ) is assumed to follow l(τ)=l0(1+δ0τ) (section 6). For the scaling estimates we pick hml0=50m (consistent with Barkan et al. 2017a), l0 = 4000 m, and the constants of proportionality used to match the initial values are 0.10 and 0.18 for the strong and weak fronts, respectively. The vertical advection terms have been multiplied by 1011 and 1012 so that they can be plotted on the same figure with the buoyancy flux terms. For information about the variability between events and release experiments refer to Figs. 9 and D1.

Fig. 11.

The model-based Lagrangian evolution of the vertical buoyancy fluxes wb (W kg−1), the vertical advection [Fvertb in Eq. (1); s−5], and the scaling [Eq. (22); W kg−1]. Quantities are averaged over the hundred frontogenetic events ( appendix D) with the largest and weakest final |hb| magnitudes (strong and weak, respectively). The δ(τ) and |hb(τ)| are interpolated directly on the particles, and l(τ) is assumed to follow l(τ)=l0(1+δ0τ) (section 6). For the scaling estimates we pick hml0=50m (consistent with Barkan et al. 2017a), l0 = 4000 m, and the constants of proportionality used to match the initial values are 0.10 and 0.18 for the strong and weak fronts, respectively. The vertical advection terms have been multiplied by 1011 and 1012 so that they can be plotted on the same figure with the buoyancy flux terms. For information about the variability between events and release experiments refer to Figs. 9 and D1.

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 (~kh2) 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 [δ~ζ; 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 |hb|2,|hu|2,ζ,δ, 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 (δζ) 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

 
Ro2(ut+uux+εRoυuy+εwuz)υRo=ΦxRo,
(A1a)
 
Ro(vt+uυx+εRoυυy+εwυz)+u=εRoΦy,
(A1b)
 
Ro2λ2(wt+uwx+εRoυwy+wwz)=1Ro(Φz+b),
(A1c)
 
Ro(bt+ubx)+ευby+Roεwbz=0,and
(A1d)
 
Ro(ux+wz)+ευy=0,
(A1e)

where all fields are nondimensional; (u,υ,w) are the velocities in the (x^,y^,z^) directions; b is the buoyancy; Φ is the geopotential; subscripts denote derivatives; and the nondimensional numbers Ro,ε, and λ are described in Eq. (8).

The buoyancy frontogenetic tendency equation Eq. (1) is derived by taking the horizontal gradient h of Eq. (A1d) and then multiplying it by hb, 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 Ωmj, namely,

 
um,j=Smj+Ωmj,Smj=12(um,j+uj,m),Ωmj=12(um,juj,m),
(A2)

the buoyancy frontogenetic equation can be simplified to

 
DDt12(b,jb,j)=Skjb,kb,jFhorbN2w,jb,jFvertb,
(A3)

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

 
Skj=[Roux12(υx+εRouy)12(υx+εRouy)ευy],b,kb,j=(bx2εbxbyεbxbyε2by2),
(A4)

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 ±(uxυy). The corresponding nondimensional Fhorb and Fvertb are

 
Fhorb=Ro[uxbx2+εRoυxbxby+ε2(uybxby+εRoυyby2)],
(A5)

and

 
Fvertb=RoεN2(wxbx+ε2wyby),
(A6)

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

 
DDt12(ui,jui,j)=Skjui,kui,jFhoruΛiw,jui,jFvertuSijΦ,ijFpresu,
(A7)

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

 
ui,kui,j=(Ro2ux2+υx2Ro2εuxuy+ευxυyRo2εuxuy+ευxυyRo2ε2uy2+ε2υy2).
(A8)

The corresponding nondimensional Fhoru, Fvertu, and Fpresu become

 
Fhoru=Ro[(ux+εRoυy)υx2+(Ro2ux)ux2+ε(Roux+ευy)uyυx+Roε2(Roux+ευy)uy2+(ε3Roυy)υy2],
(A9)
 
Fvertu=Roε[υz(wxυx+ε2wyυy)+Ro2uz(wxux+ε2wyuy)],
(A10)

and

 
Fpresu=uxΦxxεRo[ε2υyΦyy+(υx+Roεuy)Φxy].
(A11)

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ζ,Fvertζ,Fhorδ,Fvertδ, and pressureδ are

 
Fhorζ=Ro[ux(υxεRouy)+υy(εRoυxε2uy)],
(A12a)
 
Fvertζ=Roε(εRouzwyυzwx),
(A12b)
 
Fhorδ=Ro2(ux2+2εRoυxuy+ε2Ro2υy2),
(A12c)
 
Fvertδ=Roε(Rouzwx+ευzwy),and
(A12d)
 
Fpresδ=1Ro(Φxx+ε2Φyy).
(A12e)

APPENDIX B

The Asymptotic Model in Coordinate Invariant Form

In the asymptotic limit ε1,λ1,Ro~O(1) the various frontogenetic tendency terms (F) and the pressure term in the velocity gradient evolution equation (Fpresu) take the form

 
Fhorbuxbx2+H.O.Tδb,jb,j+b,jQ˜jH.O.T,
(B1a)
 
Fhoruux(ux2+υx2)+H.O.Tδui,jui,j+δdet(ui,j)H.O.T,
(B1b)
 
Fhorζuxυx+H.O.Tδζ,
(B1c)
 
Fhorδux2+H.O.Tδ2+2det(ui,j)H.O.T,
(B1d)
 
FpresδuxΦxx+H.O.Tδh2Φ+f[J(u,ug)+J(υ,υg)]H.O.T,
(B1e)

where H.O.T denotes higher-order terms; det denotes the determinant; fug=k^×Φ; the Jacobian is defined as J(a,b)=axbybxay; and Q˜j=[J(υ,b),J(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 δ=Roux+ευy, we decompose the nondimensional horizontal velocity field into its divergent and rotational components (i.e., Helmholtz decomposition), namely,

 
u=Roϕx+εψy,and
(B2a)
 
υ=Roεϕyψx,
(B2b)

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

 
uxϕxx,
(B3)

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

It is useful (section 2a) to contrast Eqs. (B1) with the dynamical regime λ1,Ro~ε1. In this dynamical regime

 
FhorbRouxbx2Roυxbxby+H.O.TSkjb,kb,j,
(B4a)
 
FhoruRo(ux+υy)υx2+H.O.Tδui,jui,j+δdet(ui,j)H.O.T,
(B4b)
 
FhorζRo(ux+υy)υx+H.O.Tδζ,
(B4c)
 
FhorδRo2(ux22υxuy+υy2)uk,iui,k,and
(B4d)
 
FpresδuxΦxxυxΦxyui,jΦ,ij
(B4e)

where the nondimensional tensors Skj and b,kb,j are written in Eq. (A4). Note that irrespective of the scaling assumption Fhorζ=δζ, and that only the divergent part of the velocity field contributes to Fhoru [the H.O.T in Eqs. (B4b) and (B1b) contain δ].

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

 
τδ=δ2g(τ)f2,
(C1)

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

 
wττ+wf2g(τ)=0.
(C2)

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(τ)=a2τn then Eq. (C2) can be solved using power series solutions, by setting w=k=0Akτk. After some algebra we get

 
k=0(k+2)(k+1)Ak+2τk+f2a2k=nAknτk=0.
(C3)

For illustration we solve Eq. (C3) for n=2 [i.e., Eq. (19)], whereby it becomes

 
2A2+6A3τ+k=2[(k+2)(k+1)Ak+2+a2f2Ak2]τk=0,
(C4)

with the recursion relation

 
Ak+4=a2f2Ak(k+4)(k+3),k=2,,,
(C5)

and with the initial conditions A2=w(0)/2,A3=w(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

To close the balances in Eqs. (1)(4) we carry out the following procedure:

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

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

  3. 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 (FhorLagFhorEul), where the FhorEul terms are computed directly from the model and are interpolated onto the particle locations. The FhorLag terms are computed as 
    FhorLagδ=τδhDhadvDtuh,FhorLagζ=τζh×DhadvDtuh,
    (D1)

and similarly for |huh|2 and |hb|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 δ,ζ,|hb|2, and |huh|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 τ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 τ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 δ,ζ, and |huh|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).

Fig. D1.

A representative illustration of the model-based Lagrangian evolution of seven particles (gray lines) that converged onto the same submesoscale feature. The mean convergence values and corresponding analytical solution [Eq. (16d)] are plotted in magenta and green, respectively. A similar agreement was found for the other gradient fields |hb|2,|huh|2, and ζ (not shown).

Fig. D1.

A representative illustration of the model-based Lagrangian evolution of seven particles (gray lines) that converged onto the same submesoscale feature. The mean convergence values and corresponding analytical solution [Eq. (16d)] are plotted in magenta and green, respectively. A similar agreement was found for the other gradient fields |hb|2,|huh|2, and ζ (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) δ(τ)corrected=δ(τ)×δ(τ)particle/δ(τ)triangle, where δ(τ) is the uncorrected divergence field (solid black line in Fig. E3d), and δ(τ)particle and δ(τ)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.

Fig. E1.

The model-based Lagrangian evolution of triangle aspect ratio (the ratio of the triangle height to triangle base), averaged over the 4480 frontogenetic events computed with the algorithm discussed in  appendix D. The shaded envelope denotes the spread (one standard deviation) over all selected events. The aspect ratio used in the drifter measurements is denoted by the dashed black line. The minimum allowable triangle scale (the RMS distance between the triangle particles) is 1000 m, as in the drifter measurements. An equidistant triangle has an aspect ratio of 3/2.

Fig. E1.

The model-based Lagrangian evolution of triangle aspect ratio (the ratio of the triangle height to triangle base), averaged over the 4480 frontogenetic events computed with the algorithm discussed in  appendix D. The shaded envelope denotes the spread (one standard deviation) over all selected events. The aspect ratio used in the drifter measurements is denoted by the dashed black line. The minimum allowable triangle scale (the RMS distance between the triangle particles) is 1000 m, as in the drifter measurements. An equidistant triangle has an aspect ratio of 3/2.

Fig. E2.

(a)–(c) A comparison between the model-based D/Dt and Fhor terms in Eqs. (2)(4) computed on particles (solid lines) and using the triangle method (dashed lines; section 3b), with a minimum aspect ratio of 3/4 and a minimum scale of 1000 m, as in the drifter-based estimate (Figs. 10, E3). The asymptotic predictions for Fhor [Eqs. (12b), (12d)] are in blue. All fields in (a) and (b) are normalized by f2 and in (c) are in s−3. (d)–(f) The corresponding Lagrangian evolution of divergence normalized by f, vorticity normalized by f, and velocity gradient variance (s−2) computed on particles (black lines) and using the triangle method (blue lines). The thin blue lines in (d)–(f) show the triangle-based estimates with a minimum triangle scale of 1000 m and a minimum triangle aspect ratio of 3/40. All fields are averaged over O(100) events when computed on particles and over O(1000) events when computed with the triangle method. Note that the D/Dt terms computed using the triangle method do not include the vertical advection, to be consistent with the LASER drifters. For information about the variability between events and release experiments refer to Figs. 9 and D1. The triangle method underestimates magnitudes for all terms.

Fig. E2.

(a)–(c) A comparison between the model-based D/Dt and Fhor terms in Eqs. (2)(4) computed on particles (solid lines) and using the triangle method (dashed lines; section 3b), with a minimum aspect ratio of 3/4 and a minimum scale of 1000 m, as in the drifter-based estimate (Figs. 10, E3). The asymptotic predictions for Fhor [Eqs. (12b), (12d)] are in blue. All fields in (a) and (b) are normalized by f2 and in (c) are in s−3. (d)–(f) The corresponding Lagrangian evolution of divergence normalized by f, vorticity normalized by f, and velocity gradient variance (s−2) computed on particles (black lines) and using the triangle method (blue lines). The thin blue lines in (d)–(f) show the triangle-based estimates with a minimum triangle scale of 1000 m and a minimum triangle aspect ratio of 3/40. All fields are averaged over O(100) events when computed on particles and over O(1000) events when computed with the triangle method. Note that the D/Dt terms computed using the triangle method do not include the vertical advection, to be consistent with the LASER drifters. For information about the variability between events and release experiments refer to Figs. 9 and D1. The triangle method underestimates magnitudes for all terms.

Fig. E3.

As in Fig. 10, but without the bias correction (Fig. E2). Because the triangle method underestimates all quantities, the analytical solutions [Eqs. (16b)(16d)] in (d)–(f) (dashed blue) only match the corresponding field values that are one standard deviation larger (for ζ and |huh|2) or smaller (for δ) than the mean.

Fig. E3.

As in Fig. 10, but without the bias correction (Fig. E2). Because the triangle method underestimates all quantities, the analytical solutions [Eqs. (16b)(16d)] in (d)–(f) (dashed blue) only match the corresponding field values that are one standard deviation larger (for ζ and |huh|2) or smaller (for δ) than the mean.

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

REFERENCES
Andersson
,
A.
,
K.
Fennig
,
C.
Klepp
,
S.
Bakan
,
H.
Graßl
, and
J.
Schulz
,
2010
:
The Hamburg ocean atmosphere parameters and fluxes from satellite data–HOAPS-3
.
Earth Syst. Sci. Data
,
2
,
215
234
, https://doi.org/10.5194/essd-2-215-2010.
Bachman
,
S. D.
, and
J. R.
Taylor
,
2016
:
Numerical simulations of the equilibrium between eddy-induced restratification and vertical mixing
.
J. Phys. Oceanogr.
,
46
,
919
935
, https://doi.org/10.1175/JPO-D-15-0110.1.
Barkan
,
R.
,
K. B.
Winters
, and
S. G. L.
Smith
,
2015
:
Energy cascades and loss of balance in a reentrant channel forced by wind stress and buoyancy fluxes
.
J. Phys. Oceanogr.
,
45
,
272
293
, https://doi.org/10.1175/JPO-D-14-0068.1.
Barkan
,
R.
,
J. C.
McWilliams
,
A. F.
Shchepetkin
,
M. J.
Molemaker
,
L.
Renault
,
A.
Bracco
, and
J.
Choi
,
2017a
:
Submesoscale dynamics in the northern Gulf of Mexico. Part I: Regional and seasonal characterization, and the role of river outflow
.
J. Phys. Oceanogr.
,
47
,
2325
2346
, https://doi.org/10.1175/JPO-D-17-0035.1.
Barkan
,
R.
,
J. C.
McWilliams
,
M. J.
Molemaker
,
J.
Choi
,
K.
Srinivasan
,
A. F.
Shchepetkin
, and
A.
Bracco
,
2017b
:
Submesoscale dynamics in the northern Gulf of Mexico. Part II: Temperature–salinity compensation, and cross-shelf transport processes
.
J. Phys. Oceanogr.
,
47
,
2347
2360
, https://doi.org/10.1175/JPO-D-17-0040.1.
Barkan
,
R.
,
K. B.
Winters
, and
J. C.
McWilliams
,
2017c
:
Stimulated imbalance and the enhancement of eddy kinetic energy dissipation by internal waves
.
J. Phys. Oceanogr.
,
47
,
181
198
, https://doi.org/10.1175/JPO-D-16-0117.1.
Boccaletti
,
G.
,
R.
Ferrari
, and
B.
Fox-Kemper
,
2007
:
Mixed layer instabilities and restratification
.
J. Phys. Oceanogr.
,
37
,
2228
2250
, https://doi.org/10.1175/JPO3101.1.
Callies
,
J.
, and
R.
Ferrari
,
2013
:
Interpreting energy and tracer spectra of upper-ocean turbulence in the submesoscale range (1–200 km)
.
J. Phys. Oceanogr.
,
43
,
2456
2474
, https://doi.org/10.1175/JPO-D-13-063.1.
Callies
,
J.
,
R.
Ferrari
,
J. M.
Klymak
, and
J.
Gula
,
2015
: Seasonality in submesoscale turbulence. Nat. Commun., 6, 6862, https://doi.org/10.1038/ncomms7862.
Capet
,
X.
,
J. C.
McWilliams
,
M. J.
Molemaker
, and
A. F.
Shchepetkin
,
2008a
:
Mesoscale to submesoscale transition in the California Current System. Part I: Flow structure, eddy flux, and observational tests
.
J. Phys. Oceanogr.
,
38
,
29
43
, https://doi.org/10.1175/2007JPO3671.1.
Capet
,
X.
,
J. C.
McWilliams
,
M. J.
Molemaker
, and
A. F.
Shchepetkin
,
2008b
:
Mesoscale to submesoscale transition in the California Current System. Part II: Frontal processes
.
J. Phys. Oceanogr.
,
38
,
44
64
, https://doi.org/10.1175/2007JPO3672.1.
Capet
,
X.
,
J. C.
McWilliams
,
M. J.
Molemaker
, and
A. F.
Shchepetkin
,
2008c
:
Mesoscale to submesoscale transition in the California Current System. Part III: Energy balance and flux
.
J. Phys. Oceanogr.
,
38
,
2256
2269
, https://doi.org/10.1175/2008JPO3810.1.
Choi
,
J.
,
A.
Bracco
,
R.
Barkan
,
J. C.
McWilliams
,
A. F.
Shchepetkin
, and
M. J.
Molemaker
,
2017
:
Submesoscale dynamics in the northern Gulf of Mexico. Part III: Lagrangian implications
.
J. Phys. Oceanogr.
,
47
,
2361
2376
, https://doi.org/10.1175/JPO-D-17-0036.1.
Cleveland
,
W. S.
, and
S. J.
Devlin
,
1988
:
Locally weighted regression: an approach to regression analysis by local fitting
.
J. Amer. Stat. Assoc.
,
83
,
596
610
, https://doi.org/10.1080/01621459.1988.10478639.
D’Asaro
,
E. A.
, and Coauthors
,
2018
:
Ocean convergence and the dispersion of flotsam
.
Proc. Natl. Acad. Sci. USA
,
115
,
1162
1167
, https://doi.org/10.1073/pnas.1718453115.
Fox-Kemper
,
B.
,
R.
Ferrari
, and
R.
Hallberg
,
2008
:
Parameterization of mixed layer eddies. Part I: Theory and diagnosis
.
J. Phys. Oceanogr.
,
38
,
1145
1165
, https://doi.org/10.1175/2007JPO3792.1.
Garrett
,
C.
, and
W. H.
Munk
,
1972
:
Space-time scales of internal waves
.
Geophys. Fluid Dyn.
,
3
,
225
264
, https://doi.org/10.1080/03091927208236082.
Gula
,
J.
,
M. J.
Molemaker
, and
J. C.
McWilliams
,
2014
:
Submesoscale cold filaments in the gulf stream
.
J. Phys. Oceanogr.
,
44
,
2617
2643
, https://doi.org/10.1175/JPO-D-14-0029.1.
Hamlington
,
P. E.
,
L. P.
Van Roekel
,
B.
Fox-Kemper
,
K.
Julien
, and
G. P.
Chini
,
2014
:
Langmuir–submesoscale interactions: Descriptive analysis of multiscale frontal spindown simulations
.
J. Phys. Oceanogr.
,
44
,
2249
2272
, https://doi.org/10.1175/JPO-D-13-0139.1.
Haza
,
A.
, and Coauthors
,
2018
:
Drogue-loss detection for surface drifters during the Lagrangian Submesoscale Experiment (LASER)
.
J. Atmos. Oceanic Technol.
,
35
,
705
725
, https://doi.org/10.1175/JTECH-D-17-0143.1.
Hoskins
,
B. J.
,
1975
:
The geostrophic momentum approximation and the semi-geostrophic equations
.
J. Atmos. Sci.
,
32
,
233
242
, https://doi.org/10.1175/1520-0469(1975)032<0233:TGMAAT>2.0.CO;2.
Hoskins
,
B.
,
1982
:
The mathematical theory of frontogenesis
.
Annu. Rev. Fluid Mech.
,
14
,
131
151
, https://doi.org/10.1146/annurev.fl.14.010182.001023.
Hoskins
,
B.
, and
F.
Bretherton
,
1972
:
Atmospheric frontogenesis models: Mathematical formulation and solution
.
J. Atmos. Sci.
,
29
,
11
37
, https://doi.org/10.1175/1520-0469(1972)029<0011:AFMMFA>2.0.CO;2.
Lapeyre
,
G.
, and
P.
Klein
,
2006
:
Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory
.
J. Phys. Oceanogr.
,
36
,
165
176
, https://doi.org/10.1175/JPO2840.1.
Large
,
W. G.
, and
S.
Yeager
,
2009
:
The global climatology of an interannually varying air–sea flux data set
.
Climate Dyn.
,
33
,
341
364
, https://doi.org/10.1007/s00382-008-0441-3.
Large
,
W. G.
,
J. C.
McWilliams
, and
S. C.
Doney
,
1994
:
Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization
.
Rev. Geophys.
,
32
,
363
403
, https://doi.org/10.1029/94RG01872.
Mahadevan
,
A.
,
2016
:
The impact of submesoscale physics on primary productivity of plankton
.
Annu. Rev. Mar. Sci.
,
8
,
161
184
, https://doi.org/10.1146/annurev-marine-010814-015912.
Mason
,
E.
,
J.
Molemaker
,
A. F.
Shchepetkin
,
F.
Colas
,
J. C.
McWilliams
, and
P.
Sangrà
,
2010
:
Procedures for offline grid nesting in regional ocean models
.
Ocean Modell.
,
35
,
1
15
, https://doi.org/10.1016/j.ocemod.2010.05.007.
McWilliams
,
J. C.
,
2016
:
Submesoscale currents in the ocean
.
Proc. Roy. Soc. London
,
472A
, 20160117, https://doi.org/10.1098/rspa.2016.0117.
McWilliams
,
J. C.
,
2017
:
Submesoscale surface fronts and filaments: Secondary circulation, buoyancy flux, and frontogenesis
.
J. Fluid Mech.
,
823
,
391
432
, https://doi.org/10.1017/jfm.2017.294.
McWilliams
,
J. C.
,
2018
:
Surface wave effects on submesoscale fronts and filaments
.
J. Fluid Mech.
,
843
,
479
517
, https://doi.org/10.1017/jfm.2018.158.
McWilliams
,
J. C.
, and
M. J.
Molemaker
,
2011
:
Baroclinic frontal arrest: A sequel to unstable frontogenesis
.
J. Phys. Oceanogr.
,
41
,
601
619
, https://doi.org/10.1175/2010JPO4493.1.
McWilliams
,
J. C.
, and
B.
Fox-Kemper
,
2013
:
Oceanic wave-balanced surface fronts and filaments
.
J. Fluid Mech.
,
730
,
464
490
, https://doi.org/10.1017/jfm.2013.348.
McWilliams
,
J. C.
,
M. J.
Molemaker
, and
E. I.
Olafsdottir
,
2009
:
Linear fluctuation growth during frontogenesis
.
J. Phys. Oceanogr.
,
39
,
3111
3129
, https://doi.org/10.1175/2009JPO4186.1.
McWilliams
,
J. C.
,
J.
Gula
,
M. J.
Molemaker
,
L.
Renault
, and
A. F.
Shchepetkin
,
2015
:
Filament frontogenesis by boundary layer turbulence
.
J. Phys. Oceanogr.
,
45
,
1988
2005
, https://doi.org/10.1175/JPO-D-14-0211.1.
Molemaker
,
M. J.
,
J. C.
McWilliams
, and
X.
Capet
,
2010
:
Balanced and unbalanced routes to dissipation in an equilibrated Eady flow
.
J. Fluid Mech.
,
654
,
35
63
, https://doi.org/10.1017/S0022112009993272.
Novelli
,
G.
,
C. M.
Guigand
,
C.
Cousin
,
E. H.
Ryan
,
N. J.
Laxague
,
H.
Dai
,
B. K.
Haus
, and
T. M.
Özgökmen
,
2017
:
A biodegradable surface drifter for ocean sampling on a massive scale
.
J. Atmos. Oceanic Technol.
,
34
,
2509
2532
, https://doi.org/10.1175/JTECH-D-17-0055.1.
Oey
,
L.-Y.
,
T.
Ezer
, and
H.-C.
Lee
,
2005
: Loop current, rings and related circulation in the Gulf of Mexico: A review of numerical models and future challenges. Circulation in the Gulf of Mexico: Observations and Models, Geophys. Monogr., Vol. 161, Amer. Geophys. Union, 31–56.
Ohlmann
,
J.
,
M.
Molemaker
,
B.
Baschek
,
B.
Holt
,
G.
Marmorino
, and
G.
Smith
,
2017
:
Drifter observations of submesoscale flow kinematics in the coastal ocean
.
Geophys. Res. Lett.
,
44
,
330
337
, https://doi.org/10.1002/2016GL071537.
Poje
,
A. C.
, and Coauthors
,
2014
:
Submesoscale dispersion in the vicinity of the Deepwater Horizon spill
.
Proc. Natl. Acad. Sci. USA
,
111
,
12 693
12 698
, https://doi.org/10.1073/pnas.1402452111.
Risien
,
C. M.
, and
D. B.
Chelton
,
2008
:
A global climatology of surface wind and wind stress fields from eight years of QuikSCAT scatterometer data
.
J. Phys. Oceanogr.
,
38
,
2379
2413
, https://doi.org/10.1175/2008JPO3881.1.
Salmon
,
R.
,
1980
:
Baroclinic instability and geostrophic turbulence
.
Geophys. Astrophys. Fluid Dyn.
,
15
,
167
211
, https://doi.org/10.1080/03091928008241178.
Shchepetkin
,
A. F.
, and
J. C.
McWilliams
,
2005
:
The Regional Oceanic Modeling System: A split-explicit, free-surface, topography-following-coordinate oceanic model
.
Ocean Modell.
,
9
,
347
404
, https://doi.org/10.1016/j.ocemod.2004.08.002.
Srinivasan
,
K.
,
J. C.
McWilliams
,
L.
Renault
,
H. G.
Hristova
,
J.
Molemaker
, and
W. S.
Kessler
,
2017
:
Topographic and mixed layer submesoscale currents in the near-surface southwestern tropical pacific
.
J. Phys. Oceanogr.
,
47
,
1221
1242
, https://doi.org/10.1175/JPO-D-16-0216.1.
Sullivan
,
P. P.
, and
J. C.
McWilliams
,
2018
:
Frontogenesis and frontal arrest of a dense filament in the oceanic surface boundary layer
.
J. Fluid Mech.
,
837
,
341
380
, https://doi.org/10.1017/jfm.2017.833.
Suzuki
,
N.
,
B.
Fox-Kemper
,
P. E.
Hamlington
, and
L. P. V.
Roekel
,
2016
:
Surface waves affect frontogenesis
.
J. Geophys. Res. Oceans
,
121
,
3597
3624
, https://doi.org/10.1002/2015JC011563.
Thomas
,
L. N.
,
2017
:
On the modifications of near-inertial waves at fronts: implications for energy transfer across scales
.
Ocean Dyn.
,
67
,
1335
1350
, https://doi.org/10.1007/s10236-017-1088-6.
Thomas
,
L. N.
,
A.
Tandon
, and
A.
Mahadevan
,
2008
: Submesoscale processes and dynamics. Ocean Modeling in an Eddying Regime, Geophys. Monogr., Vol. 177, Amer. Geophys. Union, 17–38.
Whitt
,
D. B.
, and
L. N.
Thomas
,
2015
:
Resonant generation and energetics of wind-forced near-inertial motions in a geostrophic flow
.
J. Phys. Oceanogr.
,
45
,
181
208
, https://doi.org/10.1175/JPO-D-14-0168.1.
Winters
,
K. B.
,
P.
Lombard
,
J.
Riley
, and
E. A.
D’Asaro
,
1995
:
Available potential energy and mixing in density stratified fluids
.
J. Fluid Mech.
,
289
,
115
128
, https://doi.org/10.1017/S002211209500125X.

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 kh2 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˜[huh(kh)]=E˜[ζ(kh)]=|kh|2E˜[u(kh)]~kh1, in agreement with the enstrophy cascade in geostrophic turbulence (Salmon 1980).

3

Choosing another limit would mean that the Fvert terms in Eqs. (1)(4) would appear at different higher orders of the asymptotic expansion but would not change any of the results presented in this paper, which are concerned with the leading-order balances.

4

In this scaling regime |huh|2~ζ2, and so Tu=Tζ=δ by construction.

5

When ζag=+a2f, 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).