# The Vertical Structure of Open-Ocean Submesoscale Variability during a Full Seasonal Cycle

Andrew F. Thompson, Jörn Callies, Xiaolong Yu, Alberto Naveira Garabato, Patrice Klein

► To cite this version:

Andrew F. Thompson, Jörn Callies, Xiaolong Yu, Alberto Naveira Garabato, Patrice Klein. The Vertical Structure of Open-Ocean Submesoscale Variability during a Full Seasonal Cycle. JOURNAL OF PHYSICAL OCEANOGRAPHY, 2020, 50, pp.145-160. 10.1175/JPO-D-19-0030.1. insu-03683241

HAL Id: insu-03683241

https://insu.hal.science/insu-03683241

Submitted on 31 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Vertical Structure of Open-Ocean Submesoscale Variability during a Full Seasonal Cycle

ZACHARY K. ERICKSON,ª ANDREW F. THOMPSON, AND JÖRN CALLIES

Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California

XIAOLONG YU

Ifremer, Université de Brest, CNRS, IRD, Laboratoire d'Océanographie Physique et Spatiale, IUEM, Brest, France

ALBERTO NAVEIRA GARABATO

Ocean and Earth Science, University of Southampton, Southampton, United Kingdom

PATRICE KLEIN

Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California

(Manuscript received 5 February 2019, in final form 23 October 2019)

## ABSTRACT

Submesoscale dynamics are typically intensified at boundaries and assumed to weaken below the mixed layer in the open ocean. Here, we assess both the seasonality and the vertical distribution of submesoscale motions in an open-ocean region of the northeast Atlantic. Second-order structure functions, or variance in properties separated by distance, are calculated from submesoscale-resolving ocean glider and mooring observations, as well as a 1/48° numerical ocean model. This dataset combines a temporal coverage that extends through a full seasonal cycle, a horizontal resolution that captures spatial scales as small as 1 km, and vertical sampling that provides near-continuous coverage over the upper 1000 m. While kinetic and potential energies undergo a seasonal cycle, being largest during the winter, structure function slopes, influenced by dynamical characteristics, do not exhibit a strong seasonality. Furthermore, structure function slopes show weak vertical variations; there is not a strong change in properties across the base of the mixed layer. Additionally, we compare the observations to output from a high-resolution numerical model. The model does not represent variability associated with superinertial motions and does not capture an observed reduction in submesoscale kinetic energy that occurs throughout the water column in spring. Overall, these results suggest that the transfer of mixed layer submesoscale variability down to depths below the traditionally defined mixed layer is important throughout the weakly stratified subpolar mode waters.

## 1. Introduction

Most of the energy in the ocean resides at scales of hundreds of kilometers (the mesoscale), where the ocean is primarily in geostrophic balance (Ferrari and Wunsch 2009). However, much of the vertical transport of oceanic tracers such as heat, carbon, and nutrients is thought to be accomplished at scales from hundreds of meters to kilometers (the submesoscale), where

rotation and advection components of the momentum budget can be of equal importance (Lévy et al. 2012). Submesoscale motions are largely catalyzed by the transfer of energy from the mesoscale through mixed layer instabilities (McWilliams 2016; Callies et al. 2016). Spatial and temporal patterns of submesoscale phenomena, and their impact on large-scale distributions of passive tracers, remain an active area of research. However, the vertical structure of submesoscale motions has remained particularly elusive because of the challenges in maintaining persistent in situ observations at the requisite spatial scales.

A key finding of recent observational studies is that submesoscale dynamics typically exhibit seasonality driven by annual variations in mixed layer depth (MLD) (Callies et al. 2015; Thompson et al. 2016; Buckingham

<footnote>ª Current affiliation: Ocean Ecology Laboratory, NASA Goddard Space Flight Center, Greenbelt, Maryland.</footnote>

<footnote>Corresponding author. Zachary K Erickson, zachary.k.erickson@nasa.gov</footnote>
et al. 2016; Erickson and Thompson 2018). Deep mixed layers contain high potential energy (PE) that can be released through instabilities that tend to restratify the mixed layer (Haine and Marshall 1998; Fox-Kemper et al. 2008). Modeling studies also show increased submesoscale activity during winter, largely diagnosed in terms of changes in energetic characteristics (Capet et al. 2008; Mensa et al. 2013; Su et al. 2018). Brannigan et al. (2015) find an enhancement in surface kinetic energy (KE) in the northeast Atlantic at submesoscales as model resolution is increased due to sharper fronts, stronger mixed layer baroclinic instabilities, and more frequent instances of symmetric instability. Sasaki et al. (2014) also find an increase in submesoscale activity in winter in the North Pacific, characterized by a flattening of the KE spectral slope from $k^{-3}$ during summer to $k^{-2}$ in winter, where k is an isotropic wavenumber. Following winter, energy is transferred to larger scales, resulting in a temporal shift of about 100 days between the maximum KE at scales of 200–300 km compared with scales of 10–100 km.

The depth to which this seasonality is evident is less clear. A limiting case is a mixed layer bounded by a strong pycnocline: submesoscale motions are isolated within the mixed layer and do not penetrate into the interior (Boccaletti et al. 2007; Fox-Kemper et al. 2008). Discussion of submesoscale dynamics in the literature has largely been confined to mixed layer properties and assumed to be negligible deeper in the water column (Klein et al. 2008), although recent studies have challenged this paradigm (Siegelman et al. 2019; Yu et al. 2019a). In wintertime conditions, however, many regions of the ocean have a weak pycnocline at the base of the mixed layer; in these locations submesoscale instabilities and submesoscale advection have been proposed to extend beyond the depth of a traditionally defined mixed layer (Erickson and Thompson 2018). Additionally, features within the ocean interior, such as subthermocline eddies (McWilliams 1985), can also induce significant small-scale features at depth (Hua et al. 2013). Balwada et al. (2018) found that submesoscale fluxes across the base of the mixed layer increased with finer horizontal model resolution even though the vertical stratification at the base of the mixed layer also increased, pointing to the increase of vertical velocities at the submesoscale.

In this study, we address the vertical structure of submesoscale dynamics using the framework of spatial tracer distributions that can be related to turbulent stirring properties. Quasigeostrophic (QG) theory predicts that the spectral slope of KE should scale as $k^{-3}$ far from boundaries and potential vorticity discontinuities (Charney 1971). For turbulent flow, a KE spectrum that

conforms to a power-law $k^{-n}$ can be used to predict that a passive tracer distribution will have a spectral slope of $k^{-(n-5)/2}$, for $1 < k < 3$ (Vallis 2006). This relationship, then, implies a passive tracer spectral slope of $k^{-1}$ in the interior ocean from QG theory. Near boundaries, sharp gradients in passive tracers can emerge, leading to steeper slopes of $k^{-2}$ (Klein et al. 1998). Therefore, passive tracer slopes are predicted to have a vertical structure that flattens from $k^{-2}$ near the surface to $k^{-1}$ at depth.

While these theoretical predictions are obeyed in QG models (Smith and Ferrari 2009), they are seldom borne out in field studies. Specifically, a flattening of the slope of tracer variance spectra is rarely seen. Cole and Rudnick (2012) found a $k^{-2}$ structure of spice variability with minimal seasonality throughout the upper 1000 m of the water column in the subtropical North Pacific Ocean. In the northern Pacific, Schönau and Rudnick (2015) found a $k^{-2}$ structure of spice variability irrespective of depth. Salinity gradient spectra along isopycnals in the California Current System were found to obey $k^0$, or $k^{-2}$ for salinity variance, irrespective of season (Itoh and Rudnick 2017). Kunze et al. (2015) also found a $k^0$ passive tracer gradient slope down to 100 m, and suggested non-QG stirring and internal wave/horizontal strain as possible mechanisms. Klymak et al. (2015) found an agreement with QG theory between the flow field and passive tracers near the surface, but a reddening (steepling) of passive tracer spectra with depth in the Gulf of Alaska (near Station Papa), inconsistent with a surface-intensified frontal structure— however, this reddening with depth was not seen in more open-ocean areas of the northern subtropical Pacific Ocean. Long probability-distribution-function tails of spice indicated sharp spice contrasts in both the Gulf of Alaska and the North Pacific subtropical open ocean down to 350-m depth. As shown below, relatively weak changes in tracer spectral slope in the vertical exist also in the North Atlantic, despite a significant seasonal cycle in upper-ocean KE and PE.

The mechanisms that cause passive tracer spectral slopes to steepen to $k^{-2}$ in more quiescent open ocean regions remain unclear (Callies and Ferrari 2013), although they are indicative of submesoscale activity not predicted by standard theories (Kunze et al. 2015). It also remains undetermined as to whether differences in KE/PE and tracer spectra across different studies can be explained by regional flow characteristics, temporal variability, or even measurement technique. Therefore, observational studies spanning at least a full year are important to understand these phenomena.

We use glider and mooring observations from one such observational campaign, the Ocean Surface Mixing, Ocean
Submesoscale Interaction Study (OSMOSIS), to consider the seasonality of variance in horizontal velocity (KE), buoyancy (PE), and spice (a passive tracer), and compare our results to data from a high-resolution numerical model. The studies mentioned above have largely relied on repeat hydrographic sections to calculate spectral properties of the flow field. Here, we take advantage of the multiple observational platforms that comprised the OSMOSIS study to characterize turbulent properties of the OSMOSIS region using second order structure functions, or variance in properties binned by separation distance As many ocean datasets, especially from quasi-Lagrangian instruments, are unstructured, this technique is a particularly useful method for extending regional analyses of submesoscale turbulence (e.g., Balwada et al. 2016).

The manuscript is organized as follows. The observations and modeling output are introduced in section 2, where we also give an overview of the region. Section 3 describes the structure function technique. The results of our structure function analysis are presented in section 4, and in section 5 we discuss differences between the model and observations, the effects of internal waves, theoretical models of ocean mixing, seasonality, and implications for tracer fluxes out of the surface mixed layer.

## 2. Data

### a. Glider observations

Five Seagliders (gliders) were deployed in a 20 km × 20 km region of the northeast Atlantic Ocean over a full year as part of the OSMOSIS project (Figs. 1a,b) (Thompson et al. 2016; Damerell et al. 2016; Buckingham et al. 2016). Staggered glider deployments ensured that the region was always sampled by at least two gliders, although sensor complications on one glider during November–December 2012 rendered some of the data unusable (Fig. 1c). Glider data processing, including thermal lag and salinity corrections, is described by Damerell et al. (2016). Glider conductivity–temperature–depth (CTD) measurements were made at approximately 1-m depth intervals, with a precision of 0.0003 S m⁻¹, 0.001°C, and 0.001 kg m⁻¹ for conductivity, temperature, and salinity (derived from conductivity and temperature), respectively. CTDs were calibrated with ship measurements made during deployment and recovery of each glider. A subsequent filter that removed any profile with an average salinity of less than 35.1 psu or temperature less than 9°C was also found necessary to remove erroneous measurements.

The gliders were piloted in bowtie patterns, with approximately five dives per leg, within the OSMOSIS region (Fig. 1b). Each V-shaped dive lasted approximately

5 h and the horizontal spacing between dives was generally 2–4 km; each leg of the bowtie pattern lasted approximately one day. The glider location is transmitted before and after each dive, and the horizontal glider position during the dive is linearly interpolated with respect to time between these two points. Location error produced by this interpolation is estimated using the model developed in Frajka-Williams et al. (2011) (see also the UEA Seaglider Toolbox; http://www.byqueste.com/toolbox.html) to be on average less than 1 km at the bottom of each dive. We have also completed the analysis in this paper using the glider locations given by the Frajka-Williams et al. (2011) model and found no significant differences, lending confidence to our results (not shown). Occasionally the gliders were advected out of the area shown in Fig. 1b; these data were not used in our analysis.

### b. Mooring observations

In addition to gliders, nine moorings were arrayed in two concentric quadrilaterals with side lengths of 2–3 and ~13 km around a central mooring (Fig. 1b). The moorings were instrumented with CTDs and acoustic current meters (ACMs) at 20–200-m intervals within the upper 600 m [Fig. 1d; see Buckingham et al. (2016) or Yu et al. (2019a) for more details]. ACMs recorded velocity data at 10-min intervals and CTDs at 5-min intervals; for this study CTD data were subsampled to ACM temporal resolution. Nine moorings resolve 36 different separation distances, which range from 1.2 to 18.8 km (Fig. 1e, open circles).

The moorings were subject to currents in the area, and pressure sensors on each CTD and ACM recorded deviations in depth of up to 150 m. These vertical deviations introduce error into the horizontal distance between moorings; however, in a separate analysis, Buckingham et al. (2016) found that the buoyancy applied to the mooring cables restricts their lateral movement and creates an effective pivot point near 600-m depth; stochastic modeling predicted horizontal displacements rarely exceeding 500 m. As we expect these motions to be part of a larger-scale field acting upon all of the moorings, the change in separation distance between each combination of moorings is likely to be considerably smaller than 500 m.

### c. High-resolution model

For comparison with the in situ observations, we analyze a region of the llc4320 model, a high-resolution 1/48° global MITgcm simulation. The model is initialized from Estimating the Circulation and Climate of the Ocean, Phase II (ECCO2) output (Menemenlis et al. 2008), after which the resolution is increased sequentially to 1/12°, 1/24°, and finally 1/48° (Wang et al. 2018;
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Bathymetry from ETOPO in the northeast Atlantic. The OSMOSIS region is shown in the white box at 48.7°N, 16.2°W. The larger region from the model is shown as the dotted white box. (b) Highlight on the OSMOSIS region, showing a histogram of glider surface locations (colors) and the positions of the nine moorings (x marks). (c) Periods of time in which the gliders were active. (d) Depth placement of the ACMs (x marks) and CTDs (dash marks) for each mooring. Dotted horizontal lines denote the depths over which mooring structure functions were calculated. (e) Histogram of structure function pairings from glider measurements at 200-, 400-, and 800-m depth. Bins are equally spaced logarithmically. Circle markings at the top axis show the separations between moorings.</figcaption></figure>

Torres et al. 2018). The name represents the domain configuration (latitude-longitude-polar cap) and the number of grid cells in the polar cap (4320 × 4320). The llc4320 is forced by the ECMWF Operational Model Analysis atmospheric parameters and 16 tidal components. The ECMWF forcing is interpolated from its native 0.14° and 6-hourly resolution to the higher llc4320 spatial and temporal grid using linear interpolation in time and bilinear interpolation in space (bicubic for winds; D. Menemenlis 2019, personal communication). Here we use one year (10 September 2011–9 September 2012) of model output from an approximately 120 km × 120 km box centered on the OSMOSIS location (Fig. 1a, dotted outline) and extending

from the surface to 1-km depth. The model has a horizontal resolution of approximately 1.5 km and 52 vertical levels ranging in thickness from 1 m at the surface to almost 50 m at 1-km depth. The model time step is 25 s, and data are saved as snapshots every hour. The effective spatial resolution can be estimated as 4 times the grid spacing, or approximately 6 km.

The llc4320 output has previously been compared with Argo data in the Kuroshio Extension, and showed reasonable vertical density stratification and seasonal variability (Rocha et al. 2016). Globally, Qiu et al. (2018) found consistent surface eddy KE distributions between those inferred from AVISO and llc4320 sea surface height after the latter was coarse-grained to AVISO resolution.
Comparisons with submesoscale-permitting observations are limited, but Viglione et al. (2018) found instances of surface instabilities at submesoscales that were temporally and spatially consistent with glider observations in Drake Passage.

Although previous work has validated a number of aspects of this model, care must be taken in making a comparison to in situ data. The model does not assimilate data, and in particular, does not reproduce discrete events such as the occurrence of an eddy within a domain at a specific time. While the external forcing is from reanalysis data, the llc4320 output (September 2011–September 2012) does not match the timeframe of the in situ data (September 2012–September 2013), and without observational constraints the model would decorrelate from reality even if forced by the same time period. However, statistics of ocean properties calculated over suitably large time intervals still permit a useful comparison between in situ observations and model results.

### d. Site characterization

The OSMOSIS site was chosen because of the lack of major bathymetric features or mean geostrophic currents; the eddy PE and KE are moderate (Roullet et al. 2014; Rieck et al. 2015), and this study therefore provides a complement to recent experiments concerning submesoscale dynamics in more active boundary current regions (Rocha et al. 2016; Thomas et al. 2016). The site experiences a strong seasonal cycle, which is primarily seen in annual variation of the MLD (Thompson et al. 2016; Damerell et al. 2016; Erickson and Thompson 2018), calculated from gliders and model output as the depth at which the potential density reaches $0.03 \text{ kg m}^{-3}$ above the potential density at 10 m (Fig. 2a) (de Boyer Montégut et al. 2004). This density threshold was previously found to agree with expected MLDs from optical glider measurements of chlorophyll fluorescence and backscatter in the OSMOSIS dataset (Erickson and Thompson 2018). During autumn (October–December) the MLD steadily increases, with highly variable wintertime MLDs (January–March) reaching 400 m. These seasonally deep mixed layers lead to the production of subpolar mode water in this region (McCartney and Talley 1982). The model accurately captures the autumnal deepening, wintertime variability (note that only domain-averaged values are shown), and shallow summertime values, but with a deeper mean depth during winter.

Both potential density and vertical stratification profiles show seasonality near the surface (Figs. 2b,c), with a sharp pycnocline in the summer (green) absent during winter (blue). The main pycnocline, characterized by

$N^2 \sim 10^{-5} \text{ s}^{-2}$, is located at about 800 m. We note that the model is lighter than the observations throughout the year, but this does not influence the structure function results below. The model also does not fully capture the strength of the summertime pycnocline (inset to Fig. 2b).

We treat spice (II), the component of temperature and salinity that does not contribute to density (Veronis 1972; Munk 1981), as a passive tracer. This approximation is valid in the absence of nonlinearities in the equation of state of seawater, which are small in this region (Damerell et al. 2016). Spice is calculated using the algorithm from McDougall and Krzysik (2015), and in the given region of temperature–salinity space has an estimated precision of $0.001 \text{ kg m}^{-3}$. Time series of spice during the winter-spring transition from glider measurements show remarkable small-scale features extending to the deepest glider measurements at 1000 m (Fig. 3a). Here we eliminate the heaving effects of internal waves by considering spice along potential density rather than pressure surfaces. These features are consistent over many dives (roughly five per day) and over a range of potential density surfaces. The model also shows high spice variability at depth in coherent subductive events, such as one seen in day 115–118 (Fig. 3b; note that the alignment with an observed subductive event is coincidental). Submesoscale spice variance is pervasive throughout the year during the OSMOSIS time series, even though only a short period is displayed in Fig. 3.

Snapshots of spice at potential densities of 27.03, 27.07, and $27.31 \text{ kg m}^{-3}$, corresponding to average depths of 200, 400, and 800 m, respectively, showcase the processes captured by the model (Fig. 4). A low-spice, mesoscale eddy in the northwest corner of the domain stirs water masses into filaments. A larger, anticyclonic eddy to the southeast of the domain has a high spice anomaly, and T-S characteristics suggest it is sourced from Mediterranean water outflow. A filament stretches from this eddy into the center of the region studied here (dotted black boxes). The solid black boxes represent the size of the observational OSMOSIS domain (Fig. 1b). The high-spice anomaly associated with this filament is strongest at depth. However, although the filament width is narrow compared with the size of the eddy, it is large compared with the OSMOSIS domain and the separations resolved by glider and mooring measurements. The exceptionally sharp features in the glider data are therefore not captured in the model.

Tidal influences, and especially the $M_2$ tide, are pronounced in this region, as seen in the sharp peak at the $M_2$ frequency in KE, PE, and spice power spectral density (Fig. 5). The subinertial component of these
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) MLD from glider (gray line; black line is filtered through a Gaussian window with standard deviation of 1 day) and model (black dotted line, as an average over the model region marked in Fig. 1a). For the gliders, the date is in reference to 1 Jan 2013; the model reference is to 1 Jan 2012. Summer and winter times are indicated by green and blue shading, respectively. Average (b) potential density and (c) vertical stratification $N^2 = b_z$ for summer (green) and winter (blue) from glider measurements (solid line), moorings (dashed line), and the model (dotted line). Shading indicates the 50% confidence interval for glider measurements. Inset in (b) highlights the upper 100 m of the water column. Bars in (b) represent the MLD of 90% of the measurements for winter (blue) and summer (green) from gliders (SG) and the model (LLC).</figcaption></figure>

variables agrees well between the moorings and the model; however, the model is missing considerable energy in the superinertial range of the spectra. In the mooring data, the superinertial range closely follows the Garrett and Munk (1975) (GM) spectrum for internal waves.

### 3. Methods

Wavenumber spectra are traditionally used to assess tracer variance as a function of scale. However, not all datasets are amenable to spectral decomposition. Structure functions (SFs), defined below, are a useful technique when observations, such as from surface drifters (Balwada et al. 2016) or Argo floats (McCaffrey et al. 2015), do not follow defined transects.

The nth-order SF of a scalar tracer $\theta$ is

$$D_{\theta}^{n}(\mathbf{s}) = \overline{[\theta(\mathbf{x}) - \theta(\mathbf{x} + \mathbf{s})]^{n}}, \quad (1)$$

where $\mathbf{x}$ is a position, $\mathbf{s}$ a separation distance, and boldface variables are vectors. In general $\mathbf{x}$ and $\mathbf{s}$ can be multidimensional, but for this study $\mathbf{x}$ represents a latitude/longitude position, $s = |\mathbf{s}|$ denotes a horizontal distance and an implied temporal constraint on time differences between measurements (see below), and the overbar is an average over all $\mathbf{x}$ in a given time window.

SFs provide information on how variance (or skewness, kurtosis, etc. for $n > 2$) changes as a function of separation distance, without requiring that all regions be uniformly sampled in a grid, as may be the case for power spectra decompositions. The second-order SF is related to the variance spectrum $E_{\theta}(k)$ as (Webb 1964; Babiano et al. 1985)

$$D_{\theta}^{2}(\mathbf{s}) = 2 \int_{0}^{\infty} E_{\theta}(k) [1 - \cos(ks)] dk. \quad (2)$$

Assuming $E_{\theta}(k)$ is represented by $ak^{-\lambda}$, for $1 < \lambda < 3$ the associated shape of $D_{\theta}^{2}(s)$ is $a's^{\gamma}$, where (McCaffrey et al. 2015)
<figure><img src="image_4.png" /><figcaption>FIG. 3. Spice along isopycnals for 15 days during April-May from (a) a glider and a (b) single horizontal position in the model. Average depths for the isopycnals are indicated on the right of each panel.</figcaption></figure>

$$ \gamma = \lambda - 1. \qquad (3) $$

In analogy to spectra of KE and PE, we define

$$ D_{KE} = \frac{1}{2}(D_{u}^{2} + D_{v}^{2}), \qquad (4a) $$

and

$$ D_{PE} = \frac{D_{b}^{2}}{2b_{z}}, \qquad (4b) $$

where b is buoyancy, u, v are the velocities in the x, y direction, and the z subscript denotes partial differentiation

in the vertical. Although horizontal velocities from gliders can be estimated using depth-average current calculations and assumptions of thermal wind shear, there is no reliable way to estimate the vertical structure of along-track velocities. We therefore do not attempt a calculation of $D_{KE}$ using glider data. The terms $D_{KE}$ and $D_{PE}$ are calculated on constant-depth surfaces, whereas because of the influence of isopycnal heaving on the spice field, $D_{\Pi}$ is appropriately calculated on isopycnal surfaces. The term $D_{KE}$ was also calculated on isopycnal surfaces, with only negligible differences. By definition, $D_{b} \equiv 0$ when calculated along an isopycnal surface.

<figure><img src="image_5.png" /><figcaption>FIG. 4. Model snapshots of spice on day 117 (28 Apr 2012) for isopycnals 27.03, 27.07, and 27.31 kg m⁻³, corresponding to average depths of 200, 400, and 800 m. Spice is shown as an anomaly from the average value in each panel domain. Black boxes give the size of OSMOSIS region (solid white box in Fig. 1a), and the dashed black box is the domain over which the model SFs are calculated (the dashed white box in Fig. 1a). Figure 3b is taken from the center of these boxes.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Frequency spectra of (a) KE, (b) PE, and (c) spice power at 525 m from the model (black) and the moorings (gray). Black dotted lines give the GM spectra, using the formula from Garrett and Munk (1975). Dashed vertical gray lines give the local planetary vorticity f and M2 tidal frequencies, and representative slopes of -2 and -3 are shown.</figcaption></figure>

Applying the SF framework to varied datasets in a consistent way is challenging due to differences in temporal and spatial sampling patterns between gliders, moorings, and the numerical model. In particular, a decision must be made concerning how to average measurement pairings in terms of horizontal separation distance, temporal separation between measurements, and vertical coordinate resolution. For this manuscript, we average measurement pairings using logarithmically spaced bins for s, and we make the following choices for each dataset:

Gliders—Consider temporal separations up to 3 h as contemporaneous; calculate properties vertically every 25 m (b) or 0.01 kg m⁻³ (Π); do not calculate D²u or D²v.

Moorings—Allow only near-simultaneous (<5 min) measurements when the pressure sensor is within 10 m of one of six target depths (see Fig. 1d) for u, v, and b; do not calculate D²II.

Model—Calculate simultaneous differences between 750 randomly generated points (over 500 000 pairings) within an 80 × 80 pixel (40 million possible pairings) domain; calculate u, v, b, σ (density), and Π at each modeled depth, and vertically interpolate Π to every 0.025 kg m⁻³.

We also calculate superinertial SFs from the mooring and model datasets. This is achieved by first filtering out signals less than a cutoff frequency, which is set to the local inertial period of 16 h. We tested the resilience of this filtering method by applying a cutoff of twice the inertial period to account for the width of the inertial peak; the qualitative results were not significantly affected.

We note that in our methodology, the pressure criterion for the mooring may cause mooring-based SFs to undersample knockdown periods, which strongly correlate

with high kinetic energy. However, time periods of mooring knockdown are correlated among the different moorings, suggesting that knockdown is a function of larger-scale ocean currents that are unlikely to affect the small-scale structure that is the focus of this study.

The SFs can generally be well summarized by the power law

$$D_{II}^2 = \alpha(s - s_0)^\gamma, \quad (5)$$

where α is the variance at a separation of s₀ = 20 km. We fit the SFs to Eq. (5) using a Levenberg-Marquardt minimization algorithm, with weights for the in situ calculations as inverse standard deviations, estimated as the spread of the 90% confidence interval divided by 3.29 (i.e., assuming Gaussian distributions). For the gliders, we smooth the final α and γ results vertically with a Gaussian window with a standard deviation of 75 m. For fits to modeled data we use only SF calculations for s between 4 and 20 km, as this s domain has relatively uniform slope.

## 4. Results

We construct SFs of KE [DKE; Eq. (4a)], PE [DPE; Eq. (4b)], and spice [(1/2)D²II; Eq. (1)] from the gliders, moorings, and model during winter and summer. Representative calculations for winter at 525-m depth are shown in Fig. 6 (blue lines). Slopes for DKE and DPE for in situ data (solid and dashed) are between 1/2 and 1, corresponding to a spectral slope of k⁻¹.⁵ to k⁻². The passive tracer SF, D²II, is somewhat shallower, at a slope of close to 1/2, or k⁻¹.⁵. SF slopes of KE and PE show good agreement between the model (dotted) and in situ measurements at scales of 4–20 km (Figs. 6a,b). At scales smaller than about 4 km, the in situ measurements from
<figure><img src="image_7.png" /><figcaption>FIG. 6. Structure functions (SFs) for (a) KE, (b) PE, and (c) spice (II) for winter (see Fig. 2a) at 525 m from gliders (solid), moorings (dashed), and the model (dotted). Blue lines correspond to the standard SF calculation; gray lines are only using superinertial frequencies as described in the text. The 90% confidence interval from a bootstrap analysis is given in light shading for the glider results. Confidence intervals for the mooring and model are not shown, but are small compared to those for the gliders. Representative slopes of $s^1$ and $s^{1/2}$ are shown in each panel.</figcaption></figure>

gliders and moorings are both considerably flatter than the model. Spice SFs are significantly larger in magnitude in the observations than the model results, indicating larger spice variance especially at small scales (as seen also in Figs. 3 and 5). At larger scales, the model slope decreases, representing a saturation of variance at scales approaching 100 km.

The mooring superinertial SFs (gray dashed lines) are spectrally flat, and can even be dominant over the sub-inertial component for $D_{KE}$ at scales less than 5 km (subinertial results not shown, but are equivalent to the subinertial SFs subtracted from the full SF). The model results, in contrast, show little variance associated with superinertial motions (gray dotted lines), meaning the subinertial SFs are similar to the full SF at all scales. This is an important result from this in situ dataset, as theoretical and numerical models of stirring in the ocean do not typically account for superinertial motions (e.g., Smith and Ferrari 2009).

SFs can be calculated at any depth to provide a full vertical structure of variance for a given property. An example for wintertime spice variance from gliders is shown in Fig. 7a. This calculation reveals larger spice variance with increasing depth and increasing separation.

It is convenient to approximate the resulting SFs as a power law [Eq. (5); example of the fit is shown in Fig. 7b]. The best-fit slopes $\gamma$ for winter and summer are shown in Fig. 8. Shading and error bars give the standard deviation on the fit; however, this is an incomplete description of the full uncertainty. In some instances slopes (and magnitudes, discussed below) can vary depending on the precise time range chosen for winter and summer, and since we only have one year of data, deconvolving

seasonal effects with either interannual oscillations or chance occurrences in one year, such as an eddy drifting into the region, is not possible.

The slopes for $D_{KE}$ using mooring data reveal little seasonal variance, and are approximately constant with depth, at about 0.7 (Fig. 8a). The model, in contrast, shows a striking seasonality, with slopes increasing dramatically, indicating less variance at small scales, during the summer. For $D_{PE}$, the slope varies with depth in the upper 200 m in the summer, being small (greater variance at small scales) near the surface and increasing throughout the summertime pycnocline (Fig. 8b). The slopes in spice variance are more difficult to interpret, and seasonal differences below 600 m probably rely more on the presence of eddies and small-scale filaments moving through the region (Fig. 8c). However, the model shows an increased slope with respect to the observational results, indicating less variance at small scales. During the winter, evidence of a transition across the base of the mixed layer (~300 m) in either $D_{KE}$, $D_{PE}$, or $D_{II}^2$ slopes is muted or nonexistent.

Observational $D_{KE}$ magnitudes $\alpha$ show a seasonal cycle in the mooring measurements, with higher variance in the winter (Fig. 9a). Modeled results show a slight increase at depths greater than 350 m in summer; however, if values from the beginning of the model run (the previous summer) are included, this relationship reverses itself. This seasonal difference is likely due to the model not being sufficiently “spun up” to its current horizontal resolution.

The seasonality and vertical dependence in the magnitudes for $D_{PE}$ (Fig. 9b) largely follow changes in the vertical stratification $b_z$. However, both the observational and model results show seasonality in the $D_{PE}$
<figure><img src="image_8.png" /><figcaption>FIG. 7. (a) Structure function of spice power at all isopycnals from the glider using wintertime data, expressed in terms of equivalent depth (average depth of isopycnal during winter). (b) Example of calculating best-fit slopes for a representative SF taken from 27.17 kg m⁻³ (525-m equivalent depth) in (a). The black line gives the best fit to the data, and the gray horizontal bar shows the range of values of the magnitude α, defined as the variance at 20-km separation.</figcaption></figure>

magnitude below the base of the diagnosed wintertime mixed layer, and indeed down to the permanent thermocline below 800 m, despite a lack of significant seasonality in $b_z$ below about 400 m (Fig. 2c). The magnitudes of spice variance show little seasonality for the observational data, but significantly decrease during winter in the model output. Together with the increase in summertime SF slope, this indicates a significant loss of spice variance in the simulation of the OSMOSIS region, which is inconsistent with both the glider and mooring data.

Superinertial slopes and magnitudes for $D_{KE}$ and $D_{PE}$ calculated from the moorings (Figs. 8 and 9, light dots) are relatively uniform with respect to depth and season when compared with the full results, although there is an indication that the superinertial magnitudes increase during winter within the wintertime mixed layer. This suggests a constant superinertial field reminiscent of the universal GM spectrum for internal waves, although we emphasize again that this technique is not able to

<figure><img src="image_9.png" /><figcaption>FIG. 8. Best-fit slopes of (a) KE, (b) PE, and (c) spice structure functions from gliders (solid), moorings (dots), and model (dotted line) during winter (blue) and summer (green). Light, smaller dots indicate superinertial structure functions from the moorings. Shading and error bars show the standard deviation of the fits.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 9. Best-fit magnitudes of (a) KE, (b) PE, and (c) spice structure functions from gliders (solid), moorings (dots), and model (dotted line) during winter (blue) and summer (green). Light, smaller dots indicate superinertial structure functions from the moorings. Shading and error bars show the standard deviation of the fits.</figcaption></figure>

decompose internal waves from balanced submesoscale motions.

Calculating best-fit magnitudes for the SFs at various time periods throughout the year reveals differences in seasonality between observational and model results (Fig. 10; differences in slopes throughout the year are not significant and are not shown here). KE in the model peaks roughly a month later than the observational data; the model peaks in the spring while the mooring time series has KE maximum in winter. For PE, each of the three datasets—gliders, moorings, and model results—peak in early winter at 50 m and later in winter for deeper depths. However, the model also has a twofold increase in the PE magnitude at 20 km in late spring which is not present in the observations. The cause of this increase is unclear, but is likely related to the dynamics of the spring restratification, which involves small-scale instabilities such as baroclinic mixed layer and symmetric instabilities (Erickson and Thompson 2018)

that will not be fully represented in the model. However, this difference between the model and in situ data may also be due to differences in surface forcing from different years.

## 5. Discussion

### a. Comparison between observations and high-resolution model

Because of their applicability across a wide range of datasets, it is important to assess the fidelity of modeled SFs with respect to observations, even though spectra can be directly computed from the model output. For all of the SFs analyzed in this study (KE, PE, and spice variance), the model output had steeper slopes as compared to SFs constructed from the observations. The SF magnitudes at 20 km, however, are more comparable, indicating that this increase in slope for the modeled output represents too little variance at small scales.

<figure><img src="image_11.png" /><figcaption>FIG. 10. Best-fit magnitudes of (a) KE and (b) PE applied to moving 30-day windows at 50 and 350 m (colors) for gliders (solid), moorings (dashed), and model results (dotted). Data have been smoothed in time.</figcaption></figure>
This is an expected characteristic of models, for which variance on scales smaller than about 4–5 model grid spacings (about 6 km here) is expected to be artificially low, and can also be seen in spectral decompositions (e.g., Fig. 5).

We suggest that the lack of high-frequency, super-
inertial motions in the model, excepting tidal frequen-
cies that operate at larger scales, is responsible for the
lack of variance at small scales (Fig. 6). Although the
numerical time step of the model, at 25 s, is sufficiently
fine to resolve these motions, the model is forced with
6-hourly reanalysis winds, which do not input energy at
sufficiently high frequencies. Recent work has sug-
gested that models that are subject to surface forcing at
superinertial frequencies develop greater variance at
superinterial frequencies, which from Fig. 6 we suggest
will translate into more realistic properties at small
spatial scales (Rimac et al. 2013).

b. Submesoscale motions and internal waves

Internal waves can affect the structure function slopes and magnitudes presented here. The effect of these waves is suppressed when calculating spice along isopycnals rather than depth surfaces, but will still be present for velocity and buoyancy calculations. The superinertial component of the velocity and buoyancy variance (*D*<sub>KE</sub> and *D*<sub>PE</sub>) was found to have a universally shallower spectrum and to be comparable or larger in magnitude than the subinertial component at scales less than about 5 km. This is, however, an imperfect method of disentangling the internal wave signal from the submesoscale dynamics, since balanced submesoscale motions may extend into superinertial frequencies (Torres et al. 2018). It is beyond the scope of this study to attempt a more sophisticated filtering of the internal waves and balanced submesoscale portion of the velocity or buoyancy signal.

c. Implications for theoretical models of vertical structure of passive tracers

As reviewed in the introduction, there exists con-
siderable uncertainty over the slopes of KE, PE, and
passive tracer (spice) variance in the ocean interior. In
particular, frontogenesis can lead to passive tracer
slopes of k⁻² at the surface. Under QG theory, far
from the surface, KE and PE spectra are predicted to
have slopes similar to k⁻³, which predicts passive
tracers' slopes that flatten from k⁻² to k⁻¹. A spice SF
slope of 0.5, as observed in this study, corresponds to a
KE spectral slope of k⁻¹.⁵, which is shallower than that
shown in previous studies in other parts of the ocean
(Cole and Rudnick 2012; Schönau and Rudnick 2015;
Itoh and Rudnick 2017; Kunze et al. 2015). Some of

these other studies are in regions with significant frontal
structures, such as the Pacific Ocean North Equatorial
Current (Schönau and Rudnick 2015) or the California
Current System (Itoh and Rudnick 2017). In contrast,
the OSMOSIS site is in a region of moderate eddy ki-
netic energy without a large sustained mean flow. The
Kunze et al. (2015) study spanned scales from 5 m to
50 km. At scales that match our observations they could
not distinguish between a spice gradient slope of k⁰ and
k¹/³; this latter corresponds to a SF slope of 2/3, which is
close to our estimated value. Cole and Rudnick (2012)
interpolate their glider data in the horizontal to find
passive tracer slopes of k⁻². This interpolation, though
required to perform horizontal wavenumber spectral
decompositions, can significantly steepen calculated
slopes (Callies and Ferrari 2013; Klymak et al. 2015).
Finally, Klymak et al. (2015) found spice spectral
slopes of k⁻¹.⁵, which is similar to our results, near the
surface, but these became significantly steeper with
depth, unlike the uniform slope we find here. Sum-
marizing these results, there appears to be a tendency
for the k⁻² slope in passive tracer spectra to occur in
regions with strong frontal structure, at least at scales
of 1–20 km. Our results further show that the spectral
slopes in the OSMOSIS region, which lacks the strong
frontal nature of many of the previous studies, do not
exhibit significant vertical structure. Adaptation of
SFs to other datasets may help to further explore the
regional variability of tracer distributions and the
mechanisms that give rise to them.

In the context of structure functions, for a KE SF slope of *n* the predicted passive tracer SF slope is 1 − *n*/2. KE and PE slopes estimated from observational data at the OSMOSIS site are generally between 0.6 and 1, which implies a passive tracer slope of around 0.5, close to the slopes found here. However, this theoretical prediction is not expected to hold for high-frequency motions such as internal waves, which were shown to be important here (Fig. 6).

d. Seasonality in submesoscale variance

As is typical in the extratropics, the OSMOSIS region
exhibits pronounced seasonality in MLD and the
strength of the pycnocline at the base of the mixed layer
(Fig. 2; see also Erickson and Thompson 2018). Our SF
observations show that this seasonality entails a seasonal
cycle in submesoscale energy levels, diagnosed as a
wintertime increase in KE and PE magnitude (Fig. 9).
While the seasonal cycle in MLD is faithfully repro-
duced in the model, the seasonality of the upper
ocean energetics differs: the model submesoscale KE
(α) peaks in spring rather than winter, and it remains
high even into summer (Fig. 10a).
A key feature of the winter-to-spring transition is the restratification of the upper ocean, which can be caused by direct surface forcing but also via baroclinic mixed layer instabilities (Haine and Marshall 1998; Fox-Kemper et al. 2008; Mahadevan et al. 2012; du Plessis et al. 2017). These instabilities occur at scales of 1–20 km and are therefore only marginally resolved in the model. This may partially explain the delay in spring restratification in the model, as compared with the observations (Fig. 2a). However, resolution issues are unlikely to explain the late and pronounced peak in modeled submesoscale KE—if anything, damped mixed layer instabilities should lead to a more muted seasonality.

If the difference between the model and the observations cannot be attributed to interannual variability, it suggests that the model is missing a sink of submesoscale energy. One possibility is that interactions between mixed layer turbulence and submesoscale eddies are misrepresented by the K-profile parameterization (KPP) (Large et al. 1994) used in the model (Menemenlis et al. 2008). Another possibility is that internal waves drain a significant amount of energy from submesoscale balanced flow (e.g., Taylor and Straub 2016; Barkan et al. 2017; Xie and Vanneste 2015; Rocha et al. 2018), an effect that would be underestimated in the model because the internal-wave field—particularly the near-inertial component—is unrealistically weak (Figs. 5 and 6), probably as a result from the relatively coarse (6 h) wind variability (Yu et al. 2019b; Flexas et al. 2019).

No pronounced seasonality is evident in the SF slopes estimated from the observations (Fig. 8). If anything, KE and PE slopes tend to be somewhat smaller in summer than in winter, in apparent conflict with previous shipboard ADCP observations from the western subtropical North Atlantic that showed a steeping of submesoscale wavenumber spectra in summer (Callies et al. 2015). It should be remembered, however, that the slopes considered in this study are related to scales smaller than those where Callies et al. (2015) observed steepening, and that high-frequency motion contributes significantly to the small-scale end of the observations (Figs. 6a,b). Similar to the western subtropical North Atlantic, relatively flat SFs across an entire year might emerge due to different processes in winter and summer; namely, energetic balanced submesoscales in winter and internal waves in summer. This highlights a limitation of the SF (or spectral) slopes approach, which is often insufficient to discriminate different dynamical regions when employed in isolation.

The model, in contrast to the observations, shows much steeper structure functions in summer than in winter

(Figs. 8a,b). We attribute this to the weakness of the internal-wave field in the model, which fails to provide enough small-scale variance to maintain flat SFs in summer.

### e. Submesoscale tracer fluxes out of the surface mixed layer

If submesoscale motions were largely confined to the mixed layer, a transition in passive tracer properties, such as the SF magnitudes and slopes calculated from Eq. (5), might be expected to occur at the base of the mixed layer in response to differences in stirring across the pycnocline. In fact, there is a slight change in properties associated with the shallow summertime pycnocline, but very little change across the base of the wintertime mixed layer for SFs calculated from the observational data. The model, however, shows a stronger change across the wintertime mixed layer base, which is most apparent in the changes of slope of $D_{PE}$ (Fig. 8b).

The effect of model resolution can be seen in a series of recent simulations performed by Balwada et al. (2018), where an idealized MITgcm was run at horizontal resolutions between 20 and 1 km. As the resolution increased to finer scales, stratification at the base of the mixed layer also increased, but so did the vertical transport of tracers across the mixed layer base. The authors attributed this increase in vertical transport to increased strength of baroclinic mixed layer instabilities at higher resolution. However, even at horizontal resolution of 1 km, vertical transport estimates did not converge, implying that scales less than about 5 km (the effective resolution being 4–5 times the gridscale resolution) play an important role in transporting water properties across the mixed layer base.

The high spatial resolution of horizontal velocities by the mooring array also allows a calculation of the vertical velocity through the density conservation equation. Yu et al. (2019a) used this technique to estimate strong vertical velocities, of up to 100 m day⁻¹, extending 200 m below the mixed layer during winter and spring. Vertical velocities from the llc4320 also support this observation (see Yu et al. 2019a, appendix D). These vertical velocities can lead to injection of mixed layer water with small passive tracer variance below the mixed layer, as seen in the lack of vertical structure in the winter to depths exceeding the mixed layer depth.

In this study, we used spice as a passive tracer, and found that there is not a clear sign of spice variance changing across the base of the mixed layer during winter, either in slope or magnitude. Other passive tracers, such as oxygen, nutrients, and dissolved organic carbon, can also be transported across the base
of the wintertime mixed layer in a similar fashion. Erickson and Thompson (2018) found evidence of rapid downward transport of high-oxygen waters below the wintertime mixed layer using OSMOSIS glider data, and attributed this downward flux to intermittent submesoscale instabilities, with a particular focus on symmetric instability. These instabilities provide a mechanism both for upward transport of nutrients as well as downward transport of neutrally buoyant carbon, especially during winter.

## 6. Conclusions

The OSMOSIS project provides a useful testing ground to consider seasonality in submesoscale dynamics. Second-order structure functions (SFs) of spice, buoyancy, and horizontal velocities give similar information to spectral decompositions but are better suited to data from arrays of moorings and gliders that are not aligned in a single transect.

This observational dataset is the first to enable a statistical description of turbulence properties down to scales of 1 km over a full seasonal cycle to 1-km depth. Thus we were able to probe scales that would be expected to be only marginally resolved even in relatively high-resolution global models, for example the 1/48° llc4320 analyzed in this study.

We find a difference between the observational and model results at small scales; the model output contains less variance than the observations. This is expected, as the effective model resolution is approximately a factor of 4 larger than the grid scale resolution. However, the model did compare well with the mooring data using only the subinertial components of velocity, buoyancy, and spice variables, even down to very small separations (2 km). This points to the lack of high-frequency motions, potentially related to the low-frequency surface forcing (Rimac et al. 2013), as being a crucial factor for missing small-scale variance in the model simulation.

The model SFs agree more closely with the observational results during winter than in summer. This is likely to occur because the model more accurately resolves mixed layer instabilities in the winter than the summer, as surface-enhanced submesoscale motions are larger in size when the mixed layer is deeper. In particular, the model retains a high level of submesoscale KE throughout spring and into summer, suggesting that suppressing the resolution of submesoscale motions leads to an inefficient transfer of energy to dissipative scales.

We observe a seasonality of submesoscale energy, with higher KE and PE in winter than in summer. Despite strong seasonality in KE and PE, the statistical

representation of submesoscale motions with the SFs does not exhibit a significant seasonality. This suggests either that the physical processes that are shaping the SF slope are insensitive to the partitioning of KE and PE, or that similar SF slopes can arise from different physical processes. This behavior has been highlighted previously when analyzing spectral slopes of submesoscale turbulence (Callies and Ferrari 2013).

Many previous studies have found, as in this study, that the spectral slope of tracers is relatively uniform with depth in the upper ocean. However, most of these results have reported spectral slopes of $k^{-2}$, steeper than QG theory. Here we find a uniform slope that is closer to $k^{-3/2}$, and present this dataset as a counterexample to other studies suggesting a universal $k^{-2}$ passive tracer spectral slope at depth. The passive tracer slopes we find are still far from the $k^{-1}$ expected by QG theory. However, we note that our observations to 1000 m have not penetrated below the permanent pycnocline starting near 800 m.

The observed tracer variance in the ocean's upper 1000 m in this and other studies suggest the importance of motions that can remove tracer variance below the base of the mixed layer or can generate smaller-scale stirring than is typically associated with the interior ocean. The dynamics supporting these tracer distributions will be especially important for biological properties such as particulate or dissolved carbon, where vertical exchange across the mixed layer base may lead to sequestration and carbon export on long time scales. Our observations suggest that the boundary between the ocean-atmosphere interface and the ocean interior—the pycnocline at the base of the mixed layer—is not well developed during the wintertime, meaning the conventional distinction between well-mixed surface waters and an ocean interior out of contact with the atmosphere may not apply during all seasons.

*Acknowledgments.* We are grateful for the efforts of the entire OSMOSIS team and the crews of the Royal Research Ship (RRS) *Discovery*, the R/V *Celtic Explorer*, and the RRS *James Cook*. Glider data are held at the British Oceanographic Data Centre and are accessible from https://doi.org/10/cqc6. Instructions for obtaining LLC4320 model output can be found at http://ecco2.org/llc_hires. We thank four anonymous reviewers for their comments on this manuscript. The OSMOSIS project was funded by NERC Grant NE/I019905/1 and NSF OCE-1155676. ZKE and AFT acknowledge funding from the David and Lucille Packard Foundation and from NASA Award 80NSSC17K0663. XY is supported by the French National Agency for Research (ANR) under Grant 17-CE01-0006-01. PK is supported by a NASA senior fellowship.
REFERENCES

Babiano, A., C. Basdevant, and R. Sadourny, 1985: Structure functions and dispersion laws in two-dimensional turbulence. *J. Atmos. Sci.*, **42**, 941–949, https://doi.org/10.1175/1520-0469(1985)042<0941:SFADLI>2.0.CO;2

Balwada, D., J. H. LaCasce, and K. G. Speer, 2016: Scale-dependent distribution of kinetic energy from surface drifters in the Gulf of Mexico. *Geophys. Res. Lett.*, **43**, 10 856–10 863, https://doi.org/10.1002/2016GL069405.

—, K. S. Smith, and R. Abernathey, 2018: Submesoscale vertical velocities enhance tracer subduction in an idealized Antarctic Circumpolar Current. *Geophys. Res. Lett.*, **45**, 9790–9802, https://doi.org/10.1029/2018GL079244.

Barkan, R., K. B. Winters, and J. C. McWilliams, 2017: 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.

Brannigan, L., D. P. Marshall, A. Naveira-Garabato, and A. J. G. Nurser, 2015: The seasonal cycle of submesoscale flows. *Ocean Modell.*, **92**, 69–84, https://doi.org/10.1016/j.ocemod.2015.05.002.

Buckingham, C. E., and Coauthors, 2016: Seasonality of submesoscale flows in the ocean surface boundary layer. *Geophys. Res. Lett.*, **43**, 2118–2126, https://doi.org/10.1002/2016GL068009.

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.

—, —, J. M. Klymak, and J. Gula, 2015: Seasonality in submesoscale turbulence. *Nat. Commun.*, **6**, 6862, https://doi.org/10.1038/ncomms7862.

—, G. Flierl, R. Ferrari, and B. Fox-Kemper, 2016: The role of mixed-layer instabilities in submesoscale turbulence. *J. Fluid Mech.*, **788**, 5–41, https://doi.org/10.1017/jfm.2015.700.

Capet, X., J. C. McWilliams, M. J. Molemaker, and A. Shchepetkin, 2008: 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.

Charney, J. G., 1971: Geostrophic turbulence. *J. Atmos. Sci.*, **28**, 1087–1095, https://doi.org/10.1175/1520-0469(1971)028<1087:GT>2.0.CO;2.

Cole, S., and D. Rudnick, 2012: The spatial distribution and annual cycle of upper ocean thermohaline structure. *J. Geophys. Res.*, **115**, C04012, https://doi.org/10.1029/2009JC005654.

Damerell, G. M., K. J. Heywood, A. F. Thompson, U. Binetti, and J. Kaiser, 2016: The vertical structure of upper ocean variability at the Porcupine Abyssal Plain during 2012–2013. *J. Geophys. Res. Oceans*, **121**, 3075–3089, https://doi.org/10.1002/2015JC011423.

de Boyer Montégut, C., G. Madec, A. S. Fischer, A. Lazar, and D. Iudicone, 2004: Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. *J. Geophys. Res.*, **109**, C12003, https://doi.org/10.1029/2004JC002378.

du Plessis, M., S. Swart, I. Ansorge, and A. Mahadevan, 2017: Submesoscale processes promote seasonal restratification in the Subantarctic Ocean. *J. Geophys. Res. Oceans*, **122**, 2960–2975, https://doi.org/10.1002/2016JC012494.

Erickson, Z. K., and A. Thompson, 2018: The seasonality of physically-driven export at submesoscales in the northeast Atlantic Ocean. *Global Biogeochem. Cycles*, **32**, 1144–1162, https://doi.org/10.1029/2018GB005927.

Ferrari, R., and C. Wunsch, 2009: Ocean circulation kinetic energy: Reservoirs, sources, and sinks. *Annu. Rev. Fluid Mech.*, **41**, 253–282, https://doi.org/10.1146/annurev.fluid.40.111406.102139.

Flexas, M. M., A. F. Thompson, H. S. Torres, P. Klein, J. T. Farrar, H. Zhang, and D. Menemenlis, 2019: Global estimates of the energy transfer from the wind to the ocean, with emphasis on near-inertial oscillations. *J. Geophys. Res. Oceans*, **124**, 5723–5746, https://doi.org/10.1029/2018JC014453.

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.

Frajka-Williams, E., C. C. Eriksen, P. B. Rhines, and R. R. Harcourt, 2011: Determining vertical water velocities from Seaglider. *J. Atmos. Oceanic Technol.*, **28**, 1641–1656, https://doi.org/10.1175/2011JTECHO830.1.

Garrett, C., and W. Munk, 1975: Space-time scales of internal waves: A progress report. *J. Geophys. Res.*, **80**, 291–297, https://doi.org/10.1029/JC080i003p00291.

Haine, T. W., and J. Marshall, 1998: Gravitational, symmetric, and baroclinic instability of the ocean mixed layer. *J. Phys. Oceanogr.*, **28**, 634–658, https://doi.org/10.1175/1520-0485(1998)028<0634:GSABIO>2.0.CO;2.

Hua, B. L., C. Ménesguen, S. Le Gentil, R. Schopp, B. Marsset, and H. Aiki, 2013: Layering and turbulence surrounding an anticlonic oceanic vortex: In situ observations and quasi-geostrophic numerical simulations. *J. Fluid Mech.*, **731**, 418–442, https://doi.org/10.1017/jfm.2013.369.

Itoh, S., and D. L. Rudnick, 2017: Fine-scale variability of isopycnal salinity in the California Current System. *J. Geophys. Res. Oceans*, **122**, 7066–7081, https://doi.org/10.1002/2017JC013080.

Klein, P., A.-M. Treguier, and B. L. Hua, 1998: Three-dimensional stirring of thermohaline fronts. *J. Mar. Res.*, **56**, 589–612, https://doi.org/10.1357/002224098765213595.

—, B. L. Hua, G. Lapeyre, X. Capet, S. Le Gentil, and H. Sasaki, 2008: Upper ocean turbulence from high-resolution 3D simulations. *J. Phys. Oceanogr.*, **38**, 1748–1763, https://doi.org/10.1175/2007JPO3773.1.

Klymak, J. M., W. Crawford, M. H. Alford, J. A. MacKinnon, and R. Pinkel, 2015: Along-isopycnal variability of spice in the North Pacific. *J. Geophys. Res. Oceans*, **120**, 2287–2307, https://doi.org/10.1002/2013JC009421.

Kunze, E., J. Klymak, R.-C. Lien, R. Ferrari, C. Lee, M. Sundermeyer, and L. Goodman, 2015: Submesoscale water-mass spectra in the Sargasso Sea. *J. Phys. Oceanogr.*, **45**, 1325–1338, https://doi.org/10.1175/JPO-D-14-0108.1.

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.

Lévy, M., R. Ferrari, P. J. Franks, A. P. Martin, and P. Rivière, 2012: Bringing physics to life at the submesoscale. *Geophys. Res. Lett.*, **39**, L14602, https://doi.org/10.1029/2012GL052756.

Mahadevan, A., E. D'asaro, C. Lee, and M. J. Perry, 2012: Eddy-driven stratification initiates North Atlantic spring phytoplankton blooms. *Science*, **337**, 54–58, https://doi.org/10.1126/science.1218740.

McCaffrey, K., B. Fox-Kemper, and G. Forget, 2015: Estimates of ocean macroturbulence: Structure function and spectral slope
from Argo profiling floats. *J. Phys. Oceanogr.*, **45**, 1773-1793, https://doi.org/10.1175/JPO-D-14-0023.1.

McCartney, M. S., and L. D. Talley, 1982: The subpolar mode water of the North Atlantic Ocean. *J. Phys. Oceanogr.*, **12**, 1169-1188, https://doi.org/10.1175/1520-0485(1982)012<1169: TSMWOT>2.0.CO;2.

McDougall, T. J., and O. A. Krzysik, 2015: Spiciness. *J. Mar. Res.*, **73**, 141–152, https://doi.org/10.1357/002224015816665589.

McWilliams, J. C., 1985: Submesoscale, coherent vortices in the ocean. *Rev. Geophys.*, **23**, 165–182, https://doi.org/10.1029/ RG023i002p00165.

—, 2016: Submesoscale currents in the ocean. *Proc. Roy. Soc. London*, **472A**, 20160117, 10.1098/rspa.2016.011710.1098/ rspa.2016.0117.

Menemenlis, D., J. Campin, P. Heimbach, C. Hill, T. Lee, A. Nguyen, M. Schodlok, and H. Zhang, 2008: ECCO2: High resolution global ocean and sea ice data synthesis. Mercator Ocean Quarterly Newsletter, No. 31, Mercator-Ocean, Rarmonville-Saint-Agne, France, 13–21.

Mensa, J. A., Z. Garraffo, A. Griffa, T. M. Özgökmen, A. Haza, and M. Veneziani, 2013: Seasonality of the submesoscale dynamics in the Gulf Stream region. *Ocean Dyn.*, **63**, 923–941, https://doi.org/10.1007/s10236-013-0633-1.

Munk, W., 1981: Internal waves and small-scale processes. *Evolution of Physical Oceanography*, B. A. Warren, and C. Wunsch, Eds., MIT Press, 264–291.

Qiu, B., S. Chen, P. Klein, J. Wang, H. Torres, L.-L. Fu, and D. Menemenlis, 2018: Seasonality in transition scale from balanced to unbalanced motions in the world ocean. *J. Phys. Oceanogr.*, **48**, 591–605, https://doi.org/10.1175/JPO-D-17-0169.1.

Rieck, J. K., C. W. Böning, R. J. Greatbatch, and M. Scheinert, 2015: Seasonal variability of eddy kinetic energy in a global high-resolution ocean model. *Geophys. Res. Lett.*, **42**, 9379–9386, https://doi.org/10.1002/2015GL066152.

Rimac, A., J.-S. von Storch, C. Eden, and H. Haak, 2013: The influence of high-resolution wind stress field on the power input to near-inertial motions in the ocean. *Geophys. Res. Lett.*, **40**, 4882–4886, https://doi.org/10.1002/grl.50929.

Rocha, C. B., S. T. Gille, T. K. Chereskin, and D. Menemenlis, 2016: Seasonality of submesoscale dynamics in the Kuroshio Extension. *Geophys. Res. Lett.*, **43**, 11304–11311, https://doi.org/10.1002/2016GL071349.

—, G. L. Wagner, and W. R. Young, 2018: Stimulated generation: extraction of energy from balanced flow by near-inertial waves. *J. Fluid Mech.*, **847**, 417–451, https://doi.org/10.1017/jfm.2018.308.

Roullet, G., X. Capet, and G. Maze, 2014: Global interior eddy available potential energy diagnosed from Argo floats. *Geophys. Res. Lett.*, **41**, 1651–1656, https://doi.org/10.1002/2013GL059004.

Sasaki, H., P. Klein, B. Qiu, and Y. Sasai, 2014: Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere. *Nat. Commun.*, **5**, 5636, https://doi.org/10.1038/ncomms6636.

Schönau, M. C., and D. L. Rudnick, 2015: Glider observations of the North Equatorial Current in the western tropical Pacific. *J. Geophys. Res. Oceans*, **120**, 3586–3605, https://doi.org/10.1002/2014JC010595.

Siegelman, L., M. O'Toole, M. Flexas, P. Rivière, and P. Klein, 2019: Submesoscale ocean fronts act as biological hotspot for southern elephant seal. *Sci. Rep.*, **9**, 5588, https://doi.org/10.1038/s41598-019-42117-w.

Smith, K. S., and R. Ferrari, 2009: The production and dissipation of compensated thermohaline variance by mesoscale stirring. *J. Phys. Oceanogr.*, **39**, 2477–2501, https://doi.org/10.1175/2009JPO4103.1.

Su, Z., J. Wang, P. Klein, A. F. Thompson, and D. Menemenlis, 2018: Ocean submesoscales as a key component of the global heat budget. *Nat. Commun.*, **9**, 775, https://doi.org/10.1038/s41467-018-02983-w.

Taylor, S., and D. Straub, 2016: Forced near-inertial motion and dissipation of low-frequency kinetic energy in a wind-driven channel flow. *J. Phys. Oceanogr.*, **46**, 79–93, https://doi.org/10.1175/JPO-D-15-0060.1.

Thomas, L. N., J. R. Taylor, E. A. D'Asaro, C. M. Lee, J. M. Klymak, and A. Shcherbina, 2016: Symmetric instability, inertial oscillations, and turbulence at the Gulf Stream front. *J. Phys. Oceanogr.*, **46**, 197–217, https://doi.org/10.1175/JPO-D-15-0008.1.

Thompson, A. F., A. Lazar, C. Buckingham, A. C. Naveira Garabato, G. M. Damerell, and K. J. Heywood, 2016: Open-ocean submesoscale motions: A full seasonal cycle of mixed layer instabilities from gliders. *J. Phys. Oceanogr.*, **46**, 1285–1307, https://doi.org/10.1175/JPO-D-15-0170.1.

Torres, H. S., P. Klein, D. Menemenlis, B. Qiu, Z. Su, J. Wang, S. Chen, and L.-L. Fu, 2018: Partitioning ocean motions into balanced motions and internal gravity waves: A modeling study in anticipation of future space missions. *J. Geophys. Res. Oceans*, **123**, 8084–8105, https://doi.org/10.1029/2018JC014438.

Vallis, G. K., 2006: *Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation*. Cambridge University Press, 769 pp.

Veronis, G., 1972: Properties of seawater defined by temperature, salinity, and pressure. *J. Mar. Res.*, **30** (2), 227.

Viglione, G. A., A. F. Thompson, M. M. Flexas, J. Sprintall, and S. Swart, 2018: Abrupt transitions in submesoscale structure in Southern Drake Passage: Glider observations and model results. *J. Phys. Oceanogr.*, **48**, 2011–2027, https://doi.org/10.1175/JPO-D-17-0192.1.

Wang, J., L.-L. Fu, B. Qiu, D. Menemenlis, J. T. Farrar, Y. Chao, A. F. Thompson, and M. M. Flexas, 2018: An observing system simulation experiment for the calibration and validation of the surface water ocean topography sea surface height measurement using in situ platforms. *J. Atmos. Oceanic Technol.*, **35**, 281–297, https://doi.org/10.1175/JTECH-D-17-0076.1.

Webb, E., 1964: Ratio of spectrum and structure-function constants in the inertial subrange. *Quart. J. Roy. Meteor. Soc.*, **90**, 344–346, https://doi.org/10.1002/qj.49709038520.

Xie, J.-H., and J. Vanneste, 2015: A generalised-lagrangian-mean model of the interactions between near-inertial waves and mean flow. *J. Fluid Mech.*, **774**, 143–169, https://doi.org/10.1017/jfm.2015.251.

Yu, X., A. C. Naveira Garabato, A. P. Martin, C. E. Buckingham, L. Brannigan, and Z. Su, 2019a: An annual cycle of submesoscale vertical flow and restratification in the upper ocean. *J. Phys. Oceanogr.*, **49**, 1439–1461, https://doi.org/10.1175/JPO-D-18-0253.1.

—, A. L. Ponte, S. Elipot, D. Menemenlis, E. D. Zaron, and R. Abernathey, 2019b: Surface kinetic energy distributions in the global oceans from a high-resolution numerical model and surface drifter observations. *Geophys. Res. Lett.*, **46**, 9757–9766, https://doi.org/10.1029/2019GL083074.