The Contribution of Surface and Submesoscale Processes to Turbulence in the Open Ocean Surface Boundary Layer

The ocean surface boundary layer is a critical interface across which momentum, heat, and trace gases are exchanged between the oceans and atmosphere. Surface processes (winds, waves, and buoyancy forcing) are known to contribute significantly to fluxes within this layer. Recently, studies have suggested that submesoscale processes, which occur at small scales (0.1–10 km, hours to days) and therefore are not yet represented in most ocean models, may play critical roles in these turbulent exchanges. While observational support for such phenomena has been demonstrated in the vicinity of strong current systems and littoral regions, relatively few observations exist in the open‐ocean environment to warrant representation in Earth system models. We use novel observations and simulations to quantify the contributions of surface and submesoscale processes to turbulent kinetic energy (TKE) dissipation in the open‐ocean surface boundary layer. Our observations are derived from moorings in the North Atlantic, December 2012 to April 2013, and are complemented by atmospheric reanalysis. We develop a conceptual framework for dissipation rates due to surface and submesoscale processes. Using this framework and comparing with observed dissipation rates, we find that surface processes dominate TKE dissipation. A parameterization for symmetric instability is consistent with this result. We next employ simulations from an ocean front‐resolving model to reestablish that dissipation due to surface processes exceeds that of submesoscale processes by 1–2 orders of magnitude. Together, these results suggest submesoscale processes do not dramatically modify vertical TKE budgets, though such dynamics may be climatically important owing to their ability to remove energy from the ocean.


Introduction
Mechanisms that control the exchange of momentum, heat, and trace gases between the oceans and atmosphere are critical to Earth's climate. Surface processes, including breaking waves, buoyancy loss, current shear, and Langmuir motions, generate turbulence, and this mixes the water-column vertically (Belcher et al., 2012;Large et al., 1994). More recently, however, processes taking place at horizontal and temporal scales of 0.1-10 km and hours to days scales and not yet represented in climate models have been identified as influencing turbulence and buoyancy budgets within the ocean surface boundary layer (OSBL) (Boccaletti et al., 2007;Callies et al., 2016;Fox-Kemper et al., 2008Grisouard, 2018;Haine & Marshall, 1998;Thomas, 2005;Thomas & Lee, 2005;Taylor & Ferrari, 2009Thomas & Taylor, 2010). These have collectively been referred to as submesoscale processes (McWilliams, 2016;Thomas et al., 2008). Both (a) surface and (b) submesoscale processes believed to be the leading sources of turbulence in the OSBL are illustrated in Figure 1.
Much of our understanding of the energetics at the oceanic submesoscale is based on high-resolution numerical simulations (Boccaletti et al., 2007;Brannigan et al., 2015;Capet et al., 2008;Fox-Kemper et al., 2008Grisouard, 2018;Gula et al., 2014Gula et al., , 2016Haine & Marshall, 1998;Hamlington et al., 2014;Mahadevan, 2016;Smith et al., 2016;Taylor & Ferrari, 2009Thomas & Taylor, 2010). The chief reason is that present-day observing strategies do not sufficiently resolve submesoscale motions owing to their quickly Similarly, buoyancy loss at the ocean surface reduces vertical stratification and permits upright convection (i.e., gravitational instability). (b) Stirring and straining by the mesoscale eddy field generates pronounced lateral gradients in density. Winds oriented downfront (i.e., in the direction of the geostrophic shear at the front) transport dense water into areas of less dense (more buoyant) waters. Termed the Ekman buoyancy flux, B e , this flux reduces vertical stratification and admits symmetric instability (SI) within the front. This eventually leads to turbulence dissipation through Kelvin-Helmholtz instability (Taylor & Ferrari, 2009). Similarly, buoyancy loss at the ocean surface, B o , reduces vertical stratification and permits SI and enhanced turbulence at the front. Note that wind-forced and buoyancy-forced SI approximately occur at all depths within the ML but are dominant over buoyancy convection only for depths h ≤ |z| ≤ H. evolving nature; complex sampling strategies have therefore been devised that make use of ships, satellites, and autonomous instruments (Adams et al., 2017;D'Asaro et al., 2011;Shcherbina et al., 2015; see also http://www.uhrwerk-ozean.de/index.html.en) and, in some cases, new instruments have been developed (Novelli et al., 2017). While these approaches have improved our understanding of energetics at the oceanic submesoscale, they have thus far provided only short-term records sufficient for characterizing individual oceanic events. A notable exception is the analysis of repeat transects of a merchant vessel over the Gulf Stream (Callies et al., 2015) and analysis of repeat glider transects in the Southern Ocean (du Plessis et al., 2019). In almost all cases of documented turbulence, these studies were conducted in the vicinity of strong current systems where submesoscale dynamics are expected to be intense. An outstanding question is therefore, "What are the relative contributions of surface and submesoscale processes to energy dissipation in regions distant to strong current systems-i.e., in the open ocean?" In an effort to bridge this observation-model gap, an extensive field campaign was conducted in the North Atlantic (Figure 2a). Referred to as the Ocean Surface Mixing, Ocean Submesoscale Interaction Study (OSMOSIS), observations consisted of upper ocean temperature, salinity, and currents from nine moorings (50-500 m) and two ocean gliders (0-1,000 m) for an entire year (September 2012 to September 2013) Damerell et al., 2016;Erickson & Thompson, 2018;Thompson et al., 2016;Yu et al., 2019). Several other process-oriented studies exist, as well (Buckingham et al., 2017;Lucas et al., 2019;Rumyantseva et al., 2015). Important to the present study, the inner moorings were spaced approximately 2 Journal of Advances in Modeling Earth Systems 10.1029/2019MS001801 km from one another, while outer moorings were spaced 14 km from each other, permitting horizontal gradients to be calculated at fine-scale resolution (Figure 2b). An animation is provided online as supporting information (ms01.mp4) that depicts the horizontal scales resolved at the OSMOSIS mooring site and contrasts these with mesoscale motions, as seen by satellite altimetry. Additionally, the central mooring housed an acoustic instrument capable of measuring dissipative turbulence at approximately 50-m depth and that remained within the OSBL from December 2012 to April 2013. These data are complemented by reanalyses of winds, waves, and heat and precipitation fluxes, or hereafter buoyancy fluxes. Here, we use these data in concert with existing and in some cases new turbulence parameterizations to help understand the relative contributions of submesoscale processes to energy dissipation in the open-ocean environment.
The study is outlined as follows. In section 2, we present the data sets used in the study and outline the conceptual and analytical model for both surface-and submesoscale-forced turbulent dissipation. In section 3, we present the results and examine more closely the interplay between surface and submesoscale processes. We investigate implications for upper ocean energy budgets in section 4. Finally, in sections 5 and 6, we qualify these results and summarize our conclusions, respectively.

Data Sets 2.1.1. OSMOSIS Observations
The OSMOSIS moored array consisted of nine moorings: four inner, four outer, and a single central mooring. Inner moorings were spaced approximately 2 km apart, while outer moorings were spaced approximately 14 km apart ( Figure 2b). Additionally, inner moorings were more heavily instrumented, with paired conductivity-temperature-depth (CTD) sensors and acoustic current meters located every 50-100 m, starting at approximately 50-m depth. Finally, a 600-kHz acoustic Doppler current profiler (ADCP) mounted on the central mooring enabled us to estimate dissipative (i.e., meter-scale) turbulence for an entire year. This approach assumes a functional form for the structure function of velocities measured by the instrument and inverts a regression of the data to this functional form to obtain dissipation rate (supporting information) (Lucas et al., 2014;Wiles et al., 2006). The instrument was located at a nominal depth of 45 m.
Owing to seasonal variations in mixed layer (ML) depth and given the nominal depth of the instrument, we found that the dissipation instrument was not within the OSBL for the full seasonal cycle (i.e., September 2012 to September 2013). As dissipation within this boundary layer is the topic of this manuscript, we have restricted our analysis to time periods when the instrument resided within the ML: early December through late April. It is noteworthy that this is one of the few times a dissipation time series of such long duration (5 months) has been documented for the open-ocean environment. Unique to the present study, these data are complemented by lateral buoyancy gradients for the same time period. The processing of these moored measurements are described in detail below. Finally, while one can compute ML depth from the moorings, we have elected instead to use an estimate of ML depth from ocean gliders in order to be consistent with earlier studies (e.g., Buckingham et al., 2016;Damerell et al., 2016;Thompson et al., 2016, andYu et al., 2019).
The primary observations (lateral buoyancy gradients and dissipation rates) are derived from measurements made on all nine moorings. Currents were measured by Nortek acoustic current meters, were corrected for magnetic deviations, and were interpolated to an evenly spaced grid. CTDs were checked for biases/drifts by cross-calibrating these with casts made with the ship-CTD. Temperature (T), salinity (S), and pressure (P) were converted to potential density using the equation of state of seawater and buoyancy was computed from density as b = g(1 − ∕ o ), where g = 9.81 m/s 2 is the acceleration due to gravity and o = 1, 025 kg/m 3 is a reference density.
Lateral velocity and buoyancy gradients were estimated from first-order differences in measurements made from inner moorings. Unfortunately, the CTD mounted nearest to the ocean surface on the northeast inner mooring was damaged, preventing lateral buoyancy gradients across the inner moored domain until approximately 125-m depth. As this was not suitable for our study, gradients were instead computed within the inner mooring region using only northwest (NW), southwest (SW), and southeast (SE) inner moorings. Gradients within the outer mooring domain made use of all four outer moorings. For example, we obtained differences in buoyancy between NE and SW and between NW and SE mooring, rotating these gradients by a fixed amount in order to approximate gradients in eastern and northern directions. In both cases, we approximated ∇ h b from the vertical average of ∇ h b from near the surface (i.e., 50 m) to |z| = 140 m; this suffices during winter when frontal scales penetrate well below these depths.
The relative motion of moorings introduces uncertainty into lateral gradients. One way to handle such uncertainty is by modeling relative mooring motion . As this itself has assumptions (e.g., the distribution of underlying motion) and results in excessive bounds on the true value, we have opted for a different approach in this manuscript. The magnitude-squared coherence between horizontal velocities on inner moorings was shown to be high for periods exceeding the inertial period ( Figure 3). This suggests that all moorings move together in time for periods longer than the inertial period, while they are perturbed relative to one another for shorter periods. For example, estimates of this motion suggest horizontal displacements of 300-500 m are not uncommon during winter . Thus, we have smoothed measurements of buoyancy in time, here using a 24-hr filtering window. Gradients from outer moorings are less affected by such motion owing to the larger spacing between the moorings (supporting information). The animation mentioned earlier gives some measure of confidence in this approach (ms01.mp4), as illustrated currents were averaged over 24 hr. Additionally, we wished to be more consistent with theories of Ekman transport (see section 2.5.2), which requires spin up over time scales comparable to or greater than the inertial period, that is, 16 hr at this location. In summary, this approach produces results with less error and is more consistent with Ekman dynamics.

Ocean-Atmosphere Reanalysis: ERA5
While meteorological observations exist for a portion of this year-long deployment period, we have limited knowledge of the marine boundary layer except from that given by atmospheric and wave reanalysis. Here, we have opted to supplement ocean observations with reanalyses from the European Centre for Medium-Range Weather Forecasts referred to as ERA5 (https://www.ecmwf.int/en/ forecasts/datasets/archive-datasets/reanalysis-datasets/era5). These data are expected to provide improvements over existing reanalysis data (Dee et al., 2011) (https://www.ecmwf.int/en/newsletter/150/news/ era5-aids-forecast-performance-monitoring). The atmospheric products consist of hourly winds, waves, and ocean-atmosphere buoyancy fluxes (i.e., heat and precipitation fluxes) from which surface and Ekman buoyancy fluxes, B o and B e respectively, were estimated (Cronin & Sprintall, 2001;Stull, 1988;Thomas & Lee, 2005). These terms are defined below but briefly, B o , pertains to a loss of buoyancy from the ocean surface and that causes destabilization of the upper ocean, while B e pertains to the destratification that takes place as a result of Ekman-induced transport at an ocean front, moving dense waters into areas of less dense waters. (Of course, stratifying versions of the above fluxes occur, as well, and we account for these below.) Surface wind stress and buoyancy fluxes were subsequently smoothed with a 24-hr filter in order to provide consistency with the above ocean observations. This, too, has the benefit of being consistent with Ekman theory when compared to the hourly output. It is worth mentioning that winds averaged in this manner also compare favorably with those of a nearby meteorological station about 30 km away from the OSMO-SIS site (A. Forryan, personal communication, May 2015) but which were unfortunately unavailable for the presently examined period. (Data accessible from the Porcupine Abyssal Plain website, a EuroSITES observatory [http://projects.noc.ac.uk/pap/]) Finally, we also obtained historical reanalysis of wave forecasts from ERA5. These data are at coarser horizontal resolution but were used to parameterize dissipation due to both breaking waves and Langmuir turbulence.

Realistic, High-Resolution Ocean Simulations
Lastly, we employ high-resolution numerical simulations of the upper ocean in order to further establish the relative magnitudes of surface-and submesoscale-forced turbulence and to obtain relative magnitudes as a function of depth within the OSBL. Moreover, this is a step in the direction of better understanding and representing turbulent processes in climate-scale ocean models by showing what effect submesoscale processes might have if small-scale fronts were resolved, as well as what impact this might have on the extraction of energy from ocean currents. The model used in this analysis is referred to as NATL60 (CMJ-165).
NATL60 is an innovative experiment aimed at modeling the ocean down to kilometer-scale resolution as realistically as possible for the purpose of supporting the highly anticipated Surface Water and Ocean Topography (SWOT) satellite altimeter mission (Fu & Ferrari, 2008;Morrow et al., 2019). The vertical and horizontal grid resolution, as well as model parameters, have been specifically chosen to simulate oceanic submesoscales in the North Atlantic. The domain is bounded by the 26.5 • N parallel in the south and by 66.5 • N (polar circle) in the north. In ORCA configuration, the size of the mesh is proportional to the cosine of latitude. As the horizontal resolution at the equator is (1/60) • , the horizontal resolution varies from 1,600 m at the southern boundary to 900 m at the northern boundary. To guarantee representation of surface-intensified flows (e.g., submesoscale dynamics), the NATL60 model configuration uses 300 vertical levels. The vertical resolution is 1 m at the ocean surface and increases nonlinearly to 50 m in the deepest regions. The initial and open boundary conditions are built from the GLORYS2-V3 reanalysis provided by Mercator Ocean (Garric et al., 2017) (also https://www.mercator-ocean.fr/en/science-publications/glorys/). Boundary conditions at the ocean-atmosphere interface use the Drakkar Forcing Set DFS5.2 (https://sextant. ifremer.fr/record/c837b2f8-4152-41fd-8592-d4cd887d0b51/), which is based on ERA-interim reanalysis (Dee et al., 2011). The version of NATL60 used in our study is not forced by tides.
In our study, we seek to better understand dissipation rates in the open ocean resulting from surface and submesoscale processes. We obtained 1 year of high-resolution model output from the NATL60. Using ERA5 reanalysis, we computed the same forcing terms as from the observations. For example, we used the ERA5 data to provide winds, waves, and surface fluxes. From these, we calculated the dissipation due to winds, waves, and buoyancy-forced convection. We solved for both the surface buoyancy flux, B o , and Ekman-driven buoyancy flux, B e , at each location within this model (see equations (5) and (15)). The latter, B e , requires (1) knowledge of the OSBL thickness, H, (2) a term referred to as the convective layer (CL) depth, h (Taylor & Ferrari, 2010), and (3) the orientation of surface winds and buoyancy gradients (i.e., ocean fronts) in the model. All of these terms are described below. These calculations were made for a subdomain of the model centered on the OSMOSIS site but having surface area of 330 km × 330 km and calculations were made to a depth of 500 m. The ML depth is shallower than this depth at all times. Because of the volume of data, we used daily averaged density (or buoyancy) fields from the model, which is consistent with smoothing applied to the moored data. In summary, the model output is used to provide realistic frontal structure in the horizontal and vertical for submesoscale parameterizations-as well as realistic ML depths-while the conceptual scaling model presented below is used to estimate dissipation rates. That is, we harness the ability of the model to accurately depict frontal structure expected in the ocean and apply our conceptual framework to these gradients. We quantify the degree of realism of the NATL60 by comparing relevant parameters with the OSMOSIS observations. We also note that a recent comparison of buoyancy gradients in the NATL60 and those from ship-mounted thermosalinograph measurements displays good agreement across the North Atlantic basin (G. Serazin, personal communication, January 2019). Finally, velocity gradients in the model display good agreement with the OSMOSIS moorings, suggesting that frontal structures are well resolved at our location (supporting information). In summary, ERA5 winds, waves, and ocean-atmosphere fluxes serve as the majority of inputs to our parameterizations of the turbulent dissipation rate, while the model provides h, H, and lateral buoyancy gradients. The parameterizations are described in sections 2.4.1, 2.4.2, 2.4.3, and 2.5.2.

Spatial Resolving Capability of OSMOSIS and NATL60 Data Sets
The spacing of OSMOSIS moorings was specifically designed to resolve submesoscale fronts and eddies within the OSBL while simultaneously preventing entanglement of moorings over the year-long deployment. Similarly, as stated above, the model was designed to resolve submesoscale dynamics in anticipation for the SWOT satellite altimeter mission. While not the focus of this study, it is worth considering the resolving capability of these data sets to determine to what extent the observation and model data can resolve the relevant horizontal scale. The relevant horizontal scale is here argued to be the ML deformation radius.
Submesoscale fronts and eddies have been shown to be characterized by horizontal e-folding scales consistent with the ML baroclinic Rossby deformation radius. One can estimate this magnitude from eddy scales corresponding to the maximum growth rate of a modified Eady model applied to the OSBL (Boccaletti et al., 2007;Haine & Marshall, 1998;Samelson, 1993;Stone, 1966). (As shown by Boccaletti et al., 2007, the presence of a nonzero bottom boundary condition, w ≠ 0, to this model tends to slightly increase the lateral scale of resulting eddies.) Here, the deformation radius is given as ML = ⟨NH∕ ⟩, where N is the buoyancy frequency, f is the Coriolis parameter, and H is the thickness of the OSBL. Brackets denote averaging across the OSBL. Using typical wintertime values of ⟨N∕ ⟩ = 10 and H ∼ 100-200 m, one obtains ML = 1-2 km. Doubling this value, therefore, results in a representative diameter for submesoscale eddies immediately following their formation: 2-4 km. As the OSMOSIS inner moorings have a minimum separation distance of ∼2 km, we conclude that the inner moorings resolve these eddy scales.
For similar reasons, we argue that NATL60 is able to represent the relevant horizontal scale. That is, it is able to resolve eddy-like structures down to scales of ∼4 km at midlatitudes. The effective horizontal resolution of numerical ocean models is generally considered be about 5-7 times the grid size provided the discretization uses high-order numerical schemes (Sommer et al., 2018;Soufflet et al., 2016). As stated above, NATL60's model grid size varies from 0.9 km at the northern boundary to 1.6 km close to the southern boundary. NATL60 also uses a dense vertical grid with 300 levels which allows the model to better represent surface-intensified motions that are typically observed at these small scales. Given a grid spacing of 1 km at the site location, we therefore expect that the effective resolution of NATL60 is about 8 km (wavelength). However, rationalizing one can resolve two eddies of opposite polarity (e.g., sea surface temperature or sea surface height anomaly) within a single wavelength, we estimate that NATL60 model resolves eddy-like structures with diameters near 4 km. This is slightly larger than the diameter of submesoscale eddies observed during winter. Thus, relative to the OSMOSIS moorings, we expect NATL60 resolves scales larger than those resolved by the inner moorings but smaller than those resolved by the outer moorings. This is consistent with statistical comparisons of lateral buoyancy and velocity gradients at the OSMOSIS mooring site.
The aforementioned discussion concerns the horizontal resolving capability of these data sets. However, one must also consider their vertical resolving capability. It is well known that models and observations must both have aspect ratios, = z∕ x comparable to (N∕f) −1 in order to capture the necessary dynamics (e.g., Ménesguen et al., 2018). As mentioned, ⟨N∕ ⟩ = 10 is a representative winter stratification when averaged over the ML. However, increased values are present near within the ML near submesoscale fronts and eddies. Employing a more realistic value of N∕f = 50 and using ML = x, we obtain a requirement that data sets capable of resolving the submesoscale during winter at the OSMOSIS site must have vertical resolutions of z = x = ML ∕50, or 20-40 m. This requirement is met within the model but, unfortunately, not within the observations. (This resulted from financial limitations, as well as mooring design constraints. For example, moorings experience greater drag when increased numbers of sensors are placed on the cable owing to a need to compensate weight increases with a greater number of buoyant floats.) This latter finding motivates use of a depth-averaged lateral buoyancy gradient when estimating dissipation rates due to submesoscale processes from moorings (section 2.5).

Scaling the Dissipation Rate of Turbulent Kinetic Energy
There are two classes of forcing for turbulence within the OSBL: surface-forced turbulence and submesoscale-forced turbulence. While both processes are influenced by winds and buoyancy loss, submesoscale-forced turbulence extracts its energy from the energy stored in ocean fronts (D'Asaro et al., 2011;Thomas & Taylor, 2010), whereas surface-forced turbulence extracts its energy from the atmosphere-in the form of winds, waves, and buoyancy fluxes. Below, we examine these competing classes of forcing and use the OSMOSIS observation record to place these dissipation estimates into context.
Assuming steady state, neglecting advection by the mean flow and arguing that vertical gradients of mean quantities are considerably larger than horizontal gradients of the same (these assumptions can be rationalized by scaling the turbulent kinetic energy [TKE] equation; supporting information), the equation describing the evolution of TKE simplifies as Here, u = (u, v, w) = u + u ′ is the Eulerian velocity (overbar denoting a time mean and primes turbulent fluctuations), u s is the Stokes drift, b ′ = −g ′ ∕ o is the turbulent buoyancy fluctuation (g is gravity, ′ is the turbulent density fluctuation, and o is the background density), and e = (u The first term on the right-hand side represents shear production, the second production by Stokes drift, the third buoyancy production, the fourth redistribution of TKE across the vertical through transport and pressure-work, and the fifth dissipation rate. If we further partition the mean vertical shear of the Eulerian velocity into geostrophic and ageostrophic components (Thomas & Taylor, 2010;Thomas et al., 2013), we obtain an expression that can be used to relate the TKE dissipation rate, , in terms of TKE production by surface (SF) and submesoscale (SM) processes: ( This allows us to write ≈ SF + SM . This is simply a mathematical convenience, and we emphasize that SF and SM are actually TKE production terms used to scale the dissipation resulting from each class of forcing. Also note that the above expression assumes a scale separation exists between turbulent and mean quantities. Defining mean quantities by their time-averaged value over several minutes and turbulent fluctuations as deviations from these values, we see that all primed variables are ageostrophic in nature. We simply have excluded the subscript "a" for brevity. These are perturbations due to surface waves and/or internal waves, for example. These fluctuations can then correlate with each other and extract energy from the mean shear. For example, in the case of symmetric instability (SI), which falls under the last term in equation (2), velocity perturbations caused by internal waves extract energy from the geostrophic flow at ocean fronts, producing TKE. Similarly, current shear in the absence of lateral buoyancy gradients, which falls under the first term, creates perturbation velocities which correlate and produce TKE.
For geostrophically balanced flows, the vertical derivative of the horizontal flow is related to the lateral buoyancy gradient: where subscripts x and y denote partial differentiation.
In the work that follows, we estimate the magnitude of SM (termed geostrophic shear production, GSP) from a recent parameterization (Thomas & Taylor, 2010) in combination with lateral buoyancy gradients between moorings and model grid cells.
In the above expression, we have neglected the turbulent transport term, arguing that internal waves, while present within the observations, do not systematically redistribute energy vertically. Also, we assume that processes do not interact with one another so that, for example, submesoscale turbulence does not substantially modify Langmuir turbulence, or vice versa. We revisit this assumption in section 5. Each of the processes thought to give rise to enhanced turbulence within the OSBL is illustrated in Figure 1, and we explore their contributions to turbulent dissipation rate in each of the sections below, using a combination of existing and new parameterizations and synthesizing them within a single framework.

Scaling Dissipation Rates Due to Surface Processes
The surface-forced turbulence is itself driven by three processes: wind, waves, and buoyancy loss at the air-sea interface. Each process provides a source in the TKE equation, which helps balance dissipation (cf. equation (2)). Following Belcher et al. (2012), we write this as a linear combination of dissipative terms: where s is the dissipation due to current shear, w dissipation due to wave-induced turbulence which can include Langmuir turbulence, and b is the dissipation due to surface buoyancy forcing. We ignore turbulence generated by near inertial shear (Alford et al., 2016), as well as wind shear alignment at the base of the OSBL (Brannigan et al., 2013;Burchard & Rippeth, 2009). We also ignore any mixing produced by marine organisms (Dewar et al., 2006;Rousseau et al., 2010).
The expression for the surface forcing term, SF , is somewhat complicated by the fact that, in our data set, the wind and wave-forced turbulence have similar magnitudes. That is, estimates of the ratio of the production of turbulence by current shear to the production of turbulence by Stokes shear (as summarized by the turbulent Langmuir number) yield values close to 0.3 for the duration of our record (supporting information), suggesting that while Langmuir turbulence likely dominates, it is difficult to separate these two sources of turbulence. For this reason, here we present two complementary models of surface-forced turbulence, one with and without Langmuir turbulence. For further information about the relationship between these types of surface forcing (i.e., due to buoyancy fluxes, wind, and waves) we direct the reader to a regime diagram reflective of Model B of our study (Belcher et al., 2012, Figure 3). Common to both models is the surface buoyancy forcing, which we present first. Each of the processes giving rise to surface-forced turbulence is illustrated in Figure 1a.

Buoyancy-Forced Convection: b
Traditional models of buoyancy-forced convection represent the production of turbulence due to surface buoyancy loss as constant over the depth of the ML: where B o is the surface buoyancy flux, is an empirically determined constant, and angled brackets denote a vertical average. Typical values for are 0.4-0.7 (Imberger, 1985;Shay & Gregg, 1985). However, an equally valid model of convection is one that predicts both upright convection and SI to a depth h and only SI below to a depth H (Taylor & Ferrari, 2010). The upper layer is referred to as the CL and is where upright convection dominates. The reason for adopting this more general framework for convective turbulence will become clearer when we discuss submesoscale-forced turbulence, below.
We model the turbulence generated by buoyancy production as a linear function of depth within the CL: We then approximate the dissipation rate due to buoyancy loss at the surface as A scale for the CL thickness was obtained by solving a fourth-order polynomial in z ′ = h∕H (Taylor & Ferrari, 2010;Thomas et al., 2013) (Appendix A). In all of our analysis, we approximated the OSBL thickness, H, by the ML depth, defined as the depth at which the potential density exceeds its value at 10 m by 0.03 kg/m 3 (de Boyer Montégut et al., 2004). This has the distinct advantage of being consistent with our representation of surface-and submesoscale-forced turbulence in the models below. (We also considered a definition of OSBL thickness based on the gradient in potential vorticity [PV] but found it made little difference in our results-at least in winter.) The surface buoyancy flux from the ocean to the atmosphere B o was estimated in the usual manner, namely, Here, o is the upper ocean density, S o is the surface salinity, g is the acceleration due to gravity, is the thermal expansion coefficient, is the haline contraction coefficient, c p is the specific heat capacity of seawater, Q o is the net surface heat flux from the ocean to the atmosphere, E is the evaporation rate, and P is the precipitation rate (Cronin & Sprintall, 2001;Stull, 1988). Q o was estimated as the sum of sensible and latent heat flux, as well as shortwave minus longwave radiation. Equation (5)   We also display the observation depth, |z o | (i.e., the depth corresponding to the dissipation rate, o ), whose variability reflects displacements of the central mooring. Finally, as described below, two turbulent regimes exist in which the wind-and wave-forced layer is divided: (a) breaking wave and, beneath this, (b) shear-or Langmuir-forced regimes. The transitions between these layers are z s and z L , respectively, and are also shown in this figure. During the observation window, the depth of the dissipation-recording instrument is found completely within the ML and comfortably . The observation depth corresponds to the depth of the high-frequency ADCP while depth scales include the convective layer (CL) depth, h, mixed layer (ML) depth, H, transition depth between breaking and current shear layers, z s , and transition depth between breaking and Langmuir turbulence, z L . The observation record of the observed turbulence, obs , is expected to include contributions from both surface (i.e., buoyancy-forced convection, winds, and waves) and submesoscale processes (here parameterized as only due to symmetric instability).
within the CL. Furthermore, the instrument sits almost at the transition between the breaking-dominated and current/Langmuir-dominated depths. From this point of view, we can therefore expect to see contributions from all three surface processes within our observed dissipation time series. Submesoscale turbulence extends throughout the water column to the base of the ML and therefore should also be contained in the observed dissipation rate, obs . As this tends to occur when the CL depth is smaller than the ML depth, or h∕H < 1, it is much more infrequent. How each of these sources of turbulence scale will be determined by the expressions given below.

Winds and Waves: s + w Without Langmuir (Model A)
The rate of dissipation due to winds and waves is here represented by two components: wave breaking and current shear. According to scale analysis of the upper ocean (Terray et al., 1996), dissipation below the ocean surface can nominally be represented by three layers. In the near-surface layer, which is on the order of one significant wave height, H s , dissipation is dominated by the injection of TKE by the surface wave field (Thomson, 2012). Because our moored observations do not extend to the ocean surface, we ignore this near-surface layer. Below this, but still confined near the sea surface, H s < |z| < |z s |, the dissipation of energy is dominated by wave breaking. Lastly, at greater depths, the dissipation is characteristic of that found within the bulk of the OSBL. Hence, we write the dissipation due to winds and waves as Here, w * b = ( u 2 * wc ) 1∕3 serves as the velocity scale for the turbulence generated by wave breaking [u *w is the waterside friction velocity,c = 0.1c p is the effective wave speed, c p = gT∕(2 ) is the phase speed associated with the peak waves, and T is their period], f b defines its functional form with depth, = 0.4 is the Von Kármán constant, and |z s | is the transition depth.
The dependency of dissipation rate with depth within the breaking layer is given by an inverse square law (Terray et al., 1996), namely, where again H s is the significant wave height. The transition depth between the two layers is obtained by matching dissipation in the breaking-dominated layer with that in the lower later, giving where u *a ∕u *w ≈ 31.6 is the ratio of airside to the waterside friction velocities and c p ∕u *a is the wave age (i.e., a measure of wave speed relative to wind speed). We considered a wave-age-dependent scaling for the dissipation rate (Wang & Huang, 2004), but it made little difference in the solution as the wave age in the reanalysis data exceeds 20 for most of the observation period (not shown).

Winds and Waves: s + w With Langmuir (Model B)
This surface-forced model of winds and waves closely follows that of "model A" above but with one important difference. Here, the rate of dissipation due to current shear is replaced with that resulting from Langmuir turbulence. On the basis of an ensemble of large eddy simulations (LES) (Belcher et al., 2012;Grant & Belcher, 2009), we define Here, w * L = ( u 2 * w u s0 ) 1∕3 is the velocity scale associated with Langmuir turbulence [u s0 is the surface Stokes drift velocity], H serves as the depth scale and has been found to be equal to the ML depth (Belcher et al., 2012;Grant & Belcher, 2009), and its functional form with depth, f L , has been obtained from LES of Langmuir turbulence ( Figure 5): Again, we match dissipation rates to obtain the appropriate transition depth, giving which has the following relevant (i.e., z L < 0) solution: We note that the ratio of transition depths, z s ∕z L , is approximately 1.5-2.0 over most of the observational record, supporting use of the Langmuir-based model. Nevertheless, we include both models for completion.

Scaling Dissipation Rates Due to Submesoscale Processes
The suite of submesoscale processes that give rise to enhanced dissipation is here considered for the open-ocean environment. The literature on submesoscale processes in the ocean has become vast, and here we provide only a simple introduction to that subject. Complete descriptions are found in Thomas et al. (2008) and McWilliams (2016). The representation of turbulence generated by submesoscale processes is discussed in sections 2.5.1.1 and 2.5.2, such that the reader familiar with the oceanic submesoscale regime can skip to these sections without much loss of content.

A Primer to the Oceanic Submesoscale Regime
At large horizontal scales in the ocean (e.g., >30 km at midlatitudes), fluid dynamics are well described by quasi-geostrophic theory. Here, flows are approximately in hydrostatic and geostrophic balance, and-by consequence-only minimal vertical motion is allowed. By contrast, submesoscale processes reside at smaller scales and within a dynamical regime in which both (a) vertical and horizontal shears are enhanced and (b) stratification low enough that increased vertical motion is possible. In particular, submesoscale processes are characterized by order-one gradient Rossby and gradient Richardson numbers, Ro = ∕f and Ri = N 2 ∕| z u| 2 , respectively. At midlatitudes and for flows with velocities near 0.1-1.0 m/s, typical horizontal scales are 0.1-10 km and these flows evolve with time scales near the inertial period (i.e., hours to days). These scales are larger near the equator and smaller near the poles. Because these conditions of enhanced shear and reduced stratification are met within boundary layers, we expect an active submesoscale within the OSBL. (For similar reasons, we expect submesoscale processes to be active within the bottom boundary layer; Callies, 2018;Gula et al., 2016;.) The connection of flows characterized by order-one gradient Rossby and Richardson numbers to the Ertel PV has a long history. It can be traced to conversations in the atmospheric community in which scientists examined the stability of vortices to atmospheric perturbations in order to explain motions encountered at atmospheric fronts (e.g., papers by Fjortoff, Eliassen, and Charney). While these instabilities were often considered in terms of static stability and angular momentum, they can also be conveniently summarized by the Ertel PV (Hoskins, 1974): instabilities can occur when fq < 0, where f is the Coriolis parameter and q = (2 + ∇ × u) · ∇b is the Ertel PV (Ertel, 1942). Dynamically, what this implies is that instabilities are possible when the absolute angular momentum of the fluid (as summarized by PV) approaches a marginally stable state (i.e., q = 0), or even changes sign from its background value (q < 0 in the Northern Hemisphere, q > 0 in the Southern Hemisphere).
Away from the equator, one can approximate the PV as q ≈ (f + )N 2 + h · ∇ h b, where h is the horizontal component of relative vorticity and ∇ h b is the vertical component. For flows in thermal wind balance, one can then write q ≈ (f+ )N 2 −f| z u| 2 . Assuming N 2 > 0 and f ≠ 0, we can then divide the instability criterion by f 2 N 2 to obtain a nondimensional criterion: 1 + Ro − Ri −1 < 0, where Ro and Ri are defined above. Since the last term is negative definite, we see that PV is reduced at fronts and, in the presence of forcing, is more susceptible to dissipative instabilities.
While a complete discussion of the instability criterion, fq < 0, and its implications is beyond the scope of the present study, we nonetheless highlight several outcomes of the above instability criterion. Three instability types are possible when fq < 0: (a) gravitational instability (GI) (N 2 < 0), (b) inertial or centrifugal instability (CI) (f + < 0 or Ro < −1), and (c) SI. SI is a unique form of instability present in ocean and atmospheric fronts in which fluid parcels move unbounded along lines inclined to the horizontal but the flow can be shown to be neither gravitationally unstable nor inertially unstable. Since stratification is positive, the former expression applies and we see that the gradient Richardson number (i.e., enhanced vertical shear and low stratification) plays the primary role in reducing PV.
An additional instability not mentioned above but that nevertheless occurs in the OSBL is baroclinic instability (BCI). Here, a meridional (i.e., cross-frontal) motion results in stretching or compression of fluid parcels and, because vorticity is largely conserved during this process, vorticity is generated and produces eddies. The resulting disturbances scale approximately with the ML baroclinic Rossby radius of deformation, NH∕f, where N∕f is normalized stratification and H is the depth scale of the OSBL (here, taken to be the ML depth). Typical values result in eddies with scales on the order of kilometers. Hereafter, we refer to this instability type as submesoscale BCI.

What Are the Dominant Submesoscale Processes Producing and Dissipating Energy
Given the suite of instabilities described above, we must now assess how to best represent their effect on turbulence at the OSMOSIS site. We have accounted for GI in our parameterization of buoyancy-forced convection. We are also able to exclude BCI as a mechanism for dissipation in our environment since BCI generally fluxes energy to the large scale (inverse cascade) rather than to smaller scales where dissipation takes place. (Note, however, there is evidence of a forward cascade of energy-i.e., to dissipative scales-occurring during submesoscale BCI in high-resolution numerical simulations; Molemaker et al., 2010;Skyllingstad & Samelson, 2012.) We argue that CI can be neglected since vorticity values at our moored site are much smaller than f in our domain. The standard deviation of the gradient Rossby number as measured from the inner moorings, for example, was found to be 0.3-0.5 over the observation period, even during winter when vorticity is elevated owing to the production of eddies by submesoscale BCI . More importantly, using the method of Thomas et al. (2013), we have categorized the instabilities and find that, of the observed instabilities, less than 2% have Rossby and Richardson numbers consistent with CI, while >40% have those values more consistent with SI (supporting information). Such methods depend upon thermal wind balance and assume the state of Rossby and Richardson numbers observed at a given time period are indicative of the onset of instability, which may not be true. Thus, we caution that we may simply be seeing the result of a previous instability pushing Rossby and Richardson numbers to more stable values. Additionally, it is possible that a hybrid CI-SI more akin to those examined by Grisouard (2018) would be observed-likely places include intense mean flow conditions such as the Antarctic Circumpolar Current (ACC). However, as this hybrid instability, CI/SI, has not yet been parameterized and since realistic ocean fronts are almost always characterized by lateral buoyancy gradients (Eliassen, 1983), we argue that dissipation due to CI will be represented by our parameterization in some measure (given below).
In summary, while we cannot be certain that we are accounting for all the submesoscale-forced turbulence, we elect to focus on SI in its pure form and proceed to examine the dissipation rates expected from such a process. Below, we use the symbol * to denote dissipation rate due to SI.

SI: *
SI is a generalization of CI that occurs at ocean fronts. Here, the fluid is neither gravitationally unstable nor centrifugally unstable and yet results in instability. In its purest form, SI results in unconstrained parcel motion that takes place along isopycnals and is, therefore, slanted with respect to the horizontal. This motion is eventually arrested by Kelvin-Helmholtz instability and converted into turbulence at fine horizontal and vertical scales (Taylor & Ferrari, 2009). The term "SI" comes the atmospheric community in which meteorologists were seeking to understand motions of fluid parcels within vortices. For simplicity, they considered azimuthal symmetry, examining only velocity perturbations in the radial direction, and the term SI was born. The instability therefore has a logical name given its origin but it does confuse scientists today.
There are two types of forced SI encountered in the open ocean: (a) wind-forced SI and (b) buoyancy-forced SI (cf. Figure 1b). Both occur at ocean fronts when the PV (Ertel, 1942) of a fluid approaches a marginally stable state, eventually taking on sign opposite to the Coriolis parameter (Hoskins, 1974). Also, see section 2.5.2, above. In the first case, winds blowing over the ocean surface cause a lateral transport of fluid that is to the right of the prevailing wind direction in the Northern Hemisphere and left of the wind direction in the Southern Hemisphere (Ekman, 1905). In the presence of a horizontal buoyancy gradient, ∇ h b, and if the surface winds are oriented with a component in the direction of the geostrophic shear, this transport of fluid has the potential to force dense water into regions of light water, resulting in an effective vertical buoyancy flux (Thomas, 2005;Thomas & Lee, 2005). This flux has been referred to as an Ekman-driven buoyancy flux, B e , and results in SI within the front. Similarly, in the second case (i.e., buoyancy-forced SI), a destabilizing surface buoyancy flux, B o , reduces upper ocean stratification. This in turn reduces the gradient Richardson number, Ri, in the water column and-combined with the already existing vertical shear of a front-can force SI to occur at depths below the upright convection previously described in section 2.4.1. That is, it is analogous to wind-forced SI, but here the decreased stratification is catalyzed by surface buoyancy loss rather than Ekman buoyancy flux via wind-front alignment. Of note, the source of energy for these two processes is the geostrophic shear within the front rather than the atmosphere.
Adopting the aforementioned generalized framework for convection and utilizing output from high-resolution numerical simulations of Ekman and surface buoyancy forcing, Taylor and Ferrari (2010), Thomas and Taylor (2010), and Thomas et al. (2013) proposed that the last term in equation (2), termed

10.1029/2019MS001801
GSP, could be sufficiently approximated as where, again, B e and B o are the Ekman and surface buoyancy fluxes, respectively, h and H have previously been defined, and w ′ b ′ is a linear function of depth within the CL (cf. equation (3)). Finally, the parameterization for the dissipation rate due to submesoscale processes, SM ≈ * , is essentially given by equation (13) subject to a constraint on the integral of shear production and buoyancy flux over the OSBL (supporting information). In summary, for the open-ocean, where dissipation due to other submesoscale instabilities (i.e., CI) are less frequent, we write This can further be partitioned into turbulence production due to wind-and buoyancy-forced SI by setting B o = 0 and B e = 0, respectively (supporting information).
In the above expressions, B o is the surface buoyancy flux and has been defined above. The Ekman-driven buoyancy flux has been defined as (D'Asaro et al., 2011;Thomas, 2005;Thomas & Lee, 2005) where M e = ( × f) ∕( o 2 ) is the Ekman transport, is the surface wind stress, ∇ h b is the horizontal buoyancy gradient across the front, and f =k is a vector having magnitude equal to the Coriolis parameter and oriented positive upward. As with B o , positive values of B e correspond to destabilizing conditions (e.g., heat flux out of the ocean).
We estimated B e from the along-front wind direction by using ERA5 data, multiplying (in an inner-product sense) the wind stress vector and a vector aligned with the downfront direction (i.e., more buoyant waters to the right). We have elected not to represent restratifying events in terms of their effect on turbulence except that we have "turned off" the turbulence dissipation in our parametization if the sum B o + B e is negative, consistent with equation (14).
It is worth noting that almost all of the TKE production during SI arises from shear production; very little turbulence is produced by the buoyancy flux, w ′ b ′ (L. Thomas, personal communication, October 2018), except in more sophisticated numerical simulations (Arobone & Sarkar, 2015;Skyllingstad et al., 2017). In fact, this is why the term SI is preferred over slantwise convection. This is not to say that restratification does not occur indirectly, however, since the loss of geostrophic shear results in a slumping of the front. But, in terms of TKE production, it is reasonable to neglect buoyancy flux produced by SI. This is the reason such a term is absent from equation (2). Figure 6 illustrates observed and parameterized dissipation rates due to surface and submesoscale processes at the OSMOSIS mooring site. In particular, we compare (a) our observed dissipation rates obs with those resulting from parameterizations of (b) buoyancy-forced convection b , (c) winds and waves s + w (two models), and (d) submesoscale-forced turbulence, * , here represented as due solely to SI. The reason two dissipation rates are displayed for the submesoscale portion, * I and * O , is owed to the fact that lateral buoyancy gradients were estimated in two different ways and these then feed into the parameterization for * . Finally, in (e), we display the sum of surface-and submesoscale-forced dissipation rates which we refer to as "predicted" dissipation rates and denote this quantity by = b + s + w + * .

Surface-and Submesoscale-Forced Turbulence in the Observations
First, note that there is good agreement between observed and predicted dissipation rates (cf. Figure 6e). This is further documented in the scatterplot/histogram of the two time series (Figure 7), which we describe below. This provides some confidence in our ability to model the observed dissipation rates. Second, surface processes are of leading order in terms of their effect on turbulent dissipation rates. In particular, and smoothed̃o bs using a 24-hr filter (solid line); these data are replicated in each panel for comparison; (b) dissipation rates due to buoyancy-force (upright) convection b ; (c) dissipation rate due to winds and waves ( s + w ) using a model with and without Langmuir turbulence (see text); (d) dissipation rate due to symmetric instability (SI) * ; here, we use buoyancy gradients estimated from both inner (blue) and outer (orange) moorings; (e) summed dissipation rate determined as the superposition of surface-and submesoscale-forced turbulence: b + s + w + * , where the contribution from wind and waves is given by the Langmuir-forced model. In the right-hand column, PDFs are displayed along with the most probable value (color coded to match that of their distributions), assuming a log-normal distribution for each variable; where three variables are displayed, a square denotes the third variable's most probable value.  Figure 6). Black dots indicate observation-prediction pairs falling in each bin while colors indicate their likelihood. Here, we display that prediction resulting from use of buoyancy gradients from the outer moorings in the determination of the CL depth, h, and the parameterization for SI; both inner and outer moorings yield nearly the same histogram owing to the dominance of surface-forced turbulence over submesoscale-forced turbulence.
winds and waves appear to drive most of the dissipation estimated at the observation depths of 40-to 50-m depth (cf. Figure 6c). So long as the parameterizations for surface-forced dissipation are correct, this effectively bounds the loss of energy by submesoscale processes. Lastly, the relationship between surfaceand submesoscale-forced turbulence holds even for periods when submesoscale-forced turbulence is most pronounced (e.g., February and April 2013); that is, surface forcing remains the leading-order source of turbulent dissipation rates. That the dissipation due to SI falls below the observed dissipation rates demonstrates that the parameterization is at least consistent with the observation record. In summary, from these unique observations, we conclude that surface processes dominate TKE budgets at this open-ocean site in the North Atlantic. Figure 7 displays a two-dimensional histogram of observed and predicted dissipation rates, where again predicted dissipation rates correspond to the sum of SF and SM (i.e., b , s , w , and * ). Here, we have represented s + w using Model B (i.e., with Langmuir turbulence), which we expect is more probable than Model A (supporting information). In any case, both result in similar dissipation rates. For each dissipation rate, a slight bias toward higher predicted dissipation rates is revealed. This bias could relate to simplifications made in deriving equation (2), or it could result from errors in the atmospheric reanalysis product, a subject beyond the present study. Nevertheless, the large fluctuations in obs nearly match our predicted dissipation rates (cf. Figure 6e), which together with Figure 6c suggest winds and waves drive the majority of the observed turbulence.
To better understand what specifically drives the observed dissipation rate, obs , we examined the forcing inputs to several of these parameterizations, focusing on the inputs to submesoscale-forced dissipation rates; drivers for the dissipation due to wind, wave, and buoyancy-forced convection can be inferred from the winds and B o , respectively. Figure 8 depicts these forcings (e.g., lateral buoyancy gradients, winds, B o , and B e ), as well as the predicted dissipation rates, * . Again, these are plotted alongside observed dissipation rates for comparison. Figure 9 magnifies a portion of this plot but for a time period in April 2013.
Several features are notable from these plots. First, few fronts of appreciable magnitude pass through the moored domain. This is seen both in the inner and outer mooring gradients (cf. Figures 8a and 8b). Second, though the wind stress is considerable, with magnitudes reaching as high as 0.6-0.7 N/m 2 (higher when examining the hourly winds), these winds coincide with strong lateral buoyancy gradients (i.e., fronts) and produce considerable buoyancy fluxes only a few times. That is, B e exceeds 5 × 10 −8 m 2 /s 3 only three times during the 5-month period. Note, this frequency is raised when examining B e for the inner moorings, but,   Figure 8 but magnified for 5 to 13 April 2013: (a, b) lateral buoyancy gradients from the inner and outer moorings, respectively, along with CL and ML depths, h and H; (c) mean of (a) and (b) between |z| = 50 and 140 m; (d) hourly wind stress (ERA5) and that obtained by smoothing over a 24-hr period; (e) surface buoyancy flux, B o , and Ekman buoyancy flux, B e (from inner and outer moorings); (f) observed and predicted dissipation rates, with the latter predicted from submeoscale forcing and surface forcing, respectively. In all panels, a region of enhanced dissipation levels above that predicted by winds, waves, and buoyancy forcing is highlighted in shaded boxes. As this coincides with a period in which h diverges from H, and as B e from inner moorings is pronounced, it is possible that SI is taking place. Figure 10. Comparison of mixed layer (ML) depths from the NATL60 high-resolution (1-km) ocean model and at the OSMOSIS mooring site (September 2012 to September 2013). The blue line depicts ML depths derived from OSMOSIS glider observations, while the red line and shaded region denote the NATL60 mean and ±1 standard deviation from the mean. Both were estimated as the depth at which the potential density exceeds its value at 10 m by 0.03 kg/m 3 (de Boyer Montégut et al., 2004) and OSMOSIS measurements were derived from ocean gliders that measured to within 10 m of the ocean surface  even still, only one of these events is associated with values > 1 × 10 −7 m 2 /s 3 (i.e., April 2013). Similarly, surface buoyancy fluxes B o of note are those above 2×10 −7 m 2 /s 3 . These instances of forcing then cascade down to the calculation for dissipation rates. At the observed depth, |z o |, dissipation levels are never found above 1 × 10 −7 m 2 /s 3 for outer moorings and only three times from inner mooring (not shown). It is noteworthy that, in one of the above cases at least, the sign of the Ekman buoyancy flux is opposite (up-front winds) to that of the surface buoyancy flux, reducing the dissipation. This is consistent with our parameterization for dissipation rate from the combined effects of B o and B e on SI.
It is instructive to reexamine some of these inputs, focusing on a period when * is pronounced. One such example is shown in Figure 9. Here, we display a graphic similar to that found in Figure 8 but magnified for the period 5 to 13 April 2013 (slight differences exist where we depict quantities from both inner and outer moorings). While we do not report observations of SI from the moorings (this is the topic of a separate study), Figure 11. Histograms of the magnitude of lateral buoyancy gradient, |∇ h b|, from the NATL60 high-resolution (1-km) ocean model and from the OSMOSIS outer moorings (December 2012 to April 2013).
we nonetheless indicate the potential for SI to elevate dissipation levels above that predicted by winds, waves, and buoyancy forcing (i.e., b + s + w ).

Surface-and Submesoscale-Forced Turbulence Implied by Realistic, High-Resolution Simulations and ERA5 Reanalysis
The observations presented above are necessarily finite in time and space. The moorings measure dissipation and lateral gradients of a 15-km× 15-km × 500-m volume of ocean in the North Atlantic for 5 months of the year. Thus, it is challenging to draw more widespread conclusions about the relative magnitudes of surface-and submesoscale-forced turbulence from such observations, particularly for the submesoscale-forced turbulence contributions owing to intermittent nature of ocean fronts transiting through the moored domain. Recently, however, the ocean modeling community has made significant strides in modeling ocean dynamics at small horizontal scale (e.g., 1 km) but for large fractions of the ocean (e.g., Gula et al., 2014) [also Le Sommer et al., "NATL60: A North Atlantic ocean circulation model data set based on NEMO for preparing SWOT altimeter mission." In preparation for Geophysical Model Development]. Assuming that oceanic frontal structure and ML depths in these models are reasonable, it is useful to consider the contributions of surface-and submesoscale-forced turbulence implied by such an ocean.
Before doing so, we compare ML depths and lateral buoyancy gradients in the model simulations and within the OSMOSIS observations. Figure 10 displays model and observed ML depths, H, as a function of September 2012 to September 2013 period. The agreement is good for the December-April period examined during this study, though the comparison appears to be less favorable during other portions of the year. Similarly, Figure 11 displays histograms of model and observed buoyancy gradients, |∇ h b|. A good rule of thumb is that ocean models can resolve features with horizontal scale approximately 8 times the grid resolution (Soufflet et al., 2016). Using this metric, this implies that the model resolves features similar to the scales resolved by the outer moorings. Indeed, we find remarkably good agreement between these two distributions, with peak values occurring at approximately the same value (i.e., 7 × 10 −8 s −2 ) and spread in values fairly close-although we do see slightly larger values in NATL60, as evidenced by the greater likelihood of observation in the tail of the distribution. This, in concert with a similar comparison of NATL60 and OSMOSIS velocity gradients (not shown), leads us to conclude that NATL60 is reasonably reproducing the observed gradients but with slightly higher values than those observed during winter at the OSMOSIS site. Figure 12 displays the domain-averaged turbulent dissipation rates due to surface and submesoscale processes, SF and SM , 2012-2013. Here, dissipation rates were obtained from applying our conceptual framework (idealized scaling model) to winds, waves, and buoyancy fluxes from ERA5 in combination with frontal gradients furnished by the NATL60 model. We also display in this figure both the CL and ML depths. Solving for the CL depth is computationally nontrivial even for this modest model domain size, restricting our region to that displayed in the solid black box of Figure 2a. Additionally, we note that surface-forced dissipation rates extend slightly beneath the domain-averaged ML depth, which results because there is significant spread in our ML depths, particularly during winter as a result of fronts and eddies that populate the simulations. We attribute this to the seasonality of submesoscale turbulence (Brannigan et al., 2015;Callies et al., 2015;Mensa et al., 2013).
Consistent with the observations, we see that TKE dissipation is dominated by surface-forced turbulence at all depths. Moreover, there is an onset of dissipation in late autumn and early winter that can be attributed to submesoscale turbulence, indicative of destratification by positive surface buoyancy fluxes and the onset of SI (Taylor & Ferrari, 2010;Thompson et al., 2016). (Median values depict similar characteristics as mean values but with a reduced amount of energy being dissipated by SI in median fields, perhaps revealing its infrequent or localized nature.) Canonical profiles of turbulent dissipation rate can also be obtained by first normalizing depths, |z|, by the local ML depth and then averaging in time and space. In Figure 13, we illustrate the winter-averaged (December-April) profiles of turbulent dissipation rate. For reference, and to illustrate the amount of scatter, we display the observed dissipation rates, as well as their mean. One observes several features within this plot. First, dissipation rates due to surface processes, SF , are 2 orders of magnitude greater than those of submesoscale processes. Moreover, the contribution of wind and waves, s + w , account for almost all of the dissipation, and for all depths. These two features are consistent with the moored observations. Second, dissipation profiles with and without Langmuir are slightly different but of the same magnitude. This may be why we find it difficult to distinguish between dissipation rates forced by Langmuir and current shear in our observations. (For reference, the average normalized depth of the observed dissipation record is |z|∕H ∼ 0.3, although there is substantial variation in this value.) Second, dissipation due to buoyancy-forced convection ( b ) exceeds that due to submesoscale processes above a nominal depth, say |z|∕H ∼ 0.45, suggesting that SI should take place below this. Indeed, this is observed, as dissipation due to submesoscale processes ( * ) becomes more important. It is also interesting to note that dissipation due to *w exceeds that due to *b for most depths (i.e., |z|∕H < 0.65), but below this the opposite is true. Whether this is intended as part of the underlying parameterization (Thomas & Taylor, 2010) is not known.

Discussion: Implications for Upper Ocean Energy Budgets
At the OSMOSIS site, surface processes appear to dominate over submesoscale processes in terms of their effect on turbulence dissipation rates. Here, we found that surface processes result in dissipation rates approximately 1-2 orders of magnitude higher than those of submesoscale processes. This was additionally verified using the NATL60 lateral buoyancy gradients and ML depths in concert with ocean-atmosphere reanalysis. Therefore, to the extent that the OSMOSIS site is representative of the open ocean-that is, away Figure 12. The mean energy dissipation rate due to (a) surface and (b) submesoscale processes in the North Atlantic. These time series correspond to domain-averaged turbulent dissipation rates obtained from applying the conceptual framework (i.e., scaling model) of section 2 to winds, waves, and buoyancy fluxes from ERA5 in combination with frontal gradients and ML depths provided by the 1-km ocean model (NATL60). Computations from the model were made in the exact same manner as was done for the observations, except that statistics were obtained over large domain size (330 km× 330 km × 500 m) and period (September 2012 to October 2013).
from coasts and strong current systems-this implies that neglecting submesoscale-forced turbulence in vertical budgets of TKE will not have an adverse effect on TKE budgets in Earth system models (ESM).
It is worth noting that a second-order turbulence closure model (Umlauf & Burchard, 2003) forced by the ERA forcing reproduces much of the observed variability in dissipation rates at the OSMOSIS site (not shown). Dissipation rates are biased slightly low compared to the observations, but the agreement is good. The reason such an approach was not taken, here, is that the results depend upon the particulars of its turbulence closure scheme. Instead, we felt the conceptual scaling model presented here was more transparent. Similarly, a recent study comparing several mixing models [Damerell et al. 2019, "A comparison of five surface ML models with a year of observations in the North Atlantic", in review], all of which neglect submesoscale processes, revealed that ML depths are well represented by surface processes. The fact that these models do not implement submesoscale-forced motions and yet still represent the recorded variability is therefore consistent with our results.
Exceptions to these results will be found within strong current systems. Indeed, western boundary currents such as the Kuroshio and Gulf Stream (D'Asaro et al., 2011;Thomas et al., 2013), equatorial regions (Marchesiello et al., 2011;Holmes et al., 2013), and ACC of the Southern Ocean (Adams et al., 2017; du Plessis et al., Figure 13. Winter-averaged profiles (December-April) of turbulent dissipation rate within the OSBL as functions of normalized depth, |z|∕H. Here, we display both surface-and submesoscale-forced turbulence dissipation rates, as well as their forcing terms: b + s + w and * . We have further partitioned submesoscale dissipation into dissipation due to wind-and buoyancy-forced SI, *w and *b . Gray dots denote the observed dissipation rates, obs , while the black square denotes the median value of these rates. We also illustrate the mean value of dissipation rates predicted by each process (triangles); depths assigned to these locations is the median and mean normalized depths, |z o |∕H, of the ADCP. The black triangle is not visible as it resides under red and blue triangles. The reason for the large mean value under buoyancy-forced convection, b , relative to the model-mean has not yet been determined. 2019; Viglione et al., 2018) are associated with strong mean flows that generate ocean fronts, filaments, and eddies characterized by considerable lateral density gradients. In such cases, wind-and buoyancy-forced SI will be particularly intense and the dissipation due to submesoscales, SM , might rival that of surface processes, SF . Indeed, that was found in the aforementioned studies of D' Asaro et al. (2011) and Thomas et al. (2013). The ocean front documented by D' Asaro et al. (2011), for example, was associated with a cross-frontal gradient in buoyancy of |∇ h b| ≈ 5-10×10 −6 s −2 and an Ekman-driven buoyancy flux as large as B e = 1.5×10 −5 m 2 /s 3 . For comparison, our largest documented buoyancy gradient (estimated from the inner mooring array) is approximately 2 × 10 −7 s −2 and the maximum observed value of the wind-driven buoyancy flux is B e = 4 × 10 −7 m 2 /s 3 . These values are 1-2 orders of magnitude smaller and occur infrequently during our observation period. In conclusion, we suggest that our results are not characteristic of oceans with such intense frontal dynamics but are more characteristic of areas where the ambient mesoscale eddy and/or Rossby wave field (Chelton & Schlax, 1996;Chelton et al., 2011) gives rise to ocean frontal gradients.
We have found it difficult to differentiate between current shear and Langmuir turbulence. In our case, this is most obviously due to the depth of the acoustic Doppler current profiler in the water column, |z o | ≈ 50 m. The vertical profiles of dissipation rate shown in Figure 13 are mean values for the winter season. However, instantaneous profiles show similar relationships. Together, these mean profiles suggest that, given an accurate dissipation estimate, one might be able to distinguish between Langmuir and current shear at normalized depths, close to |z|∕H ∼ 0.85. The profiles also suggest that in the absence of winds and waves it might be possible to observe a difference between surface buoyancy-forced convection and SI. Additionally, progress might be made by examining scaling laws for tracers in the presence of these varied turbulent mechanisms and then compare observations and predictions of upper ocean tracer concentrations as functions of depth. The argument being that the different types of turbulence will have very different effects on tracer concentrations. For example, Langmuir turbulence and current shear-induced turbulence have different signatures, with Langmuir circulations enhancing vertical motion, bringing parcels to depths greater than by shear-driven turbulence, alone (Thorpe, 1992). Here, the predictions could be informed by LES, such as those of Kukulka and Brunner (2015) and Brunner et al. (2015). It is probable that one could find scaling laws for tracer concentrations under SI (Smith et al., 2016) that could then be tested against observations. Indeed, repeat oxygen profiles were obtained at the OSMOSIS mooring site over the 2012-2013 period and show interesting dynamics in winter (U. Binetti and K. Heywood, personal communication, January 2017). See also Erickson and Thompson (2018).
An important difference between surface-and submesoscale-forced turbulence is that surface-forced turbulence extracts energy from the atmosphere, while submesoscale-forced turbulence draws energy from the kinetic and potential energy stored in ocean currents (D'Asaro et al., 2011;Thomas & Taylor, 2010). Winds inject energy into the general ocean circulation mostly through surface waves (Wunsch & Ferrari, 2004, p. 298), whereas submesoscale flows directly remove energy from the general ocean circulation (Müller et al., 2005). For this reason, the diminutive nature of the dissipation due to submesoscale processes (found at our moored location) does not necessarily imply that submesoscale processes are negligible. On the contrary, they can have an integrated effect upon upper ocean circulation owing to the lateral scale of the oceans and the ubiquity of ocean fronts. To illustrate this possibility, below, we consider the potential for submesoscale processes to modify upper ocean energy budgets, using SI as a particular case.

Implied Energy Loss Due to SI
Using the aforementioned model simulations, we integrated the total energy dissipated by SI (cf. equation (14)) within the NATL60 from the surface to 500 m. (The ML depth, H, never exceeds 500 m during the simulations.) This was done for the same 300-km × 300-km × 500-m volume of ocean, where we obtained an annual-averaged power loss per unit surface area of 420 W/km 2 . Assuming the entire ocean is characteristic of that found in the OSMOSIS region of the ocean-which we know is not correct but nonetheless will give us an approximate number-we multiplied the aforementioned value by the surface area of the global ocean, obtaining 1.42 × 10 11 W, or 0.14 TW. In making this latter calculation, we have approximated the ocean surface area as the sum of the areas in the following basins: North Atlantic (41,490,000 km 2 ), South Atlantic (40,270,000 km 2 ), North Pacific (77,010,000 km 2 ), South Pacific (84,750,000 km 2 ), Indian (70,560,000 km 2 ), and Southern Oceans (21,960,000 km 2 ).
To place this number into context, it is helpful to compare this value with the input into the ocean circulation by the wind, which is expected to be 1.1 TW (Ferrari & Wunsch, 2009;Wunsch & Ferrari, 2004). Thus, approximately 10-15% of the total rate of energy input into the oceans by the winds (i.e., "wind work") might be lost due to SI. Given the approximations and dependency of these results on our model, we emphasize that this estimate is provided merely as a scale value. In these calculations, we have not accounted for the reduction in wind work resulting from the differential motion of the ocean relative to the surface wind stress (Duhaut & Straub, 2006). This is believed to modify wind stress in key areas where the oceanic mean flow follows prevailing winds (e.g., in the ACC and Gulf Stream). Current estimates suggest the wind input is reduced by about 15-20% (i.e., 0.17-0.22 TW) (Hughes & Wilson, 2008;Zhai & Greatbatch, 2007), leading to a revised estimate of about 17-18% of the wind work. While this is simply a "back of the envelope calculation," it nonetheless suggests that, were the effect of SI able to be properly represented within ESM, the position of major ocean currents might change. Since these carry with them large reservoirs of heat, this would affect ocean-atmospheric heat fluxes and atmospheric circulation.
Submesoscale processes have other important effects not focused on here, including restratification of the surface ocean (Boccaletti et al., 2007;Fox-Kemper et al., 2011). SI, for example, might do this through the dissipation of geostrophic shear and subsequent relaxation of vertical buoyancy gradients into horizontal ones (Taylor & Ferrari, 2010). Furthermore, SI may be important for the supply nutrients to the euphotic zone (Smith et al., 2016). Thus, despite our main conclusion that they appear negligible for vertical turbulence budgets, they still may be relevant for ESM.

Limitations and Sources of Error
The interaction between ocean processes must be considered in any energy budget. Here, we have neglected the interaction between surface and submesoscale processes in order to simplify the TKE equation (cf. equation (2)). Moreover, several parameterizations utilized here have been developed under idealized ocean conditions, further simplifying the parameterization process. Lastly, internal waves can redistribute energy throughout the water column through the pressure-work term in the TKE equation, thereby potentially modifying dissipation rates from those presented. While we cannot comment on all of these without ourselves performing LES simulations capable of resolving these dynamics, some information might be gleaned from already existing high-resolution numerical studies in concert with some thought about the nature of the observations themselves.
From theoretical considerations, the nonlinearity of ocean surface waves and the associated Stokes shear is expected to modify frontal and submesoscale dynamics (Haney et al., 2015;Suzuki et al., 2016). For example, in a challenging numerical study involving both SI and Langmuir turbulence, Hamlington et al. (2014) found that Langmuir dynamics shift TKE production from vertical shear to one where convection plays a greater role. In another, wind-driven study using LES, Skyllingstad et al. (2017) found wind forcing at a lateral buoyancy gradient extracts energy from an ocean front in two forms: ageostrophic current shear and frontal geostrophic shear. Thus, two sources of turbulence exist: wind-forced current shear and GSP associated with SI. The coexistence of these two processes might, therefore, invalidate the parameterizations for SI utilized above. With respect to redistribution of energy due to internal waves, we argue that this is minimal owing to the duration over which the data were collected (i.e., 5 months). That is, the statistics are most likely robust to internal waves since these should not systematically redistribute energy in the vertical. In summary, while the number of these LES studies is few, we believe that the interaction between surface and submesoscale processes could affect our results significantly. Further research is therefore needed to better understand the interactions between processes within the OSBL.
An important limitation lies in our use of ERA5 winds, waves, and buoyancy fluxes. Recall that in the absence of coincident in situ meteorological observations, we elected to use ERA5 reanalysis to provide surface forcing for our dissipation parameterizations. However, coupling of the ocean and atmosphere might modify the results presented here. It is now well established, for example, that the mesoscale ocean has a measurable effect on the atmosphere above (Chelton & Xie, 2010;Foussard et al., 2019;Park & Cornillon, 2002;Small et al., 2008). These mesoscale processes are typically captured in the ERA5 product owing to the assimilation of satellite measurements in European Centre for Medium-Range Weather Forecasts's processing stream (Dee et al., 2011). However, considerably less is known about the impact of small-scale fronts on atmospheric conditions within the marine boundary layer (Renault et al., 2018;Wenegrat & Arthur, 2018) and these fine scales would be absent from ERA5. Thus, the effect of smaller-scale frontal structures that are absent from the ERA5 data may impact our results. This is one of the weaknesses of this study.
A significant source of error or uncertainty exists when calculating lateral gradients from the moored data. For one, there can be substantial errors when computing lateral gradients of velocity and density between inner moorings, which are spaced approximately 2 km apart and reside in 4,800 m of water (cf. Figure 2). Buckingham et al. (2016) found that the coherence between horizontal velocities on inner moorings was high (>0.7) for periods exceeding the inertial period (approximately 16 hr at this location), and this was the impetus for averaging for 1 day. Nevertheless, it must be acknowledged that the horizontal position of the moorings is uncertain and that this affects inner mooring gradients. The outer moorings do not suffer from this error as the uncertainty is at most 3% (supporting information).
Another potential error is that the surface ocean appears to have horizontal buoyancy gradients that extend to small spatial scales. That is, there is evidence from numerical models (J. Gula, personal communication, September 2019) and observations (D'Asaro et al., 2018) that gradients encountered at the submesoscale may become even more pronounced with increased resolution-for example, on the order of hundreds of meters. Since the inner moorings sample at 2-km resolution, this suggests estimates of B e reported here would be biased low. One possibility for checking if this is true (as suggested by a Reviewer) would be to simulate flows of the OSMOSIS region at ultrahigh resolution (e.g., dx = 75-100 m), perhaps for a month, and systematically explore what effect increasing the mooring separation has on the implied B e and * . We leave this as a future area of research.

Conclusions
We have demonstrated with this simple scaling framework applied to both observations that, relative to surface processes, submesoscale processes dissipate minimal TKE dissipation within the OSBL. This is only true (i) during winter, (ii) at the OSMOSIS site, and (iii) at approximately 50-m depth. However, we applied the same framework to the high-resolution numerical simulations (NATL60) and obtained comparable results throughout the entirety of the OSBL. The surface processes considered in this study included winds, waves, and buoyancy-forced (upright) convection, and the implied surface-forced dissipation rates were reproduced in two complementary models of surface-forced turbulence: with and without Langmuir. In contrast, the submesoscale processes considered here have been restricted only to SI. (An argument was made on account of the open-ocean nature of the environment and based on estimated categorized instabilities at the OSMO-SIS site that these are likely the more relevant submesoscale instabilities. Close to coasts and in the vicinity of strong current systems, horizontal shear should play a more important role; Gula et al., 2016;Jiao & Dewar, 2015.) Another notable finding is the challenge of distinguishing between turbulence production by Langmuir motions and current shear-where the latter that follows law-of-the-wall scaling with distance from the surface. This is chiefly because of the location of the instrument in the water column (i.e., |z| ∼ 40-50 m). It is noteworthy that there might be a greater ability to distinguish between these processes with increased depth, as suggested by averaged vertical profiles from the NATL60 model. Despite this finding, on account of the smallness of the turbulent Langmuir number (supporting information), we believe Langmuir turbulence is the dominant mechanism.
Additional to that documented above, we have highlighted a period in April 2013 when the probability of SI is particularly high. This was determined on account of the departure of z ′ = h∕H from unity, as well as the observed dissipation rates relative to surface-forced dissipation rates; that is, the observed rates appeared higher than expected from a prediction of winds and waves, alone. We refrain from examining this in detail, here, and leave this for a future study.
A final outcome of our analysis with the NATL60 model is that SI may play an important role in large-scale energy budgets, removing a nonnegligible fraction of geostrophic energy from the wind-driven ocean circulation. While dependent on the particular model simulations and resolutions used here, this nevertheless helps support the postulate that SI may serve as an important sink of energy otherwise available to the ocean circulation Thomas and Taylor (2010). Here, we estimate this value to be 10-20% of the wind work but we caution use of these values in any absolute sense; an additional study with two or more different models is necessary.
In summary, though turbulent dissipation rates due to submesoscale processes are small relative to those of surface processes in this open-ocean environment, we find that submesoscale motions may still be important for energy loss from the gyre-scale ocean circulation when integrated over larger domains. Finally, while we have not focused on this in our study, SI (and other submesoscale processes) may have important impacts on tracers (e.g., of carbon or nutrients) (Smith et al., 2016), may have an important restratifying effects through the slumping of ocean fronts (Taylor & Ferrari, 2010), and may be responsible for locally enhanced turbulence dissipation in instances of strong fronts and winds (D'Asaro et al., 2011;Thomas et al., 2013). These are active areas of study. For these reasons, we remain encouraged by early efforts to efficiently represent submesoscale processes within climate-scale ocean models . It will be interesting to know if the position and/or intensity of large-scale ocean currents are altered when examining models with the energetic effects of SI represented. It will additionally be interesting to know if this has an impact upon the atmosphere through large-scale changes in ocean-atmosphere heat exchange. can be solved for h. We solved this for the duration of the OSMOSIS observational record, 2012-2013. It should be noted that this was first solved for at the OSMOSIS site by Thompson et al. (2016). Figure A1 displays a portion of this time series, December 2012 to April 2013. Additionally, the likelihood of an observed value of z ′ falling into each of the bins is displayed to the right. Note that one can interpret z ′ = 1.0 as a baseline value such as in the absence of ocean fronts or in the presence of fronts but in the absence of atmospheric forcing, B e and B o . From these time series, we conclude that SI occurs only rarely during winter, occupying the OSBL (movement from z ′ = 1.0 toward z ′ = 0) rather infrequently. This maybe be due to the separation distance of outer moorings, as inner moorings result in a larger number of occurrences. However, we feel less confident about these statistics owing to mooring motion.