5.1. Length scale ratio for $\varGamma _{\mathcal {M}}$

First, we assume that the largest energy containing eddies of length $L_i$ inject energy into the system and that such an input is ultimately balanced by an energy sink at the viscous dissipation scale after a (net) forward cascade. This leads to the classic turbulent (and inherently unstratified) integral scale,

(5.1)\begin{equation} L_i \propto \frac{k^{3/2}}{\epsilon}, \end{equation}

where $k$ is the TKE. This scaling has at its heart that the turbulence being considered is not strongly affected by stratification.

Next, we invoke classic mixing length theory (Taylor 1915; Prandtl 1925), which gives the turbulent diffusivity in terms of a mixing length $L_{m}$ and $k$

(5.2)\begin{equation} \kappa_T\propto L_{m}\ k^{1/2}. \end{equation}

Here, the (possibly significant) dynamical influence of ambient stratification may be embedded within the mixing length, but it also contains an implicit assumption that the mixing of the buoyancy is assumed to be set by (and so in a sense locked to) the mixing of momentum, as quantified by the right-hand side of the expression. Equivalently, it is assumed that the turbulent Prandtl number $Pr_T \sim O(1)$, defined as

(5.3a,b)\begin{equation} Pr_T \equiv \frac{\nu_T}{\kappa_T};\quad \nu_T \equiv \frac{\left \langle u' w' \right \rangle}{ S}\simeq \frac{P}{S^{2}}, \end{equation}

where ${\nu }_T$ is the eddy diffusivity (of momentum), and we assume that the turbulent production $P$ is dominated by Reynolds stress extraction from the mean shear $S$. Therefore, it is most definitely not appropriate to think of the flow as ‘strongly’ stratified in any meaningful sense.

Also, by its very nature ‘overturning’ mixing must be occurring in (at least locally) a subcritical stratification as observed by Alford & Pinkel (Reference Alford and Pinkel2000). Such overturning mixing is qualitatively different from the ‘scouring’ mixing in the vicinity of ‘sharp’ and strong density interfaces, a classification distinction drawn by Woods et al. (Reference Woods, Caulfield, Landel and Kuesters2010) (see Caulfield (Reference Caulfield2020,  Reference Caulfield2021) for more discussion).

Finally, we invoke the Osborn balance presented previously in (2.6):

(5.4)\begin{equation} \kappa_T\approx \varGamma_{\mathcal{M}} \frac{\epsilon}{N^{2}}. \end{equation}

These three ingredients can form the basis for a simple parameterization that fits all the data introduced earlier, provided various data are thought of as being sampled at different stages in the time evolution of a particular mixing event. To achieve this, however, it is also important to distinguish between the two characteristic length scales defined in (5.1) and (5.2). Remembering that $\epsilon$ and $N$ can in turn be related to the Ozmidov length scale since $L_O^{2}=\epsilon /N^{3}$, it is possible to use these equations to arrive at

(5.5)\begin{equation} \varGamma_{\mathcal{M}}\propto \left[\frac{L_i L_{m}^{3}}{L_O^{4}} \right]^{{1}/{3}}. \end{equation}

To develop a parameterization, it is now necessary to identify appropriate candidate lengths for $L_m$ and $L_i$ that can actually be quantified. Since here we are interested in overturns far from boundaries, we expect that the distance to the boundary is not a relevant candidate. Motivated by the phenomenology discussed in § 3, we argue that there are different appropriate candidates for these scales at different stages of the flow evolution of an overturn in a sufficiently weakly stratified flow.

5.2. ‘Young’ turbulence scaling for $\varGamma _{\mathcal {M}}$

During the first growing phase of the turbulence, $L_T \gg L_O$. During that phase, it is appropriate to identify the mixing length $L_m \sim L_T$, as the Thorpe scale is the relevant scale over which the fluid gets stirred and mixed against the background stratification. On the other hand, the actual turbulent motions (and the associated inertial range for the forward cascade of energy) are expected to be constrained to lengths bounded above by $L_O$, and so a reasonable estimate for $L_i$ should be $L_O$. Therefore, (5.5) reduces to

(5.6)\begin{equation} \varGamma_{\mathcal{M}} \propto R_{OT}^{{-}1} \end{equation}

for this first ‘young’ turbulence phase.

5.3. ‘Fossilization’ scaling for $\varGamma _{\mathcal {M}}$

Conversely, during the final decaying or ‘fossilization’ phase when $L_T \ll L_O$, it seems reasonable that the mixing length can still be identified with the overturning scale, so $L_m \sim L_T$. In this stage, $L_T$ can also be closely related to the integral scale of the turbulence, $L_i$, as the energy for turbulence and mixing is injected by the remaining, and smaller scale than earlier in the flow evolution, overturnings; $L_O$ may now be interpreted as the largest eddies which could in principle be present within the given ambient stratification. However, for such eddies actually to arise, energy has to be injected in to the system at a sufficient rate, and since $L_T$ has dropped below $L_O$, there is no longer a mechanism by which that injection can occur. Crucially, this should not be interpreted as ‘strong’ stratification suppressing (in particular) vertical motions inducing strongly stratified anisotropic flow, but rather that turbulent dissipation has converted a large proportion of the KE, and there is (at the moment) no forcing or energy injection mechanism occurring. Such a gap between the actual eddy length scale (scaling with $L_T$) and the notionally possible eddy length scale (scaling with $L_O$, if only the flow were energetic enough) is one of the hallmarks of ‘fossil turbulence’ (Gibson Reference Gibson1987). It is important to remember that both $L_T$ and $L_O$ decrease with time, but $L_T$ decreases more rapidly, and so the $L_O$ is ‘left behind’ by the more rapid decay of the energy injection scale associated with $L_T$.

Thus, in the final fossilization phase, (5.5) reduces to

(5.7)\begin{equation} \varGamma_{\mathcal{M}} \propto R_{OT}^{{-}4/3}, \end{equation}

where once again it is important to remember that this argument is fundamentally associated with the stratification always being below some critical, threshold strength, not least because of the implicit assumption that buoyancy and momentum are mixed on the same scales and by the same processes. In summary, we thus suppose that (5.6) holds in the $R_{OT} \ll 1$ limit, corresponding to the first growing ‘young turbulence phase’, while (5.7) holds in the $R_{OT} \gg 1$ limit, corresponding to the final decaying or ‘fossilization phase’, with the variation in the dependence of the flux coefficient with $R_{OT}$ being associated with different (temporal) phases in the evolution of the overturning stratified mixing event, where buoyancy forces never dominate the dynamics.

5.4. Goldilocks scaling for $\varGamma _{\mathcal {M}}$

As previously reported in observations (e.g. see Dillon Reference Dillon1982; Ferron et al. Reference Ferron, Mercier, Speer, Gargett and Polzin1998), and further confirmed by our analyses of these oceanographic datasets, most of the observed overturns actually exist around $R_{OT}\sim 1$, which, from § 3, occurs during the intermediate, energetically mixing phase. As pointed out in Mashayek et al. (Reference Mashayek, Caulfield and Peltier2017a), this also actually turns out to be a regime of efficient (indeed apparently optimal) mixing, which we refer to as Goldilocks mixing. To achieve a parameterization for this regime, we simply combine (5.6) and (5.7) using the following formula:

(5.8)\begin{equation} \varGamma_{\mathcal{M}} =A\frac{R_{OT}^{{-}1}}{1+R_{OT}^{{1/3}}}, \end{equation}

which is constructed to have the correct asymptotic properties at small and large $R_{OT}$.

Before discussing the physical interpretation of (5.8) and attempting to determine the scaling coefficient ‘A’ on physical grounds, it is very important to appreciate that many of the arguments presented above are not new. In fact, employing $R_{OT}$ as a turbulence age proxy goes back several decades as reviewed by Smyth & Moum (Reference Smyth and Moum2000), who also proposed that it be used to quantify $\varGamma _{\mathcal {M}}$. Numerical simulations and observations have also shown clear dependence of mixing on $R_{OT}$ (Smyth et al. Reference Smyth, Moum and Caldwell2001; Mashayek et al. Reference Mashayek, Caulfield and Peltier2017a; Ijichi & Hibiya Reference Ijichi and Hibiya2018; Ijichi et al. Reference Ijichi, St. Laurent, Polzin and Toole2020; Smith Reference Smith2020). Equation (5.7), in particular, has been previously proposed by many authors (Weinstock Reference Weinstock1992; Schumann & Gerz Reference Schumann and Gerz1995; Baumert & Peters Reference Baumert and Peters2000; Ijichi & Hibiya Reference Ijichi and Hibiya2018).

Indeed, these previous authors typically invoked turbulence properties largely unaffected by ‘weak’ stratification. We go a step further by arguing that (5.7) is appropriate only late in the life cycle of a time-varying overturning mixing event, where the stratification is only required not to be so strong that it always remains sub-dominant in the flow evolution. Furthermore, $L_T$ was argued by Ivey & Imberger (Reference Ivey and Imberger1991) as representative for the actual mixing length scales (once again without explicit consideration of time dependence), an idea which was verified later via numerical simulations (see Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005), for example). However, we are not aware of a previous proposal of (5.6) for ‘young’, yet vigorously overturning mixing events, and as we shall show, its combination with (5.7) leads to an accurate parametric modelling of $\varGamma _{\mathcal {M}}$ through (5.8) (or indeed alternative definitions like $\varGamma _{\chi }$) in a way that would not be possible merely based on (5.7).

5.5. Context in terms of previous studies

It should also be noted that while we use the paradigm of a single (with three temporal phases) mixing life cycle for arguments herein, our arguments also hold in a system composed of a multitude of interacting, co-existing overturns, at different stages of their individual life cycle. In such systems, the young phase simply refers to the phase where large quasi-laminar and relatively recently formed overturns exist and provide APE to small scale overturns to grow, while the fossilization phase corresponds to when the APE source has been largely depleted and the actual still-occurring overturns are smaller than the notionally possible size scale as marked by $L_O$.

Mater et al. (Reference Mater, Venayagamoorthy, St. Laurent and Moum2015) showed, using several oceanic datasets, that the evolution of $R_{OT}$ in data does in fact resemble that based on the shear-driven KHI induced mixing flow evolution shown in figure 1. In a companion paper, Scotti (Reference Scotti2015) argued that in shear-driven mixing $L_T\sim L_O$, whereas in convectively driven mixing (where the underlying source of TKE is APE), $L_T$ can actually overestimate mixing. However, we would argue such a clear and binary distinction between the two types of mixing is not entirely appropriate. A simple mixing event due to a shear instability can be dominated by a cascade of ‘secondary’ shear instabilities growing on shear instabilities, for example as in an energetic estuary as shown in Geyer et al. (Reference Geyer, Lavery, Scully and Trowbridge2010), or it can be more appropriately characterized as being convectively driven once the primary billow has rolled up, as observed in the thermocline by Woods (Reference Woods1968). We argue that if a distinction is allowed in principle between the mixing length and the turbulent integral scale, as done here, at least some ‘convective’ mixing can fit within the same parameterization as ‘shear-driven’ overturnings; oceanic datasets appear to support this argument. Needless to say, our argument does not extend to truly convection-driven mixing such as that in deep convection zones, but rather to patches of mixing where it is reasonable to think of relatively isolated time-varying overturning mixing events where the stratification is sufficiently subcritical so as not to suppress such overturning.

Caldwell (Reference Caldwell1983) provided a nice description of the observed data on $R_{OT}$ on physical grounds and summarized the opposing views that became a source of controversy: on one hand, Gibson (Reference Gibson1987) argued that most observations were what he referred to as fossilized turbulence, while Gregg (Reference Gregg1987) argued that they were of active turbulence. Caldwell proposed that the largest eddies (of scale $L_O$) feed on the original overturn (of scale $L_T$). He provided scaling arguments, observational support, and other references, to show that $L_O\sim L_T$ in this growth phase. This was also pointed out by Itsweire et al. (Reference Itsweire, Koseff, Briggs and Ferziger1993). This is consistent with the idealized picture presented in figure 1, as well as with simulations of other non-shear-instability-like turbulent overturns (e.g. see Fritts et al. Reference Fritts, Bizon, Werne and Meyer2003; Chalamalla & Sarkar Reference Chalamalla and Sarkar2015). Caldwell argued that while the fossilization phase is also captured in the data, most of the observed overturns are in a state of ‘continuous production’.

This continuous production corresponds to the intermediate energetically mixing phase identified in § 3. As the mixing during this phase seems both to be most efficient and also to be associated with $L_O \sim L_T$, we refer to this phase as being the Goldilocks mixing phase, where the turbulent overturning is neither too ‘hot’ ($L_T \gg L_O$) nor too ‘cold’ ($L_T \ll L_O$) but ‘just right’ ($L_O \sim L_T$). As argued in Mashayek et al. (Reference Mashayek, Caulfield and Peltier2017a), in this Goldilocks phase, the energy is most efficiently supplied to turbulence as the external driving (associated with the primary KHI billow roll up) is at a scale which essentially matches the upper bound of an assumed inertial subrange largely unaffected by the stratification for which the (dynamic) scalar mixing of buoyancy is subordinate to the turbulence cascade.

5.6. Quasi-equilibrium and self-organized criticality

It has been suggested (Turner Reference Turner1973; Sherman, Imberger & Corcos Reference Sherman, Imberger and Corcos1978; Salehipour, Peltier & Caulfield Reference Salehipour, Peltier and Caulfield2018; Smyth, Nash & Moum Reference Smyth, Nash and Moum2019; Smyth Reference Smyth2020) that various kinds of shear-induced geophysically relevant turbulence are in a state of self-organized criticality (SOC) in which the system is marginally unstable and overturns spontaneously emerge when the flow properties intermittently drop below a stability criterion with respect to a critical Richardson number, i.e. in the sense discussed in § 1, the stratification intermittently becomes (just) subcritical, self-organizing and relaxing back towards criticality. Although Howland, Taylor & Caulfield (Reference Howland, Taylor and Caulfield2018) demonstrated that KHI are exceptionally weak when the minimum Richardson number approaches the classical Miles–Howard (Howard Reference Howard1961; Miles Reference Miles1961) criterion of $1/4$, calling into question the appropriateness of describing an energetically turbulent flow as being ‘marginally unstable’ with respect to that specific value, as proposed by Thorpe & Liu (Reference Thorpe and Liu2009), a body of evidence suggests that stratified turbulent flows can indeed adjust towards a state where characteristic values of $Ri \simeq 0.17\text {--}0.25$ (Zhou, Taylor & Caulfield Reference Zhou, Taylor and Caulfield2017; Portwood, de Bruyn Kops & Caulfield Reference Portwood, de Bruyn Kops and Caulfield2019; van Reeuwijk, Holzner & Caulfield Reference van Reeuwijk, Holzner and Caulfield2019). As originally argued by Turner (Reference Turner1973) and Sherman et al. (Reference Sherman, Imberger and Corcos1978), this evidence is at least suggestive that shear-driven turbulent flows adjust to a ‘kind of equilibrium’ with $Pr_T \sim O(1)$, and characteristic values of $Ri$ close to the linear stability critical value, perhaps fortuitously.

We combine these ideas to propose an appropriate value for the scaling coefficient ‘$A$’ in (5.8). Keeping in mind that the $R_{OT} \ll 1$ and $R_{OT} \gg 1$ limits of (5.8) represent the young and decaying or fossilizing limiting phases of a turbulent mixing event, $R_{OT}\sim 1$ thus represents the intermediate Goldilocks mixing phase. Let us assume that such a mixing phase is in the ‘kind of equilibrium’ proposed by Turner, which in the modern language of physics is essentially equivalent to a state of SOC. However, it might be described, emergence of the balance $R_{OT}$ on physical grounds is consistent with the empirical fact that the majority of observed overturns appear to be in this phase of efficient, and indeed ‘optimal’ mixing, as has been discussed previously (Dillon Reference Dillon1982; Ferron et al. Reference Ferron, Mercier, Speer, Gargett and Polzin1998; Mater et al. Reference Mater, Venayagamoorthy, St. Laurent and Moum2015; Mashayek et al. Reference Mashayek, Caulfield and Peltier2017a).

Therefore, we assume that the flow has a characteristic critical Richardson number $Ri_{cr}$ close to, but perhaps somewhat less than $1/4$, i.e. that the stratification of the flow is just slightly subcritical. Although it is perhaps stating the obvious, it important to appreciate that such a Richardson number should be considered as being significantly less than one, and to be associated with a flow where it is not to be expected that the turbulence will be dominated by the stratification. For such a subcritically stratified flow at $R_{OT}\approx 1$, (5.8) reduces to

(5.9)\begin{equation} \varGamma_{gold}=\frac{A}{2}. \end{equation}

Combining (2.5) and (5.3a,b),

(5.10)\begin{equation} Pr_T = \frac{P}{[{-}B]} \frac{N^{2}}{S^{2}} = \frac{Ri}{Ri_f}, \end{equation}

defining the ‘flux’ Richardson number $Ri_f$ (Turner Reference Turner1973). Over sufficiently long averages and ignoring transport terms, we can assume $Ri_f$ is equivalent to a particular definition of a mixing efficiency, since $P \simeq {\mathcal {M}} +\epsilon$ and $-B \simeq {\mathcal {M}}$. Therefore

(5.11)\begin{equation} Ri_f \simeq \frac{\varGamma_{\mathcal{M}}}{1+\varGamma_{\mathcal{M}}} , \end{equation}

(see Ivey & Imberger Reference Ivey and Imberger1991; Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Salehipour & Peltier Reference Salehipour and Peltier2015) Therefore,

(5.12)\begin{equation} \varGamma_{\mathcal{M}} = \frac{\dfrac{Ri}{Pr_T}}{1- \dfrac{Ri}{Pr_T}}. \end{equation}

Since we are always assuming that the mixing of scalar is locked to the mixing of momentum so that $Pr_T\sim 1$, $Ri_f \sim Ri \lesssim 1/4$, and so $\varGamma _{\mathcal {M}} \lesssim 1/3$. As we shall see, the underlying assumption that the flow is ‘subcritically stratified’ but time evolving, leading to these various scalings, must always be remembered when comparing with previous studies, particularly those of Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2019) and Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) based around the turbulent Froude number $Fr_T = {\epsilon }{/[N k]}$.

If we further assume that the Goldilocks optimally efficient mixing occurs at some specific $Ri_{cr}$ rather than over some range bounded above by $Ri_{cr}$, consistently with the flow being in a self-organized critical equilibrium, we can then express the scaling parameter ‘$A$’ as

(5.13)\begin{equation} A=\frac{2\dfrac{Ri_{cr}}{Pr_t}}{1-\dfrac{Ri_{cr}}{Pr_t}}. \end{equation}

Evidence from a range of flows suggests that vigorous turbulence in a subcritically stratified flow (i.e. in the specific sense that the flow has sufficiently small $Ri$ that the turbulence, characterized by overturnings, can be sustained throughout the flow) has $Pr_T \simeq 1$ (Zhou et al. Reference Zhou, Taylor and Caulfield2017; Portwood et al. Reference Portwood, de Bruyn Kops and Caulfield2019; van Reeuwijk et al. Reference van Reeuwijk, Holzner and Caulfield2019). There is also evidence that appropriate values of $Ri_{cr}$ range from $0.16$ (as reported in Portwood et al. Reference Portwood, de Bruyn Kops and Caulfield2019) through $0.21$ (Zhou et al. Reference Zhou, Taylor and Caulfield2017; van Reeuwijk et al. Reference van Reeuwijk, Holzner and Caulfield2019) to the canonical ‘marginally stable’ value of $1/4$ (Salehipour et al. Reference Salehipour, Peltier and Caulfield2018; Smyth et al. Reference Kaminski and Smyth2019; Smyth Reference Smyth2020). Indeed, due to a variety of reasons, in particular associated with the presence of ambient turbulence, the critical value of the Richardson number associated with linear instability may be expected to be less than this canonical value of $1/4$ (Smyth et al. Reference Smyth, Moum and Caldwell2001; Thorpe, Smyth & Li Reference Thorpe, Smyth and Li2013; Kaminski & Smyth Reference Kaminski and Smyth2019).

Using these expressions $Pr_T \simeq 1$ and $1/6 \lesssim Ri_{cr} \lesssim 1/4$ suggests that $2/5 \lesssim A \lesssim 2/3$, i.e. $1/5 \lesssim \varGamma _{\mathcal {M}} \lesssim 1/3$, with the canonical lower value being associated with $Ri_{cr} =1/6$ and $Pr_T =1$, consistently with the direct calculations of Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019) in a continuously forced flow, and also the (apparently) most efficient mixing in an evolving Kelvin–Helmholtz unstable shear flow (Mashayek et al. Reference Mashayek, Caulfield and Peltier2013). Alternatively $\varGamma _{\mathcal {M}}=0.2$ can also be identified with $Ri_{cr}=1/4$, $Pr_T=5/4$, which is also largely consistent with available data, as one might argue that a continuously forced flow should correspond to the intermediate mixing phase of an evolving overturning event. As we will soon show, for the data considered in this work, ‘$A$’ lies within $2/5 \lesssim A \lesssim 2/3$ for individual and combined datasets. This implies that the above assumptions about $Pr_T$, $Ri_{cr}$ are reasonable.

At the risk of belabouring the point, it is also important to appreciate two fundamental aspects of our argument concerning this optimal mixing phase with intermediate values of $R_{OT} \sim O(1)$ between the two asymptotic regimes of ‘young’ turbulence and the decaying ‘fossilization’ phase. First, although this regime is indeed intermediate, it still occurs when the flow has relatively small characteristic Richardson numbers significantly less than one, with crucially $Pr_T \simeq 1$, and so also $Ri_f \simeq Ri$. Therefore, the mixing of buoyancy is set by the mixing of momentum, and the density is mixed essentially as a passive scalar, with buoyancy playing an at most sub-dominant dynamical role. Specifically, and critically for our modelling, the properties of the turbulence itself, in particular characteristic length scales and time scales, are largely unaffected by stratification in this intermediate mixing phase within our framework. Second, the time dependence of the flow evolution is central to our argument, and this phase naturally follows the ‘young’ (and crucially vigorously overturning) phase of the flow evolution, and precedes the later decaying or fossilizing phase. Fundamentally, at no stage in the flow evolution do buoyancy effects dominate any aspect of the flow evolution.

To illustrate our proposed parameterization, figure 2 shows the parameterization in (5.8) plotted for $A=2/3$ along with the asymptotic limits given by (5.6) and (5.7). Also shown is the naturally equivalent definition for an (instantaneous) mixing efficiency, i.e.

(5.14)\begin{equation} \mathcal{E}_i \equiv \frac{\varGamma_{\mathcal{M}}}{1+\varGamma_{\mathcal{M}}}.\end{equation}

This definition of mixing efficiency tends to one in the limit of a non-turbulent flow due to non-zero molecular diffusion (at $Pr=1$) but the complete absence of turbulent dissipation, hence yielding a perfectly efficient (laminar) mixing, but importantly, no turbulent mixing. Yet again, it is important to remember that the flows of interest should never be thought of as strongly stratified, and so $R_{OT} \rightarrow 0$ because $\epsilon \rightarrow 0$ in the numerator of the Ozmidov length $L_O$ as there is no turbulence at all, not because $N$ is very large in the denominator, thus suppressing the turbulence due to buoyancy forces. In the limit of $R_{OT}\rightarrow \infty$, $\mathcal {E}_i$ tends to zero as would be expected in decaying turbulence, and $R_{OT}\rightarrow\infty$ since the Thorpe length $L_T$ is decaying more rapidly than $L_O$. Finally, at the heart of the intermediate Goldilocks phase, $R_{OT}=1$, $\mathcal {E}_i=A/(2+A)$.

Figure 2. Variation with ${R_{OT}}$ of the parameterization (5.8) (with $A=2/3$, implying $\varGamma _{\mathcal {M}}=1/3$ and $\mathcal {E}_i=1/4$ at $R_{OT}=1$). Also shown are the two asymptotic scalings for growing turbulence (5.6) and for decaying turbulence (5.7) with thin grey lines, and the corresponding mixing efficiency, with a dashed red line. The light blue box in the middle of the plot represents the bounding box within which the ocean turbulence data considered in this work lie, as will be shown in figure 3. Conceptually, time for a particular shear-driven mixing event runs from left to right as $R_{OT}$ increases.

Figure 3(a–f) shows an excellent agreement between the parameterized form of $\varGamma _{param}=\varGamma _{\mathcal {M}}$ using (5.8) and measurements of $R_{OT}$ and the observationally inferred $\varGamma _{data}=\varGamma _\chi = \chi /\epsilon$ constructed from measurements of $\chi$ and $\epsilon$ ‘directly’ from the six datasets introduced earlier. Each panel also includes the $\varGamma =A\,R_{OT}^{-1}$ and $\varGamma =A\ R_{OT}^{-4/3}$ limits. These two are expected to be relevant only in the limits of $R_{OT} \ll 1$ and $R_{OT}\gg 1$, as was shown in figure 2. However, the datasets in figure 3 lie near the Goldilocks limit $R_{OT}\sim 1$ (the blue box in the middle of figure 2). Thus, the two limiting cases are, perhaps unsurprisingly, inaccurate in the range relevant to the data whereas their joint scaling, as formulated by (5.8), much better represents the data. This is particularly clear when all the data are combined in (g) where a total of $\sim$50 000 turbulent patches provide sufficient statistics to provide a fair comparison between the data and the parameterization (note that at least a few tens of thousands of samples are needed to capture the distribution statistically significantly, as discussed in Cael & Mashayek (Reference Cael and Mashayek2020)).

Figure 3. (a–g) Observationally inferred turbulent flux coefficient $\varGamma _{data}=\varGamma _\chi$ (constructed from measurements of $\chi$ and $\epsilon$) as a function of $R_{OT}$ for six oceanic datasets introduced in § 4, as well as for the combined data. In each panel, blue scattered dots in the background represent actual patches, the solid blue line with filled circles and error bars represents the data binned (in log($\varGamma _\chi$) space) and the solid red line with square symbols represents the parameterized $\varGamma _{param}=\varGamma _{\mathcal {M}}$ using (5.8). Error bars represent $\pm 1$ standard deviation. The value of ‘$A$’ for each panel is inferred based on linear regression of the corresponding data to log (5.8). The right side inset compares probability density functions (p.d.f.s) of $\varGamma _{data}$ vs $\varGamma _{param}$ and the top inset shows the p.d.f. of data in log($R_{OT}$). The solid straight black lines represent $\varGamma =A\,R_{OT}^{-1}$ and $\varGamma =A\,R_{OT}^{-4/3}$. Panel (h) shows p.d.f.s of $\varGamma _{data}/\varGamma _{param}$ for the various datasets.

As discussed earlier in connection with (5.13), for $Pr_T \simeq 1$ and $1/6 \lesssim Ri_{cr} \lesssim 1/4$, we obtain $2/5 \lesssim A \lesssim 2/3$. As shown in the captions of individual panels in figure 3, the range of values for ‘$A$’ obtained based on regression of individual datasets to the parameterization (5.8) falls within this range except for TIWE which is data-limited from a statistical standpoint but nevertheless gives an ‘$A$’ very close to the upper bound. The corresponding ‘$A$’ for the combined data is almost exactly $2/3$ which corresponds to $Pr_T \simeq 1$ and $Ri_{cr} \simeq 1/4$. Of course, this close agreement might well be somewhat fortuitous considering we only use six oceanic datasets. But given the wide range of processes, geographical locations and oceanic depths spanned by these datasets, it is pleasing that their mixing can be captured so well using a parameterization entirely based on physical grounds: (5.8) depends on the Kolmogorov theory, the mixing length theory and the TKE budget, while ‘$A$’ is determined based on a classical theoretical prediction for the critical Richardson number necessary for shear-induced instability. No empiricism, tuning, or filtering/modification of the data was required to obtain the agreement. Indeed, allowing some empiricism associated with numerical simulations such as those reported in Deusebio, Caulfield & Taylor (Reference Deusebio, Caulfield and Taylor2015), Zhou et al. (Reference Zhou, Taylor and Caulfield2017), Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019) and van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019) so that the critical Richardson number $Ri_{cr} \sim 0.2$ would nevertheless only modify the result slightly.

Crucially, this agreement also suggests that ocean mixing likely occurs at scales sufficiently small that the local gradient Richardson number falls below a critical value so that mixing in such a state remains close to and perhaps slightly below marginally unstable, and fundamentally, scalings for the properties of the turbulence are never dominated by the ambient stratification. The turbulence itself may be generated at a scale close to the shear instability, through for example mixing in the thermocline through the breakdown of KHI. Alternatively, the initial generation may be at a much larger scale through (for example) tidal generation, which through linear and nonlinear processes, leads to generation of shear at sufficiently small scales at which local $Ri$ can become sufficiently small (Nikurashin & Legg Reference Nikurashin and Legg2011). In other words, we hypothesize that the classic paradigm of a parallel shear flow, as idealized as it undoubtedly is, may plausibly lie at the heart of ocean interior mixing.

To compare the data and parameterization more quantitatively, figure 3(h) shows the histograms of the ratio of $\varGamma$ from data and the parameterized $\varGamma$. The agreement is remarkable given the range of dominant driving mechanisms for the different datasets: shear instability triggered by tropical instability waves in the equatorial undercurrent (TIWE), near-inertial-shear (FLX91), internal tides (BBTRE) and hydraulically controlled flow over abyssal canyons and sills (DoMORE and Samoan Passage). For the latter two cases, the agreement is somewhat surprising considering that the measured turbulence is close to the seafloor and possibly partially convectively driven, although even such convectively breaking overturns are shear driven in the sense that the background shear is a main source of energy. Nevertheless, it is reasonable to expect an influence of another length scale, which is absent in the construction of (5.8) for near boundary turbulence. Perhaps this explains the larger spread in the data for the Samoan Passage than predicted by (5.8). Nevertheless, the agreement remains good at the peak of the p.d.f. of the data where $R_{OT}\sim 1$.

5.7. Comparison to $Fr_T$-based parameterizations

Here, we have been focussed on mixing events which can be related to shear-driven processes with subcritical stratification in the specific sense that vigorous overturning events can be observed. As cogently argued by Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2019), building on the insights of Ivey & Imberger (Reference Ivey and Imberger1991), there are several reasons why it is appropriate to attempt to construct parameterizations in terms of the turbulent Froude number $Fr_T$ defined (using our notation) as

(5.15)\begin{equation} Fr_T \equiv \frac{\epsilon}{N k}. \end{equation}

They observed (from the results of body-forced numerical simulations) that $\varGamma \propto Fr_T^{-2}$ for large $Fr_T$, which is also consistent with gravity current observations as discussed in Wells, Cenedese & Caulfield (Reference Wells, Cenedese and Caulfield2010).

Of course, large $Fr_T$ corresponds to ‘weak’ stratification, and so it is entirely expected that the unstratified scaling (5.1) applies. Therefore, remembering also the definition of $L_O$ (2.1a–c), appropriately large $Fr_T$ is proportional to

(5.16)\begin{equation} Fr_T \propto \left ( \frac{L_O}{L_i} \right) ^{2/3}.\end{equation}

Considering this weakly stratified class of flows, Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2019) (effectively) assumed that $Pr_T \sim 1$, $\varGamma _\chi = \varGamma _{\mathcal {M}}$ and that the mixing length $L_m \sim L_i$ (using our notation), so that (2.1a–c) becomes

(5.17)\begin{equation} \varGamma_{\mathcal{M}} \propto \left ( \frac{L_i}{L_O} \right) ^{4/3} \propto Fr_T^{{-}2}. \end{equation}

Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) extended the discussion of Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2019), constructing and interpreting various Froude number scalings in terms of time scales and length scales, motivated by constructing a parameterization that would be practically useful in terms of oceanic measurements. In the high $Fr_T$ limit, they also argued (in our notation) that $L_i \sim L_T$, and so that $\varGamma _{\mathcal {M}} \propto Fr_T^{-2} \propto R_{OT}^{-4/3}$. This corresponds to the final decaying fossilization phase described above, though neither Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2019) nor Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) considered time dependence of mixing events (by design). In contrast to our inherently time-dependent viewpoint, their argument was that $L_T < L_O$ was a generic characteristic of weakly stratified turbulent flows, whereas we postulate that this scaling only occurs late in the flow evolution of a shear-driven overturning mixing event, and in point of fact, our arguments only require the stratification to be subcritical (and thus allow shear-driven overturnings) with sufficiently small but not necessarily negligible characteristic Richardson number.

Furthermore, although not couched in terms of shear flows by Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2019), the $Fr_T^{-2}$ scaling can also be interpreted in terms of (sufficiently small) Richardson numbers, as the natural equivalent scaling is that

(5.18)\begin{equation} Fr_T \sim \frac{U_0}{L_i N} \propto Ri^{{-}1/2} ,\end{equation}

and so

(5.19)\begin{equation} \varGamma_{\mathcal{M}} \propto Fr_T^{{-}2} \propto Ri . \end{equation}

This is of course entirely as expected for flows with $Pr_T \sim 1$, as then $\varGamma _{\mathcal {M}}\lesssim Ri_f \simeq Ri \lesssim 1/4$.

Interestingly, Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) also presented scaling arguments for two other (effectively assumed to be quasi-steady) classes of flows, neither of which should be construed as being equivalent to the (inherently time-dependent) phases of a subcritically stratified flow as considered here. Building on the observations of Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2019), which in turn were motivated by experimental observations dating back (at least) to Kato & Phillips (Reference Kato and Phillips1969), (see also Park, Whitehead & Gnanadeskian Reference Park, Whitehead and Gnanadeskian1994; Oglethorpe, Caulfield & Woods Reference Oglethorpe, Caulfield and Woods2013; Olsthoorn & Dalziel Reference Olsthoorn and Dalziel2015) Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) presented a scaling argument for the observed constant flux coefficient $\varGamma _{\chi } \propto Fr_T^{0}$ in (at least apparently) strongly stratified flows as $Fr_T \rightarrow 0$. They argued that the strong stratification completely dominates the dynamics, so that both the turbulent dissipation rate $\epsilon$ and the buoyancy flux $B$ (and hence $\chi$ as the flow is quasi-steady) scale with $w^{2} N$, where $w$ is the (perturbation) vertical velocity. Therefore $\varGamma _\chi$ is constant in very strong stratification. Clearly, this class of flows is not considered accessible by the (vertical) shear-driven mixing events with subcritical stratification allowing overturnings to reach significant amplitude on which we are focusing in this paper.

Furthermore, as argued in more detail in Caulfield (Reference Caulfield2021), care must be taken in inferring the existence of this class of strongly stratified and yet vigorously turbulent flows. The motivating experimental data were typically extracted from flows with emergent density staircases, where relatively well-mixed and deep ‘layers’ are separated by relatively thin ‘interfaces’ of enhanced density gradient. Indeed, as demonstrated by Portwood et al. (Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016), flows with low $Fr_T$ (when averaged over relatively large scales) typically exhibit significant spatio-temporal variability, with vigorous turbulence (characterized by large local values of $\epsilon$) being restricted to patches of locally significantly weaker stratification. It is also necessary to be cautious when comparing scaling arguments from different studies, as there are often significant differences in the way in which superficially equivalent non-dimensional parameters are defined.

In particular, as discussed in detail in Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), very strongly stratified turbulent flows at sufficiently high Reynolds number can be profoundly anisotropic, with characteristic horizontal length scales significantly larger than characteristic vertical scales of turbulent patches, which in turn are larger than the Ozmidov scale (see also Caulfield (Reference Caulfield2020) for further discussion). Therefore, the particular length and velocity scales which are used to construct Froude numbers must be chosen with care and consistency before direct comparisons between different flows (and indeed scaling arguments) are made.

Connecting this (assumed) strongly stratified class of flows as $Fr_T \rightarrow 0$ and the weakly stratified class with $Fr_T \gg 1$, Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) argued for the existence of an intermediate ‘moderately stratified’ class of flows, where the properties are effectively hybrid between the two end members. They argued that the buoyancy variance destruction rate should scale as in the strongly stratified limit so that $\chi \propto w^{2} N$, while the turbulent kinetic energy (and its dissipation rate) are both (largely) unaffected by the stratification, so $k \sim w^{2}$, and so, in this class of flows, they argued that

(5.20)\begin{equation} \varGamma_\chi \propto \frac{kN}{\epsilon} = Fr_T^{{-}1} .\end{equation}

Furthermore, they argued that this moderately stratified class is also characterized by $Fr_T \sim O(1)$, which finally allowed them to argue that this intermediate class should have the length scaling $\varGamma _\chi \propto R_{OT}^{-1}$ as when $Fr_T \sim O(1)$ it is expected that $R_{OT} \sim O(1)$ as well.

The presented evidence of this intermediate class is perhaps not as convincing as for the two end members. Whatever the quality of the data fit, it is important to appreciate that the arguments leading to this scaling are fundamentally different from our arguments leading to the (superficially same) scaling (5.6) for the early-time ‘young’ phase of an inherently time-dependent and subcritically stratified shear-driven mixing event. Furthermore, our ‘subcritically stratified’ intermediate Goldilocks phase is also in marked conceptual contrast to the intermediate, moderately stratified class of flows parameterized by Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019), in the sense that the physical processes and balances should be thought of as arising for fundamentally different reasons.

In their argument, the stratification (and hence the characteristic values of the buoyancy frequency) in their intermediate class has an order-one effect on the mixing. That can be seen perhaps most clearly in terms of the implications of their scaling arguments for $Pr_T$. Since they argue that $\chi \propto w^{2} N$,

(5.21)\begin{equation} \kappa_T \propto \frac{\chi}{N^{2}} \propto \frac{w^{2}}{N }.\end{equation}

On the other hand, since they argue that the TKE $k \propto w^{2}$ is largely unaffected by stratification, it is natural to suppose that the mixing length for momentum is given by (5.1), and so

(5.22)\begin{equation} \nu_T \propto \left (\frac{w^{3}}{\epsilon} \right ) w \rightarrow Pr_T \propto \frac{N w^{2}}{\epsilon} \propto Fr_T^{{-}1} ,\end{equation}

using the scaling (5.21) for $\kappa _T$.

It could perhaps be argued that since in this regime $Fr_T \sim 1$, $Pr_T \sim 1$ still using this scaling. Nevertheless, a consistent interpretation of this moderately stratified class of flows leads to the prediction that the dynamics of the mixing of the buoyancy (as described by $Pr_T$) does depend on the stratification, with in particular $Pr_T$ increasing as the stratification becomes stronger, and hence $Fr_T$ decreases. Importantly such a variation in $Fr_T$ can occur straightforwardly for different quasi-steady forced flows within their assumed class of flows.

This is a qualitatively different mixing behaviour from the mixing behaviour we argue occurs during the intermediate (in time) Goldilocks mixing phase. Our argument implicitly assumes that the stratification is sufficiently subcritical that $Pr_T \simeq 1$ throughout the mixing life cycle, and in particular during the intermediate Goldilocks phase. During this phase, it appears that the stratification is in some sense self-organized close to a critical value of the Richardson number, and so there is not expected to be any parameter dependence of $Pr_T$: the flow is after all neither too hot nor too cold, but just right. In summary, their argument assumes an intermediate quasi-steady class of flows with hybrid dynamical balance between buoyancy and turbulence, while our intermediate (in time) phase is still (and always) subcritically stratified. We argue that it is the middle phase of an inherently time-dependent flow evolution, that happens after a primary shear-driven overturning breaks down, but before it enters a final fossilization decay phase.

Indeed, caution must always be shown in drawing any connections between descriptions in terms of $Fr_T$ and $Ri$, particularly as $Ri \simeq 1$ is actually really a quite strongly stratified flow in the context of shear-driven mixing events. Of course, there is no necessity for equality in the various scalings used in the construction of the parameterizations, and the presence of an order-one constant can always allow $Ri$ to be sufficiently small to be subcritical in the sense used in this paper, i.e. for a shear-driven overturning to develop to sufficiently large amplitude to undergo vigorous turbulent breakdown, which we believe is actually the dominant process far from boundaries.

In summary, while (5.7) was proposed previously, to our knowledge (5.5), (5.6) and thereby (5.8) are new. Furthermore, the previously proposed scaling relation only strictly applies for $R_{OT}\gg 1$ which corresponds to a small fraction of the data. Through proposing a scaling for the $R_{OT}\ll 1$ limit and combining the two limits, we have managed to propose a comprehensive parameterization and also to determine the coefficient of proportionality on physical grounds. Note that earlier works used (5.7) with an adjustable coefficient that could be tuned to fit the data better. Figures 2 and 3 suggest that such a tuning is incorrect, and the data ‘feel’ both (5.6) and (5.7). In the next section, we will further discuss the significance of (5.6) for bulk estimates of mixing on regional and global scales.

6.1. The devil is in the tails

An important question, one mostly left out of the literature discussing properties of various definitions of the turbulent flux coefficient $\varGamma$, is the extent to which the nuances of variations in $\varGamma$ with other parameters matter in practice. There is no simple answer to this question, as it depends on the application for which the rate of mixing is being measured. Furthermore, a key missing piece is how one would connect the fluid-dynamical understanding of mixing of turbulent patches (as studied herein and in articles cited throughout), to an appropriate time- and space-averaged rate of mixing. The latter is out of the scope of our paper (but some discussion is provided in Cael & Mashayek Reference Cael and Mashayek2020). Since our primary motivation behind this work is to relate efficiency of mixing to transformation of water masses in the ocean interior, which depends on the vertical divergence of the buoyancy flux (see Mashayek et al. Reference Mashayek, Ferrari, Nikurashin and Peltier2015; Ferrari et al. Reference Ferrari, Mashayek, McDougall, Nikurashin and Campin2016), we attempt to answer the question by simply quantifying the extent to which various measures of $\varGamma$ will affect the total buoyancy flux obtained by summing the flux due to all the overturns in the datasets that we consider here.

Taking $\varGamma \epsilon$ (for various definitions of $\varGamma$) as a proxy for the buoyancy flux $\kappa N^{2}$ through the Osborn balance (2.6), in figure 4(a) we plot the flux based on $\varGamma _{data}=\varGamma _{\chi }$ from the data, parameterizations and as constants. While we of course appreciate the many questionable assumptions underlying the approximation of actual flux with $\varGamma \epsilon$, as discussed in Mashayek & Peltier (Reference Mashayek and Peltier2013) and Mashayek et al. (Reference Mashayek, Caulfield and Peltier2013) for example, they do not affect our arguments here in any significant way. As expected based on figure 3, flux parameterized using (5.8) agrees well with the observed flux, while flux parameterized using (5.7) is biased high. Most importantly, if we consider the total flux due to all overturns, given the several orders of magnitude that $\epsilon$ spans, the tail of the histograms makes a significant contribution and so the extent to which $\varGamma$ is accurately estimated for the tails is key. As shown in panel (b), it is the higher $\epsilon$ tail that leads to an overestimate when using parameterization (5.7). It is also notable that, while $\varGamma =0.2$ underestimates the total flux, $\varGamma =1/3$ actually leads to a much better agreement. Since $\varGamma =1/3$ is linked (at least within our modelling approach) to assuming that the critical Richardson number $Ri_{cr} \simeq 1/4$, as was discussed above, it seems a more appropriate choice in the simplest approach of using a single constant value of $\varGamma$.

Figure 4. (a) Value of $\varGamma \times \epsilon$, considered a proxy for the buoyancy flux through the Osborn relation (2.6), calculated based on $\varGamma _{data}$, $\varGamma _{param}$ using our full parameterization (5.8) and $\varGamma _{-4/3}$ using only the asymptotic scaling (5.7), as well as based on use of constant values $\varGamma =0.2$ and $\varGamma =1/3$. The last choice is based on (5.13) if we consider $Ri_{cr}=1/4$. (b) Cumulative flux for the cases of panel (a), all normalized by the total $\varGamma _{data} \epsilon _{data}$ both obtained directly from data. (c–e) Bivariate histogram plots for the data over $\varGamma -\epsilon$, $\varGamma _{data}/\varGamma _{param}-\epsilon$ and $R_{OT}-\epsilon$ parameter spaces. All plots are made for the data combined for all six experiments in table 1.

Figure 4 plots all six datasets combined. The disagreement between various estimates, however, can be much larger when looked at case by case. To do so, in table 1 we compare the total flux based on various estimates (i.e. the high end of the curves in figure 4b) for individual datasets. It is clear that parameterization (5.8) is quite accurate in capturing the total flux for five of the six cases (within 10 %). Even for the Samoan Passage case, which is somewhat of an outlier, it is still within a factor of 3. On the other hand, parameterization (5.7) is less accurate. $\varGamma =0.2$ underestimates the total flux by as much as a factor of four, and fixing $\varGamma =1/3$ does a much better job.

Table 1. Comparison of the inaccuracies between the total buoyancy fluxes using the Osborn approximation $\varGamma \epsilon$ estimated based on constant values and various other parameterizations of $\varGamma$.

It is useful to interpret figure 4(a,b) and the agreement of data with parameterization (5.8) by looking at the distribution of data over various parameter spaces, as shown in (c–e) of the figure. While the data are spread around $R_{OT}\sim 1$, the spread is much larger at low $\epsilon$, representing both the early stages of turbulence growth and the later stages of the turbulence after significant decay, as one would expect from consideration of figure 1. For energetic patches that correspond to the optimal Goldilocks mixing phase, $R_{OT}$ nicely collapses around the value $R_{OT} \simeq 1$ as shown in (e). For this Goldilocks phase, $\varGamma \sim 1/3$ as shown in (c) and in agreement with the assumption of $Ri_{cr}\sim 1/4$, while at early and later stages of turbulence life cycles, $\varGamma$ can be very large or very small, respectively. Panel (d) shows that the parameterization (5.8) manages to collapse this spread on $R_{OT}$, and that importantly it works well for the tail of the histograms which make the largest contribution to the total flux. In summary, given the accuracy, simplicity, and physically justified nature of the parameterization (5.8), it seems like a natural way to model the turbulent flux coefficient $\varGamma$ associated with overturns in shear-induced turbulence.

6.2. Bulk $\varGamma$: the importance of young turbulence events

In ocean and climate general circulation models, small scale wave breaking is not resolved but the combined power available for mixing from tides, winds and other sources is available as a bulk value per grid cell. For example, an energetic turbulent zone in the ocean, such as that shown in figure 5(a) for the Samoan Passage, is roughly represented by 1–2 grid cells (in each direction) in a typical climate model. A bulk $\varGamma$ is thus required to divide the power available (from tides, winds, currents, etc.) into bulk measures of mixing and dissipation for a given grid cell (see Cimoli et al. (Reference Cimoli, Caulfield, Johnson, Marshall, Mashayek, Naveira Garabato and Vic2019), for a thorough discussion). In practice, this $\varGamma _{Bulk}$ is often taken to be 0.2. Even without all the variations and caveats associated with $\varGamma _i$ for individual turbulent patches, $\varGamma _{Bulk}$ inevitably depends on the statistics of turbulent patches within each grid cell. Several studies have applied parameterizations for $\varGamma$, derived based on patch data, to coarse grid-based calculations to establish the leading-order importance of variations in $\varGamma$ for the larger scale ocean circulation (De Lavergne et al. Reference De Lavergne, Madec, Le Sommer, Nurser and Garabato2016; Mashayek et al. Reference Mashayek, Salehipour, Bouffard, Caulfield, Ferrari, Nikurashin, Peltier and Smyth2017c; Cimoli et al. Reference Cimoli, Caulfield, Johnson, Marshall, Mashayek, Naveira Garabato and Vic2019). That approach, while illuminating, is also not correct since it applies a patch-based quantity to grid cells of $O(100\ \mathrm {km})\times O(100\ \mathrm {km}) \times O(100\ \mathrm {m})$ in size. In this subsection we show how $\varGamma _{Bulk}$ depends on statistics of turbulent patches and that not only does it depend on Goldilocks mixing of energetic patches, but equally on the young patches which possess large $\varGamma$ and decaying patches with vanishing $\varGamma$.

Figure 5. (a) An example of an energetic oceanic turbulence zone, in the Samoan Passage, to demonstrate the intermittency of turbulence; reproduced from Alford et al. (Reference Alford, Girton, Voet, Carter, Mickett and Klymak2013): northward velocity (colours), potential temperature (black contours) and dissipation rate measured by the velocity microstructure profiles (shaded black profiles) and from Thorpe scales (blue profiles). The scale for $\epsilon$ is given at the lower left and redrawn above each profile. (b) A snapshot of vertical shear from an observationally forced and tuned wave-resolving simulation of the Drake Passage, borrowed from the work of Mashayek et al. (Reference Mashayek, Ferrari, Merrifield, Ledwell, St. Laurent and Garabato2017b). (c–f) Bivariate histograms showing the distribution of patches over the $\epsilon -\varGamma$ parameter space for four of the datasets considered herein. The values of bulk $\varGamma$ in each panel's title is obtained from (6.1). Note that the data in (e) are associated with (a).

Imagining a typical energetic zone, such as the example in figure 5(a), we assume the domain to contain $n$ number of turbulent patches (within a grid cell of a coarse resolution model which would represent the region) and over a period associated with a coarse resolution model time stepping. Panel (b) shows a snapshot from a wave-resolving high resolution realistic simulation to demonstrate the intermittency of vertical shear generated by surface winds and flow–topography interaction in the very energetic Drake Passage in the Southern Ocean where high flow speeds occur due to the constriction of the Antarctic Circumpolar Currents and the overlying strong westerly winds known as the roaring forties. Even in such an energetic region, the flow is mostly non-turbulent (from a small scale stratified turbulence perspective – a rich mesoscale and submesoscale turbulence field exists throughout the domain). However, in the ‘quiet’ regions, background weak mixing almost always exists: $97\,\%$ of the data acquired from $\sim$750 full-depth microstructure profiles from 14 field experiments considered in Cael & Mashayek (Reference Cael and Mashayek2020) possessed a diffusivity larger than $10^{-6} \ \textrm {m}^{2}\ \textrm {s}^{-1}$ (also see Waterhouse et al. Reference Waterhouse2014). In the limit of completely laminar flow, $\kappa \rightarrow O(10^{-7})\ \textrm {m}^{2}\ \textrm {s}^{-1}$ and so the flux $\mathcal {M}\approx \kappa N^{2}$ stays finite while $\epsilon$ tends to zero. Thus, while for energetic turbulent patches $\varGamma \approx \varGamma _{gold}$, for quieter regions it tends to $O(10)$ and larger values (see histograms in figure 3).

To show how the energetic, young and decaying turbulent patches collectively set the bulk $\varGamma$, we define

(6.1)\begin{equation} \varGamma_{Bulk}=\frac{\mathcal{M}_{tot}}{\epsilon_{tot}}\approx \frac{\sum\limits_{i=1}^{n} \varGamma_i \times \epsilon_i}{\sum\limits_{i=1}^{n} \epsilon_i} , \end{equation}

where $n$ represents the number of patches in a region of interest. Figure 5(c–f) shows the estimated $\varGamma _{Bulk}$ for four deep energetic turbulent datasets considered herein. As discussed in figure 4, for energetic patches with high $\epsilon$, $\varGamma _i \rightarrow \varGamma _{gold}$ while small $\epsilon$ corresponds to young and decaying patches which possess large and small $\varGamma _i$, respectively. For BBTRE and DoMORE, abundance of young patches result in a large $\varGamma _{Bulk}$ even though most of the data correspond to Goldilocks mixing. For IH18, a larger fraction of the data correspond to Goldilocks mixing and so does $\varGamma _{Bulk}$. A simple thought experiment can help illustrate the impact of young turbulent patches. Imagine a box with three turbulent patches: a young patch ($\epsilon =0.001$, $\varGamma =100$), a moderately turbulent patch ($\epsilon =0.1$, $\varGamma =10$) and an energetic patch ($\epsilon =1$, $\varGamma =1/3$). This results in $\varGamma _{Bulk}=0.44$.

Our key message is that while most of the collective community efforts have focused on finding the right parameterization for energetic turbulent patches, the net mixing also critically depends on young patches and the large fraction of the water column where only weak background turbulence occurs. This highlights the importance of the new scalings presented in this work in (5.5), (5.6). It also highlights the need for investing efforts in understanding better the statistics of turbulence to quantify its intermittency. Such statistics are key to connecting our understanding of physics of small scale stratified turbulence, as discussed in this work and the ones cited herein, to the larger scale applications. This statistical bridge is arguably highly underdeveloped – see Cael & Mashayek (Reference Cael and Mashayek2020) for a discussion. Finally, while we used the patch data for four experiments in figure 5, it was merely to illustrate the simple point that $\varGamma _{Bulk}$ can be large even though most patches have $\varGamma \sim \varGamma _{gold}$. We have no reason to believe that the number of patches within each experiment was statistically sufficient to quantify mixing accurately in the geographical region of each experiment. Also, the data we used are only for the already identified patches. The actual microstructure profiles also include regions that are not identified as patches, but should contribute to $\varGamma _{Bulk}$.

6.3. Statistics of $L_O/L_T$: active vs fossil turbulence

Thus far we have repeated on several occasions and on empirical grounds (based on the data used in this study and other works cited herein) that ‘most of the data are centred around $R_{OT}\sim 1$’. This, however, lies at the heart of a historically significant debate. As discussed in the introduction, Gibson (Reference Gibson1987) argued that most observations were of fossilized turbulence, while Gregg (Reference Gregg1987) argued that they were of active turbulence. The latter, referred to as ‘continuous production’ by Caldwell (Reference Caldwell1983) and as Goldilocks mixing herein, is supported by a seemingly universal concentration of data around ${R_{OT}\sim 1}$.

Figure 6(a) shows the histogram of the combined data considered in this work to highlight the clustering around $R_{OT}\sim 1$. Importantly, the clustering improves for more energetic turbulent patches with larger $\epsilon$. Panels (b,c) of the figure demonstrate, based on data from direct numerical simulations of shear instabilities, that once the parameter ranges approach those of oceanic turbulent patches (i.e. large $Re$ and subcritical $Ri$), a larger fraction of the total overturning life cycle corresponds to sustained energetic turbulence (for which $R_{OT} \sim 1$). This explains the statistical distribution of $R_{OT}$ in panel (a). In contrast, less energetic turbulence events, such as those created in laboratory or numerical simulations, or those observed in lakes, are expected to have a more broadly distributed (over $R_{OT}$) distributions. Once again, we argue that $R_{OT} \sim 1$ is associated with an intermediate phase of an energetic time-dependent subcritically stratified mixing event, not an intermediate, hybrid class of moderately stratified quasi-steady mixing events as argued by Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019).

Figure 6. (a) Probability density function of $R_{OT}$, grouped into quartiles of $\epsilon$, from the combined six oceanic datasets considered in this work. (b,c) Temporal fraction of turbulence life cycle as well as the ratio of $R_{OT}$ during turbulent phase of the flow to its mean value over the whole life cycle. Each symbol represents a life-cycle-averaged quantity from a direct numerical simulation for a given pair of Reynolds and Richardson numbers. All cases in (b) are for $Ri=0.12$ while all cases in (c) are for $Re=6000$. Turbulent phase of the life cycles are defined as the times when $Re_b>20$. Panels (b,c) are produced from analysis of simulations originally discussed in Mashayek & Peltier (Reference Mashayek and Peltier2011a,Reference Mashayek and Peltierb) and Mashayek & Peltier (Reference Mashayek and Peltier2013); Mashayek et al. (Reference Mashayek, Caulfield and Peltier2013).