Characterizing the Role of Non-Linear Interactions in the Transition to Submesoscale Dynamics at a Dense Filament

Ocean dynamics at the submesoscale play a key role in mediating upper-ocean energy dissipation and dispersion of tracers. Observations of ocean currents from synoptic mesoscale surveys at submesoscale resolution (250 m–100 km) from a novel airborne instrument (MASS DoppVis) reveal that the kinetic energy spectrum in the California Current System is nearly continuous from 100 km to sub-kilometer scales, with a k −2 spectral slope. Although there is not a transition in the kinetic energy spectral slope, there is a transition in the dynamics to non-linear ageostrophic interactions at scales of 𝐴𝐴  (1 km). Kinetic energy transfer across spatial scales is enabled by interactions between the rotational and divergent components of the flow field at the submesoscale. Kinetic energy flux is patchy and localized at submesoscale fronts. Kinetic energy is transferred both downscale and upscale from 1 km in the observations of a cold filament.

break down at the submesoscale. At the larger end of the submesoscale, the surface quasigeostrophy framework presupposes that surface density fronts modify geostrophic balance (Klein & Lapeyre, 2009), while other theoretical results emphasize the role of non-linear advection in submesoscale dynamics (Barkan et al., 2019).
In this work we characterize the transition to submesoscale dynamics at scales smaller than 10 km and provide observational analysis of the kinetic energy cascade that has been hypothesized from models and theory. We observe submesoscale ocean surface velocity using remote sensing from airplanes during the submesoscale ocean dynamics experiment (S-MODE) field campaign (Farrar et al., 2020). We find substantial kinetic energy at the submesoscale, with a kinetic energy spectral slope that is nearly continuous from 100 to 1 km spatial scales. The dynamics that result in the spatial distribution of kinetic energy at the submesoscale are diagnosed through analysis of velocity cross spectra. These reveal that non-linear interactions between balanced and unbalanced dynamics contribute to submesoscale energy and illuminate the dynamics influencing upper-ocean velocity gradient distributions.

Remote Sensing
The observations used in this study were collected by the DoppVis instrument (Lenain et al., 2023), a new sensor that is part of the Modular Aerial Sensing System (MASS; Melville et al., 2016), that infers currents from optical observations of the spatio-temporal evolution, that is, dispersion relationship, of surface waves. This method infers the depth-resolved Lagrangian current in the upper ocean. Here, we use the depth averaged current over the upper 2 m. Details about the DoppVis instrument are available in Lenain et al. (2023). The instrument package was installed on a Twin Otter DH-6 aircraft, flying at constant altitude above mean sea level (hereafter, altitude), with a flight profile consisting of repeated reciprocal straight tracks. Consistency between the reciprocal passes is used to validate velocity measurements. Velocity observations are binned to 256 m or 500 m prior to analysis.
Observations from two field campaigns are considered in this study. The first field campaign sampled across a cold filament approximately 70 nautical miles offshore of California, as part of the NASA S-MODE program (Farrar et al., 2020). This region is subsequently referred to as the "filament region" and is the focus of this study. These observations occurred on 3 November 2021 from 18:23 to 23:33 UTC while flying at approximately 500 m altitude and on 5 November 2021 from 22:40 to 23:00 UTC while flying at 940 m altitude (Figures 1a-1c).
The higher altitude flight on November 5 enables collection of multiple data points in the cross-swath direction, resulting in a 1.5 km wide swath that is used to compute velocity gradients using central differences. The observations from November 5 are binned at 256 m prior to analysis.
The second field campaign collected observations across two counter-rotating eddies approximately 45 nautical miles offshore of San Diego on 19 May 2021 from 20:56 to 23:26 UTC (Lenain et al., 2023). This region is referred to as the "eddy region." Observations using a vessel mounted acoustic Doppler current profiler (ADCP) were collected under the long northwest-to-southeast leg of the DoppVis observations from 19 May 2021 10:00 to 20 May 2021 13:00 UTC.
Confidence intervals of velocity gradient results including kinetic energy flux and fronotgenesis are estimated using a bootstrapped confidence interval and a velocity error of 0.05 m s −1 (Lenain et al., 2023).
The data analyzed to generate each figure is noted in the figure caption and in Table S1 in Supporting Information S1.

Spectra
We analyze both the kinetic energy spectrum (̂( ) ) and the cross spectrum (̂( ) ) with 95% confidence intervals calculated following Bendat and Piersol (2011). Both the kinetic energy spectra and the cross-spectrum between along-track and across-track velocity are computed using Welch's method with Hanning windows.

Kinetic Energy Spectrum
The multi-scale nature of the flow is quantified using energy spectra, which can also be used to make inferences about the dominant dynamics governing the flow (Callies & Ferrari, 2013). The filament region is more energetic than the eddy region ( Figure 1d), with approximately twice the amount of energy at nearly all spatial scales sampled. The kinetic energy spectra of the DoppVis observations have slopes that are approximately k −2 ( Figure 1d). The observed kinetic energy spectrum crossing the eddies has magnitude and spectral slope similar to that of the spectra from currents (15 m depth) taken with a vessel mounted ADCP on a nearby transect on the same day for 5-100 km scales.
This analysis extends the observations to smaller spatial scales than have been observed previously. Notably, these scales are smaller than those resolved by state-of-the-art global and regional models. As an example, we show the kinetic energy spectrum from a 2 km grid spacing MITgcm regional model of the California current system (Mazloff et al., 2020) (Figure 1d, red line). This model is forced with ERA5 atmospheric state, Hybrid Coordinate Ocean Model + Navy Coupled Ocean Data Assimilation boundary conditions, and both local and remote tides. The effective resolution of this model is 20 km with the velocity spectrum falling off steeply below that scale due to grid scale dissipation. Even at larger scales, both regions are more energetic than the 2 km grid spacing ocean model of the same region (the eddy region is 5 times more energetic). The discrepancy between model and observations at lower wavenumbers is likely caused by an inverse cascade of submesoscale energy energizing surface mesoscale features in ways that are not represented in the model (Lévy et al., 2001;Mahadevan & Tandon, 2006) and by biased observational sampling toward more energetic features. It is important to note that only the larger end of submesoscale dynamics are resolved by 2 km models (Sinha et al., 2022;Su et al., 2018). This is especially important to keep in mind when considering cross-scale energy fluxes that may be modified by dynamics at small spatial scales.
The observed kinetic energy spectral slopes are consistent with previous observations from this region: a comprehensive analysis of surface velocities measured from vessel-mounted ADCPs in the California Current region from 1993 to 2004 found that the kinetic energy spectral slope in this region is approximately k −2 to k −5/3 at scales of 10-200 km (Chereskin et al., 2019). This is in contrast to the steeper spectral slope (k −3 ) in more energetic regions such as the Antarctic Circumpolar Current, which implies geostrophic dynamics (Rocha et al., 2016). Modeling studies in the California Current System have found the kinetic energy spectrum to be continuous from the mesoscale to submesoscale, with a slope of approximately k −2 (Capet et al., 2008a).
A range of dynamics could result in the observed spectral slope including internal gravity waves (k −2 ), surface quasigeostrophy (k −5/3 ), and fronts (k −2 ) (Boyd, 1992;Lapeyre & Klein, 2006). The observations available in this study cannot sufficiently distinguish these spectral slopes, nor can we identify whether the observed slope has transitions in the submesoscale regime from existing methods. We therefore rely on further analysis to infer the dynamics in this region.

Distributions of Vorticity, Divergence, and Strain Rate
One of the implications of a kinetic energy spectrum E(k) with a k −2 slope is that the velocity derivative spectrum V(k) is flat because the spectra are linked through the relationship V(k) = k 2 E(k). Velocity derivatives are essential to understanding and defining the submesoscale, which is defined as Ro∼(1) and where the flow becomes more fully three dimensional (velocity divergence δ ∼ U/L). The key velocity derivative quantities, divergence (δ = u x + v y ), vorticity (ζ = v x − u y ), and strain are related to each other through a system of coupled non-linear ordinary differential equations (cf. Barkan et al., 2019) and at the submesoscales these non-linear terms become more important. For example, the rate of change of vorticity in an adiabatic system is Only the first term on the right hand side, which does not involve a feedback, is present in a quasigeostrophic system. In addition, at the submesoscale the inertial term in the equations of motion (u ⋅∇u ∼ U 2 /L) is of the same order as the Coriolis term (uf ∼ Uf), facilitating cross-scale kinetic energy transfers (Johnson, 2020;Vela-Martín, 2022).
In the observations that allow for computing gradient across the track, which are only available for one track in the filament case (Figures 1b and 1c (Rudnick, 2001;Shcherbina et al., 2015). This skewness can arise from conservation of potential vorticity at fronts. Strain-driven frontogenesis at the sea surface, in the absence of dissipation, results in an infinitely sharp front in finite time with ageostrophic flow that has skewed distributions of divergence (negative) and vorticity (positive) (Barkan et al., 2019;Hoskins & Bretherton, 1972). Compared with anticyclonic fronts, cyclonic fronts are thought to progress more slowly to singularities during frontogenesis, which could result in longer lived cyclonic fronts (Shakespeare, 2016). In addition, the dynamical feedbacks are such that large negative relative vorticity is typically unstable to symmetric and centrifugal instabilities but positive relative vorticity stabilizes the flow to these instabilities, and therefore a skewed distribution develops (Buckingham et al., 2016;Rudnick, 2001). However, it is notable that strain-driven frontogenesis can suppress the growth of symmetric instability (Thomas, 2012), which could reduce the efficiency of this mechanism for generating skewness. Regardless, in boundary layers, negative potential vorticity can arise from frictional and atmospheric forcing, which can trigger symmetric instability.
We next examine the relationship between divergence, vorticity, and strain. The strain is composed of shear strain (σ s = v x + u y ) and normal strain (σ n = u x − v y ). In the filament observations studied here, vorticity is strongly correlated with shear strain (Figure 2d). Vorticity is not passively advected by shear strain and is not directly forced by shear strain (Equation 1), but in an integrated sense vorticity and shear strain become aligned. This correlation can persist due to the relative stability of cyclonic vorticity at straight fronts (Buckingham et al., 2021). This provides an explanation for the strain-vorticity relationship that has been observed in high-resolution simulations (Balwada et al., 2021). However, there is not a strong correlation between divergence and normal strain (σ n = u x − v y ). The lack of correlation is likely due to the short timescales of submesoscale dynamics; while non-zero vorticity can be maintained in an adiabatic system in the absence of divergence and vertical motion, a similar balance does not exist for divergence (Figure 2e).

Interactions Between Rotational and Divergent Flow
The approximately k −2 spectral slope in both the filament and eddy regions is informative but inconclusive about the dominant dynamics operating in these regions. The nearly uniform slope across the observed spatial scales leaves open questions about the scales at which a transition to submesoscale dynamics may occur.
The submesoscale feedback between vorticity and divergence (Equation 1) results in a correlation between the geostrophically balanced rotational (streamfunction) flow and the divergent (potential) component of the velocity. These components are not expected to be correlated if waves are not modified by geostrophic or non-wave ageostrophic dynamics (Bôas & Young, 2020). We diagnose when the streamfunction and potential become correlated using the cross spectrum (̂( ) ) between the along-track (u) and cross track (v) velocity components.
When the streamfunction and potential are uncorrelated, as is typically true at the mesoscale and larger, the cross spectrum of the u and v velocity components is a superposition of the spectra of the streamfunction and velocity potential. In this case, since spectra are real, the cross spectrum is real (Bühler et al., 2017). This is a key assumption of the "wave-vortex" decomposition introduced by Bühler et al. (2014). However, when the rotational and divergent flow components correlate, the cross spectrum between the along-track and across-track velocity is complex (̂( ) =̂( ) +̂( ) where ̂≠ 0 ). We are therefore able to diagnose the spatial scale where a shift to submesoscale dynamics occurs as the scale at which the cross spectrum becomes complex. We use the cross spectral phase (tan(ϕ uv (f))) to summarize this relationship ) . ( The cross spectral phase should only be interpreted in this manner if the coherence, which is the normalized cross spectrum between the along-track and across track velocity, is above the significance threshold (Text S2 in Supporting Information S1).
The squared coherence has contrasting dependence on spatial scale in the two regions studied here (Figure 3a). In the eddy region (blue lines in Figure 3), the squared coherence is large at the largest spatial scales sampled (∼100 km) and decreases toward smaller spatial scales, but remains significantly different from zero. In the filament region (green lines), the squared coherence is also large at the largest spatial scales sampled (∼10 km), decreases at scales larger than 6 km, and then increases again toward the smallest spatial scales sampled (∼1 km). Fronts and filaments are expected to be anisotropic at the scale of the feature, as is observed. In all observations considered here, the coherence is large enough to be statistically significant, allowing for analysis of the cross-spectral phase. The only exception is in the eddy region between 0.1 and 0.6 cpkm where the coherence falls below the significance threshold. The phase is noisy in this wavenumber range but is not significant and so should not be interpreted.
The cross-spectral phase summarizes the relationship between the real and imaginary parts of the cross spectrum. When the cross spectrum is purely real, the phase is 0° or 180°; when it is purely imaginary, the phase is ±90°. We find abrupt transitions at a scale slightly smaller than 10 km in the eddy region and 6 km in the filament region, where the imaginary part of the cross spectrum becomes larger than the real part ( Figure 3b). This 6 km spatial scale is the same scale where the coherence increases in the filament region (suggesting increased anisotropy), providing consistent evidence of a change to increasingly non-linear frontal dynamics at these scales. By contrast,  Figure S1 in Supporting Information S1) using 10 km windows. in a surface quasigeostrophic model (Abernathey et al., 2022), which neglects ageostrophic advection, the real part of the cross spectrum dominates at all spatial scales ( Figures S2 and S3 in Supporting Information S1). During the eddy observations, the mixed-layer depth was 40-55 m with some regions as shallow as 15 m. In contrast, for the filament observations, the mixed-layer depth was approximately 35 m or shallower and the stratification was approximately 3 × 10 −5 s −2 . Therefore, the mixed-layer deformation radius was 2-4 km for these locations, implying that the fastest growing baroclinic mode is around 8-24 km (Dong et al., 2020). Thus, the transition to non-linear ageostrophic dynamics observed here occurs in the approximate range of the scale of mixed-layer baroclinic instability.
There are a number of mechanisms that could be responsible for the interaction between rotational and divergent velocity. In the filament case, the interaction of the ageostrophic frontal divergence and larger scale geostrophic flow is likely the dominant mechanism. Here we find that the shift to a mostly imaginary cross spectrum is localized in the regions of largest velocity gradient ( Figure S4 in Supporting Information S1). The eddy case likely encompasses a larger range of dynamics, including near-inertial oscillations modified by the vorticity of the observed features, frontal dynamics, and submesoscale vortices.

Spectral Energy Transfers
The distribution of kinetic energy across spatial scales reflects dynamics that are local in wavenumber, but importantly also reflects energy transfers across scales. At the submesoscale, major open questions remain regarding the direction of the energy cascade, the mechanisms that lead to a forward energy, and the rate of the forward energy cascade (McWilliams, 2016;Müller et al., 2005). Forward energy flux precedes dissipation at small spatial scales by turbulent processes.
The energy transfer across scales can be quantified using coarse graining (Aluie et al., 2018;Eyink, 2005;Germano, 1992). The kinetic energy flux is defined here as where = − and ⋅ is a top hat filter. Positive (negative) values indicate a flux of energy toward smaller (larger) spatial scales. We use velocity observed on a 256 m grid and a top hat filter with a scale of 1 km to compute an instantaneous energy flux across the observed transect.
In the frontal regions in this flow, there is a strong forward energy flux localized in a 1 km region at the frontal outcrop, where there is also a peak in frontogenesis (Figures 4a and 4b). The energy flux to smaller spatial scales is driven by the first term in Equation 3 ( Figure S5 in Supporting Information S1). This term involves the shear strain multiplied by the scale-dependent covariance between the along-track and cross-track velocity. In fact, over the entire 60 km section, there is a strong correlation between the shear strain and the kinetic energy flux (Figure 4c).
The influence of the shear strain on the kinetic energy flux is modulated by the covariance between the u and v velocity components (or, equivalently, the anistropy of the flow), which becomes large below scales of 6 km in this filament region (Figure 3a). Barotropic shear instabilities extract kinetic energy from sheared mean flows when smaller scale features lean into the shear, resulting in a forward energy cascade.
The observed kinetic energy flux is patchy (Figure 4c), with the largest flux concentrated in small spatial scales even within the 60 km filament region observed here. On this whole transect, the kinetic energy flux varies over three orders of magnitude (Figure 4c). The typical kinetic energy flux across 1 km in the filament region is  ( 10 −6 m 2 s −3 ). This rate is about an order of magnitude larger than the kinetic energy flux obtained from mooring based observations using a filter scale of 5 days (Naveira Garabato et al., 2022), drifter based observations at 1 km scale that average over a large region of the northern Gulf of Mexico (Balwada et al., 2022), and a model at 500 m spatial resolution (Srinivasan et al., 2023). Given that we present direct observations of the kinetic energy flux terms, this suggests that the magnitude of instantaneous kinetic energy flux has been underestimated by previous modeling and observational work.

Discussion and Conclusion
The airborne observations presented here reveal large and patchy kinetic energy flux localized at submesoscale fronts and advance observational characterization of submesoscale dynamics. The synoptic sampling from submesoscale to mesoscale allows us to extend an observational kinetic energy spectrum to scales below 1 km.
Dense filaments such as the one observed here have an important role in the energetics of upwelling systems with submesoscale dynamics influencing the fate of upwelled waters.
We demonstrate that although there is not a clear change in the kinetic energy spectral slope, there is a transition in the dynamics to non-linear interactions that characterize submesoscales at scales of 6-10 km. In particular, this transition is characterized by the interaction between divergent and rotational velocity components. This transition would not occur with surface quasigeostrophic dynamics, and we attribute it instead to a dominance of ageostrophic dynamics in the observations.
The observed transition to non-linearity has important implications for observations of ocean velocity from remote sensing. For example, the SWOT mission aims to infer mesoscale to submesoscale velocities through observation of sea surface height. These velocities are computed through geostrophic balance, which only accounts for the rotational component of the flow (but does not include higher order rotational velocity components). Not only do we find that a significant amount of the kinetic energy is likely in the divergent component of the flow at scales below 10 km in this region-and potentially at larger scales in more energetic regions (Callies et al., 2015)-but also that the rotational and divergent flows interact such that filtering of the divergent processes (e.g., waves) to recover the rotational component of the flow may be more challenging.
These observations are also the first direct observations of snapshots of kinetic energy flux and frontogenesis in the ocean. This allows us to investigate the relationships between the kinetic energy flux and hydrographic features. We find that kinetic energy flux is patchy but can be large (10 −6 m 2 s −3 ) at submesoscale fronts. The patchiness of kinetic energy flux has important implications for resolving the dynamics that contribute to an energy cascade. Due to the difficulty resolving scales ranging from mesoscale straining to turbulent dissipation in models, these observationswhere that challenge is observationally addressed using a novel remote sensing platform-are particularly valuable. These aircraft measurements provide a preview of what might be possible from future satellite-based radar snapshots from platforms such as Harmony and SEASTAR (Gommenginger et al., 2019;López-Dekker et al., 2019). In these observations, kinetic energy is transferred both downscale and upscale from 1 km.
Recent modeling work has suggested that resolving frontogenesis is essential to accurate representation of submesoscale kinetic energy transfers (Naveira Garabato et al., 2022;Srinivasan et al., 2023). The observations analyzed here demonstrate a large forward energy transfer localized at fronts, although not exclusively during active large-scale frontogenesis. Recent work in the Gulf of Mexico, another region with an active submesoscale, has hinted that a forward cascade of kinetic energy occurs at scales of 500 m-5 km (Balwada et al., 2022) in observations (with smaller scales during the summer) and at scales of 5 km in models (Srinivasan et al., 2023).
The modeling study of Sullivan and McWilliams (2018), which simulated a dense filament, also found an important role for the horizontal Reynolds stress term (u′v′v x ) during the frontal arrest phase of a dense filament, which is consistent with our observation that the shear strain term dominated kinetic energy flux. This relationship may arise from certain aspects of the feature studied here and may not generalize to all fronts. For example, Srinivasan et al. (2023) analyzed kinetic energy fluxes in 500 m and 2 km resolution ocean models, which resolve dynamics at larger scales than those that are the focus of our study. They find an equipartition between strain-driven and convergence-driven forward energy cascade at submesoscale scales (Srinivasan et al., 2023). While we observe that the forward energy transfer is strain-driven in our observations, it is important to note that we have only one snapshot of a filament that appears to be partially restratifying, so this does not invalidate the role of convergence in forward energy flux.
These results suggest an out-sized role for fronts and filaments as hotspots of surface kinetic energy flux. Barotropic energy transfer is enabled by interactions between the rotational and divergent components of the flow field at submesoscale fronts. Fronts are spatially inhomogeneously distributed in the ocean and vary seasonally (Drushka et al., 2019;Mauzole et al., 2020), but the distributions of fronts are distinct from the distributions of mesoscale kinetic energy (Busecke & Abernathey, 2019). Surface kinetic energy dissipation may similarly vary substantially in space and time, but understanding how it varies relies on increased mechanistic understanding of the energetics of submesoscale features. Disentangling these would require more observations to establish the effect of particular submesoscale features on the regional statistics.

Data Availability Statement
All presented data are available at UCSD Library Digital Collection, https://doi.org/10.6075/J0F76CRK.