# The Annual Cycle of Upper-Ocean Potential Vorticity and Its Relationship to Submesoscale Instabilities

Xiaolong Yu, Alberto C. Naveira Garabato, Adrian P. Martin, David P.

Marshall

► To cite this version:

Xiaolong Yu, Alberto C. Naveira Garabato, Adrian P. Martin, David P. Marshall. The Annual Cycle of Upper-Ocean Potential Vorticity and Its Relationship to Submesoscale Instabilities. Journal of Physical Oceanography, 2021, 51 (2), pp.385-402. 10.1175/JPO-D-20-0099.1. hal-04203489

HAL Id: hal-04203489

https://hal.science/hal-04203489v1

Submitted on 15 Sep 2023

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 CC BY 4.0 - Attribution - International License
The Annual Cycle of Upper-Ocean Potential Vorticity and Its Relationship to Submesoscale Instabilities

XIAOLONG YU,ª ALBERTO C. NAVEIRA GARABATO,ᵇ.c ADRIAN P. MARTIN,ᶜ AND DAVID P. MARSHALLᵈ

(Manuscript received 11 May 2020, in final form 27 October 2020)

**ABSTRACT:** The evolution of upper-ocean potential vorticity (PV) over a full year in a typical midocean area of the northeast Atlantic is examined using submesoscale- and mesoscale-resolving hydrographic and velocity measurements from a mooring array. A PV budget framework is applied to quantitatively document the competing physical processes responsible for deepening and shoaling the mixed layer. The observations reveal a distinct seasonal cycle in upper-ocean PV, characterized by frequent occurrences of negative PV within deep (up to about 350 m) mixed layers from winter to mid-spring, and positive PV beneath shallow (mostly less than 50 m) mixed layers during the remainder of the year. The cumulative positive and negative subinertial changes in the mixed layer depth, which are largely unaccounted for by advective contributions, exceed the deepest mixed layer by one order of magnitude, suggesting that mixed layer depth is shaped by the competing effects of destratifying and restratifying processes. Deep mixed layers are attributed to persistent atmospheric cooling from winter to mid-spring, which triggers gravitational instability leading to mixed layer deepening. However, on shorter time scales of days, conditions favorable to symmetric instability often occur as winds intermittently align with transient frontal flows. The ensuing submesoscale frontal instabilities are found to fundamentally alter upper-ocean turbulent convection, and limit the deepening of the mixed layer in the winter-to-mid-spring period. These results emphasize the key role of submesoscale frontal instabilities in determining the seasonal evolution of the mixed layer in the open ocean.

**KEYWORDS:** North Atlantic Ocean; Instability; Potential vorticity; Mixed layer; In situ oceanic observations; Time series

**1. Introduction**

The ocean's surface mixed layer is a vital component of the global climate system, as it mediates the exchanges of physical and biogeochemical tracers between the atmosphere and the ocean interior. As a result of its direct contact with the atmosphere, the mixed layer is continuously influenced by solar radiation; air-sea heat; and freshwater transfers, winds, and waves at the atmosphere-ocean boundary. This causes the mixed layer to typically host active turbulent mixing, weak vertical stratification, and nearly uniform vertical tracer (e.g., temperature and salinity) distributions. The temporal variability of the mixed layer thickness is important in determining the rates of water mass formation, the vertical structure of oceanic properties, and biological productivity. However, it remains challenging to accurately reproduce the deepening and shoaling of the mixed layer in climate-scale ocean models on a range of time scales (e.g., Fox-Kemper et al. 2011; Sallée et al. 2013). A primary reason for this difficulty is the fundamentally multiscale character of the dynamics shaping mixed layer evolution. Recent in situ and satellite observations reveal the occurrence of a wide spectrum of mixed layer-controlling processes spanning horizontal scales from millimeters to many hundreds of kilometers (Ferrari 2011; Belcher et al. 2012; Klein et al. 2019). Thus, an essential step in advancing model

representations of the mixed layer is to parameterize the key cross-scale processes governing upper-ocean mixing and restratification. These are simplistically approximated as one-dimensional in the vertical by most current parameterizations (D'Asaro 2014).

Submesoscales provide the most glaring example of a class of mixed layer-controlling processes that are not yet recognized by the mixed layer parameterization schemes used in most ocean models. Submesoscale flows are ubiquitous within the mixed layer throughout the global ocean (McWilliams 2016). They are manifested at horizontal scales of 0.1–10 km and time scales of several hours to several days, and are dynamically associated with Rossby and Richardson numbers of *O*(1) (Thomas et al. 2008). Importantly for mixed layer evolution, a variety of frontal instabilities may occur at the submesoscale for which vorticity, divergence and strain can be locally intense (i.e., comparable to the local Coriolis frequency). The associated ageostrophic motions may convert lateral buoyancy gradients into vertical stratification via enhanced upward buoyancy transport (Klein and Lapeyre 2009; Su et al. 2018, 2020; Siegelman et al. 2020), triggering restratification and shoaling of the mixed layer.

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

This article is licensed under a Creative Commons Attribution 4.0 license (http://creativecommons.org/licenses/by/4.0/).

**Publisher's Note:** This article was revised on March 11 to designate it as open access.

**Publisher's Note:** This article was revised on 12 May 2021 to include the CCBY license that was missing when originally published.

<footnote>ªIfremer, Univ. Brest, CNRS, IRD, Laboratoire d'Océanographie Physique et Spatiale, IUEM, Brest, France</footnote>

<footnote>ᵇOcean and Earth Science, University of Southampton, Southampton, England</footnote>

<footnote>ᶜNational Oceanography Centre, Southampton, United Kingdom</footnote>

<footnote>ᵈDepartment of Physics, University of Oxford, Oxford, United Kingdom</footnote>

<footnote>ª Denotes content that is immediately available upon publication as open access.</footnote>

<footnote>Corresponding author: Xiaolong Yu, xiaolong.yu@ifremer.fr</footnote>

<footnote>DOI: 10.1175/JPO-D-20-0099.1</footnote>
Among submesoscale frontal instabilities, submesoscale baroclinic instability is widely regarded as particularly important in restratifying the mixed layer (Boccaletti et al. 2007; Fox-Kemper et al. 2008; Callies et al. 2016). As demonstrated by these numerical studies, the instability drives frontal slumping and restratification through the release of available potential energy by an eddy-induced overturning circulation. This effect has been parameterized as an overturning streamfunction within the mixed layer (Fox-Kemper et al. 2008). Another submesoscale frontal instability that has come into focus more recently is symmetric instability (SI), which arises from the interaction between destabilizing atmospheric forcing and submesoscale fronts. Special attention has been paid to the case of downfront winds (i.e., winds oriented in the direction of the front's geostrophic shear), which destabilize the water column as Ekman flow moves water from the dense side of the front over lighter water in the light side of the front, thereby preconditioning the flow to SI (e.g., D'Asaro et al. 2011; Thomas et al. 2013). Generalization of this problem to the forcing of a submesoscale front by either downfront winds or surface buoyancy loss (Taylor and Ferrari 2010) predicts an oceanic response with two layers of distinct dynamics. In a near-surface convective layer, available potential energy is the dominant energy source for overturning motions, convective mixing develops, and the water column remains gravitationally unstable. Ageostrophic shear production may also play an important role in the turbulent kinetic energy budget of this convective layer (Skyllingstad et al. 2017). Beneath the convective layer and above the base of the surface boundary layer (hereinafter, the forced-SI layer), the dominant energy source is instead the background vertical shear, slanted overturning motions linked to SI dominate over convective mixing, and the water column is restratified. All in all, the mixed layer impacts of submesoscale frontal instabilities have been extensively assessed with theoretical approaches and in numerical models of varying complexity [see McWilliams (2019) for a recent review]. However, rigorously testing such assessments against observations has proven more problematic, due to the “snapshot” nature of many observations (e.g., Adams et al. 2017) and the few constraints on dynamical budgets afforded by the observational time series–based analyses available to date (Thompson et al. 2016; du Plessis et al. 2019).

Ertel potential vorticity (PV; Ertel 1942; Schubert et al. 2004) provides a natural framework within which to appraise the degree of realism of theoretical and numerical predictions, for two reasons. First, the dynamical behavior of geophysical fluid systems is commonly expressed in terms of PV (Hoskins et al. 1985). Second, in the ocean, PV is materially conserved along Lagrangian trajectories in the absence of forcing and dissipation. [Equivalently, in the flux form of the conservation equation (Haynes and McIntyre 1987, 1990), PV substance can only be injected or extracted through the ocean's boundaries (Marshall and Nurser 1992).] These properties of PV have been exploited by, for example, Thomas et al. (2008) and Brannigan (2016), who analyzed numerical simulations to show that SI acts as an upper-ocean PV pump, upwelling high-PV water from the pycnocline and subducting low-PV water from the mixed layer. Similarly, other authors (e.g., Thomas 2005; Czaja and Hausmann 2009; Maze et al. 2013; Deremble et al. 2014; Wenegrat et al. 2018), using both numerical models and climatologies, have focused on the diabatic and frictional modifications of upper-ocean

PV in order to gauge the relative contributions of those forcings to driving the ocean circulation, as well as the role of submesoscale frontal instabilities in shaping the oceanic response.

Here, we build on and expand this body of work by diagnosing the annual cycle of upper-ocean PV, and its regulation by submesoscale frontal instabilities, in a typical midocean region using observations. Our overarching goal is to quantitatively assess current views—largely grounded on theoretical and modeling investigations—on the processes governing midlatitude mixed layer evolution on time scales of days to seasons. Our approach is to construct a budget of PV in a 13 km × 13 km × 500 m upper-ocean volume using measurements from a mooring array deployed in the northeast Atlantic as part of the OSMOSIS (Ocean Surface Mixing, Ocean Submesoscale Interaction Study) experiment, complemented with glider observations and an atmospheric reanalysis. The mooring dataset is exceptional in that it samples a complete annual cycle of the mixed layer evolution concurrently to many of its controlling 3D processes (such as submesoscale frontal instabilities).

This work represents the culmination of a long succession of studies of submesoscale turbulence based on the OSMOSIS mooring and glider observations, and integrates insights and diagnostics generated by those investigations. These include the work of Buckingham et al. (2016), who demonstrated the seasonality of submesoscale motions, evidenced by the wintertime occurrence of positively skewed relative vorticity in mooring data, and that of Thompson et al. (2016), who characterized the PV conditions for submesoscale frontal instabilities, and the instabilities' impact on upper-ocean stratification, using the year-long glider measurements. Yu et al. (2019a) followed by diagnosing the annual cycle of upper-ocean vertical motion and restratification associated with submesoscale processes from mooring observations, and showed that submesoscale restratification events are generally triggered by mesoscale frontogenesis. Evans et al. (2018) and Buckingham et al. (2019) used the glider and mooring datasets to examine the contribution of submesoscale frontal instabilities to sustaining turbulent kinetic energy dissipation in the mixed layer, and concluded that such contribution is generally modest. Erickson et al. (2020) highlighted the vertical penetration of wintertime submesoscale motions to depths well in excess of the mixed layer, by applying a horizontal structure function approach to both glider and mooring measurements. Finally, Callies et al. (2020) developed and applied a novel frequency-resolved horizontal structure function methodology to demonstrate that submesoscale flows have largely subinertial time scales.

Our work adds to this body of work by showing that surface forcing of fronts is centrally involved in symmetric and gravitational instabilities, and can fundamentally alter upper-ocean turbulent convection in the context of a PV budget. The paper is organized as follows. Sections 2 and 3 respectively introduce the data and theoretical PV framework. The annual cycle of upper-ocean PV is described in section 4, and the diagnosed PV budget is analyzed and discussed in section 5. Conclusions are offered in section 6.

## 2. Data

### a. Mooring data

The data analyzed in this study were primarily collected from a mooring array deployed at an approximate water depth
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Geographical location of the OSMOSIS mooring array (marked as a white filled star) in the northeast Atlantic. The annual mean of climatological mixed layer depth in 2° × 2° bins (de Boyer Montégut et al. 2004) is colored as background. (b) Locations of the central, four inner, and four outer moorings.</figcaption></figure>

of 4800 m over the Porcupine Abyssal Plain (PAP; 48.63°–48.75°N, 16.09°–16.27°W) in the northeast Atlantic (Fig. 1). The array's primary purpose was to measure the detailed evolution of the mixed layer, its controlling submesoscale processes and its mesoscale context, over a complete annual cycle. The mooring area was intentionally chosen to be representative of the midlatitude open ocean far away from western boundaries and complex topography, a regime that spans a substantial fraction of the global ocean (Fig. 1a). Nine subsurface moorings were deployed for the period September 2012–September 2013, arranged in two concentric quadrilaterals with side lengths of ~13 km (outer cluster) and 1–2 km (inner cluster) around a centrally located single mooring (Fig. 1b).

The mooring array design enabled simultaneous measurements to be made of horizontal flows on spatial scales of $O(1)$ km and $O(10)$ km, from the inner and outer mooring clusters, respectively. Mooring sensors comprised a series of paired MicroCAT conductivity–temperature–depth (CTD) sensors and Nortek Aquadopp acoustic current meters installed at depths spanning the approximate range 50–520 m. The central mooring was equipped with 13 CTD/current meter pairs, and the inner and outer moorings with 7 and 5 such pairs, respectively. The present study predominantly uses data from the CTD/current meter pairs. More detailed information on mooring instrumentation is provided by Yu et al. (2019a).

The moored instrumentation returned a full annual cycle of upper-ocean temperature, salinity, pressure, and horizontal velocity. The CTDs and current meters sampled at 5- and 10-min intervals, respectively. For each mooring, we linearly interpolate the data onto a uniform depth grid with 10-m spacing between depths of 50 and 520 m, and average onto hourly bins. Potential density (referenced to the ocean surface) and depth are calculated from interpolated temperature, salinity, and pressure using the Gibbs Seawater Oceanographic Toolbox (McDougall and Barker 2011). Compressibility effects are considered to be negligible over the top 520 m (Yu et al. 2019a). A fourth-order low-pass

Butterworth filter with a cutoff of one inertial period (16 h) is applied to the hourly data to remove unbalanced motions, such as internal tides, near-inertial flows and other high-frequency motions.

A quality control of the mooring data is carried out prior to analysis. For the available measurements, missing values occasionally occurred, especially for salinity measurements (about 0.004% in total). We also delete obviously erroneous values in the year-long time series of each property. Gaps in the mooring data were addressed via linear interpolation whenever possible. The pressure record of the CTD sensor installed at a nominal depth of 262 m on the central mooring contained one distinct downward shift (~30 m) from July to September 2013. These values were corrected to the nominal sensor depth. The top CTD sensor with a nominal depth of 54 m was damaged on the northeast inner mooring, limiting the calculation of buoyancy gradients across the inner cluster above approximately 115-m depth.

b. Additional datasets

In addition to the mooring observations, the OSMOSIS domain was also continuously sampled by at least two (five in total) autonomous underwater gliders for the entire year (Damerell et al. 2016; Thompson et al. 2016; Erickson and Thompson 2018; Evans et al. 2018). The gliders navigated in a bow-tie pattern across the mooring array, measuring temperature and salinity profiles within the top 1000 m of the ocean at approximately 1-m depth intervals. The mixed layer depth $H$ is calculated from the glider data using a threshold value of potential density increase ($\Delta\rho = 0.03 \text{ kg m}^{-3}$) from a near-surface value at 10 m (Damerell et al. 2016).

Air–sea heat and freshwater fluxes and wind stress data are taken from the ECMWF (European Centre for Medium-Range Weather Forecasts) ERA-Interim reanalysis product, with a time interval of 3 h (Dee et al. 2011). Using fields with a horizontal resolution of 0.75°, data are linearly interpolated onto the OSMOSIS central mooring site. The net heat flux is calculated as the sum of shortwave radiation, longwave radiation, latent heat flux, and sensible heat flux.
### c. Definition of seasons

In this work, the seasons are defined as follows: fall (September–November), winter (December–February), spring (March–May), and summer (June–August). The four seasons during the observational period are indicated in Fig. 3a below.

## 3. Potential vorticity framework

### a. Calculation of potential vorticity

The Ertel PV is a useful diagnostic to study the evolution and stability of ocean flows, and under the Boussinesq approximation can be defined as

$$q = \omega_a \cdot \nabla b = (f\mathbf{k} + \nabla \times \mathbf{u}) \cdot \nabla b, \quad (1)$$

where $\omega_a$ is the 3D absolute vorticity, $\nabla$ is the spatial gradient operator, $b = g(1 - \rho/\rho_0)$ is buoyancy (with $g$ as the gravitational acceleration, $\rho$ as potential density, and $\rho_0 = 1025 \text{ kg m}^{-3}$ as a reference density), $f = 2\Omega \sin\phi$ is the Coriolis parameter (with $\Omega$ as Earth's angular velocity and $\phi$ as latitude), $\mathbf{k}$ is the vertical unit vector, and $\mathbf{u} = (u, v, w)$ is the 3D velocity vector.

A range of submesoscale instabilities may arise when the Ertel PV takes the opposite sign to the planetary vorticity (Hoskins 1974; Haine and Marshall 1998), which is positive in the Northern Hemisphere. Negative PV values may occur when the fluid column is unstably stratified (gravitational instability) or experiences horizontally sheared flows (centrifugal instability) or strong anticyclonic along-isopycnal shear (symmetric instability). To identify the causes of the Ertel PV becoming negative, it is useful to decompose $q$ into vertical and horizontal components. The vertical component of PV,

$$q_v = (f + \zeta)N^2, \quad (2)$$

is associated with the vertical component of the absolute vorticity, $f + \zeta$, and the vertical stratification $N^2 = \partial b/\partial z$, where $\zeta = \mathbf{k} \cdot \nabla \times \mathbf{u}$ is the vertical component of the relative vorticity.

The horizontal component of PV, $q_h = [((\partial w/\partial y) - (\partial v/\partial z)](\partial b/\partial x)] + [((\partial u/\partial z) - (\partial w/\partial x)](\partial b/\partial y)]$, is associated with the horizontal component of the absolute vorticity and the horizontal buoyancy gradient. By neglecting terms that include vertical velocity derivatives, a more compact expression for $q_h$ is obtained:

$$q_h = \frac{\partial v}{\partial z} \frac{\partial b}{\partial x} + \frac{\partial u}{\partial z} \frac{\partial b}{\partial y}. \quad (3)$$

This assumption is justified by the numerical study of Brannigan et al. (2017) in a model domain analogous to our mooring area, where they found that vertical velocity derivatives are everywhere negligible. Then the Ertel PV becomes

$$q = q_v + q_h = (f + \zeta)N^2 - \frac{\partial v}{\partial z} \frac{\partial b}{\partial x} + \frac{\partial u}{\partial z} \frac{\partial b}{\partial y}. \quad (4)$$

By further assuming that the flow is in thermal wind balance, $[(\partial u/\partial z), (\partial v/\partial z)] = (1/f)[-(\partial b/\partial y), (\partial b/\partial x)]$, we can derive geostrophic versions of the Ertel PV components, such as

$$q_{h_s} = -f \left| \frac{\partial \mathbf{u}_h}{\partial z} \right|^2 = - \frac{|\nabla_h b|^2}{f}, \quad (5)$$

and

$$q_{\text{geo}} = q_v + q_{h_s} = (f + \zeta)N^2 - f \left| \frac{\partial \mathbf{u}_h}{\partial z} \right|^2 = (f + \zeta)N^2 - \frac{|\nabla_h b|^2}{f}, \quad (6)$$

where the $h$ subscript denotes horizontal component. In the geostrophic expression, $q_{h_s}$ is always negative and reduces $q_{\text{geo}}$, whereas $q_v$ may be either positive or negative. Importantly, SI develops when $q_{h_s}$ overcomes $q_v$ with $q_v > 0$.

We assess the assumption of geostrophy in section 4b. The vertical shear estimated on the central mooring is vertically smoothed over 60 m to match with the vertical resolution of the horizontal buoyancy gradient estimated from inner-cluster measurements. Note that, as the top CTD sensor on the northeast inner mooring was damaged, we choose to estimate lateral buoyancy gradients above 110 m (the depth of the second sensor closest to the ocean surface) using only the three remaining inner moorings.

### b. Potential vorticity flux equation

To quantify the effects of competing processes deepening and shoaling the mixed layer, the flux form of the PV equation is employed (Haynes and McIntyre 1987, 1990; Marshall and Nurser 1992),

$$\frac{\partial q}{\partial t} + \nabla \cdot \mathbf{J} = 0, \quad (7)$$

where $\mathbf{J}$ represents the advective and nonconservative transport of PV in the ocean, defined as

$$\mathbf{J} = q\mathbf{u} + \nabla b \times \mathbf{F} - \omega_a \frac{Db}{Dt}. \quad (8)$$

The $\mathbf{J}$ vector has advective ($q\mathbf{u}$), frictional ($J^F = \nabla b \times \mathbf{F}$), and diabatic [$J^D = -\omega_a (Db/Dt)$] components, where $\mathbf{F}$ is the frictional force and $Db/Dt$ is the Lagrangian rate of change of buoyancy.

The impermeability theorem suggests that PV cannot be fluxed across isopycnal surfaces, but nonconservative processes can inject or extract PV substance through a boundary (Haynes and McIntyre 1987). Here we shall only consider the upper ocean with no contact with topography. In the presence of a horizontal buoyancy gradient (i.e., buoyancy surfaces that are nearly vertical at the upper boundary), $J^F$ will be dominated by its vertical component and may be approximated as a vertical flux. Likewise, $\omega_a$ is typically dominated by its vertical component, so that $J^D$ can be represented by the convergence of a vertical flux (Thomas 2008). Thus, using the continuity equation (i.e., $\nabla \cdot \mathbf{u} = 0$), the flux form of the PV equation reduces to

$$\frac{\partial q}{\partial t} + \mathbf{u} \cdot \nabla q + \frac{\partial}{\partial z} (J^F + J^D) = 0. \quad (9)$$

For a depth-integrated PV budget, the PV equation in the mixed layer becomes
<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Finite-difference configuration used to compute PV at the central mooring site. (b) Illustration of the finite-difference configuration used to estimate PV within a triangle-shaped region surrounded by the central and two outer moorings.</figcaption></figure>

$$ \left\langle \frac{\partial q}{\partial t} \right\rangle_{\text{ML}} + \langle \mathbf{u} \cdot \nabla q \rangle_{\text{ML}} + J_z^D + J_z^F = 0, \quad (10) $$

where $J_z^D$ and $J_z^F$ denote the vertical components of the diabatic and frictional forcings; the angle brackets $\langle \cdot \rangle$ indicate a depth integration; and the subscript “ML” refers to the diagnosed mixed layer. Note that the PV flux across the mixed layer base is neglected.

In this work, the PV equation [Eq. (10)] will be applied on a horizontal scale of $O(10)$ km, using the measurements from the central and four outer moorings. The inner mooring cluster is not used for the PV budget analysis for two reasons. First, as mentioned in section 2a, the top CTD sensor on the northeast inner mooring was damaged, limiting the calculation of buoyancy gradients (e.g., vertical stratification $N^2$) above approximately 115-m depth. Second, the horizontal gradients of PV estimated within the inner cluster are found to be excessively noisy, as could be expected from the computation of second derivatives of horizontal velocity and buoyancy over the reduced horizontal separation between moorings.

We also considered an alternative estimate of $H$ based on a PV criterion (Fig. 3c), that is, $H$ defined as the shallowest depth where $fq > 0$ (Bachman et al. 2017). We found our key findings to be insensitive to this choice. For instance, adopting the $q = 0$ surface as an approximation of $H$ can quantitatively alter the magnitudes of $J_z^D$ and $J_z^F$ by as much as a factor of 2, but the qualitative temporal evolution of these variables remains unchanged. Thus, in this study, we conservatively use the density-threshold definition of the mixed layer depth as an approximation to the surface boundary layer depth, to be consistent with previous OSMOSIS works (Thompson et al. 2016; Buckingham et al. 2019; Erickson et al. 2020). A caveat of our analysis is that the mooring observations do not include the uppermost 50 m, where PV dynamics can be modified by small-scale turbulent processes, including surface wave breaking and Langmuir turbulence (Canuto 2015; Bodner and Fox-Kemper

2020). Errors introduced to our diagnostics by mooring motion and instrumental noise are discussed in appendix A.

1) TEMPORAL CHANGE OF POTENTIAL VORTICITY ($\partial q/\partial t$)

The temporal derivative of PV, $\partial q/\partial t$, is calculated as a second-order centered finite difference in time. To determine the estimates of $q$ at the central mooring site, the vertical derivatives of buoyancy and horizontal velocity are computed as centered finite differences in depth. The relative vorticity at the central mooring site is estimated using the horizontal derivatives of horizontal velocity from the outer cluster (see Fig. 2a). To do so, horizontal coordinates are rotated counterclockwise by an angle of 45° to approximately match the cross shape of the mooring array.

2) HORIZONTAL ADVECTION OF POTENTIAL VORTICITY [$u(\partial q/\partial x) + v(\partial q/\partial y)$]

The horizontal advection of PV, $u(\partial q/\partial x) + v(\partial q/\partial y)$, is derived from estimates of $q_{\text{geo}} = (f + \zeta)N^2 - f\partial\mathbf{u}_h/\partial z|^2$ at four triangle-shaped areas and horizontal velocity measurements at the central mooring (Fig. 2b). It will be shown in section 4b that the subinertial flows are to leading order in thermal wind balance. Estimates of $\zeta$ are obtained from each triangle-shaped area following Stokes’ theorem, $\zeta = \mathbf{k} \cdot \nabla \times \mathbf{u} = (1/A)\int_{\mathbf{u}_h} \mathbf{u}_h \cdot d\mathbf{s}$. The vertical stratification $N^2$ and the vertical shear $\partial\mathbf{u}_h/\partial z$ of each triangle-shaped region are calculated as spatial averages of the respective estimates at the three surrounding moorings.

3) VERTICAL ADVECTION OF POTENTIAL VORTICITY [$w(\partial q/\partial z)$]

The vertical advection of PV, $w(\partial q/\partial z)$, is computed separately within the mixed layer and the ocean interior. In the mixed layer, direct estimates of $w$ are not available, so the vertical advection of PV is approximated as
$$w \frac{\partial q}{\partial z} \approx \frac{\partial(wq_v)}{\partial z} - q_v \frac{\partial w}{\partial z}, \quad (11)$$

where the vertical gradient of the horizontal component of PV has been neglected. We can eliminate $w$ from the first term on the right hand side of (11) by using the buoyancy equation,

$$\frac{\partial b}{\partial t} + u \frac{\partial b}{\partial x} + v \frac{\partial b}{\partial y} + w \frac{\partial b}{\partial z} = D, \quad (12)$$

and assuming that, to a good approximation, $D$ is small (i.e., that local changes in buoyancy are generally driven by horizontal advection). We can use the continuity equation,

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0, \quad (13)$$

to eliminate $w$ from the second term on the right hand side of (11). We thus obtain the following expression for the vertical advection of PV within the mixed layer:

$$w \frac{\partial q}{\partial z} \approx \partial \left[ \left( -\frac{\partial b}{\partial t} - u \frac{\partial b}{\partial x} - v \frac{\partial b}{\partial y} \right) (f + \xi) \right] / \partial z - \left[ (f + \xi) N^2 \right] \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right). \quad (14)$$

Below the mixed layer, $w(\partial q/\partial z)$ can be calculated directly from the vertical velocity diagnosed by Yu et al. (2019a) and $q$ computed at the central mooring site.

### 4) DIABATIC POTENTIAL VORTICITY FLUX ($J_z^D$)

Following Marshall and Nurser (1992), the diabatic component of the **J** vector at the ocean surface (hereinafter the diabatic PV flux) is approximated by

$$J_z^D = \frac{f}{H} B_0, \quad (15)$$

where $B_0 = -[g\alpha Q_{net}/\rho_0 c_p + g\beta S(P - E)]$ is the surface buoyancy flux, $Q_{net}$ is the air-sea heat flux, $c_p$ is the specific heat capacity, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients, $S$ is the surface salinity, $(E - P)$ is the freshwater flux, and $H$ is the mixed layer depth. We adopt the salinity measurements at the top CTD sensor on the central mooring as an approximation to the surface salinity. Note that $J_z^D$ is positive for surface buoyancy loss (e.g., surface cooling or salinification), and negative for surface buoyancy gain (e.g., surface heating or freshening).

### 5) FRICTIONAL POTENTIAL VORTICITY FLUX ($J_z^F$)

Following Thomas (2005), the frictional component of the **J** vector (hereinafter the frictional PV flux) at the ocean surface is approximated by

$$J_z^F = \frac{f}{H} B_e, \quad (16)$$

where $B_e = (\tau \times \mathbf{k}) \cdot \nabla_h b / \rho_0 f$ is the Ekman buoyancy flux and $\tau$ is the surface wind stress. Note that $J_z^F$ is positive if the wind stress has a downfront component (i.e., directed with the geostrophic shear), and negative if the wind stress has

an upfront component (i.e., directed against the geostrophic shear).

### c. Convective layer depth ($h$)

Following Bachman et al. (2017), a quartic equation is used to solve for the theoretical prediction of the convective layer depth $h$,

$$\left(\frac{h}{H}\right)^4 - c^3 \left(1 - \frac{h}{H}\right)^3 \left[ \frac{w_*^3}{|\Delta \mathbf{u}_g|^3} + \frac{u_*^2}{|\Delta \mathbf{u}_g|^2} \cos \theta_w \right]^2 = 0, \quad (17)$$

where $w_* = (B_0 H)^{1/3}$ is the convective velocity, $u_* = \sqrt{|\tau|/\rho_0}$ is the frictional velocity, $\mathbf{u}_g$ is the geostrophic velocity vector, $\theta_w$ is the angle between the wind vector and the geostrophic shear, and $c = 14$ is an empirical constant. Here, we adopt the mixed layer depth $H$ as an approximation to the surface boundary layer depth. The change in geostrophic velocity across the mixed layer, $\Delta \mathbf{u}_g$, is computed by assuming geostrophic balance as $\Delta \mathbf{u}_g = |\partial \mathbf{u}_g/\partial z| H = |\nabla_h b| H/f$. Note that the convective layer depth $h$ solved for here is always no greater than the mixed layer depth $H$.

### d. Categorization of instability types

To identify the types of submesoscale frontal instabilities that may potentially develop for the measured PV conditions, the balanced Richardson number angle, $\phi_{RiB} = \tan^{-1}(-f^2 N^2 / |\nabla_h b|^2)$, and the critical angle, $\phi_C = \tan^{-1}(-\xi/f)$, are computed following Thomas et al. (2013). This approach has been used by previous studies to assess the susceptibility of the flow to submesoscale instabilities (e.g., Thompson et al. 2016; Naveira Garabato et al. 2017; Ramachandran et al. 2018; Viglione et al. 2018; Naveira Garabato et al. 2019). The instability criteria for gravitational, symmetric, and centrifugal instabilities may be synthesized as follows:

(i) For unstable vertical stratification (i.e., $N^2 < 0$), gravitational instability is expected to develop when $-180^\circ < \phi_{RiB} < -135^\circ$, and hybrid gravitational/symmetric instabilities will occur when $-135^\circ < \phi_{RiB} < -90^\circ$.

(ii) For stable stratification (i.e., $N^2 > 0$) and cyclonic vertical vorticity, SI is predicted to develop when $-90^\circ < \phi_{RiB} < \phi_C$, with $\phi_C < -45^\circ$.

(iii) For stable stratification (i.e., $N^2 > 0$) and anticyclonic vertical vorticity, SI is expected for $-90^\circ < \phi_{RiB} < -45^\circ$ with $\phi_C > -45^\circ$, and hybrid symmetric/centrifugal instabilities will occur when $-45^\circ < \phi_{RiB} < \phi_C$.

Note that this analysis does not take into account submesoscale baroclinic instability (which can arise when $fq > 0$) or the modification of stratification by surface waves (Hamlington et al. 2014).

## 4. Evolution of upper-ocean potential vorticity

### a. Annual cycle of upper-ocean potential vorticity

The surface heat flux (which overwhelmingly dominates the surface buoyancy flux) and upper-ocean potential density at the central mooring site from September 2012 to September
<figure><img src="image_5.png" /><figcaption>FIG. 3. Year-long time series of (a) surface heat flux $Q_{net}$, (b) subinertial potential density $\rho$, and (c) PV $q$ at the central mooring site. The black lines in (b) are isopycnal contours at $0.05 \text{ kg m}^{-3}$ intervals, and the white line is the $27.1 \text{ kg m}^{-3}$ isopycnal. The four seasons (fall, winter, spring, and summer) are indicated by colored horizontal lines along the $x$ axis. The mixed layer depth and the PV = 0 contour are overlaid in (c) as black and magenta lines, respectively. Depths not sampled by moored instrumentation in (b) and (c) are colored in gray. The winter-to-mid-spring period is indicated in between two yellow dashed lines.</figcaption></figure>

2013 are displayed in Figs. 3a and 3b. The surface heat flux swings from sizeable surface cooling periods through fall, winter and most of spring, to moderate surface heating during late spring and summer. Intense cooling exceeds $-400 \text{ W m}^{-2}$ in winter, and the year-mean surface heat flux reaches about $-45 \text{ W m}^{-2}$, indicating a net destabilization of the mixed layer by surface buoyancy forcing. High-frequency variability, such as a diurnal cycle, is not present in the frequency spectrum of surface heat flux (not shown). The annual cycle of upper-ocean potential density is mainly driven by the local surface heat flux. In summer and fall, the upper ocean is strongly stratified, with a sharp pycnocline mostly above 100 m. Slowly evolving deep baroclinic eddies are evident below the shallow pycnocline, with steep potential density contours down to 500 m (e.g., end of November or 6–11 August). Throughout winter and spring, the most striking feature is the absence of the shallow pycnocline under persistent surface cooling, implying that that shallow pycnocline is a seasonal feature. This is endorsed by the glider observations which, with higher vertical resolution and a wider sampling range (0–1000 m), reveal the occurrence of a permanent pycnocline below 600 m throughout the year (Thompson et al. 2016).

Using Eq. (5) and inner cluster measurements, we document the year-long time series of upper-ocean Ertel PV computed on horizontal scales of $O(1)$ km (Fig. 3c).

PV exhibits substantial seasonality, typically in synchrony with the seasonal evolution of the mixed layer depth. Instances of negative PV are frequently observed within the deep mixed layers in winter and early-to-mid-spring (December 2012-late April 2013). For the remainder of the year (September-December 2012 and May-September 2013), the seasonal pycnocline is manifested in the significant enhancement of positive PV at the mixed layer base. PV is consistently positive below the mixed layer, and the maximum PV value in the seasonal pycnocline reaches $10^{-8} \text{ s}^{-3}$, exceeding typical values in the ocean interior by one order of magnitude.

Further to this seasonality, PV and mixed layer depth both display abrupt seasonal transitions. The mixed layer depth is about 20 m at the beginning of the record (September 2012), and gradually deepens to 100 m through fall. A strong convective event, caused by destabilizing surface forcing at the fall-to-winter transition (around 25 November–5 December), reduces mixed layer PV to negative values and deepens the mixed layer by over 50 m within just a few days. The transition from spring to summer occurs in mid-April, as evidenced by the rapid shoaling of the mixed layer in response to a reversal in surface forcing from cooling to heating. The mixed layer depth continues to exhibit substantial variability until July, after which the mixed layer remains shallow at about 20 m. Waters with low (but positive) PV and prominent isopycnal
<figure><img src="image_6.png" /><figcaption>FIG. 4. Cumulative time integrals of the temporal change of mixed layer depth $\int_{0}^{t}(\partial H/\partial t) dt$ in yellow, $\int_{0}^{t}(\partial H/\partial t)^{+} dt$ (i.e., mixed layer deepening) in orange, and $\int_{0}^{t}(\partial H/\partial t)^{-} dt$ (i.e., mixed layer shoaling) in blue.</figcaption></figure>

displacements are present around 6 August 2013, associated with the propagation of an anticyclonic eddy across the mooring array.

The mooring array area is also influenced by high-frequency processes. Frequency spectra of horizontal velocity and potential density at the central mooring display high-energy peak at the $M_2$ semidiurnal tidal frequency and at the inertial frequency (Yu et al. 2019a). Erickson et al. (2020) further show that the superinertial range closely follows the Garrett and Munk (1975) spectrum for internal waves. However, only near-inertial signals are detected from the glider-based mixed layer depth (not shown). To quantify the competition between processes deepening and shoaling the mixed layer, we calculate $\int_{0}^{t}(\partial H/\partial t)^{+} dt$ and $\int_{0}^{t}(\partial H/\partial t)^{-} dt$, where superscripts “+” and “-” respectively represent positive (i.e., deepening) and negative (i.e., shoaling) changes in the mixed layer depth. Cumulative integrals of the positive and negative subinertial changes in $H$ (Fig. 4) indicate that the mixed layer base undergoes frequent and large vertical excursions in winter to mid-spring, but is relatively invariable during the rest of the year. This suggests that the winter-to-mid-spring period hosts a persistent and vigorous competition between destabilizing atmospheric forcing and restratifying upper-ocean processes, such as submesoscale frontal instabilities. Integrated over the full annual cycle, destabilizing forcing and restratifying processes respectively account for over 8000 m of mixed layer deepening and shoaling, with the seasonal cycle of the mixed layer depth (peaking at about 350 m) arising as a small residual between the two. The accumulated magnitude of the mixed layer deepening or shoaling is increased by a factor of 2.2 when superinertial variability is considered, although the patterns in Fig. 4 are insensitive to the inclusion of this variability. Note that this picture of competing processes is not an artifact of horizontal advection of an inclined mixed layer base past the mooring array. Estimating this advection as $u(\partial H/\partial x) + v(\partial H/\partial y)$ (where $H$ is calculated at each outer mooring using the same density threshold method as previously used in the glider

data, and $u$ and $v$ are taken from the central mooring) reveals that horizontal advection accounts for only a modest fraction of the local variability in the mixed layer depth in the winter-to-mid-spring period (not shown). Similarly, examining the vertical velocity at the mixed layer base [$w_{-H}$, quantified as in Yu et al. (2019a)], which is associated with the divergence of horizontal flow within the mixed layer [i.e., $\int_{-H}^{0} \nabla_h \cdot \mathbf{u}_h dz = \int_{-H}^{0} -(\partial w/\partial z) dz = w_{-H}$], indicates that this effect contributes unimportantly to the rate of change of the mixed layer depth, e.g., the root-mean-square value of $w_{-H}$ is approximately 30% of that of $\partial H/\partial t$ in winter to mid-spring.

Overall, the OSMOSIS moorings capture the bulk of the mixed layer from winter to mid-spring, the shallow seasonal pycnocline for the remainder of the year, and the weakly stratified layer just above the main pycnocline all year round. Our primary focus in this work is on the subinertial frontal processes that are associated with negative PV in the mixed layer in the winter-to-mid-spring period (Fig. 3c). While negative PV is often observed when PV destruction is effected by surface cooling, a direct link between negative PV and destabilizing surface heat fluxes is not necessarily expected. This is because locally generated low-PV waters might be advected to surrounding areas, or low-PV waters generated elsewhere might be advected into the upper-ocean volume sampled by the moorings. By constructing a PV budget of this volume, we will account for advective PV transport, and thus more precisely identify the impacts of submesoscale frontal instabilities on upper-ocean stratification and mixed layer depth.

b. Vertical and horizontal components of potential vorticity

To gain deeper insight into the dynamical evolution of upper-ocean Ertel PV, its vertical and horizontal components (estimated from the inner mooring cluster) are examined next. Within the mixed layer, PV is largely determined by its vertical component (Figs. 5a,b), but the horizontal component becomes comparable in magnitude to the vertical component
<figure><img src="image_7.png" /><figcaption>FIG. 5. Time series of mixed layer–averaged (a) PV, (b) its vertical component, (c) its horizontal component, and (d) the geostrophic approximation to the horizontal component. Gray scatter points represent all estimates in the mixed layer. (e) Cumulative distribution functions of PV (black), its vertical component (blue), its horizontal component (orange), and the geostrophic approximation to the horizontal component (green).</figcaption></figure>

during relatively short-lived, intermittent frontal events asso-
ciated with weak stratification. Geostrophic shear dominates
over ageostrophic shear during the entire winter-to-mid-spring
period (Figs. 5c,d). The ageostrophic shear acts to slightly
weaken the horizontal component of PV, and thus tends to
systematically stabilize the upper ocean in the presence of
strong lateral buoyancy gradients.

The overwhelmingly geostrophic character of the sub-
inertial velocity field in the OSMOSIS area can be explicitly
illustrated by considering the vertical structure of the cur-
rents (Fig. 6). The winter-to-mid-spring period hosts in-
stances of strong flow, exceeding 0.5 m s⁻¹ near the surface
and decreasing gently with depth down to 500 m. The vertical

shear in this flow exhibits a good agreement with the lateral
buoyancy gradient estimated from the inner mooring cluster,
with a correlation coefficient of 0.75 (Yu et al. 2019a). This
indicates that subinertial currents with horizontal scales of
*O*(1) km are in thermal wind balance to leading order, in line
with the frequency-resolved horizontal structure function
results of Callies et al. (2020).

In quantitative terms, the cumulative distribution functions
of PV and its components during the winter-to-mid-spring
period show that PV in the mooring-observed mixed layer is
negative approximately 34% of the time (Fig. 5e). About
70% of these negative PV events result from a negative ver-
tical component, and the remaining 30% from the horizontal
<figure><img src="image_8.png" /><figcaption>FIG. 6. Time series of subinertial horizontal velocities, vertical shear terms, and horizontal buoyancy gradient terms in the thermal wind balance equation, at the central mooring site in the winter-to-mid-spring period: (a) $u$, (b) $v$, (c) $\partial u/\partial z$, (d) $\partial v/\partial z$, (e) $-(1/f)(\partial b/\partial y)$, and (f) $(1/f)(\partial b/\partial x)$. The mixed layer depth is superimposed as a black line in all panels. Depths not sampled by the moored instrumentation are colored in gray.</figcaption></figure>

component. This indicates that at least 30% of the negative PV events exhibit conditions favorable to the growth of SI, in accord with the analysis of OSMOSIS glider measurements (Thompson et al. 2016). Recall, however, that the uppermost 50 m is excluded from the mooring observations, where negative PV is expected to occur most often under destabilizing surface forcing. Last, we compare the PV estimates with and without thermal wind balance (i.e., $q_{geo}$ versus $q$), and find that assuming geostrophy increases the occurrence of negative PV by 1% (not shown). This suggests that the inhibition of

submesoscale frontal instabilities by the ageostrophic shear is insignificant in our study region.

Below the mixed layer, PV is almost exclusively determined by its vertical component, which is typically one order of magnitude larger than the horizontal component (Fig. 7). In winter and spring, PV is reduced to small values well below the mixed layer, as a result of the embedding of that part of the water column (down to 500 m) within the weakly stratified subpolar mode water that overlies the main thermocline (Thompson et al. 2016; Callies et al. 2020). PV estimated with

<figure><img src="image_9.png" /><figcaption>FIG. 7. Year-long time series of PV (black), and its vertical (blue) and horizontal (orange) components estimated from the inner mooring cluster. The geostrophic approximation to the horizontal component of PV is shown in green. All variables are from below the mixed layer and depth averaged.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Probability histogram of the occurrence of submesoscale instabilities (GI, gravitational instability; SI, symmetric instability; CI, centrifugal instability) in the convective and forced-SI layers, calculated using (a) the density-threshold mixed layer depth (black line in Fig. 3c) and (b) the $q = 0$ surface (magenta line in Fig. 3c) as the base of the ocean surface boundary layer.</figcaption></figure>

the thermal wind balance assumption is barely distinguishable from PV computed with inclusion of the ageostrophic shear (not shown).

### c. Submesoscale instabilities in the mixed layer

We may now draw on the estimates of upper-ocean PV and its components to classify the submesoscale frontal instabilities predicted to develop during instances of negative PV in the winter-to-mid-spring period, following the balanced Richardson angle $\phi_{RIB}$ formulation in section 3d. The point-wise instability diagnostics are synthesized in histogram form in Fig. 8, where results are grouped into a convective layer (defined in section 3c) and a forced-SI layer (extending from the lower boundary of the convective layer to the mixed layer base), respectively based on the glider-derived density-threshold mixed layer depth and on the mooring-derived $q = 0$ surface (Fig. 3c). In both cases, conditions favorable to the growth of gravitational instability are generally found within the convective layer, and become rare below it. In turn, instances of stable stratification and SI-favorable conditions are often observed in the forced-SI layer, but much less frequently in the convective layer. A limitation of this set of diagnostics is that they may not be representative of the entire mixed layer, as gravitational instability is more likely to occur than SI in the uppermost 50 m not sampled by the moorings.

In contrast to western boundary current regions, which exhibit relatively persistent frontal and wind forcing patterns, the OSMOSIS area is characterized by ephemeral and intermittent SI events, likely triggered by the occasional alignment of the winds with transient upper-ocean fronts generated by mesoscale frontogenesis. Yu et al. (2019b) reported unambiguous evidence of a SI event forced by downfront winds at one such transient front in early April 2013, using the OSMOSIS mooring and glider data, and found the instability to be associated with elevated upper-ocean kinetic energy, rapid restratification and intensified turbulent dissipation. Our results suggest that this

type of event was a relatively frequent occurrence in our study area during the winter-to-mid-spring period.

## 5. Upper-ocean potential vorticity budget

As shown in the previous section, the upper ocean exhibits numerous instances of mixed layer deepening and shoaling, the small residual of which yields the seasonal cycle in mixed layer depth. Episodes of high variability in the mixed layer depth are most obvious in the winter-to-mid-spring period, and often occur in association with negative PV events. Next, we quantify the competing processes inducing deepening and shoaling of the mixed layer during the entire year of OSMOSIS observations, through application of the PV budget framework introduced in section 3b.

### a. Annual cycle of surface potential vorticity fluxes

We commence by examining the surface forcing of PV in the mooring array area. This forcing is synthesized in Fig. 9, which displays year-long time series of the diabatic ($J_z^D$) and frictional ($J_z^F$) PV fluxes through the ocean surface. The diabatic PV flux reverses from a sizeable PV extraction ($J_z^D > 0$) in winter to a very strong PV injection ($J_z^D < 0$) in summer, with respective peak amplitudes of $1 \times 10^{-12}$ and $-3 \times 10^{-12}$ m s$^{-4}$. The diabatic PV flux $J_z^D$ is larger in magnitude in summer than in winter because summer heating confines the diabatic flux to a much shallower depth range (e.g., down to 20 m between June and September in 2013) than winter cooling, which operates on a much deeper mixed layer. The strong seasonality of $J_z^D$ is mainly determined by the air-sea heat flux, with the freshwater flux contribution being more intermittent and smaller in magnitude than the air-sea heat flux.

By contrast, the frictional PV flux is highly variable throughout the year, due to the episodicity of wind stress magnitudes and directions and the transient nature of upper-ocean
<figure><img src="image_11.png" /><figcaption>FIG. 9. Year-long time series of (a) surface diabatic PV flux $J_z^D$ and (b) surface frictional PV flux $J_z^F$ calculated from outer (orange) and inner (gray) mooring clusters. The gray dashed line in (a) shows the contribution from the air-sea heat flux, for reference. Note that PV is extracted from the ocean when $J_z^D + J_z^F > 0$.</figcaption></figure>

fronts in the area. Consistent estimates of $J_z^F$ are obtained from inner and outer mooring cluster measurements, although the inner cluster-based calculation produces larger $J_z^F$ magnitudes as a result of that cluster's ability to sample lateral buoyancy gradients at higher horizontal resolution. This is especially true in the winter-to-mid-spring period, when submesoscale flows are most active. The $J_z^F$ values are in the range from $-6 \times 10^{-13}$ to $6 \times 10^{-13}$ m s$^{-4}$, and are thus substantially smaller than $J_z^D$ estimates except in winter to mid-spring, when the magnitudes of $J_z^F$ and $J_z^D$ are comparable.

The annual-mean surface PV flux ($J_z^D + J_z^F$) in the OSMOSIS region is overwhelmingly dominated by diabatic processes, as found by previous numerical studies in the North Atlantic (e.g., Maze and Marshall 2011). The annual-mean $J_z^D$ is found to be directed into the ocean, with a value of $-2.46 \times 10^{-13}$ m s$^{-4}$, indicating that diabatic PV injection more than compensates for diabatic PV extraction over the year as a whole. In turn, $J_z^F$ averages to near-zero values over any period of several months and longer. The annual-mean $J_z^F$ is thus $-1.37 \times 10^{-15}$ or $-1.52 \times 10^{-15}$ m s$^{-4}$ (for inner and outer cluster measurements, respectively), i.e., two orders of magnitude smaller than the annual-mean $J_z^D$.

In summary, the surface forcing of PV in the mooring area is dominated by diabatic fluxes on time scales of months to a year, but frictional fluxes can contribute significantly on shorter time scales in the winter-to-mid-spring period. We now examine the upper ocean's response to this forcing in winter to mid-spring in two stages: in the mixed layer (section 5b) and below the mixed layer (section 5c).

## b. Potential vorticity budget in the mixed layer

The daily-averaged local temporal change of PV ($\langle \partial q/\partial t \rangle_{\text{ML}}$) and total (horizontal + vertical) advection of PV ($\langle \mathbf{u} \cdot \nabla q \rangle_{\text{ML}}$), integrated between the base of the mixed layer and 50 m, exhibit a statistically significant anticorrelation ($R = -0.64$, $p < 0.001$; Fig. 10a). This indicates that advective transports

of PV play an important role in determining the local PV variability within the mixed layer and, thus, that advection must be accounted for in order to unravel the upper ocean's response to surface forcing. Combining the local temporal change of PV $\langle \partial q/\partial t \rangle_{\text{ML}}$ and total advection of PV $\langle \mathbf{u} \cdot \nabla q \rangle_{\text{ML}}$ in the material change of PV, $\langle Dq/Dt \rangle_{\text{ML}}$, reveals that the PV of mixed layer water parcels is highly variable in time, and is modified by processes operating on time scales of a few days (Fig. 10b). The material change of PV is similar in pattern and magnitude, but of opposite sign, to the surface PV flux, $J_z^D + J_z^F$ (Fig. 10c). This provides an observational demonstration of the relationship between PV modification and surface forcing expressed by the PV equation (section 3b).

The closure of the mixed layer PV budget is examined in detail in Fig. 11, which displays cumulative time integrals of the material change of mixed layer PV and of the surface forcing terms. The diabatic PV flux stands out as the dominant surface forcing, and acts to extract PV from the ocean in winter to mid-spring. By contrast, the overall effect of the frictional PV flux during this period is to inject PV into the ocean, although this contribution is an order of magnitude smaller than that of the diabatic PV flux. The reason for this modest role of frictional forcing is that the impacts of upfront and downfront wind events tend to average out on time scales of months. The temporal change of PV $\langle \partial q/\partial t \rangle_{\text{ML}}$, the horizontal advection term $\langle u(\partial q/\partial x) + v(\partial q/\partial y) \rangle_{\text{ML}}$ and the vertical advection term $\langle w(\partial q/\partial t) \rangle_{\text{ML}}$ all act to oppose the effect of the diabatic forcing, with $\langle u(\partial q/\partial x) + v(\partial q/\partial y) \rangle_{\text{ML}}$ being the largest term. The evolution of the material change of mixed layer PV $-\langle Dq/Dt \rangle_{\text{ML}}$ follows that of the total surface PV flux, $J_z^D + J_z^F$, but is approximately a factor of 1.5 smaller in magnitude. This is likely a consequence of the moorings' failure to sample the top 50 m of the water column, which make up a substantial fraction of the mixed layer volume. To assess this interpretation, we re-estimate the material change of mixed layer PV by assuming that $Dq/Dt$ within $-50 \text{ m} < z < 0$ is equal to the measured
<figure><img src="image_12.png" /><figcaption>FIG. 10. (a) Scatterplot of vertically integrated (all measurements in the mixed layer) $\langle \partial q / \partial t \rangle_{\text{ML}}$ vs $\langle \mathbf{u} \cdot \nabla q \rangle_{\text{ML}}$, with the best-fit regression line in black. Time series of (b) $\langle Dq / Dt \rangle_{\text{ML}}$ and (c) $J_z^D + J_z^F$ during the winter-to-mid-spring period.</figcaption></figure>

$Dq/Dt$ at $z = -50$ m. This is, of course, only a reasonable assumption when the top 50-m layer hosts the same dynamical processes. The revised estimate of the material change of mixed layer PV is much closer to the total surface PV flux, integrating to 89% of this flux over the course of the winter-to-mid-spring period. That this revised estimate comes a little short of the integrated total surface PV flux is consistent with

the expectation that the top 50 m, being in direct contact with the atmosphere, host more pronounced changes in PV than underlying waters. Note that our revised estimate of the material change of mixed layer PV may also be affected by its (partial or complete) omission of the Ekman transport of PV, which will be largest near the ocean surface. However, an estimate of this transport based on reanalysis wind stress data

<figure><img src="image_13.png" /><figcaption>FIG. 11. Cumulative time integrals of terms in the mixed layer PV budget during the winter-to-mid-spring period. (a) The surface diabatic PV flux $J_z^D$ and surface frictional PV flux $J_z^F$ are indicated by the blue and orange lines, respectively. The temporal change of PV $\langle \partial q / \partial t \rangle_{\text{ML}}$, horizontal advection of PV $\langle u(\partial q / \partial x) + v(\partial q / \partial y) \rangle_{\text{ML}}$, and vertical advection of PV $\langle w(\partial q / \partial t) \rangle_{\text{ML}}$ integrated over the observed mixed layer are shown by the magenta, yellow, and green lines, respectively. (b) The sum of the surface diabatic and frictional PV fluxes ($J_z^D + J_z^F$) is indicated by the black line. The material change of PV $-\langle Dq / Dt \rangle_{\text{ML}}$ integrated over the observed mixed layer is shown by the red line, and an estimate of the material change of PV with the top 50 m included $-\langle (Dq / Dt)_{\text{ML}} + (Dq / Dt)_{50\text{m}} \rangle$ is shown by the red dashed line. The Ekman transport of PV $\langle u_{\text{Ek}}(\partial q / \partial x) + v_{\text{Ek}}(\partial q / \partial y) \rangle$ is indicated by the purple line. The shaded regions illustrate the 95% confidence envelope of cumulative time integrals, estimated using a Monte Carlo approach.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIG. 12. (a) Scatterplot of the vertically integrated (all measurements in the convective layer) $\langle \partial q/\partial t \rangle_{\text{CL}}$ vs $(\mathbf{u} \cdot \nabla q)_{\text{CL}}$, with the best-fit regression line in black. (b) Scatterplot of $J_z^F$ vs $J_z^D$ colored by $\langle Dq/Dt \rangle_{\text{CL}}$. (c) As in (a), but for the forced-SI layer. (d) As in (b), but for the forced-SI layer.</figcaption></figure>

and lateral buoyancy gradients measured at 50 m suggests that its contribution to the mixed layer PV budget is trivial (appendix B).

We next test the theoretically generated hypothesis that, in the presence of upper-ocean fronts, the mixed layer's dynamical response to surface forcing is distinct in a near-surface convective layer and in a deeper forced-SI layer (Taylor and Ferrari 2010). The mixed layer is notionally subdivided into convective and forced-SI layers following the approach outlined in section 3c: the convective layer (denoted by the subscript “CL”) spans $0 > z > -h$, and the forced-SI layer (denoted by the subscript “SI”) extends across $-H < z < -h$. Inspection of the PV budget terms in each layer reveals significant differences that are in accord with theoretical expectations (Fig. 12). In the convective layer, the anticorrelation between the local temporal change of PV and the total advection of PV ($R = -0.52, p < 0.001$; Fig. 12a) is weaker than for the mixed layer as a whole (Fig. 10a), though still significant. This is consistent with diabatic and frictional processes having a more prominent imprint on PV evolution within the convective layer than in underlying waters. Figure 11b, illustrating the relationship of $\langle Dq/Dt \rangle_{\text{CL}}$ to surface diabatic and frictional forcings, endorses this interpretation. For example, for downfront wind ($J_z^F > 0$) or surface cooling ($J_z^D > 0$) conditions, PV in the convective layer is generally found to

decrease (i.e., $\langle Dq/Dt \rangle_{\text{CL}} < 0$; red shading). Conversely, for upfront wind ($J_z^F < 0$) or surface heating ($J_z^D > 0$) conditions, PV in the convective layer generally increases (i.e., $\langle Dq/Dt \rangle_{\text{CL}} > 0$; blue shading). Further, $\langle Dq/Dt \rangle_{\text{CL}}$ is regularly small (white shading) in the transition between destabilizing and restratifying forcing conditions. Our analysis thus indicates that diabatic and frictional processes play a leading-order role in PV modification within the convective layer.

A very different regime is diagnosed in the forced-SI layer. There, the anticorrelation between the local temporal change of PV and the total advection of PV ($R = -0.72, p < 0.001$; Fig. 12c) is greater than for the convective layer (Fig. 12a) or the mixed layer as a whole (Fig. 10a), suggesting that advective processes exert the primary control on local PV variability within the forced-SI layer. The material change of PV is generally small within this layer (Fig. 12d), and exhibits no discernible relation to surface diabatic and frictional forcings. Thus, our analysis suggests that PV is largely conserved following the flow within the forced-SI layer.

### c. Potential vorticity budget below the mixed layer

Below the mixed layer down to 500 m (denoted by the subscript “int”), the depth-integrated local temporal change of PV $\langle \partial q/\partial t \rangle_{\text{int}}$ and the total advection of PV $(\mathbf{u} \cdot \nabla q)_{\text{int}}$ are significantly anticorrelated, with a correlation ($R = -0.81, p < 0.001$)
<figure><img src="image_15.png" /><figcaption>FIG. 13. (a) Scatterplot of vertically integrated (from 500 m to the base of the mixed layer) $\langle \partial q/\partial t \rangle_{\text{int}}$ vs $\langle u \cdot \nabla q \rangle_{\text{int}}$, with the best-fit regression line in black. (b) Magnitude-squared coherence and (c) cross-spectrum phase between $\langle \partial q/\partial t \rangle_{\text{int}}$ and $-\langle u \cdot \nabla q \rangle_{\text{int}}$.</figcaption></figure>

that exceeds that in the mixed layer (Fig. 13a). The diagnosed balance between $\langle \partial q/\partial t \rangle_{\text{int}}$ and $\langle u \cdot \nabla q \rangle_{\text{int}}$, characterized by a best-fit linear regression with a slope -0.93, indicates that the PV below the mixed layer is approximately conserved following the flow. This observational result is in agreement with the expectation that PV be conserved for adiabatic, frictionless motion below the turbulent surface boundary layer. PV conservation is manifested in high coherence between $\langle \partial q/\partial t \rangle_{\text{int}}$ and $\langle u \cdot \nabla q \rangle_{\text{int}}$ on time scales of 5–10 days (Fig. 13b), and a near-zero lag in the phase of the cross-spectrum between both variables (Fig. 13c). Horizontal advection dominates over vertical advection in balancing $\langle \partial q/\partial t \rangle_{\text{int}}$ over the frequency range of high coherence (not shown).

## 6. Conclusions

In this work, we have diagnosed the first—to our knowledge—observation-based budget of PV in the upper ocean over a full annual cycle, in order to assess the processes governing midlatitude mixed layer evolution. Our results portray a picture of the mixed layer as a highly variable and dynamically active boundary layer experiencing large changes in depth and stratification on time scales as short as days, and whose seasonal cycle arises as a small residual between much larger contributions from destratifying and restratifying processes. Further, our analysis provides quantitative observational verification for a range of theoretical predictions on the nature of the mixed layer response to surface forcing in the presence of upper-ocean fronts. Our main conclusions are summarized as follows:

1) Surface diabatic and frictional PV fluxes both contribute significantly to driving the evolution of PV in the mixed layer in the study area. The diabatic PV flux dominates on time scales of months, and underpins the local seasonal cycle in mixed layer depth and stratification (specifically, a

persistent mixed layer deepening in the fall–winter transition and an abrupt mixed layer shoaling in the spring–summer transition). However, the frictional PV flux has a substantial impact on shorter time scales of days, and can induce abrupt changes in mixed layer depth (both shoaling and deepening) during the winter-to-mid-spring period.

2) The local rate of mixed layer PV destruction from winter to mid-spring, which is associated with mixed layer deepening, is driven strongly by surface buoyancy loss but also modulated significantly by advective processes.

3) In spite of persistent atmospheric cooling in winter-to-mid-spring period, at least 30% of the negative PV events in the mixed layer show conditions conducive to the growth of SI, and those events are typically associated with rapid restratification of the mixed layer as shown in Yu et al. (2019b). This highlights the key role of submesoscale restratifying instabilities in shoaling the mixed layer in our study region.

4) The mixed layer may be conceptualized as a two-layer system, consisting of a near-surface convective layer in which gravitational instability dominates and an underlying forced-SI layer in which SI is more prevalent. This is in accord with theoretical predictions (e.g., Taylor and Ferrari 2010).

5) Below the mixed layer, PV is approximately conserved following the flow, consistent with the theoretical expectation that PV is a conservative tracer in the ocean interior (McIntyre 2015).

Our finding of the important role of submesoscale frontal instabilities in upper-ocean restratification echoes the results of Thompson et al. (2016), who characterized such instabilities from the hydrographic data collected by the OSMOSIS gliders. In our work, the direct velocity measurements provided by the moorings have enabled us to quantitatively document the instabilities' effects, exerted via the advective redistribution of
PV, in the context of the local PV budget. We have also shown that this advection is associated with geostrophic, subinertial flows and, often, with conditions of elevated vertical shear of horizontal velocity and weak vertical stratification that are conducive to the development of symmetric instability. All in all, both glider and mooring observations consistently highlight the significance of submesoscale frontal instabilities in shaping the seasonal evolution of the mixed layer in a typical midocean region.

*Acknowledgments.* The OSMOSIS experiment was funded by the U.K. Natural Environment Research Council (NERC) through Grants NE/I019999/1, NE/I019905/1, NE/I019921/1, and NE/I01993X/1. We are grateful to the engineers, scientists, captain, and crew of the RRS *Discovery*, RRS *James Cook*, and R/V *Celtic Explorer*, who participated in the deployment and recovery of the moorings and gliders. All data are archived at the British Oceanographic Data Centre (glider data are accessible at https://doi.org/10/cqc6). We thank two anonymous reviewers for their constructive comments on this manuscript. We acknowledge helpful discussions with Christian E. Buckingham, Liam Brannigan, and Florian Sévellec. XY is supported by the French National Agency for Research under Grant 17-CE01-0006-01. ACNG acknowledges the support of the Royal Society and the Wolfson Foundation.

# APPENDIX A

## Observational Uncertainty in the PV Budget

We quantify the observational uncertainties in the PV budget introduced by two sources of error: mooring motion and instrumental errors in the velocity and buoyancy measurements.

Because PV estimates involve spatial gradients computed from mooring measurements, unknown mooring motions could result in uncertainty in the locations of the measurements. Stochastic modeling predicts that horizontal displacements rarely exceed 500 m (Buckingham et al. 2016). Following this work, we model the intermooring distance perturbations associated with mooring motion as a Gaussian white noise process with zero mean and nonzero variance, estimated from the time integration of differential horizontal currents.

Instrumental error arises from the accuracy of the moored sensors, and is unavoidable during the measuring process. According to manufacturer specifications, Seabird MicroCAT CTDs have precisions of 0.001°C for temperature and 0.0003 S m⁻¹ for conductivity, combining for an estimated 0.0001 m s⁻² for buoyancy. The precision of Nortek Aquadopp current meters is 0.005 m s⁻¹.

To estimate the total error, we introduce random errors from the two sources simultaneously, and allow these errors to accumulate in the PV budget. In doing so, the errors associated with mooring motion and instrumental noise are assumed to be independent. 95% confidence intervals are finally estimated (Fig. 11) from a Monte Carlo approach by repeating this exercise 1000 times.

# APPENDIX B

## Ekman Transport of Potential Vorticity

The Ekman horizontal velocities can be expressed as

$$u_{\text{Ek}} = \frac{\tau_y}{\rho_0 \delta_e f'}, \qquad (B1)$$

$$v_{\text{Ek}} = - \frac{\tau_x}{\rho_0 \delta_e f'}, \qquad (B2)$$

where $\tau = (\tau_x, \tau_y)$ is the surface wind stress vector and $\delta_e$ is the Ekman layer depth. The Ekman layer depth is approximated by $\delta_e = 0.4u_*/f$, with $u_* = \sqrt{|\tau|/\rho_0}$, where 0.4 is an empirical constant determined from observations (Wang and Huang 2004). The contribution to the material change of PV from Ekman advection, $u_{\text{Ek}}(\partial q/\partial x) + v_{\text{Ek}}(\partial q/\partial y)$, is then estimated by assuming the Ekman flows ($u_{\text{Ek}}$, $v_{\text{Ek}}$) have the same direction at all depths within the Ekman layer, and $\partial q/\partial x$ and $\partial q/\partial y$ obtained from 50-m depth. We find that the Ekman layer depth is almost always shallower than the mixed layer depth with a mean depth of 51 m during winter and mid-spring (not shown), and the Ekman advection of PV is likely insignificant at our study region (Fig. 11).

## REFERENCES

Adams, K. A., P. Hosegood, J. R. Taylor, J. B. Sallee, S. Bachman, R. Torres, and M. Stamper, 2017: Frontal circulation and submesoscale variability during the formation of a southern ocean mesoscale eddy. *J. Phys. Oceanogr.*, **47**, 1737–1753, https://doi.org/10.1175/JPO-D-16-0266.1.

Bachman, S. D., B. Fox-Kemper, J. R. Taylor, and L. N. Thomas, 2017: Parameterization of frontal symmetric instabilities. I: Theory for resolved fronts. *Ocean Modell.*, **109**, 72–95, https://doi.org/10.1016/j.ocemod.2016.12.003.

Belcher, S. E., and Coauthors, 2012: A global perspective on Langmuir turbulence in the ocean surface boundary layer. *Geophys. Res. Lett.*, **39**, 39, https://doi.org/10.1029/2012GL052932.

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.

Bodner, A. S., and B. Fox-Kemper, 2020: A breakdown in potential vorticity estimation delineates the submesoscale-to-turbulence boundary in large eddy simulations. *J. Adv. Model. Earth Syst.*, **12**, e2020MS002049, https://doi.org/10.1029/2020MS002049.

Brannigan, L., 2016: Intense submesoscale upwelling in anticyclonic eddies. *Geophys. Res. Lett.*, **43**, 3360–3369, https://doi.org/10.1002/2016GL067926.

—, D. P. Marshall, A. C. N. Garabato, A. J. G. Nurser, and J. Kaiser, 2017: Submesoscale instabilities in mesoscale eddies. *J. Phys. Oceanogr.*, **47**, 3061–3085, https://doi.org/10.1175/JPO-D-16-0178.1.

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.

—, N. S. Lucas, S. E. Belcher, T. P. Rippeth, A. L. M. Grant, J. Le Sommer, A. O. Ajayi, and A. C. N. Garabato, 2019: The contribution of surface and submesoscale processes to turbulence in the open ocean surface boundary layer.
J. Adv. Model. Earth Syst., **11**, 4066–4094, https://doi.org/10.1029/2019MS001801.

Callies, J., 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.

—, R. Barkan, and A. Naveira Garabato, 2020: Time scales of submesoscale flow inferred from a mooring array. *J. Phys. Oceanogr.*, **50**, 1065–1086, https://doi.org/10.1175/JPO-D-19-0254.1.

Canuto, V. M., 2015: PV dynamics: The role of small-scale turbulence, submesoscales and mesoscales. *J. Geophys. Res. Oceans*, **120**, 6971–6985, https://doi.org/10.1002/2015JC011043.

Czaja, A., and U. Hausmann, 2009: Observations of entry and exit of potential vorticity at the sea surface. *J. Phys. Oceanogr.*, **39**, 2280–2294, https://doi.org/10.1175/2009JPO4024.1.

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.

D'Asaro, E. A., 2014: Turbulence in the upper-ocean mixed layer. *Annu. Rev. Mar. Sci.*, **6**, 101–115, https://doi.org/10.1146/annurev-marine-010213-135138.

—, C. Lee, L. Rainville, R. Harcourt, and L. Thomas, 2011: Enhanced turbulence and energy dissipation at ocean fronts. *Science*, **332**, 318–322, https://doi.org/10.1126/science.1201515.

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.

Dee, D. P., and Coauthors, 2011: The ERA-Interim reanalysis: Configuration and performance of the data assimilation system. *Quart. J. Roy. Meteor. Soc.*, **137**, 553–597, https://doi.org/10.1002/qj.828.

Deremble, B., N. Wienders, and W. K. Dewar, 2014: Potential vorticity budgets in the North Atlantic Ocean. *J. Phys. Oceanogr.*, **44**, 164–178, https://doi.org/10.1175/JPO-D-13-0871.1.

du Plessis, M., S. Swart, I. J. Ansorge, A. Mahadevan, and A. F. Thompson, 2019: Southern ocean seasonal restratification delayed by submesoscale wind-front interactions. *J. Phys. Oceanogr.*, **49**, 1035–1053, https://doi.org/10.1175/JPO-D-18-0136.1.

Erickson, Z. K., and A. F. 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.

—, —, J. Callies, X. L. Yu, A. N. Garabato, and P. Klein, 2020: The vertical structure of open-ocean submesoscale variability during a full seasonal cycle. *J. Phys. Oceanogr.*, **50**, 145–160, https://doi.org/10.1175/JPO-D-19-0030.1.

Ertel, H., 1942: Ein neuer hydrodynamischer wirbelsatz. *Meteor. Z.*, **59**, 277–281.

Evans, D. G., and Coauthors, 2018: Annual cycle of turbulent dissipation estimated from Seagliders. *Geophys. Res. Lett.*, **45**, 10 560–10 569, https://doi.org/10.1029/2018GL079966.

Ferrari, R., 2011: A frontal challenge for climate models. *Science*, **332**, 316–317, https://doi.org/10.1126/science.1203632.

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.

—, and Coauthors, 2011: Parameterization of mixed layer eddies. III: Implementation and impact in global ocean climate simulations. *Ocean Modell.*, **39**, 61–78, https://doi.org/10.1016/j.ocemod.2010.09.002.

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. N., 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.

Hamlington, P. E., L. P. Van Roekel, B. Fox-Kemper, K. Julien, and G. P. Chini, 2014: Langmuir-submesoscale interactions: Descriptive analysis of multiscale frontal spindown simulations. *J. Phys. Oceanogr.*, **44**, 2249–2272, https://doi.org/10.1175/JPO-D-13-0139.1.

Haynes, P. H., and M. E. McIntyre, 1987: On the evolution of vorticity and potential vorticity in the presence of diabatic heating and frictional or other forces. *J. Atmos. Sci.*, **44**, 828–841, https://doi.org/10.1175/1520-0469(1987)044<0828:OTEOVA>2.0.CO;2.

—, and —, 1990: On the conservation and impermeability theorems for potential vorticity. *J. Atmos. Sci.*, **47**, 2021–2031, https://doi.org/10.1175/1520-0469(1990)047<2021:OTCAIT>2.0.CO;2.

Hoskins, B. J., 1974: The role of potential vorticity in symmetric stability and instability. *Quart. J. Roy. Meteor. Soc.*, **100**, 480–482, https://doi.org/10.1002/qj.49710042520.

—, M. E. McIntyre, and A. W. Robertson, 1985: On the use and significance of isentropic potential vorticity maps. *Quart. J. Roy. Meteor. Soc.*, **111**, 877–946, https://doi.org/10.1002/qj.49711147002.

Klein, P., and G. Lapeyre, 2009: The oceanic vertical pump induced by mesoscale and submesoscale turbulence. *Annu. Rev. Mar. Sci.*, **1**, 351–375, https://doi.org/10.1146/annurev.marine.010908.163704.

—, and Coauthors, 2019: Ocean-scale interactions from space. *Earth Space Sci.*, **6**, 795–817, https://doi.org/10.1029/2018EA000492.

Marshall, J. C., and A. J. G. Nurser, 1992: Fluid-dynamics of oceanic thermocline ventilation. *J. Phys. Oceanogr.*, **22**, 583–595, https://doi.org/10.1175/1520-0485(1992)022<0583:FDOOTV>2.0.CO;2.

Maze, G., and J. Marshall, 2011: Diagnosing the observed seasonal cycle of Atlantic subtropical mode water using potential vorticity and its attendant theorems. *J. Phys. Oceanogr.*, **41**, 1986–1999, https://doi.org/10.1175/2011JPO4576.1.

—, J. Deshayes, J. Marshall, A. M. Treguier, A. Chronis, and L. Vollmer, 2013: Surface vertical PV fluxes and subtropical mode water formation in an eddy-resolving numerical simulation. *Deep-Sea Res. II*, **91**, 128–138, https://doi.org/10.1016/j.dsr2.2013.02.026.

McDougall, T. J., and P. M. Barker, 2011: Getting started with TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox. SCOR/IAPSO WG127, 28 pp., http://www.teos-10.org/pubs/Getting_Started.pdf.

McIntyre, M. E., 2015: Potential vorticity. *Encyclopedia of Atmospheric Sciences*, 2nd ed. G. R. North, J. Pyle, and F. Zhang, Eds., Elsevier, 375–383, https://doi.org/10.1016/B978-0-12-382225-3.00140-7.

McWilliams, J. C., 2016: Submesoscale currents in the ocean. *Proc. Roy. Soc.*, **472**, 20160117, https://doi.org/10.1098/rspa.2016.0117.
—, 2019: A survey of submesoscale currents. *Geosci. Lett.*, **6**, 3, https://doi.org/10.1186/s40562-019-0133-3.

Naveira Garabato, A. C., and Coauthors, 2017: Vigorous lateral export of the meltwater outflow from beneath an Antarctic ice shelf. *Nature*, **542**, 219–222, https://doi.org/10.1038/nature20825.

—, and Coauthors, 2019: Rapid mixing and exchange of deep-ocean waters in an abyssal boundary current. *Proc. Natl. Acad. Sci. USA*, **116**, 13 233–13 238, https://doi.org/10.1073/pnas.1904087116.

Ramachandran, S., and Coauthors, 2018: Submesoscale processes at shallow salinity fronts in the Bay of Bengal: Observations during the winter monsoon. *J. Phys. Oceanogr.*, **48**, 479–509, https://doi.org/10.1175/JPO-D-16-0283.1.

Sallée, J. B., E. Shuckburgh, N. Bruneau, A. J. S. Meijers, T. J. Bracegirdle, and Z. Wang, 2013: Assessment of southern ocean mixed-layer depths in CMIP5 models: Historical bias and forcing response. *J. Geophys. Res. Oceans*, **118**, 1845–1862, https://doi.org/10.1002/jgrc.201517.

Schubert, W., E. Ruprecht, R. Hertenstein, R. N. Ferreira, R. Taft, C. Rozoff, P. Ciesielski, and H. C. Kuo, 2004: English translations of twenty-one of Ertel's papers on geophysical fluid dynamics. *Meteor*. Z., **13**, 527–576, https://doi.org/10.1127/0941-2948/2004/0013-0527.

Siegelman, L., P. Klein, P. Rivière, A. F. Thompson, H. S. Torres, M. Flexas, and D. Menemenlis, 2020: Enhanced upward heat transport at deep submesoscale ocean fronts. *Nat. Geosci.*, **13**, 50–55, https://doi.org/10.1038/s41561-019-0489-1.

Skyllingstad, E. D., J. Duncombe, and R. M. Samelson, 2017: Baroclinic frontal instabilities and turbulent mixing in the surface boundary layer. Part II: Forced simulations. *J. Phys. Oceanogr.*, **47**, 2429–2454, https://doi.org/10.1175/JPO-D-16-0179.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.

—, H. Torres, P. Klein, A. F. Thompson, L. Siegelman, J. Wang, D. Menemenlis, and C. Hill, 2020: High-frequency submesoscale motions enhance the upward vertical heat transport in the global ocean. *J. Geophys. Res. Oceans*, **125**, e2020JC016544, https://doi.org/10.1029/2020JC016544.

Taylor, J. R., and R. Ferrari, 2010: Buoyancy and wind-driven convection at mixed layer density fronts. *J. Phys. Oceanogr.*, **40**, 1222–1242, https://doi.org/10.1175/2010JPO4365.1.

Thomas, L. N., 2005: Destruction of potential vorticity by winds. *J. Phys. Oceanogr.*, **35**, 2457–2466, https://doi.org/10.1175/JPO2830.1.

—, 2008: Formation of intrathermocline eddies at ocean fronts by wind-driven destruction of potential vorticity. *Dyn. Atmos. Oceans*, **45**, 252–273, https://doi.org/10.1016/j.dynatmoce.2008.02.002.

—, A. Tandon, and A. Mahadevan, 2008: Submesoscale processes and dynamics. *Ocean Modeling in an Eddying Regime, Geophys. Monogr.*, Vol. 117, Amer. Geophys. Union, 17–38, https://doi.org/10.1029/177GM04.

—, J. R. Taylor, R. Ferrari, and T. M. Joyce, 2013: Symmetric instability in the Gulf Stream. *Deep-Sea Res. II*, **91**, 96–110, https://doi.org/10.1016/j.dsr2.2013.02.025.

Thompson, A. F., A. Lazar, C. Buckingham, A. C. N. 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.

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, W., and R. X. Huang, 2004: Wind energy input to the Ekman layer. *J. Phys. Oceanogr.*, **34**, 1267–1275, https://doi.org/10.1175/1520-0485(2004)034<1267:WEITTE>2.0.CO;2.

Wenegrat, J. O., L. N. Thomas, J. Gula, and J. C. McWilliams, 2018: Effects of the submesoscale on the potential vorticity budget of ocean mode waters. *J. Phys. Oceanogr.*, **48**, 2141–2165, https://doi.org/10.1175/JPO-D-17-0219.1.

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.

—, —, —, D. G. Evans, and Z. Su, 2019b: Wind-forced symmetric instability at a transient mid-ocean front. *Geophys. Res. Lett.*, **46**, 11 281–11 291, https://doi.org/10.1029/2019GL084309.