## Abstract

The importance of submesoscale instabilities, particularly mixed layer baroclinic instability and symmetric instability, on upper-ocean mixing and energetics is well documented in regions of strong, persistent fronts such as the Kuroshio and the Gulf Stream. Less attention has been devoted to studying submesoscale flows in the open ocean, far from long-term, mean geostrophic fronts, characteristic of a large proportion of the global ocean. This study presents a year-long, submesoscale-resolving time series of near-surface buoyancy gradients, potential vorticity, and instability characteristics, collected by ocean gliders, that provides insight into open-ocean submesoscale dynamics over a full annual cycle. The gliders continuously sampled a 225 km^{2} region in the subtropical northeast Atlantic, measuring temperature, salinity, and pressure along 292 short (~20 km) hydrographic sections. Glider observations show a seasonal cycle in near-surface stratification. Throughout the fall (September–November), the mixed layer deepens, predominantly through gravitational instability, indicating that surface cooling dominates submesoscale restratification processes. During winter (December–March), mixed layer depths are more variable, and estimates of the balanced Richardson number, which measures the relative importance of lateral and vertical buoyancy gradients, depict conditions favorable to symmetric instability. The importance of mixed layer instabilities on the restratification of the mixed layer, as compared with surface heating and cooling, shows that submesoscale processes can reverse the sign of an equivalent heat flux up to 25% of the time during winter. These results demonstrate that the open-ocean mixed layer hosts various forced and unforced instabilities, which become more prevalent during winter, and emphasize that accurate parameterizations of submesoscale processes are needed throughout the ocean.

## 1. Introduction

One-dimensional, vertical budgets of the ocean surface boundary layer are now recognized to have limitations in reproducing observed changes in mixed layer depth (MLD), especially at subseasonal time scales (Fox-Kemper et al. 2011; Belcher et al. 2012; Sallée et al. 2013). Three-dimensional fluid motions, typically arising from strong lateral buoyancy gradients, may be one mechanism that explains this discrepancy. Lateral buoyancy gradients serve as reservoirs of both potential and kinetic energy, the latter occurring through thermal wind (Haine and Marshall 1998; Ferrari and Wunsch 2009). These gradients may also catalyze instabilities that efficiently transfer energy to smaller scales, which impacts upper-ocean mixing, stratification, ventilation, and ecosystem dynamics (Lévy et al. 2012). These instabilities are largely trapped within the mixed layer and occur at temporal and spatial scales faster and smaller than the ocean’s mesoscale, a regime known as the submesoscale.

Insight into submesoscale processes has, to date, largely been achieved through careful numerical simulations (Haine and Marshall 1998; Boccaletti et al. 2007; Thomas et al. 2008; Capet et al. 2008; Taylor and Ferrari 2009, 2010; Mahadevan et al. 2010; Rosso et al. 2014) that directly resolve the evolution of submesoscale fronts and eddies. Most, but not all of these studies, examine idealized ocean configurations, such as periodic fronts with unidirectional wind stress forcing or constant buoyancy forcing. Observations of submesoscale processes are challenging because of the constraint of capturing intermittent events that evolve rapidly on the order of 1 to 10 days. Observational programs have largely focused on regions that are preconditioned toward an active submesoscale field, for example, western boundary currents, including strong surface forcing and strong frontal gradients (D’Asaro et al. 2011; Thomas et al. 2013). The prevalence of these motions in more quiescent regions of the ocean remains unresolved. The results presented in this study provide direct observational evidence, complimentary to Buckingham et al. (2016), that submesoscale motions are active throughout the ocean, at least during certain times of the year.

The mixed layer is a turbulent environment responding to a variety of mixing processes including internal and surface wave breaking, Langmuir turbulence (Belcher et al. 2012), and convection responding to surface heat and freshwater fluxes (Taylor and Ferrari 2010). The observational component of the Ocean Surface Mixing, Ocean Submesoscale Interaction Study (OSMOSIS) project was designed to observe a full annual cycle of turbulent motions in the ocean surface boundary layer spanning small, dissipative scales to the transition between submesoscales and mesoscales. The time series presented in this study provides a census of instability processes throughout the year. The focus is on those dynamical processes that are likely to impact MLD variability and upper-ocean stratification, specifically surface buoyancy forcing, frictionally generated surface buoyancy forcing (Thomas 2005), symmetric instability (SI; Taylor and Ferrari 2009; Thomas et al. 2013), and mixed layer baroclinic instability (BCI; Boccaletti et al. 2007; Fox-Kemper et al. 2008). With the exception of the surface buoyancy forcing, all of these processes are related to motions that span a spatial range of many hundreds of meters to 10 km and time scales of hours to days. Dynamically, these motions approach a regime where the Rossby number Ro = *ζ*/*f* is *O*(1), where *ζ* is the vertical vorticity and *f* is the Coriolis frequency. Furthermore, the balanced Richardson number

where *N*^{2} is the vertical buoyancy gradient and *M*^{2} is the horizontal buoyancy gradient, is frequently less than or close to one. These values represent a dynamical definition of the ocean’s submesoscale (Thomas et al. 2008). From this characterization, it becomes clear that submesoscale processes are intensified in regions of strong lateral buoyancy gradients, strong vorticity, and weak vertical stratification. These conditions are often met in or near western boundary currents of the Pacific (D’Asaro et al. 2011) and Atlantic (Thomas et al. 2013; Shcherbina et al. 2013), in coastal oceans, and in the Southern Ocean.

Persistent, strong frontal currents, similar to the Gulf Stream, are absent in the interior of midlatitude gyres, yet these regions still host coherent mesoscale eddies either generated locally via mesoscale baroclinic instability or advected into the gyre from boundary currents. There are a number of reasons to investigate whether these features and their associated surface flow fields are sufficient to generate submesoscale variability. First, submesoscale instabilities are an efficient mechanism for generating a forward energy cascade that transfers energy from large-scale geostrophically balanced motions to smaller turbulent scales associated with dissipation (Capet et al. 2008; D’Asaro et al. 2011; Thomas et al. 2013). SI is an example of this process. Energy transfer to dissipative scale may occur if secondary (Kelvin–Helmholtz) instabilities develop (Taylor and Ferrari 2009). Additionally, submesoscale instabilities can give rise to high-frequency variations in MLD. In particular, lateral gradients of buoyancy in the surface mixed layer provide a source of potential energy, which when released through a vertical buoyancy flux may rapidly stratify the upper ocean (Boccaletti et al. 2007). Finally, Thomas et al. (2013) and Omand et al. (2015) have shown that intermittent turbulent isopycnal fluxes, associated with strong straining fields at the edge of geostrophic fronts and eddies, may inject tracers, such as particulate organic carbon (POC), into the interior at depths below the mixed layer. Omand et al. (2015) also considered a region removed from boundary currents and estimated that isopycnal injection of POC may account for as much as half of the total springtime export of POC in subpolar regions.

The OSMOSIS dataset provides a unique opportunity to resolve variations in mixed layer instability characteristics across an entire year. This offers a higher degree of granularity than recent studies that have analyzed bulk seasonal differences in flow characteristics, which are then used to infer submesoscale activity. Both numerical models (Mensa et al. 2013; Sasaki et al. 2014) and observations (Callies et al. 2015; Swart et al. 2015; Buckingham et al. 2016) provide strong evidence that there are seasonal modulations in submesoscale activity. In general, these reflect a steepening of the kinetic energy spectra during summer, typically with slopes around *k*^{−3}, and a shoaling of the spectra during winter to values close to *k*^{−2} (Callies et al. 2015). This wintertime shoaling is an indication that a larger fraction of the total kinetic energy resides at smaller scales. This seasonal cycle is consistent with an enhancement of BCI, as the potential energy stored in the upper-ocean boundary layer scales with the MLD and with the preferential generation of submesoscale eddies during winter. Furthermore, vertical stratification decreases either in response to the deepening mixed layer or increased surface wind and buoyancy forcing, which reduces Ri_{b}. It remains less clear whether seasonal modulations in lateral surface buoyancy gradients also occur. As discussed below, a moderate enhancement of the mixed layer lateral buoyancy gradients are observed during winter at the OSMOSIS site.

This study provides a more direct look at variations in submesoscale characteristics that occur throughout the year by identifying stratification and flow configurations susceptible to different instabilities based on 292 individual short (~20 km), glider-obtained sections in the North Atlantic subtropical gyre. This perspective provides an unprecedented observational record of the pervasiveness of submesoscale motions, and the results, if extended throughout midgyre environments, suggest that submesoscales are likely to play a dominant role in the global ocean energy budget as well as the subduction of heat, carbon, and other tracers.

In the following section, we describe the OSMOSIS field program and the glider observations; we also provide background information on theoretical approaches to quantifying mixed layer instabilities and their contribution to upper-ocean stratification. Section 3 contains a summary of the annual cycle of submesoscale characteristics. Section 4 combines the glider data with existing parameterizations to make additional inferences about the impact of submesoscale motions at the OSMOSIS site; a discussion of error sources is also included. We present our conclusions in section 5.

## 2. Methods

### a. Glider data

The goal of the OSMOSIS project was to constrain the potential vorticity (PV) budget over a small patch of an oceanic gyre throughout a complete annual cycle (September 2012 to September 2013). To accomplish this aim, nine full-depth moorings were deployed in a 15 km by 15 km box at the Porcupine Abyssal Plain (PAP) in the northeast Atlantic, centered at 48.7°N, 16.2°W (Fig. 1). Moored observations from the OSMOSIS region show a seasonal cycle in the vorticity skewness at scales smaller than 5 km, indicative of temporal changes in the intensity of submesoscale turbulence (Buckingham et al. 2016). Moorings were complemented by measurements from ocean gliders that provided continuous sampling from the surface to a depth of 1000 m, obtaining measurements of temperature, salinity, dissolved oxygen, fluorescence, optical backscatter, and photosynthetically available radiation (PAR). Each glider deployment (three separate deployments of two gliders each) lasted at least 4 months with one deployment of 5 months. In this study, we analyze measurements from three gliders concatenated into a single year-long time series. The results are consistent with measurements from the gliders that are not presented here.

The observational design consisted of four outer moorings (Figs. 1d,e) and five inner moorings, the latter clustered in a 5 km by 5 km box. The gliders were piloted within the array describing a “butterfly” pattern: alternating meridional (glider 1) and zonal (glider 2) sections, with diagonal sections across the array (Fig. 1e). The glider observations provide a rich, but challenging, dataset that convolves both spatial and temporal variability as the glider moves through the domain. Additionally, advective processes carry mesoscale and submesoscale features across the array throughout the year. A typical velocity scale at the OSMOSIS site is 0.1 m s^{−1}, which gives rise to an advective time scale of 2.3 days to cross the array. Our approach is to consider the glider dataset as a large number (292) of discrete hydrographic sections that span the domain. Each individual section provides information about the presence of submesoscale structure contained within the array. A compilation of these sections provides a census, or a statistical picture, of submesoscale dynamics in the region.

Each glider dive provides both downward and upward inclined profiles to a depth of 1000 m; each dive took 4 to 5 h to complete. Each glider dive describes a V-shaped pattern spanning a horizontal distance of anywhere between 2 and 5 km depending on the strength of the background flow and the desired position of the next surfacing (Fig. 2). The glider data are a fundamental component of the OSMOSIS observing system since the uppermost mooring measurements begin at 40-m depth, while MLDs can shoal to less than 30 m during summer. The strategy to sample to depths of 1000 m was largely chosen to minimize power consumption.

Two features guided the choice to conduct the OSMOSIS field program at the PAP site. First, this site hosts a multidisciplinary observing system, coordinated by the National Oceanography Centre, Southampton (NOCS),^{1} which has been sustained for over 20 yr to investigate the response of open-ocean and deep-sea ecosystems to climatic change. These long-term measurements provide a valuable historical context for the new observations. Second, the site is characteristic of a considerable fraction of the ocean’s midlatitude gyres, including (i) a deep-water location (4800 m) removed from complex topography minimizing the effect of internal tides; (ii) a MLD that undergoes a significant seasonal cycle, spanning roughly 30 to 300 m; (iii) a location at the eastern edge of the North Atlantic subtropical gyre in a region of relatively weak mean flow and moderate eddy kinetic energy (Fig. 1b; Jia et al. 2011); and (iv) the frequent occurrence of slowly evolving, deep baroclinic eddies with radius of ~50 km and surface velocities of ~0.3 m s^{−1}.

Along-track time series of tracers (salinity) and buoyancy, *b* = *g*(1 − *ρ*/*ρ*_{0}), where *g* is gravity, *ρ* is density, and *ρ*_{0} = 1025 kg m^{−3} is a reference density, reveal variability covering a range of scales (Fig. 2). The glider sampling is designed to have a higher frequency in the upper 300 m than in the lower 700 m. Tracer anomalies are resolved by multiple vertical profiles, such that the structure in Figs. 2c and 2d is not related to aliasing of the along-track data. Temperature and salinity largely compensate in the upper ocean; however, lateral buoyancy gradients persist in the mixed layer with a typical horizontal scale of 5 to 10 km. The water column is rife with compensated thermohaline interleaving features at depths between 600 and 1000 m throughout the entire year. For the analysis described below, in particular the calculation of vertical and horizontal buoyancy gradients, the glider data are optimally interpolated onto a regular grid in the time domain, with a vertical grid resolution of 10 m and a temporal grid resolution of 0.1 days (Figs. 2c,d). The horizontal glider position is then interpolated onto the time series to produce a monotonically increasing along-track distance, which is then used to calculate horizontal spatial gradients. It remains challenging to distinguish variations in the interpolated property fields as along-track variability, temporal variability, a change in the glider orientation after completion of a transect, or some combination of these.

Submesoscale processes will, from the definition Ro ≈ 1, be composed of both geostrophic and ageostrophic components. We are limited in our velocity calculations to using thermal wind balance along the path of the glider with the additional information obtained from the glider-derived depth-averaged current. Below we identify periods when PV is negative that may lead to instabilities, such as gravitational and symmetric instability. It is unlikely that the gliders resolve the scale of individual vertical or slantwise convective plumes. The separation between dives, ~4 km, combined with processing of the glider data likely produces water mass properties that represent an average over multiple plumes. On the other hand, it is clear that the gliders, on a number of occasions, sample a mixed layer that is actively convecting and records PV values less than zero. We emphasize here, and throughout our analysis, that this glider data provide a statistical view of the frequency at which submesoscale processes occur in the OSMOSIS region.

### b. Potential vorticity calculation

We use the Ertel PV as a diagnostic to indicate the potential for flow instabilities. For the first stage of this study, we focus on instabilities that arise when PV takes the opposite sign of the Coriolis parameter or, for our Northern Hemisphere study, when PV is less than zero (Hoskins 1974). This is an approach that has been used recently by Thomas et al. (2013) for observations collected near the Gulf Stream. The generation of negative PV values may occur when the fluid column is unstably stratified (gravitational instability) or experiences horizontally sheared flows (inertial or centrifugal instability) or strong lateral buoyancy gradients (SI). The Ertel PV is defined as

where *ω*_{a} is the absolute vorticity, Ω is Earth’s angular velocity, **u** = (*u*, *υ*, *w*) is the three-dimensional velocity vector, and *b* is the buoyancy as defined above. If lateral variations in velocity and buoyancy are small, PV may be negative if the flow is unstably stratified or *N*^{2} ≡ ∂*b*/∂*z* < 0, which gives rise to gravitational instability. As horizontal buoyancy gradients and horizontal vorticity become stronger, additional terms must be considered, which gives rise to a baroclinic component of *q*_{Ertel}:

where *f* = 2Ω sin*ϕ*, *ϕ* is latitude, *ζ* = *υ*_{x} − *u*_{y} is the vertical relative vorticity, subscripts indicate partial derivatives, and *q*_{nt} contains terms related to the nontraditional component of the Coriolis frequency.

We rely on certain assumptions in our calculation of the PV fields, some relating to the dynamical regime of the flow and others imposed by the glider sampling capabilities. First, we neglect terms in *q*_{Ertel} that involve the vertical velocity *w*. Preliminary estimates of the vertical velocity from the glider data confirm that this is a reasonable assumption. Second, we neglect *q*_{nt}, which also makes a small correction to *q*_{Ertel}. Finally, we will assume the flow is in thermal wind balance, such that the vertically sheared horizontal velocities can be directly related to lateral buoyancy gradients (*u*_{z}, *υ*_{z}) = *f*^{−1}(−*b*_{y}, *b*_{x}). This approximation is justified because the horizontal length scale resolved by the glider is on the order of 4 km, and recent analyses of high-resolution numerical models (Brannigan et al. 2015) show that the flow is in geostrophic balance down to scales of 5 km. Furthermore, Buckingham et al. (2016) show from the mooring data that the Rossby number associated with these scales, calculated from the inner moorings, has values |Ro| < 0.5. The application of these approximations allows *q*_{bc} to be written in a more compact form:

Following Thomas et al. (2013), we identify instability types by introducing the balanced Richardson angle, given by

We refer the reader to Fig. 1 of Thomas et al. (2013), which shows schematically the range of angles associated with gravitational instability (−*π*/2 < < −*π*/4), a mixed gravitational and symmetric instability (−*π*/4 < < 0), and SI (0 < < *ϕ*_{c}). The criterion for SI depends on a critical angle *ϕ*_{c} that is a function of the sign of the relative vorticity:

SI is of particular interest as this process results in vigorous along-isopycnal motions, sometimes referred to as slantwise convection (Taylor and Ferrari 2009). These motions may give rise to a secondary, Kelvin–Helmholtz instability that induces diapycnal mixing and rapidly restratifies the water column. SI, through both restratification and a relaxation of lateral buoyancy gradients, pushes Ri_{b} back toward marginal stability or a value of 1 (Bachman and Taylor 2014).

Ideally, the calculation of PV would employ observations with high vertical and horizontal resolution. The latter is required to fully resolve *q*_{bc} as well as *ζ*. The OSMOSIS gliders provide a high spatial resolution in the vertical but in only one horizontal direction, which additionally has a time-varying orientation; thus, further approximations are necessary to estimate the PV field. Most importantly, geostrophic velocities can only be calculated in the direction perpendicular to the glider’s orientation. Here, we refer to the along-track direction as *x*, regardless of orientation, and the cross-track direction as *y*. This unidirectional velocity field is not a serious constraint on the analysis if the glider orientation during a given transect is perpendicular to the front orientation (Shcherbina et al. 2013). Our choice of sampling pattern results in a glider orientation that is arbitrary with respect to front orientation within the OSMOSIS domain. The error associated with glider orientation is discussed in detail in section 4c, where it is shown that on average the glider underestimates both the magnitude of the horizontal buoyancy gradient and the strength of mechanically driven surface buoyancy forcing.

With the available glider data, an “observational” expression for the PV becomes

where *M*^{2}(*x*, *z*, *t*) = *b*_{x} is now the along-track horizontal buoyancy gradient. This expression makes clear that the contribution of lateral buoyancy gradients to periods of negative PV, as *M*^{4} is a positive definite quantity. A complicating term here is *ζN*^{2}, since *υ*_{x} cannot be estimated from the thermal wind relationship alone. We calculate the barotropic component of *υ* using the depth-averaged current from the glider dives. Thus, *ζ* is an approximation of the vertical relative vorticity, given by the along-track gradient of the cross-track velocity. We find that this term tends to be large in regions where *N*^{2} is weak and *M*^{4} is large and does not tend to modify regions where *q*_{obs} < 0. Finally, the critical angle that determines situations conducive to SI reduces to

### c. Wind-driven and mixed layer instabilities

While SI modifies the near-surface buoyancy field toward a state where Ri_{b} = 1, horizontal buoyancy gradients may persist once the system becomes marginally stable to SI. These lateral gradients provide a reservoir of potential energy that may be released via BCI (Tandon and Garrett 1994; Haine and Marshall 1998; Boccaletti et al. 2007). The relaxation of isopycnals back toward the horizontal occurs over a period of many days, during which time Ri_{b} continues to increase. This relaxation is slower than the equilibration to marginal stability that happens during SI, which is on the order of a day (Bachman and Taylor 2014) but considerably faster than mesoscale variations in the OSMOSIS domain. Fox-Kemper et al. (2008) derive a parameterization for the restratification process in terms of an overturning streamfunction *ψ*_{BCI}, involving the MLD *H* and the amplitude of lateral buoyancy gradients |∇*b*|. We reproduce the parameterization replacing |∇*b*| with the along-stream buoyancy gradient *b*_{x}:

Fox-Kemper et al. (2008) provide a structure function *μ*(*z*) that describes the vertical distribution of the restratification process, which peaks at the middle of the MLD and tapers to zero at the surface and at the base of the mixed layer. Since we will treat (9) as a scaling for the effect of BCI, we will take *μ* = 1 for simplicity, although it is straightforward to retain the depth dependence here. The coefficient 0.06 is determined empirically from numerical simulations and we acknowledge that directly applying this value may not be valid for this regional study. For this reason, variations in restratification due to BCI are emphasized and direct comparison to the surface heat flux (section 3d) below should be viewed with caution. Mahadevan et al. (2012) showed that *ψ*_{BCI} can be represented in terms of a buoyancy flux, or an equivalent heat flux, in order to compare its effects to surface heating and surface freshwater fluxes. The expression for an equivalent heat flux is

where *C*_{p} is the specific heat and *α*(*T*) is the thermal expansion coefficient. Ideally, it would be useful to consider a similar expression for the restratification effect that arises via SI. An expression for this process has not been previously published, although it is an active area of research (J. Taylor 2015, personal communication). Derivation of such an expression is beyond the scope of the present work, although we address the relative importance of forced symmetric instability to surface forcing in section 4a.

Thomas (2005) showed that overturning circulations related to surface Ekman divergence may also have a significant impact on upper-ocean stratification. In particular, a surface wind stress oriented downfront, or in the direction of the geostrophic shear, produces an Ekman transport that advects less buoyant fluid over more buoyant fluid, a condition that is gravitationally unstable. An upfront wind stress, on the other hand, acts to restratify the mixed layer. As was done for BCI, the impact of these processes on the stratification can be expressed in terms of an overturning circulation, which involves the component of the wind stress aligned with the front (Thomas 2005). Again, we alter this parameterization to account for the limitations of the glider data and introduce an overturning streamfunction:

where *τ*^{y} is the component of the wind stress perpendicular to the glider section or perpendicular to the observed lateral buoyancy gradient. We discuss the deviation of *ψ*_{Ek} from the true Ekman overturning in section 4c. With the additional knowledge of the surface buoyancy gradient, the impact on the stratification can be expressed as an equivalent heat flux:

Here and throughout, a negative heat flux reduces the stratification of the mixed layer. Importantly, *Q*_{BCI} is positive definite, whereas *Q*_{Ek} may be either positive or negative.

### d. Additional datasets

In the following section, the impact of submesoscale motions on MLD variability and upper-ocean buoyancy (heat) forcing will be considered. The base of the mixed layer is calculated from the glider data using the criterion that the density must exceed the density at 10 m by a value of Δ*ρ* = 0.03 kg m^{−3}. Additional criteria based on density gradients, temperature differences, and temperature gradients give qualitatively similar mixed layer variations. The surface heat and freshwater fluxes as well as the surface wind stress are taken from the ECMWF ERA-Interim reanalysis product (Dee et al. 2011) and are averaged over a period of 1 day. Sea surface temperature (SST) data are obtained from the GHRSST level 4 Multiscale Ultrahigh Resolution (MUR) global analysis that incorporates measurements from multiple instruments including Advanced Microwave Scanning Radiometer 2 (AMSR2) and MODIS.^{2}

## 3. Results

### a. Characterization of the PAP site

The year-long atmospheric and hydrographic observations collected at the OSMOSIS site from September 2012 to September 2013 reveal annual cycles in surface-observable properties (Fig. 3). SST steadily decreases from early September to the end of December, after which SST remains uniformly cool around 12°C. These cool temperatures persist throughout the winter, abruptly warming in May with the most rapid increase in SST occurring in late June, consistent with the period of the strongest surface heat flux. The surface velocities are characteristic of a region of the ocean with moderate eddy kinetic energy and weak mean flow. Remotely sensed products such as SST and sea surface height show the clear signature of multiple coherent mesoscale vortices that pass through the OSMOSIS site throughout the year. The surface eddy kinetic energy, derived from mesoscale measurements of sea surface height, varies on relatively slow temporal scales, peaking at 0.025 m^{2} s^{−2} in May. Depth-averaged currents calculated from the glider trajectories show more variability in upper-ocean velocities than is apparent from the 0.25° gridded sea level anomaly product from AVISO (Damerell et al. 2016, manuscript submitted to *J. Geophys. Res.*). The surface heat flux cools the ocean surface throughout most of the year, with the transition to a steady positive heat flux occurring in late May (Fig. 3c). The surface heat flux is strongest in June and July and weakens in August. The freshwater flux at the OSMOSIS site is more intermittent, with episodic rainfall exceeding 5 mm day^{−1}. However, as the site is located along a boundary of net precipitation to the north (over the subpolar gyre) and net evaporation to the south (over the subtropical gyre), the annual-averaged freshwater flux (precipitation minus evaporation, *P* − *E*) is small: 0.014 mm day^{−1}.

Surface wind forcing at the OSMOSIS site is highly variable, although the array is small enough that the wind stress is uniform across the OSMOSIS site (Fig. 4). The year-long mean wind stress is 0.09 N m^{−2}, but occasional wind stress values approach 1 N m^{−2}. Unlike wind forcing over the western boundary currents, the orientation of the wind stress is highly variable. Figure 4b provides a time series of wind stress orientation over the duration of the OSMOSIS field program in a manner that retains information about the wind stress rotation, for example, over the fall the rotation is largely clockwise, accumulating a negative angle. The color indicates the modulus of the angle between −*π* and *π*. Although the orientation is variable, the decorrelation time scale of the wind orientation over the entire year is 1.30 days, where we have estimated the decorrelation scale from the *e*-folding scale of the autocorrelation function. Since, on average, a glider leg across the OSMOSIS site is completed in a little under a day, below we will use the approximation that the wind orientation is constant during a single glider transect.

### b. Upper-ocean potential vorticity

Using along-track definitions provided in section 2b, it is possible to generate a year-long times series of potential vorticity (Fig. 5). The dominant signal throughout the year is the strong stratification at the base of the mixed layer, which becomes weaker and considerably more variable throughout winter; multiple instances of rapid surface restratification occur throughout the winter. The restratification of the surface boundary layer occurs in two stages across the spring–summer transition. The first occurs around day 245 (3 May), caused by the appearance of a cyclonic intrathermocline eddy that displaces the isopycnals toward the surface. The permanent stratification of the surface coincides with a transition to a persistent positive heat flux around day 260 (18 May). Over the full depth of the water column, there is a stronger vertical stratification in the depth range that hosts the pronounced lateral salinity gradients extending from the outflow across the Straits of Gibraltar, 600 to 1000 m (Koltermann et al. 2011).

Our focus in this study is on instances where PV becomes negative because of changes in vertical or horizontal buoyancy gradients. Figure 6 summarizes changes in the lateral buoyancy gradient *M*^{2} throughout the year. During summer (September 2012 and July and August 2013) the mixed layer is shallow, and changes in *M*^{2} are strongly correlated with changes in the MLD. Variations in MLD across consecutive dives, which are likely related to internal wave fluctuations, can produce strong lateral buoyancy gradients in the interpolated fields. These are neglected since they do not represent reservoirs of potential energy. In the fall, day 85 (Fig. 6b) marks a transition between a regime where the largest *M*^{2} values occur along the base of the mixed layer to a regime where the largest *M*^{2} values are found within the mixed layer itself. This behavior persists throughout the winter (Fig. 6c). As density gradients at the base of the mixed layer weaken, strong lateral buoyancy gradients are able to extend to depths of 500 m or more. Strong lateral buoyancy gradients in the interior may also develop along the flanks of intrathermocline eddies (Fig. 6d, day 250); the core of the eddy is identified by minima in PV and *M*^{2}. The strength of *M*^{2} at both the middle of the MLD *H*/2 and at a fixed depth of 25 m is summarized in Fig. 6a. Values of *M*^{2} exceeding 10^{−7} s^{−2} occur less than 5% of the time. The large *M*^{2} values at 25 m during spring are a result of this depth sitting at the base of the mixed layer. Both winter and fall periods show an elevated probability of observing lateral buoyancy gradients in the range 10^{−8} to 10^{−7} s^{−2}, while close to 90% of lateral buoyancy gradients during the spring–summer period are less than 10^{−8} s^{−2}. This provides evidence that there is a seasonal cycle in the strength of horizontal buoyancy gradients in the mixed layer.

A closer look at an individual hydrographic section reveals further details about submesoscale structure in this region. Figure 7 shows a section (day 114, 23 December) in which the glider track is oriented perpendicular to both the surface SST gradient and, more importantly, to the mixed layer buoyancy gradient. SST shows a nearly 1°C change in surface temperature over less than 10 km. Figures 7b–e show depth and along-track variations in potential temperature (*θ*) and salinity (*S*), both from the glider profiles themselves and the interpolated values for comparison. Subsequent panels only show interpolated values. The ocean is predominantly temperature stratified, but strong, compensated variations in *θ* and *S* are found throughout the water column. Finally, since *S* acts more like a passive tracer, we can see that lateral gradients here are also associated with high salinity variance and high spiciness variance (not shown) throughout the water column (Ferrari and Rudnick 2000; Damerell et al. 2016, manuscript submitted to *J. Geophys. Res.*). Although the temperature and salinity distributions are compensating, a lateral buoyancy gradient persists in the mixed layer (Fig. 7f). For this particular front,

Figures 7f and 7g show that the base of the mixed layer deepens from west to east across most of the section but shoals again over the easternmost 5 km. Overall, the MLD varies by almost 100 m over a 20-km section.

Calculating the Ertel PV *q*_{obs} [(7)] from the buoyancy field (Fig. 7g), the base of the mixed layer is a region of enhanced PV, dominated by the vertical stratification. However, within the mixed layer, the white curve indicates a region where *q*_{obs} < 0. The plunging of this contour throughout the depth of the mixed layer, where the lateral buoyancy gradient *M*^{2} is strongest, indicates the presence of SI. Indeed, in Fig. 7h, the ratio of the square of the lateral buoyancy gradient to the vertical stratification is given by the nondimensional term *M*^{4}/*f*^{2}*N*^{2} = . The weak vertical stratification allows the lateral buoyancy gradient contribution to the PV to be as much as two orders of magnitude larger, conditions that could support SI. The banded structure in Fig. 7h may arise from different stages of SI having modified Ri_{b} back toward *O*(1) or neutral stability.

### c. Instabilities in the mixed layer: Seasonal variations

Using the year-long time series from the glider data, we can now document the frequency and classification of instabilities where PV < 0. An important assumption here is that the cross-track velocities, or the along-track density variations, make the dominant contribution to *q*. This is, of course, only a reasonable assumption when the glider track is perpendicular to a front. Thus, any individual submesoscale instability may be poorly represented by the glider-derived estimate of PV. However, it is likely that the gliders underestimate the intensity of submesoscale fronts and thus underestimate the frequency of SI in the OSMOSIS region (section 4c). Our approach assumes that the gliders’ ability to observe many submesoscale features provides a robust interpretation of the instabilities that occur in the mixed layer throughout the year.

Following Thomas et al. (2013), we use [(5)] to classify instabilities along the glider track. Figure 8 collapses the depth–along-track calculation of the into a single time series by showing the fraction of the mixed layer where *q*_{obs} < 0. The time series is divided into three panels (three different deployments), which naturally partitions the glider time series between fall and early winter (upper panel), winter and early spring (middle panel), and summer (bottom panel). Between day 20 and 80, the mixed layer is shallow, and only a few grid points are classified as being within the mixed layer. This produces the spiky structure of the instabilities during this period. This early period is dominated by gravitational instability, which is consistent with the strong cooling that occurs between September and December. As the mixed layer deepens, for example, around day 85 (Fig. 6), the diversity of the instabilities increases. In particular, a strong convective event occurring at day 83 coincides with a rapid deepening of the mixed layer (~40 m day^{−1}). This is immediately followed by a period of enhanced mixed layer variability as well as an appearance of conditions conducive to SI (Fig. 8). The mixed layer shoals again at day 100, coincident with a reduction in the occurrence of negative PV events.

Around day 105, a second regime starts in which negative PV events appear intermittently. This coincides with the deepening of the mixed layer to values close to the wintertime maximum (~300 m). This second regime extends throughout the winter, during which we frequently observe conditions conducive to SI that persist for a day or more. These instances also tend to occupy a larger fraction of the MLD. The period between days 135 and 150 is particularly active, while the period between days 165 and 180 is more quiescent. This behavior is closely related to surface forcing as discussed in the following section. This analysis does not follow the evolution of any particular instability event. The combination of advection of fronts across the OSMOSIS domain and the constant reorientation of the glider makes this difficult to achieve. The emphasis is on documenting the frequency and classification of instabilities in regions of negative PV.

Instances of negative PV stop abruptly around day 235. This behavior is closely linked to the characteristics of the mixed layer, with a layer of very high PV and strong vertical stratification developing at the surface. While day 235 does not correspond to the first zero crossing of the surface buoyancy flux in the OSMOSIS region, it does correspond to the first time that the surface heat flux exceeds 20 W m^{−2}. The mixed layer is observed to deepen again between days 250 and 260, but there are no significant periods with PV < 0.

To summarize, at the onset of cooling in the fall, instabilities in regions of negative PV are dominated by gravitational instability, consistent with a deepening of the mixed layer throughout this period by convective mixing. As the mixed layer deepens, negative PV events become more often associated with lateral buoyancy gradients, indicative of SI. Lateral buoyancy gradients may also lead to active BCI, as discussed further below. Finally, these instabilities cease abruptly with the onset of a positive heat flux in the spring. During the restratification period, remnants of deep mixed layers remain, but there is little evidence of negative PV events, emphasizing the importance of surface processes for preconditioning instabilities.

### d. Restratification processes

We now modify our focus from only those instabilities that occur for PV < 0 and address other physical processes that may modify the stratification of the mixed layer: restratification due to surface fluxes, surface wind forcing, and the release of potential energy by BCI (Stone 1966). In the previous section, a state of zero PV represents a neutrally stable state, even with *N*^{2} > 0, such that any destabilizing force will immediately trigger small-scale diapycnal mixing. Surface fluxes and BCI, on the other hand, do not simply restore PV to its neutrally stable state but can modify PV to be positive through a restratification of the mixed layer. Our analysis will not account for situations where a period of sustained restratification generates conditions where a destabilizing forcing must reduce PV to zero prior to the onset of small-scale diapycnal mixing. This approach is adopted partially for simplicity but also because the glider data do not sample in a Lagrangian sense, making it difficult to know the time history of a fluid parcel before and after it is observed by the glider.

Similar to the calculation of potential vorticity along a single transect, there are complications with attempting to estimate the size of the Ekman buoyancy forcing term and BCI when the glider orientation is changing frequently. In general, the glider path is oriented randomly with respect to mixed layer gradients and the surface wind stress (Fig. 9). Yet, both the glider and wind stress orientations are relatively stable across a single transect (Fig. 10); the time to complete a single transect is *O*(1) day. Thus, for the following analysis, we determine the glider orientation, the wind stress orientation, and amplitude as an average over the dives of a given transect. An example is shown in Fig. 10. Conversely, the lateral buoyancy gradient and MLD can vary significantly over a single transect as shown in Figs. 10b and 10c. Therefore, the calculation of *Q*_{BCI} [(10)] and *Q*_{Ek} [(12)] is continuous in time. Calculation of *Q*_{BCI} and *Q*_{Ek} using continuous glider and wind orientations is noisier but qualitatively similar.

During winter, the heat flux *Q*_{surf} is largely negative (Fig. 11a), with isolated instances of surface heating occurring after day 160. The decorrelation time scale of the surface heating is roughly a week. The equivalent heat flux coming from the freshwater forcing is small compared with the surface heat flux during this period. In contrast, changes in the surface stratification arising from submesoscale processes are more intermittent in time (Fig. 11b). The equivalent heat flux from the surface wind forcing may take both signs, reflecting downfront (negative) and upfront (positive) wind orientations, whereas restratification due to BCI implies an equivalent heat flux that is always positive. The decorrelation time scale for both processes is approximately 6 h. In addition to the abruptness of these features, their magnitudes are considerably larger than the surface heat flux. Winds generate surface forcing amplitudes in excess of 500 W m^{−2} more than 10 times during a 3-month period. BCI tends to generate an even larger equivalent heat flux than the wind forcing. There is also evidence, for example, day 218 to 222, of a strong mixed layer restratification immediately following the destruction of stratification due to wind forcing. While the exact value of the restratification due to BCI relies on empirical constants not tested in the field, the dependence of *Q*_{BCI} on *H*^{2} suggests that strong restratification occurring during a period of deep mixed layer and strong lateral buoyancy gradients may limit the extent to which *H* can grow.

Summing *Q*_{surf}, *Q*_{BCI}, and *Q*_{Ek} to arrive at a total surface buoyancy forcing *Q*_{tot} (Fig. 11c) indicates a persistent negative buoyancy forcing during the winter punctuated by occasional instances of strong negative forcing from winds but more frequent restratification events mainly due to BCI. Convective events, identified by *N*^{2} < 0, happen sporadically throughout the winter season. The gray blocks in Fig. 11c indicate periods when at least 15% of the mixed layer has a buoyancy distribution that is statically unstable; results are qualitatively similar for different fractions of the mixed layer between 5% and 25%. Considering the sampling pattern of the gliders, linking the surface forcing to the dynamical evolution of these instability events must be done cautiously. Nevertheless, there is a striking correlation between instances of strong positive buoyancy forcing arising from BCI and stable stratification in the mixed layer. Gravitational instability is more frequent during periods of persistent negative buoyancy forcing, whereas all instances of the total surface forcing exceeding 500 W m^{−2}, except for day 152, are associated with stable stratifications throughout the mixed layer. These results, occurring in a region without high levels of eddy kinetic energy, suggest that Ekman buoyancy forcing and BCI, extracting energy from submesoscale fronts, are likely to play a critical role in setting wintertime MLDs globally and not just in persistent western boundary currents.

## 4. Discussion

### a. The impact of submesoscales on surface restratification

In the previous section, observations of buoyancy distributions, MLDs, and surface forcing were used to examine the impact of submesoscale processes on the restratification of the ocean surface boundary layer. Critically, the methods used to arrive at these equivalent heat fluxes, in particular BCI, derive from parameterizations that summarize an ageostrophic overturning that arises from correlations between eddy fluctuations of the vertical velocity and the buoyancy . Therefore, the evolution of any particular submesoscale front observed by the gliders may undergo quite a different evolution than that predicted by the parameterized overturning. However, over the course of the year, and especially throughout late fall and winter when the mixed layer is deep, the general pattern emerges of a moderate but persistent negative heat flux that is frequently spiked by more intense forcing. BCI has the greatest impact on the mixed layer stratification, as compared to the wind forcing and surface buoyancy forcing, at least for the parameters used in (9), and is single signed, always acting to restratify. As the strength of this parameterized overturning depends on the square of the MLD, it provides a strong constraint on wintertime MLDs.

The wintertime series (day 130 to 230) shown in Fig. 11 can be binned by heat flux or equivalent heat flux to produce a cumulative distribution function (CDF; Fig. 12). Surface forcing alone (heat and freshwater fluxes) *Q*_{surf} (black curve) acquires almost exclusively negative values and is confined to fluxes with amplitudes less than 500 W m^{−2}. Similar CDFs for *Q*_{BCI} and *Q*_{Ek} have much longer tails, and the asymmetry in *Q*_{BCI} produces a skewness in the distribution of *Q*_{tot}, where

The difference between the CDF of *Q*_{tot} and *Q*_{surf} (dashed curve) indicates that forcing due to submesoscale fronts can reverse the sign of the equivalent surface buoyancy forcing up to 25% of the time during the winter.

Using these binned equivalent heat flux values, the mean MLD for a given surface forcing is displayed in Fig. 12b; black and red refer to *Q*_{surf} and *Q*_{tot} in Fig. 12a. Darker colors indicate that at least three observations went into calculating the mean MLD for a given heat flux bin. During winter, the mean mixed layer depth is only weakly dependent on the surface heat flux (black squares), with a mean value of ~200 m. In contrast, analysis of the full surface forcing *Q*_{tot} (red circles) shows that shallow MLDs are associated with the strongest fluxes, both positive and negative. For positive values of *Q*_{tot}, this can be interpreted as a period when *Q*_{surf} has been weak. For negative values of *Q*_{tot}, the largest values are dominated by *Q*_{Ek}, which is required to achieve negative forcing with amplitudes greater than 500 W m^{−2}. In this regime, restratification due to BCI is unable to oppose the wind-driven deepening because the mixed layers are shallow and contain relatively little potential energy.

This result suggests that in a region where persistent submesoscale fronts exist, a rough estimate of the maximum MLD may be achieved by assuming that wind-driven mixed layer deepening is balanced by restratification due to BCI. Indeed, this balance may exist around day 220 (8 April), in Fig. 11b, when a strong negative *Q*_{Ek} is rapidly followed by a strong restratification arising from *Q*_{BCI}. This is accompanied by a rapid deepening and subsequent shoaling of the mixed layer at the same time (Fig. 8b). Equating *Q*_{BCI} and *Q*_{Ek} in (10) and (12) gives the relationship

This expression is related to the “restratification ratio” defined by Mahadevan et al. (2010), where *H*_{max} is the value that makes this ratio equivalent to 1. However, our interpretation here is slightly different in that this ratio may be interpreted as bounding winter MLDs. Since the winter buoyancy flux is persistently negative, (15) can be modified to include the surface buoyancy flux:

Note that the term associated with the surface buoyancy forcing has a stronger sensitivity to the strength of the lateral buoyancy gradients. Neglecting *Q*_{surf} and using characteristic values observed by the OSMOSIS gliders during the period of strong wind forcing in early winter, *τ*^{y} = 0.4 N m^{−2}, *ρ*_{0} = 1025 kg m^{−3}, and *b*_{x} = 1 × 10^{−7} s^{−2}, predicts a MLD of 255 m, which is consistent with observations of the deepest wintertime MLDs around 300 m. The MLD exhibits substantial variability throughout the winter, and therefore *H*_{max} should be interpreted as a bound on the mean winter MLD. Furthermore, the amplitude of *b*_{x} may be constrained by the resolution of the glider data, and smaller-scale structures that yield larger lateral buoyancy gradients would reduce the value of *H*_{max}. Nevertheless, this is an intriguing relationship that could be exploited in parameterizations representing the impact of submesoscale dynamics on MLDs.

The analysis resulting in (15) should be viewed with some circumspection. First, our approach assumes the estimated fluxes are uniformly distributed throughout the mixed layer, or at least that they are all evaluated at a single depth. However, each of these flux contributions has a distinct vertical structure, of which some have been explored in numerical simulations but remain unvalidated by observations. The ability of the glider to resolve the vertical structure within the mixed layer is marginal, especially when the mixed layer is shallow, and thus trying to describe the full vertical structure of the instabilities in the mixed layer is beyond the scope of this study. Nevertheless, this complexity is reflected in the tendency in late winter and early spring for there to be a separation between the depth of an actively mixing layer and the MLD as defined by hydrographic criteria (A. Rumyantseva et al. 2016, manuscript submitted to *Limnol. Oceanogr*.).

To further link the restratification analysis to the observed SI-favorable conditions at the OSMOSIS site, it would be useful to apply an overturning relationship to the SI events, similar to (12) and (10). Without this parameterization, though, the importance of SI, as compared to wind and buoyancy forcing, can be assessed using the convective layer depth *h* proposed by Taylor and Ferrari (2010). Appealing again to the vertical structure of processes occurring in the mixed layer, restratification by SI is only effective below the convective depth. In regimes of strong (negative) surface forcing, convection extends throughout the mixed layer and dominates SI. Taylor and Ferrari (2010) used the ratio of the convective depth to the MLD *h*/*H* to predict when SI would be more important than convection. The ratio *h*/*H*, derived from frontal strength, mixed layer depth, and surface fluxes, is a prediction of SI conditions that is independent of PV in the mixed layer. We use the observed wintertime series of surface buoyancy forcing and Ekman buoyancy flux (Fig. 12) to calculate *h* following equation (24) of Thomas et al. (2013), with *α* and *β* set to 0. We find that most of the significant SI events, as diagnosed from calculating (Fig. 8), coincide with low values of *h*/*H* (Fig. 13). This is especially true in late February through April (days 170 to 230); earlier in the winter the surface buoyancy flux is more intense. Figure 13b probes the relationship between *h*/*H* and MLD and shows that the reduction in *h*/*H* is *not* caused by an increase in *H*. Most SI-favorable events are associated with moderate values of MLD (~150 m). This result is consistent with the deep mixed layers being associated with periods of strong surface forcing and convection. Since the OSMOSIS site frequently hosts submesoscale fronts, when the surface forcing weakens, SI is able to take over (a low *h*/*H* regime). Again, the gliders do not sample in a Lagrangian sense, so after observing a deep mixed layer, the restratification process is not necessarily captured. In summary, the time series of *h*/*H* confirms that SI will also have an active role in intermittently modifying the stratification of the mixed layer, especially in the late winter and early spring before the permanent restratification of the mixed layer.

### b. Mixed layer instability wavenumber and growth rates

We have used parameterizations to infer the strength of submesoscale activity. While these parameterizations, some of which involve linear theory, have been successfully tested against numerical simulations, in situ validations have been difficult to achieve. It is beyond the scope of the present study to follow the evolution of one or more submesoscale fronts or filaments and determine whether its evolution is in agreement with predictions. This task is complicated by the small size of the domain and lateral advection of features through the site. Nevertheless, the availability of this high-fidelity time series allows us to diagnose characteristics of BCI.

We now consider the growth rate and fastest-growing mode following the analysis by Stone (1970) and reproduced by Fox-Kemper et al. (2008). The linear growth rate ^{−1} is a function of wavenumber, Ri, and the mixed layer baroclinic vertical shear. The wavelength of the fastest-growing mode ℓ_{max} and the inverse growth rate of the fastest-growing mode _{max} are given by

Here, the vertical shear, represented by a velocity difference *U*, is determined by multiplying *H* by the geostrophic shear calculated from the glider-derived lateral buoyancy gradient. The mean values of *b*_{x} and *H* are calculated for each glider section, leading to 292 estimates of ℓ_{max} and _{max}.

The year-long time series shows a distinct seasonal cycle in ℓ_{max} but not in _{max} (Fig. 14). The size of ℓ_{max} is smaller in late spring and summer compared to the winter. The increase in ℓ_{max} during winter is consistent with an increase in the MLD and an associated increase in the mixed layer deformation radius. Note that the wavelength of the fastest-growing mode during winter, ~12 km, is comparable to the size of the OSMOSIS domain. Taking the whole year into consideration, a histogram of ℓ_{max} (Fig. 14b) shows a broad peak at scales between 5 and 20 km, all considerably smaller than the (mesoscale) first baroclinic Rossby deformation radius. Consecutive glider transects may sample conditions that cause _{max} to vary by an order of magnitude or more. The _{max} values are consistent with submesoscale perturbations with inverse growth rates peaking at values between around 3 days and occasional conditions that indicate inverse growth rates less than a day. The simulations carried out by Fox-Kemper et al. (2008) have values of ℓ_{max} = 3.9 km and _{max} = 16.8 h (for Ri = 1), which are conditions that are captured by the glider but are, respectively, smaller and more rapid than most values inferred here.

In summary, there is evidence that the eddy scale associated with the fastest-growing mode exhibits seasonal changes, whereas the growth rate of these perturbations is more uniform, with most probable values of a few days.

### c. Observational bias

The challenges in interpreting both potential vorticity and mixed layer overturning circulations from the along-track glider data are significant. The identification of property distributions that are consistent with SI and BCI events depend strongly on capturing mixed layer horizontal buoyancy gradients. The observed buoyancy gradient from the gliders will always underestimate this effect as crossing a submesoscale front at any angle different from a perpendicular orientation increases the distance over which a given buoyancy change is observed (assuming the front is uniform in one direction). Averaging over all possible angles leads to a root-mean-square error (underestimate) of the buoyancy gradient amplitudes by a factor of 1/.

Estimating the representativeness of buoyancy forcing resulting from surface winds is more difficult because it involves the orientation of three vector quantities: the glider path, the horizontal buoyancy gradient, and the wind direction. Figure 15 explores the glider interpretation of the Ekman buoyancy flux (EBF). The inset panels show three idealized scenarios of a surface buoyancy gradient *b*_{x} > 0 sampled with a glider orientation shown in black and an imposed wind stress shown in orange. Observationally, the ideal orientation is a glider path that is perpendicular to both the front and the wind orientation. In Fig. 15a, the glider observations would capture most of the amplitude of the Ekman buoyancy flux. However, it is possible to grossly underestimate the EBF if the glider path is aligned along the front as in Fig. 15b. The glider data may also grossly overestimate the EBF when the buoyancy observations are interpreted to have a lateral gradient orthogonal to the wind stress, when in reality the wind stress is aligned perpendicular to the front, as in Fig. 15c. The curve showing the ratio of the observed EBF to the actual EBF is constructed by taking all of the possible angles of the glider and wind orientation with respect to a fixed buoyancy gradient and plotting the histogram.^{3} The most probable values of the ratio (Fig. 15) correspond to capturing almost all or little of the true EBF. However, there are long tails in the distribution, which can lead to errors not only in the magnitude of the EBF but also the sign. Ultimately, though, the observed EBF is smaller than the actual EBF. For instance, the ratio of the root-mean-square of EBF_{obs} to EBF_{actual} is 0.71. This gives us confidence that over the duration of an entire year, or even over a season, the glider sampling is likely to underestimate the occurrence of large wind-driven buoyancy forcing.

### d. Implications

The results of this study provide evidence that submesoscale instabilities, and the seasonal variability in these processes, are active in a region that is characteristic of a large portion of the ocean’s surface. While the OSMOSIS project represents the most comprehensive field program to date aimed at studying submesoscale motions over a full year, it also provides support for earlier studies that have made use of autonomous vehicles in similar regions. In particular, our calculation of mixed layer overturning circulations are consistent with the suggestion by Omand et al. (2015) that eddy-driven vertical transport, generated by ageostrophic overturning near filamentary features, *O*(1–10) km, make a significant contribution to the total vertical flux of POC, up to 50% in the North Atlantic and parts of the Southern Ocean. Both Omand et al. (2015) and the present study use the same parameterization for BCI, and observational verification of the turbulent fluxes arising from these processes are still required.

Additionally, in this study, the direct calculation of an Ekman-generated mixed layer overturning [(11)] indicates that the interaction between surface winds and submesoscale fronts can induce changes, either restratifying or destroying stratification, that are comparable to BCI. In the case where the surface wind forcing acts to remove stratification and deepen the mixed layer, this may provide a mechanism for creating larger reservoirs of potential energy and more vigorous vertical velocities. This sensitivity to the orientation of local winds and finescale fronts (section 4c) is a particularly difficult feature to capture as a parameterization in a GCM.

The purpose of this study is to provide a first overview of the mesoscale and submesoscale structure captured in the OSMOSIS domain. This has allowed us to estimate the impact of these features on the upper-ocean stratification using existing parameterizations. We acknowledge, however, that there have been few, if any, in situ validations of parameterizations, which in some cases have been developed from more idealized process model configurations. Using the OSMOSIS data to document the temporal evolution of one or more of these features is a tantalizing target. The static nature of the OSMOSIS domain make the attribution of observed changes in the mixed layer challenging as features are changing rapidly due to surface forcing, advection, and local dynamics.

## 5. Conclusions

Between September 2012 and September 2013, ocean gliders were used to observe mixed layer variability at ocean submesoscales at the Porcupine Abyssal Plain site in the North Atlantic. This provides the first year-long observational record of submesoscale processes in a region removed from a strong, persistent, geostrophically balanced mean flow. The major results include the following:

The mixed layer is shown to host strong submesoscale fronts, with horizontal scales of 5 to 10 km, throughout the year. Mixed layer buoyancy gradients show a seasonal cycle in amplitude, elevated in late fall and winter.

A distinct seasonal cycle in the type and frequency of instabilities occurring in regions of negative PV is revealed in a region of moderate eddy kinetic energy. Conditions conducive to SI are rampant throughout late fall and winter. From existing parameterizations and estimates of the convective layer depth, both BCI and SI are likely to contribute to mixed layer restratification.

There is evidence that buoyancy forcing induced by a wind-driven Ekman transport may be a precursor to SI and then BCI in cases where the Ekman buoyancy flux reduces the mixed layer stratification or equivalently destroys potential vorticity.

Parameterization of surface buoyancy forcing, cast in terms of an equivalent heat flux, arising from Ekman processes and mixed layer BCI suggests that the OSMOSIS region intermittently undergoes a forcing that is significantly larger than the surface heat and freshwater fluxes. Including the effects of these processes leads to a total surface buoyancy flux that is positive roughly 25% of the winter period.

Using linear theory, the fastest-growing wavelength of BCI in the mixed layer exhibits seasonal variability, whereas the inverse growth rate of this mode is more uniformly distributed throughout the year, with typical values between 3 and 5 days.

The sampling pattern of the gliders is likely to underestimate the amplitude of mixed layer buoyancy gradients as well as the inferred strength of mixed layer overturning circulations. Thus, submesoscale instabilities may be even more prevalent than is documented here, although we caution that any single event may be misinterpreted because of limitations in the glider data.

This work provides evidence that submesoscale motions have a significant impact on upper-ocean stratification throughout the global ocean and need to be represented in GCMs and climate models. The intermittency of events captured by the gliders indicates that the mixed layer responds sensitively to the interaction between surface wind and submesoscale fronts on scales of only a few kilometers. This will be a particularly difficult aspect of these flows to incorporate into parameterizations.

## Acknowledgments

We are grateful for the efforts of the OSMOSIS team and the crews of the RRS *Discovery*, the R/V *Celtic Explorer*, and the RRS *James Cook*. The data presented in the study are the result of five research cruises and a full year of glider piloting and involved many contributions not represented in the author list. The authors acknowledge constructive reviews from John Taylor and Jörn Callies as well as helpful conversations with Liam Brannigan, Baylor Fox Kemper, Amala Mahadevan, Leif Thomas, Patrice Klein, and Eric D’Asaro. AFT was supported by NSF-OCE 1155676, AL was supported by NSF-OCE 1235488, CB and ACNG were supported by OSMOSIS NERC Grant NE/I019999/1, and GMD and KJH were supported by OSMOSIS NERC Grant NE/I019905/1.

## REFERENCES

*Geophys. Res. Lett.*,

**39**, L18605, doi:.

*Atlantic Ocean*. Vol. 3,

*Hydrographic Atlas of the World Ocean Circulation Experiment (WOCE)*, International WOCE Project Office. [Available online at http://www-pord.ucsd.edu/whp_atlas/atlantic_index.html.]

*Geophys. Res. Lett.*,

**39**, L14602, doi:.

*J. Geophys. Res.*,

**115**, C03017, doi:.

*Ocean Modeling in an Eddying Regime*,

*Geophys. Monogr.*, Vol. 177, Amer. Geophys. Union, 17–38.

## Footnotes

^{1}

Information on the PAP site observing system can be found at http://www.noc.soton.ac.uk/pap/index.php.

^{3}

We thank J. Callies for pointing out the analytical expression to this curve: [1/(2*π*^{2})]{[log*r*^{2} − log(*r* −1)^{2}]/[*r* − (1/2)]}.