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

: 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 frequentoccurrences of negative PV within deep(up to about 350m) mixedlayers fromwinter to mid- spring, and positive PV beneath shallow (mostly less than 50m) 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, exceedthedeepestmixedlayer byone orderof magnitude, suggestingthat mixedlayer depthisshapedby 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 ﬂows. 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 determiningthe seasonalevolution of the mixed layer in the open


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 onedimensional 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. 2018Su et al. , 2020Siegelman et al. 2020), triggering restratification and shoaling of the mixed layer.
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 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 andMcIntyre 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 3 13 km 3 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.

a. Mooring data
The data analyzed in this study were primarily collected from a mooring array deployed at an approximate water depth of 4800 m over the Porcupine Abyssal Plain (PAP; 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 highfrequency 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 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 (Dr 5 0.03 kg m 23 ) from a near-surface value at 10 m .
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.758, 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.

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 where v a is the 3D absolute vorticity, = is the spatial gradient operator, b 5 g(1 2 r/r 0 ) is buoyancy (with g as the gravitational acceleration, r as potential density, and r 0 5 1025 kg m 23 as a reference density), f 5 2V sinf is the Coriolis parameter (with V as Earth's angular velocity and f as latitude), k is the vertical unit vector, and u 5 (u, y, 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, is associated with the vertical component of the absolute vorticity, f 1 z, and the vertical stratification N 2 5 ›b/›z, where z 5 k Á = 3 u is the vertical component of the relative vorticity. The horizontal component of PV, q h 5 {[(›w/›y) 2 (›y/›z)](›b/›x)} 1 {[(›u/›z) 2 (›w/›x)](›b/›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: 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 By further assuming that the flow is in thermal wind balance, [(›u/›z), (›y/›z)] 5 (1/f)[2(›b/›y), (›b/›x)], we can derive geostrophic versions of the Ertel PV components, such as and where the h subscript denotes horizontal component. In the geostrophic expression, q hg is always negative and reduces q geo , whereas q y may be either positive or negative. Importantly, SI develops when q hg overcomes q y with q y . 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 andMcIntyre 1987, 1990;Marshall and Nurser 1992), where J represents the advective and nonconservative transport of PV in the ocean, defined as The J vector has advective (qu), frictional (J F 5 =b 3 F), and diabatic [J D 5 2v a (Db/Dt)] components, where 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, v 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., = Á u 5 0), the flux form of the PV equation reduces to For a depth-integrated PV budget, the PV equation in the mixed layer becomes ›q ›t where J D z and J F z denote the vertical components of the diabatic and frictional forcings; the angle brackets hÁi 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 . We found our key findings to be insensitive to this choice. For instance, adopting the q 5 0 surface as an approximation of H can quantitatively alter the magnitudes of J D z and J F z 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 smallscale 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 (›q/›t)
The temporal derivative of PV, ›q/›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 458 to approximately match the cross shape of the mooring array.
2) HORIZONTAL ADVECTION OF POTENTIAL VORTICITY [u(›q/›x) 1 y(›q/›y)] The horizontal advection of PV, u(›q/›x) 1 y(›q/›y), is derived from estimates of q geo 5 (f 1 z)N 2 2 fj›u h /›zj 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 z are obtained from each triangle-shaped area following Stokes' theorem, The vertical stratification N 2 and the vertical shear ›u h /›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(›q/›z)]
The vertical advection of PV, w(›q/›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 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, 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, 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: Below the mixed layer, w(›q/›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 D z )
Following Marshall and Nurser (1992), the diabatic component of the J vector at the ocean surface (hereinafter the diabatic PV flux) is approximated by where B 0 5 2[gaQ net /r 0 c p 1 gbS(P 2 E)] is the surface buoyancy flux, Q net is the air-sea heat flux, c p is the specific heat capacity, a and b are the thermal and haline expansion coefficients, S is the surface salinity, (E 2 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 D z 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 F z )
Following Thomas (2005), the frictional component of the J vector (hereinafter the frictional PV flux) at the ocean surface is approximated by where B e 5 (t 3 k) Á = h b/r 0 f is the Ekman buoyancy flux and t is the surface wind stress. Note that J F z 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, h H where w * 5 (B 0 H) 1/3 is the convective velocity, u * 5 ffiffiffiffiffiffiffiffiffiffiffi jtj/r 0 p is the frictional velocity, u g is the geostrophic velocity vector, u w is the angle between the wind vector and the geostrophic shear, and c 5 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, Du g , is computed by assuming geostrophic balance as Du g 5 j›u g /›zjH 5 j= h bjH/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, f RiB 5 tan 21 (2f 2 N 2 /j= h bj 2 ), and the critical angle, f C 5 tan 21 (2z/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 21808 , f RiB , 21358, and hybrid gravitational/symmetric instabilities will occur when 21358 , f RiB , 2908. (ii) For stable stratification (i.e., N 2 . 0) and cyclonic vertical vorticity, SI is predicted to develop when 2908 , f RiB , f C , with f C , 2458. (iii) For stable stratification (i.e., N 2 . 0) and anticyclonic vertical vorticity, SI is expected for 2908 , f RiB , 2458 with f C . 2458, and hybrid symmetric/centrifugal instabilities will occur when 2458 , f RiB , f 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).

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 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 2400 W m 22 in winter, and the year-mean surface heat flux reaches about 245 W m 22 , 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 upperocean 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 ).
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 28 s 23 , 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 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 nearinertial 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 Ð t 0 (›H/›t) 1 dt and Ð t 0 (›H/›t) 2 dt, where superscripts ''1'' and ''2'' 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 midspring, 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(›H/›x) 1 y(›H/›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 y 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 2H , quantified as in Yu et al. (2019a)], which is associated with the divergence of horizontal flow within the mixed layer [i.e., 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 2H is approximately 30% of that of ›H/›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 FIG. 4. Cumulative time integrals of the temporal change of mixed layer depth Ð t 0 (›H/›t) dt in yellow, Ð t 0 (›H/›t) 1 dt (i.e., mixed layer deepening) in orange, and Ð t 0 (›H/›t) 2 dt (i.e., mixed layer shoaling) in blue.
during relatively short-lived, intermittent frontal events associated 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 subinertial velocity field in the OSMOSIS area can be explicitly illustrated by considering the vertical structure of the currents (Fig. 6). The winter-to-mid-spring period hosts instances of strong flow, exceeding 0.5 m s 21 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 vertical component, and the remaining 30% from the horizontal  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 . 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 Callies et al. 2020). PV estimated with 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 f 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 5 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.

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 D z ) and frictional (J F z ) PV fluxes through the ocean surface. The diabatic PV flux reverses from a sizeable PV extraction (J D z . 0) in winter to a very strong PV injection (J D z , 0) in summer, with respective peak amplitudes of 1 3 10 212 and 23 3 10 212 m s 24 . The diabatic PV flux J D z 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 D z 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 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 5 0 surface (magenta line in Fig. 3c) as the base of the ocean surface boundary layer.

FEBRUARY 2021
Y U E T A L .
fronts in the area. Consistent estimates of J F z are obtained from inner and outer mooring cluster measurements, although the inner cluster-based calculation produces larger J F z 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 F z values are in the range from 26 3 10 213 to 6 3 10 213 m s 24 , and are thus substantially smaller than J D z estimates except in winter to mid-spring, when the magnitudes of J F z and J D z are comparable. The annual-mean surface PV flux (J D z 1 J F z ) 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 D z is found to be directed into the ocean, with a value of 22.46 3 10 213 m s 24 , indicating that diabatic PV injection more than compensates for diabatic PV extraction over the year as a whole. In turn, J F z averages to near-zero values over any period of several months and longer. The annual-mean J F z is thus 21.37 3 10 215 or 21.52 3 10 215 m s 24 (for inner and outer cluster measurements, respectively), i.e., two orders of magnitude smaller than the annual-mean J D z . 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 (h›q/›ti ML ) and total (horizontal 1 vertical) advection of PV (hu Á =qi ML ), integrated between the base of the mixed layer and 50 m, exhibit a statistically significant anticorrelation (R 5 20.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 h›q/›ti ML and total advection of PV hu Á =qi ML in the material change of PV, hDq/Dti 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 D z 1 J F z (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 midspring. 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 h›q/›ti ML , the horizontal advection term hu(›q/›x) 1 y(›q/›y)i ML and the vertical advection term hw(›q/›t)i ML all act to oppose the effect of the diabatic forcing, with hu(›q/›x) 1 y(›q/›y)i ML being the largest term. The evolution of the material change of mixed layer PV 2hDq/Dti ML follows that of the total surface PV flux, J D z 1 J F z , 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 reestimate the material change of mixed layer PV by assuming that Dq/Dt within 250 m , z , 0 is equal to the measured Dq/Dt at z 5 250 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-tomid-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 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 D z and surface frictional PV flux J F z are indicated by the blue and orange lines, respectively. The temporal change of PV h›q/›ti ML , horizontal advection of PV hu(›q/›x) 1 y(›q/›y)i ML , and vertical advection of PV hw(›q/›t)i 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 D z 1 J F z ) is indicated by the black line. The material change of PV 2hDq/Dti 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 2(hDq/Dti ML 1 hDq/Dti 50m ) is shown by the red dashed line. The Ekman transport of PV hu Ek (›q/›x) 1 y Ek (›q/›y)i is indicated by the purple line. The shaded regions illustrate the 95% confidence envelope of cumulative time integrals, estimated using a Monte Carlo approach.

FEBRUARY 2021
Y U E T A L . 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 . 2h, and the forced-SI layer (denoted by the subscript ''SI'') extends across 2H , z , 2h. 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 5 20.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 hDq/Dti CL to surface diabatic and frictional forcings, endorses this interpretation. For example, for downfront wind (J F z . 0) or surface cooling (J D z . 0) conditions, PV in the convective layer is generally found to decrease (i.e., hDq/Dti CL , 0; red shading). Conversely, for upfront wind (J F z , 0) or surface heating (J D z , 0) conditions, PV in the convective layer generally increases (i.e., hDq/Dti CL . 0; blue shading). Further, hDq/Dti 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 5 20.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 h›q/›ti int and the total advection of PV hu Á =qi int are significantly anticorrelated, with a correlation (R 5 20.81, p , 0.001) that exceeds that in the mixed layer (Fig. 13a). The diagnosed balance between h›q/›ti int and hu Á =qi int , characterized by a best-fit linear regression with a slope 20.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 h›q/›ti int and hu Á =qi int on time scales of 5-10 days (Fig. 13b), and a nearzero lag in the phase of the cross-spectrum between both variables (Fig. 13c). Horizontal advection dominates over vertical advection in balancing h›q/›ti int over the frequency range of high coherence (not shown).

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 springsummer 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-midspring 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.