A Marginal Stability Paradigm for Shear‐Induced Diapycnal Turbulent Mixing in the Ocean

Turbulent mixing induced by breaking internal waves is key to the ocean circulation and global tracer budgets. While the classic marginal shear instability of Richardson number ∼1/4 has been considered as potentially relevant to turbulent wave breaking, its relevance to flows that are not steady parallel shear flows has been suspect. We show that shear instability is indeed relevant in the ocean interior and propose a new marginal stability paradigm that relates the stability criterion based on Richardson number to one based on the ratio of Ozmidov and Thorpe turbulence scales. The new paradigm applies to both ocean interior and boundary layer flows. This allows for accurate quantification of the transition from downwelling to upwelling zones in a recently emerged paradigm of ocean circulation. Our results help climate models more accurately calculate the mixing‐driven deep ocean circulation and fluxes of tracers in the ocean interior.


Introduction
Breaking internal waves induce widespread turbulent mixing at all ocean depths and across the globe, thereby playing an important role in (a) closure of the ocean circulation by upwelling of dense waters that form in polar regions and sink to the abyssal ocean basins and (b) regulating the budgets of heat, carbon, nutrients and other tracers important to the climate system (Talley et al., 2016;Wunsch & Ferrari, 2004).Among other sources, internal waves are forced at both surface and bottom ocean boundaries through winds, tides and ocean currents and eddies (Alford, 2020;Garrett & Kunze, 2007;Legg, 2021;Nikurashin & Ferrari, 2013).
Dating back to the pioneering work of G. I. Taylor (originally described in the Adams Prize Essay of 1915, and then in modified form in Taylor, 1931), it has long been known that under certain circumstances, locally enhanced vertical shear in a stratified fluid leads to an array of flow instabilities that can trigger the transition to turbulence and (ultimately) irreversible mixing.A key controlling parameter for stratified shear flows is the (gradient) Richardson number Ri (z, t), the ratio of the square of the buoyancy frequency N and the background vertical shear S: where, g is the acceleration due to gravity, ρ 0 is a reference density, and   ( ) and ( , ) are appropriate averages of the density and horizontal velocity.
In two classic papers, Miles (1961) and Howard (1961) established that a necessary condition for flow instability in a laminar steady inviscid shear layer is that     = 1 4 somewhere within the layer, where the flow is linearly marginally stable.Although this result is only proven for a very idealized flow situation, many studies (and indeed parameterizations) have been based around heuristic energetic arguments leading to the criterion that Ri (z, t) ≲ Ri L somewhere is a necessary condition for turbulence to occur within a flow.In practice, the specific theoretical value of   =  = 1 4 associated with the linear stability theory of laminar flow is identified as being a marginal Richardson number Ri m for the occurrence of turbulence.
Abstract Turbulent mixing induced by breaking internal waves is key to the ocean circulation and global tracer budgets.While the classic marginal shear instability of Richardson number ∼1/4 has been considered as potentially relevant to turbulent wave breaking, its relevance to flows that are not steady parallel shear flows has been suspect.We show that shear instability is indeed relevant in the ocean interior and propose a new marginal stability paradigm that relates the stability criterion based on Richardson number to one based on the ratio of Ozmidov and Thorpe turbulence scales.The new paradigm applies to both ocean interior and boundary layer flows.This allows for accurate quantification of the transition from downwelling to upwelling zones in a recently emerged paradigm of ocean circulation.Our results help climate models more accurately calculate the mixing-driven deep ocean circulation and fluxes of tracers in the ocean interior.(1979) argued that stratified shear flows can naturally adjust into a "kind of equilibrium" where local (in space and/or time) intensification of shear will lead to a local drop in Ri, allowing instability, enhanced turbulent dissipation and hence reduction in the shear until the flow becomes at least close to a marginally stable state.In an analogous fashion, S. A. Thorpe and Liu (2009) argued that generically stratified shear flows will adjust towards a state of such "marginal instability", with   ≃  ≃  = 1 4 .Indeed, as originally shown by Woods (1968) by direct dye measurements in the Mediterranean, and subsequently in many observational studies (see for example van Haren & Gostiaux 2010;Geyer et al., 2010), vortical structures at least visually reminiscent of the classic shear-driven overturning "Kelvin-Helmholtz" instability appear to play an important role in turbulent mixing in the world's oceans, and observational measurements of Richardson number have also been shown to be peaked around the theoretical linear stability value of   = 1 4 in various oceanic environments (W.D. Smyth, 2020; W. D. Smyth & Moum, 2013;W. Smyth et al., 2017).
Furthermore, both numerical simulations and observational data suggest that shear-driven mixing events are not only well-characterized by being (close to) marginal linear stability, but also tend to be in another "critical" or marginal state where the length scale ratio R OT ∼ 1 (e.g., Dillon, 1982;Imberger & Ivey, 1991).R OT is defined as the ratio of the Ozmidov length scale L O to the Thorpe scale L T

∶=
; where, ϵ is the turbulent kinetic energy (density) dissipation rate, and L T is the (purely geometric) scale constructed from the root-mean square of the displacements required to sort fluid parcels in a (locally statically unstable) overturning density profile into a monotonic, statically stable distribution.L O may be interpreted as the largest vertical scale which is essentially unaffected by stratification, while the Thorpe scale may be interpreted as the energetic overturning (and hence turbulence injection) scale.Mashayek et al. (2021) (hereafter, MCA21) argued that this critical matching of length scales is characteristic of an intermediate mixing phase in a turbulent patch's life cycle where, L O is both close to its maximum value during a transient shear-driven mixing event (and so the turbulence is most energetic) and the mixing is most efficient (as the overturning scale L T ∼ L O , and so energy is injected right at the top of the dynamic range of turbulence largely unaffected by the stratification, and so with at least the potential to exhibit an isotropic cascade for all scales smaller than the overturning scale).Due to the optimally tuned nature (in that the turbulence is simultaneously most energetic and most efficient) of this phase of the flow evolution, they referred to it as the "Goldilocks" mixing phase.One of the two primary goals of this manuscript is to investigate whether there is a connection between the marginal stability paradigm based on Ri and the Goldilocks paradigm of MCA21 based on R OT , that is, is it possible to identify this emergent optimal state (i.e., most energetic and most efficient state) conceptually with some marginally stable state of the flow.It is our second primary goal herein to test whether in complex, realistic flows, specifically those not characterized by directly forced shear instability (as discussed in W. D. Smyth, 2020), a generalised marginal stability paradigm for shear-driven overturning mixing is of value for describing significant mixing events.
To investigate the above-mentioned two goals, namely (a) the potential connection between shear-driven marginal stability and the R OT -based Goldilocks mixing paradigm, and (b) the relevance of (generalized) marginal stability to dynamically complex oceanic turbulence zones, we consider three regions where it is not immediately clear that shear-driven mixing should be significant.We consider mixing in the Drake Passage, in the Brazil Basin tidal region, and in an abyssal canyon also in the Brazil basin.

Drake Passage
Figure 1 shows one of the most energetic zones in the ocean, the Drake Passage of the Southern Ocean, where strong currents and eddies pass through a narrow constriction and over rough topography.The results are from a snapshot output from an observationally forced and verified ocean model at unprecedentedly high spatial resolution.The model has been shown to reproduce hydrographic structures, mesoscale dynamics, and tracer transports in excellent agreement with observations-see Supporting Information S1 for more details.Of particular relevance to this work is the ability of the model to reproduce the small scale mixing in close agreement with observations (see top-right inset in panel a of Figure 1).Panel a shows how an energetic complex full depth internal wave field is excited by flow topography interactions at the bottom and strong westerly winds at the surface.Nonlinear wave-wave and wave-current interactions lead to localized increases in shear, here defined as where U and V are the zonal and meridional velocities contributing to the horizontal velocity as mentioned in Equation 1. S 2 , overlain by density contour lines, is shown in panel b for a longitude-depth slice of the simulation.The stability criterion marking the emergence (or lack thereof) of shear-induced turbulence is a marginal value of the Richardson number, as previously described.
Panel c shows Ri, normalized by its marginal value, that is, the value above which the flow is assumed to be unable to sustain turbulence.Given the high vertical resolution in the model, we set Ri m to  1 3 , which is close to but actually slightly higher than the canonical linear marginal value of   = 1 4 .This choice is made because there have been suggestions that existent ambient residual turbulence can bias the stability criterion slightly higher or lower (Kaminski & Smyth, 2019;Li et al., 2015;S. Thorpe et al., 2013).Furthermore, mixing layers often actually re-laminarize at   ∼ 1 3 (Pham & Sarkar, 2010;Mashayek, Caulfield, & Peltier, 2017;W. D. Smyth & Moum, 2000).Conversely, there is also numerical and observational evidence that appropriate values of Ri m (i.e., for sustained turbulence) can actually be as low as 0.16 (as reported in (Holleman et al., 2016;Portwood et al., 2019)) through 0.21 (Zhou et al., 2017;van Reeuwijk et al., 2019)  Note that the coarser a model, the higher Ri m needs to be set to account for the unresolved subgrid-scale turbulence.Indeed, climate models often use Ri m = 0.7.Panel c shows that Ri ≈ Ri m at enhanced sheared locations in the interior and along the top margin of the bottom and top boundary layers (the latter hard to see in the figure).Within the boundary layers, the flow is increasingly less stratified, leading to lower values of Ri.Boundary layer turbulence is complex, a topic of extensive active research, and parameterized heavily in ocean models including ours.Crudely speaking, at the risk of oversimplification, most such parameterizations involve "shearing" and "convective" processes within boundary layers, with the former often parameterized based on Ri.In such layers, the flow is strongly unstable to shear instability as opposed to the interior flow which is predominantly marginally stable (S.Thorpe & Liu, 2009).
It is the interior mixing, however, which is the focus of this study.In the interior of our model, mixing is parameterized in terms of Ri using the now classic "KPP" model of Large et al. (1994): where, κ max = 5 × 10 −3 m 2 s −1 , and   = 1 3 as previously mentioned, and κ is the turbulent or eddy diffusivity.As the local Ri approaches Ri m , the diffusion, and hence the turbulence, is parameterized to decrease, until it is "switched off" at the chosen input value of Ri m .When no turbulence is excited by shear in "quiet" regions, the model's mixing is set to a background value.Since the focus of our study is on a parameterization based around some kind of marginal stability criterion in the ocean interior, we will not discuss the inference of a turbulent diffusivity in the model.However, we'd like to point out that the diffusivity output of the model agrees excellently with observed profiles of microstructure turbulence in the Drake Passage (see top-right corner inset in Figure 1a; also see Supporting Information S1), implying that the dynamics and the resulting shear distributions are sufficiently resolved in the model for the confident discussion of a stability criterion.Here, we will merely use the model's diffusivity output to separate the turbulent regions from the quiet background flow where diffusivity is set to the background κ b .Thus, employing a model, we can achieve a separation which would not be feasible in the real world.

Marginal Stability and Goldilocks Mixing
Recently MCA21 proposed a unifying parameterization for the turbulent flux coefficient Γ in terms of R OT in the form of where, ρ′ is the perturbation density, w′ is the perturbation velocity (and so   is the appropriately averaged vertical (specific) density flux, or (negative) buoyancy flux).Assuming that the mixing occurs at Ri go in the range of marginal Richardson numbers  1∕6 ≤  ≤ 1 4 , MCA21 inferred (on physical grounds) that   in the R OT ≫ 1 limit, which corresponds to decaying turbulence.As reviewed in MCA21, this scaling has been proposed by many authors, dating back to the "fossil" turbulence arguments of (Gibson, 1987).In the young turbulence limit of R OT ≪ 1, however, Equation 4 reduces to  Γ ≈   −1  which was proposed by MCA21.Moreover, they argued that it is the intermediate adjustment phase between the two limits which corresponds to optimal mixing, a phase they dubbed "Goldilocks Mixing" since in parallel with the fairy-tale it exists when there is just the perfect balance between energy available to turbulence from the background shear and the local stratification which induces stratified turbulence, yet suppresses vertical turbulent motions at the same time.As noted above, this optimality manifests itself in the turbulent flow simultaneously being at its most energetic and its most efficient in mixing, which appears to be closely related to this balance and matching of the Thorpe overturning scale and the Ozmidov scale.Importantly, MCA21 showed that oceanic data of turbulent patches seem to cluster consistently around this Goldilocks balance of R OT ∼ 1, thus suggesting another as yet unexplained emergent phenomenon of flows organizing so that R OT ∼ 1.
Figure 2-top, reproduced from MCA21, again shows the excellent agreement between the Goldilocks parameterization (4) with 6 oceanic datasets comprising a total of ∼50,000 turbulent patches excited by different turbulence processes at different geographical and oceanic depths (see Supporting Information S1 for a brief description of the data and MCA21 for comparisons for individual datasets).Interestingly, inferring A from regressing (4) to the diverse collated data set yields A = 0.68 ∼ 2/3 which corresponds to   =  =  = 1 4 .This is a clear hint that marginally stable shear instability is broadly relevant to interior ocean mixing.It also highlights the power of Equation 4 in that it is entirely based on physical grounds, even up to the coefficient A (i.e., assuming that maximally efficient and energetic, i.e., "optimal" mixing occurs close to the marginal Richardson number at which shear turbulence can be maintained), and agrees well with data.
The use of appropriate definitions of a buoyancy Reynolds number Re b : where, ν is the kinematic viscosity, and L K is the Kolmogorov microscale, and/or a Richardson number Ri to quantify the flux coefficient is relatively well established (Caulfield, 2021;Gregg et al., 2018;Ivey et al., 2008;Peltier & Caulfield, 2003).However, various datasets do not overlap when mapped onto these parameters (Bouffard & Boegman, 2013; Mashayek, Salehipour, et al., 2017;Monismith et al., 2018).To highlight this tendency in the datasets employed in this article, Figure 2-bottom plots the same data set as in the top panel but against Re b .Put bluntly, the data are all over the place, and so some further analysis is required to resolve the scatter.
On dimensional more than one nondimensional parameter is required to quantify mixing (Ivey & Imberger, 1991;Mashayek & Peltier, 2011, 2013;Mater & Venayagamoorthy, 2014;Shih et al., 2005).(Re b , Ri) is perhaps an obvious pair of parameters on which the flux coefficient Γ might reasonably be assumed to depend, particularly for shear-driven stratified mixing, although other pairs have been proposed before (Ivey & Imberger, 1991;Mater & Venayagamoorthy, 2014).Through the lens of the three stage "Goldilocks mixing" life cycle, a parameterization based on Re b and Ri may be proposed in such a way that demonstrates at least reasonable agreement with the data.There is substantial empirical evidence that for sufficiently large Re b ,  Γ ∝  −1∕2  (Bouffard & Boegman, 2013;Ivey et al., 2008;Mashayek, Salehipour, et al., 2017;Monismith et al., 2018).There is also some evidence from experimental and numerical data (see review in Bouffard & Boegman 2013; also see Mashayek, Salehipour, et al., 2017) that for smaller Re b ,  Γ ∝  1∕2  .Finally, there is also substantial evidence for Γ ∝ Ri in not particularly strongly stratified flows especially when Ri is identified with the inverse square of an appropriate Froude number Fr ≡ U/NL where, U and L are characteristic velocity and length scales (Lozovatsky & Fernando, 2013;Maffioli et al., 2016;Salehipour & Peltier, 2015;Shih et al., 2005;Wells et al., 2010;Zhou et al., 2017).Combining these various observations into a simple empirical relation consistent with Equation 4, we obtain (see Supporting Information S1) where, and Ri m and Re bm are the values of Ri and Re b where, the Goldilocks mixing phase occurs, that is, where R OT ∼ 1.This generalizes (to include Re b ) the observation that the Goldilocks mixing phase appears to occur effectively at a marginal or critical value of Ri go ≃ Ri m .
By normalizing Re b and Ri with their values evaluated at R OT ∼ 1, Equation 6 may be interpreted analogously to Equation 4 as representing two limits of the young/growing turbulence (   *  ≪ 1 ) and fossilizing/decaying turbulence (   *  ≫ 1 ), smoothly connected at   *  ≈  * ≈  ∼ 1 , that is, at the special (in at least some sense marginal) values of buoyancy Reynolds number and Richardson number at which R OT ∼ 1, within the Goldilocks mixing phase.This has two significant implications.
First, it explains (at least partially) the shift of the peak of the curves in Figure 2-bottom (with variations in Ri explaining the rest); the mere existence of peaks in those curves implies the relevance of the idea of a flow case-sensitive critical or marginal Re b = Re bm .Second, it implies that the concept of marginal stability based purely on a special value of Ri = Ri L associated with the onset of linear shear instability is only a part of the picture and an appropriate two-parameter generalization is essential as implied by the observed peaks of distributions varying with both Ri and Re b .For example, W. D. Smyth (2020) subsampled the most intense turbulent patches (i.e., focused on a particular narrower range of Re b ) to show that those patches appeared to be in an apparently linearly marginal state with   ∼  = 1 4 .Using the hypothesis that turbulent flows organize and mix with R OT ∼ 1, on the other hand, seems to be a more natural generalized critical or "marginal stability" criterion, at least with respect to parameterizing the turbulent flux coefficient Γ.Furthermore, while Ri and Re b vary over a few to several orders of magnitude, 80% of data in Figure 2 lie within a factor of 3 of R OT = 1, implying that expressing and quantifying mixing in terms of values of R OT is much more suitable.

Direct Application to Ocean Data
We now show that Equation 4 and Equation 6 appear to agree reasonably well when applied to two oceanic datasets that include all the parameters involved in the two formulations.Two of the datasets used in Figure 2 were from the Brazil Basin Tracer Release Experiment (BBTRE), and the Dynamics of Mid-Ocean Ridge Experiment (DoMORE).The former sampled turbulence induced by internal tide shear in the deep Brazil Basin over the mid-Atlantic ridge (Polzin et al., 1997), while the latter sampled turbulence over a sill on a canyon floor also in the Brazil Basin (Clément et al., 2017).So, while for BBTRE it is expected that nonlinear wave-wave interactions will downscale energy to small scales where shear instabilities can ultimately kick in (e.g., see Nikurashin & Legg, 2011), the flow is expected to be somewhat hydraulically controlled for DoMORE with at least part of the sampled turbulence corresponding to boundary layer turbulence, somewhat similar to the deepest canyons in Figure 1c.The two datasets were recently analyzed by Ijichi et al. (2020) where they identified turbulence patches and calculated their corresponding L O , L T , and Ri values.This is convenient for our purposes as most often explicit shear measurements co-located with profiles of density and ϵ are lacking.
The first test is whether the connection we made between the Ri L -based linear marginal instability criterion and the proposed generalized stability criterion of R OT -based Goldilocks mixing holds, based on the heuristic equivalency of Equation 4 and Equation 6.We start by inferring appropriate values of Ri m and Re bm through appropriately averaging the observational values of these quantities associated with observations where R OT ∼ 1.We then equate the value of Γ constructed directly from the observations of R OT using Equation 4 with the value of Γ given by Equation 6.In Equation 6, we input measured Re b , and the inferred values of Ri m and Re bm to obtain a parameterization prediction for Ri param .Figures 3a and 3b show this Ri param , plotted against the directly measured Ri.The agreement is impressive for the peak of both pdfs.For BBTRE, the agreement seems to also hold nicely over the entire range of Ri, consistent with the hypothesis that the ocean interior is close to marginally stable and turbulence bursts are excited by localized increases in shear (thereby a local drop in Ri).For DoMORE, on the other hand, the data consists of patches both within and outside of the bottom boundary layer.Thus, the range of Ri extends from close to   = 1 4 to much smaller values compared to BBTRE (somewhat similar to Figures 2b and 2c).Richardson numbers inferred from equating the two expressions for Γ, Equations 4 and 6, tend to over-predict the low end of Ri in Figure 3b as Equations 4 and 6 are meant to only apply to interior mixing, "sufficiently" far from boundary layer turbulence so that shear layers are "free".According to Figure 1c, the upper end of the bottom boundary layer marks the limits of the relevance of Equation 4 and Equation 6 and coincidentally is where they are quite actively relevant as that interface is rife with shear-driven turbulence, similar to mixing at the base of the surface ocean mixed-layer.
Figures 3c and 3d further illustrate the correspondence between our proposed generalized marginal stability criterion based on R OT and a more conventional marginal stability criterion based purely on Ri m .Panel c shows that for turbulent patches sampled during BBTRE, the joint probability distribution of Ri and R OT is sharply peaked at Ri ≈ 0.2 and R OT ≈ 1.In physical terms, most turbulence patches are likely excited through shear instability (hence, with Ri slightly less than   = 1 4 ) and are in the Goldilocks phase where the APE absorbed from the mean flow is actively energizing the turbulence inertial subrange with L O as an upper bound.In contrast, panel d shows that for turbulent patches sampled during DoMORE, the distribution of Ri is much more spread, with a peak probability density at a modal Ri value of ∼0.1 .As discussed above, this is due to the near boundary nature of sampling in DoMORE which implies lower stratification and hence lower Ri.The R OT distribution, however, remains tightly bound (similar to BBTRE) with the mode , implying that non-shear-instability-like overturns likely contribute at least to some of the data.Hydraulically induced overturns are some candidates (e.g., MCA21, through analysis of data from Carter et al. (2019), show that the paradigm that R OT ∼ 1 seems to have relevance to hydraulically induced flow over sills in the abyssal Samoan Passage).Panels c and d, therefore, reinforce our argument that the R OT -based generalized marginal stability criterion (which we have demonstrated can be thought of as having an implicit 2D Re b − Ri dependence) is more general as compared to one solely based on Ri.A purely Ri-based marginal stability criterion is incomplete without consideration of an appropriate companion buoyancy Reynolds number and is not directly relevant to non-shear-instability-induced overturns which seem to abound in boundary layers.8 of 11 appears to hold while R OT ∼ 1 for most patches except for the top quantile of ϵ for which the modal value is R OT ∼ 2. This is perhaps unsurprising since the most energetic phase of turbulence often occurs simultaneously with or shortly after turbulence has grown rapidly (and so L O is maximal) at the cost of the conversion of some of the stored APE in the overturn (and so L T has fallen from its maximum).For DoMORE, as discussed above, the patches closer to the boundary correspond to relatively low stratification and Ri while those further away and in the interior are still likely prone to shear instability (noting that shear instabilities can still be relevant in deep canyons and flows over sills e.g., Alford et al., 2013).For all the quantiles of data in DoMORE, R OT ∼ 1, highlighting again the more general applicability of our proposed generalized marginal stability criterion based on R OT .

Discussion
We have investigated the relevance of a criterion based on the concept of marginally stable shear instabilities and found that in the ocean interior such a criterion is relevant even in dynamically complex regions.In such regions, energy downscales from mesoscales and submesoscale dynamics, or from the internal wave field (e.g., induced by tides) to scales sufficiently small that localized shear can induce small scale shear instability and mixing.Within the boundary layers, using a Ri-based criterion alone is not so successful.
We have shown that a generalized marginal stability criterion based on assuming that the flow adjusts towards R OT ∼ 1 holds in the interior and seemingly also within the boundary layers.We have demonstrated that the R OTbased marginal stability criterion can be related to an alternative criterion, which explicitly depends jointly on Ri and Re b .We have argued that a collection of oceanic datasets show that each field experiment (or turbulence region) has its own flow case-sensitive generalized marginal stability criterion with respect to the particular marginal values of Re bm and Ri m , whereas the R OT ∼ 1 criterion holds more universally, and hence more usefully.Our finding sheds light on the seeming significant discrepancy between parameterization of ocean mixing based on Re b alone and puts bounds on the relevance of the conventional marginal shear instability criterion based purely on linear stability arguments that   ≃  = 1 4 .The parameterization proposed in MCA21 for the flux coefficient Γ, as reviewed in Figure 2a, together with the ability of the generalized marginal stability criterion based around R OT ∼ 1 to cross over from interior turbulence to boundary layer dynamics, have the potential to be very useful.In tandem, they allow unified consideration of several different mixing regimes.For example, it is possible to consider: (a) weak mixing in the ocean interior far from the boundaries (where, there is less energy available to turbulence and the stratification is relatively strong); (b) energetic efficient mixing in the vicinity of the boundary between interior dynamics and boundary layers (where the "right" balance of shear and stratification leads to the most efficient mixing, dubbed "Goldilocks mixing" by MCA21); and (c) weakly mixing regions deep within boundary layers (where, while energy is available, stratification is weak and hence there is less to mix).
There has been a recent paradigm shift in our understanding of the role of deep ocean turbulence in the global ocean overturning circulation, thereby in the climate system.The new paradigm suggests that the turbulence in the ocean interior above rough topography can lead to densification and downwelling of water masses, while boundary layer turbulence is primarily responsible for lightening and upwelling of the dense waters that form as such plus the dense waters that form at high latitude and sink to the deep ocean (de Lavergne et al., 2016b(de Lavergne et al., , 2017;;Ferrari et al., 2016;Mashayek et al., 2015).The transition from interior downwelling to boundary upwelling is marked by a change in the sign of the vertical gradient of the effective flux of buoyancy, approximated as the multiplication of the local flux coefficient and the local rate of dissipation of kinetic energy.Thus, accurate quantification of the flux coefficient is key to an accurate calculation of the deep branch of ocean circulation.Before the work presented here, parameterizations of small scale mixing based on Ri and/or Re b have so far been incapable of fully capturing the subtlety of transitioning from interior to boundary mixing (de Lavergne et al., 2016a;Cimoli et al., 2019;Mashayek, Salehipour, et al., 2017).

Figure 1 .
Figure 1.(a)A sector of the Southern Ocean in which strong westerly winds at the surface and interaction of the Antarctic Circumpolar Current system with rough topography create an energetic internal wave field, shown by red/blue shading representing upward/downward vertical velocity.From an observationally forced and tuned high resolution numerical simulation(Mashayek, Ferrari, et al., 2017).The top-right inset shows the model diffusivity, κ, compared with microstructure observations from DIMES.(b) A longitude-depth slice of squared shear (i.e., vertical gradient of horizontal velocity) superimposed by density layers (white lines).(c) Richardson number normalized by the marginal value of 0.33 used in the model.(d) Probability density function of Richardson number for turbulent regions for the full 3D domain, for the full domain excluding the top and bottom 200 m, and for the non-turbulent background regions.Turbulent regions are where turbulence level is above the model background value of κ b = 5 × 10 −5 m 2 s −1 .The top and bottom 200 m are excluded in one case to exclude the top mixed layer and the bottom boundary layer where other parameterized processes, in addition to shear instability, can create mixing in the model.

Figure
Figure1dshows the probability density function (pdf) of Ri in the full 3D domain over turbulent regions and "quiet" regions.The pdf has a sharp peak at 0.276, marked by the dashed line, which has emerged to occur quite close to the imposed input marginal value of   = 1 3 .This implies that co-evolution and nonlinear interactions between mesoscale-submesoscale-wave processes have led to a downward cascade of energy that ultimately has

Figure 2 .
Figure 2. Distribution of the flux coefficient as a function of R OT = L O /L T (top) and Re b = ϵ/(νN 2 ) (bottom) from ∼50,000 turbulent patches from a total of six oceanic experiments covering a wide range of depth, geographic locations, and turbulence generation processes; see Supporting Information S1 for a description of the datasets.The original patch data is shown in small black dots (with their histograms shown along the right and top axes) and the experiment-binned/mean distributions are shown in thick lines with large symbols in the top panel and with thick lines in the bottom panel.Both panels are reproduced from Mashayek et al. (2021) where details of the datasets may be found.

Figures
Figures 3e and 3f further emphasize the points made above by showing how the modal values of the distributions of Ri, R OT and their joint distribution vary for patches with quantiles of ϵ.For BBTRE, for all the data (weakly turbulent and strongly turbulent patches) the marginal shear instability criterion based on a marginal value Ri m

Figure 3 .
Figure 3. Joint probability distribution for Ri from data and inferred from parameterization (6) as described in the text for Brazil Basin Tracer Release Experiment (BBTRE) (a) and DoMORE (b).White line corresponds to Ri data = Ri param. .Scatterpoint color indicates local probability density.Ri param is inferred from Equation 6 by taking Re b and Γ from data, and solving for Ri(c and d) Joint probability distribution of Ri and R OT for BBTRE and DoMORE.Green cross indicates estimated mode of probability distribution.The modal values are shown in the legend (not in log)(e and f) Modal values for probability distribution of Ri alone and for joint distribution of Ri and R OT for different ϵ quintiles of patches in BBTRE and DoMORE.