2.1. Definitions

We first split the total local specific kinetic energy of the flow $K(\boldsymbol {x},t) \equiv (1/2) \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {u} = \bar {K} + K'$ into a mean and a turbulent (or perturbation) kinetic energy, respectively,

(2.1a,b)\begin{equation} \bar{K}(y,z) \equiv \tfrac{1}{2} \bar{\boldsymbol{u}} \boldsymbol{\cdot} \bar{\boldsymbol{u}} \quad \textrm{and} \quad K'(\boldsymbol{x},t) \equiv \tfrac{1}{2} \boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{u}' , \end{equation}

where we recall that ${\boldsymbol{x}}, t, {\boldsymbol{u}}, \rho$ are respectively the `shear-layer rescaled' spatial coordinates, time, velocity, and density (see Part 1, (3.3)), that the bar averages are $\bar {\cdot } \equiv \langle \cdot \rangle _{x,t}$, and that the prime variables are perturbations with respect to these $x-t$ averages (see Part 1, (4.1)).

By analogy, we also define the total scalar density variance $K_\rho \equiv (1/2) Ri_b^s \, \rho ^2 = \bar {K}_\rho + K'_\rho$ into a mean and a turbulent (or perturbation) scalar variance, respectively,

(2.2a,b)\begin{equation} \bar{K}_\rho(y,z) \equiv \tfrac{1}{2} Ri_b^s \, \bar{\rho}^2 \quad \textrm{and} \quad K'_\rho(\boldsymbol{x},t) \equiv \tfrac{1}{2} Ri_b^s \, \rho'^2, \end{equation}

where $Ri_b^s$ is the shear-layer bulk Richardson number defined in Part 1 (3.4b). These variances are useful and more direct and convenient alternatives to potential energies when estimating mixing. In particular $\bar {K}_\rho$ is more informative in SID flows than in most canonical stratified shear layers since the average density field $\bar {\rho }$ results entirely from mixing inside the duct, rather than being set as an initial condition. No mixing, i.e. the bimodal $\pm 1$ distribution from the external reservoirs, corresponds to a maximum $(1/Ri_b^s)\langle \bar {K}_\rho \rangle = 1/2$ (and $K'_\rho =0$). By contrast, complete mixing (uniform $\bar {\rho }=0$) corresponds to a minimum $(1/Ri_b^s) \langle \bar {K}_\rho \rangle = 0$, and a linear stratification with uniform gradient across the shear layer $\partial _z \bar {\rho } =-1$ corresponds to an intermediate value of $(1/Ri_b^s)\langle \bar {K}_\rho \rangle = 1/3$ (also note that in this latter case the non-dimensional buoyancy frequency simplifies to $N^2=Ri_b^s$).

2.2. Evolution equations

The averaged equations of $\bar {K}, \bar {K}_\rho,$ and the temporal evolution equations of $K',K'_\rho$ follow from the equations of motion (3.5) in Part 1

(2.3a)$$\begin{gather} \overline{\partial_t K}(y,z) = \varPhi^{\bar{K}} - \mathcal{P} + \mathcal{F} - \bar{\epsilon}, \end{gather}$$

(2.3b)$$\begin{gather}\partial_t K'(\boldsymbol{x},t) = \varPhi^{K'} + \mathcal{P} - \mathcal{B} - \mathcal{E}, \end{gather}$$

(2.3c)$$\begin{gather}\overline{ \partial_t K_\rho}(y,z) = \varPhi^{\bar{K}_\rho} - \mathcal{P}_\rho, \end{gather}$$

(2.3d)$$\begin{gather}\partial_t K'_\rho (\boldsymbol{x},t) = \varPhi^{K'_\rho } + \mathcal{P}_\rho - \chi , \end{gather}$$

where the mean temporal gradients over the data set interval $t \in [0, L_t]$ have the form $\overline {\partial _t K}= (1/L_t)(\langle K \rangle _x(t=L_t)-\langle K \rangle _x(t=0)) \approx 0$ in quasi-steady state (similarly for $\bar {K}_\rho$). All $\varPhi$ terms are transport terms that will be discussed in § 2.4.

The mean kinetic energy equation (2.3a) has three source/sink terms: the production of turbulent kinetic energy $\mathcal {P}$ (generally positive) by interaction of the off-diagonal (deviatoric) Reynolds stresses with the mean shear, the gravitational forcing term $\mathcal {F}$ (generally positive) transferring energy from the mean potential energy (not shown here) and the viscous dissipation of the mean $\bar {\epsilon }$ (always positive)

(2.4a–c)\begin{equation} \mathcal{P} \equiv{-} \overline{u'v'} \partial_y \bar{u} - \overline{u'w'} \partial_z \bar{u}, \quad \mathcal{F} \equiv Ri_b^s \, \sin \theta \, \bar{u}\bar{\rho} , \quad \bar{\epsilon} \equiv \frac{2}{Re^s} \bar{s}_{ij} \bar{s}_{ij}>0, \end{equation}

where $Re^s$ is the shear-layer Reynolds number defined in Part 1 (3.4a), the mean strain rate tensor is $\bar {s}_{ij} \equiv (\partial _{x_i} \bar {u}_j + \partial _{x_j} \bar {u}_i)/2$ and we implicitly sum over repeated indices (unless specified otherwise).

The remaining equations (2.3b)–(2.3d) have four further volumetric terms: the turbulent buoyancy flux $\mathcal {B}$ (transferring energy from the turbulent kinetic energy, generally positive), the production of turbulent scalar variance $\mathcal {P}_\rho$ (generally positive), the turbulent dissipation $\mathcal {E}$ (always positive) and the turbulent scalar dissipation $\chi$ (always positive)

(2.5a–d)\begin{align} \mathcal{B} \equiv Ri_b^s \, \overline{w'\rho'}, \quad \mathcal{P}_\rho \equiv{-} Ri_b^s \, \overline{w' \rho'} \partial_z \bar{\rho}, \quad \mathcal{E} \equiv \frac{2}{Re^s} s'_{ij} s'_{ij}>0, \quad \chi \equiv \frac{Ri_b^s}{Re^s \, Pr} \partial_{x_j} \rho'\partial_{x_j} \rho'>0, \end{align}

where $s'_{ij} \equiv (\partial _{x_i} u'_j + \partial _{x_j} u'_i)/2$. All terms in (2.4) and (2.5) are functions of $y,z$ only, except for $\mathcal {E}$ and $\chi$, which are functions of $\boldsymbol {x},t$.

We see in (2.5) that $\mathcal {P}_\rho$ is proportional to $\mathcal {B}$ in the simple case of linear stratification. Moreover, $\mathcal {P}_\rho =\mathcal {B}$ if $\partial _z \bar {\rho }=-1$ (linear mixing layer spanning the entire shear layer), since in this case $\mathcal {B}$ is a source term for the turbulent potential energy, which is exactly equal to $K'_\rho$ (as noted by Taylor et al. Reference Taylor, de Bruyn Kops, Caulfield and Linden2019, § 3).

2.3. Approximations

A few simplifying approximations were made in (2.3)–(2.5). First, in (2.3c) we neglected the molecular scalar dissipation $Ri_b^s/(Re^s Pr) (\partial _{x_j} \bar {\rho }\partial _{x_j} \bar {\rho }) \approx 0$ (requiring $|\partial _{x_j} \bar {\rho }|\ll \sqrt {Re^s Pr/Ri_b^s}$ which is true here for $Pr=700$). Second, in the definition of $\mathcal {P}$ we assumed parallel mean flow, i.e. $\bar {v},\bar {w} \approx 0$, and in $\mathcal {P}_\rho$ we assumed $\partial _y\bar {\rho } \approx 0$ (which are good approximations). Third, in $\mathcal {F}$ we assumed no mean vertical buoyancy flux, i.e. $\bar {w} \bar {\rho } \approx 0$ (this term is key in horizontal exchange flows at $\theta =0$, but negligible in long ducts at $\theta >0$ since the mean slope of the density interface is small, as explained in LPL19, § 4.3). Fourth, in $\mathcal {B}$ we assumed $\overline {u'\rho '},\overline {v'\rho '} \approx 0$ and $\cos \theta \approx 1$.

2.4. Boundary fluxes

The transport terms $\varPhi$ in (2.3) represent the divergence of advective, pressure and viscous/molecular fluxes

(2.6)\begin{gather} \left.\begin{gathered} \varPhi^{\bar{K}} \equiv{-} \overline{\partial_x (u K)} - \overline{\partial_x (u p)} + \frac{2}{Re^s} \overline{\partial_x(u s_{11}}), \quad \varPhi^{K'} \equiv \partial_{x_i} \left( - u'_i K' - u'_i p' + \frac{2}{Re^s} u'_j s'_{ij} \right), \\ \varPhi^{\bar{K}_\rho} \equiv{-} \overline{\partial_x(uK_\rho)}, \qquad \varPhi^{K'_\rho} \equiv \partial_{x_i} \left(- u'_i K'_\rho + \frac{1}{Re^s \, Pr} \partial_{x_i} K'_\rho\right), \end{gathered}\right\} \end{gather}

respectively, where, as above, we assumed parallel mean flow and negligible molecular transport in $\varPhi ^{\bar {K}}$ and $\varPhi ^{\bar {K}_\rho }$, and where $\overline {\partial _x \phi }$ denotes the mean gradient along $x$ (non-zero if $\phi$ is non-periodic). When averaged over a volume, these divergence terms become boundary fluxes.

These boundary fluxes are typically neglected in the stratified turbulence literature, because they usually conveniently vanish in idealised geometries (e.g. for periodic boundary conditions), greatly simplifying (2.3). In the SID geometry, they are unfortunately slightly more complicated as we explain below.

In the $y$ and $z$ directions, $\varPhi ^{K'}$ and $\varPhi ^{K'_\rho }$ will not generally cancel if the volume average is done over the shear layer (as in this paper) because the boundaries do not include the duct walls (whereas LPL19, § 4.2.1 included them). In other words, turbulent fluctuations can in principle be transported freely across our shear layer ‘imaginary’ boundary ($y=\pm L_y, z=\pm 1$) to (or, more rarely, from) the near-wall region.

More importantly, in the $x$ direction, most boundary fluxes can generally be neglected when $\theta \gtrsim \arctan A^{-1}$, where $A\equiv L/H$ is the length-to-height aspect ratio of the duct (high in the long ducts of interest here, $A=30$ in our set-up). In these so-called forced flows, the mean slope of the density interface is small and the flow is approximately periodic (see LPL19, § 4.3 and their Appendix B). This applies in particular to $-\overline {\partial _x (uK)}$, which is important and ${<}0$ when $\theta \approx 0$, but unimportant and ${\approx }0$ in forced flows.

We, however, note two exceptions. First, our Part 1 results on the unexpected nature of the estimated mean pressure gradient $\varPi = -\overline {\partial _x p}$ (weakening $\bar {u}$ rather than strengthening it) suggest that the simple hydrostatic pressure assumed in LPL19's Appendix B may not be correct and, consequently, that $-\overline {\partial _x (u p)}$ may not be neglected. However, for simplicity, and due to our inability to measure it directly, we ignore it in this paper until future work sheds light on it. Second, the flux of mean scalar variance $\varPhi ^{\bar {K}_\rho }>0$ represents the continuous inflow of unmixed fluid from the reservoirs, countering the effects of mixing, and must be retained to ensure that a steady state for $\bar {K}_\rho$ is possible.

2.5. Steady-state balances

In our unsteady flows, steady state ($\partial _t = 0$) cannot be expected in the pointwise and instantaneous sense of (2.3). It can, however, be expected in a time- and volume-averaged sense, leading to the following balances:

(2.7a)$$\begin{gather} 0 \approx \langle \mathcal{F}\rangle - \langle \mathcal{P} \rangle - \langle \bar{\epsilon} \rangle, \end{gather}$$

(2.7b)$$\begin{gather}0 \approx \langle \mathcal{P}\rangle - \langle \mathcal{B} \rangle - \langle \mathcal{E} \rangle, \end{gather}$$

(2.7c)$$\begin{gather}0 \approx \langle \varPhi^{\bar{K}_\rho} \rangle - \langle \mathcal{P}_\rho \rangle , \end{gather}$$

(2.7d)$$\begin{gather}0 \approx \langle \mathcal{P}_\rho \rangle - \langle \chi \rangle, \end{gather}$$

where we recall from Part 1 that $\langle \cdot \rangle \equiv \langle \cdot \rangle _{x,y,z,t}$. In the above, we assumed for simplicity that all boundary fluxes were negligible, except the essential $\varPhi ^{\bar {K}_\rho }$ sustaining the steady-state scalar dissipation. We also assumed that all mean temporal gradients are negligible $\langle \partial _t \phi \rangle _t = (1/L_t)(\phi (L_t)-\phi (0)) \approx 0$ (verified in our data).

These above balances can alternatively be expressed as two ‘independent’ estimations of the mean turbulent dissipation rates $\langle \mathcal {E}\rangle,\langle \chi \rangle$

(2.8a)\begin{align} \langle \mathcal{E} \rangle &\approx \langle \mathcal{P}\rangle - \langle \mathcal{B} \rangle \end{align}

(2.8b)\begin{align} &\approx \langle \mathcal{F}\rangle - \langle\mathcal{B}\rangle - \langle \bar{\epsilon}\rangle , \end{align}

(2.8c)\begin{align} \langle \chi \rangle &\approx \langle \mathcal{P}_\rho \rangle \end{align}

(2.8d)\begin{align} &\approx \langle \varPhi^{\bar{K}_\rho} \rangle . \end{align}

Equation (2.8a) represents the classical balance of Osborn (Reference Osborn1980), (2.8c) represents the classical balance of Osborn & Cox (Reference Osborn and Cox1972), while (2.8b) and (2.8d) are more specific to SID flows.

3.1. Time and volume averages

3.1.1. Energy reservoirs

In figure 1 we plot the steady-state energy reservoirs in all 16 data sets, both as function of $Re^s$ (panels a–f), and as correlation plots (panels g–j).

Figure 1. Steady-state energy reservoirs in all 16 data sets. (a–c) Mean and turbulent kinetic energies, and their square-root ratio, as function of $Re^s$ (separating flow regimes). (d–f) Mean and turbulent scalar variances and their square-root ratio (rescaled from (2.1a,b) to obtain $\langle \bar {\rho }^2 \rangle$, $\langle \rho '^2 \rangle$). (g–j) Correlation between scalar variance and kinetic energies (in (i-j) we show the full $\bar {K}_\rho, K'_\rho$ to test for potential-to-kinetic energy partitions). Symbol shapes and colours follow the flow regimes, as in Part 1. All dashed lines have slope one. Dotted lines are labelled explicitly.

Note that our definition of turbulent perturbations around the $x-t$ mean flow can attribute artificially high energies to ${\rm L}$ and ${\rm H}$ flows, whose perturbations $u'\equiv u - \bar {u}$ and $\rho ' = \rho -\bar {\rho }$ can exhibit slight residual $x-t$ structure due to the nature of our exchange flow (slightly non-parallel in $x$ and/or accelerating or decelerating in $t$). Therefore, in figure 1 we removed this artefact (not due to the turbulent or wave motions or interest) by subtracting from $\langle K'\rangle$, $\langle K'_\rho \rangle$ the mean $x-t$ variance corresponding to the zero $x$-wavenumber and temporal frequency content of their respective spectra (i.e. we subtracted from $\langle u'^2\rangle _{x,y,z,t}$, the components $\langle \langle u'\rangle ^2_x\rangle _{y,z,t}$ and $\langle \langle u'\rangle ^2_t\rangle _{x,y,z}$ and similarly for $\rho '$). We verified that ${\rm I}$ and ${\rm T}$ flows are almost unaffected by this correction. We return to this in our discussion of energy spectra in § 3.3 and Appendix A.3.

The mean kinetic energy $\langle \bar {K} \rangle$ (panel a) is approximately constant around $0.2$ in all flows, with values decreasing from $0.25$ in ${\rm L}, {\rm H}$ flows to $0.15$ in ${\rm T}$ flows. The turbulent kinetic energy $\langle K' \rangle$ (panel b) increases from $0$ in ${\rm L}$ flows to around $0.01$ in ${\rm T}$ flows (e.g. T2, T3). The square-root ratio of turbulent-to-mean kinetic energies $\sqrt {\langle K' \rangle /\langle \bar {K} \rangle }$ (panel c), indicating the relative magnitude of velocity fluctuations, is 10 %–15 % in ${\rm H}$ flows, and 10 %–20 % in ${\rm I}$ flows (with significant spread) and up to 25 % in ${\rm T}$ flows.

We now turn to the scalar variance reservoirs. Although $\bar {K}_\rho, K'_\rho$ are preferred when discussing energy fluxes (as in § 2) because of their interpretation as proxy for potential energy under linear stratification, we first consider in panels (d–h) the rescaled quantities $\langle \bar {K}_\rho \rangle / Ri_b^s \equiv \langle \bar {\rho }^2 \rangle$ and $\langle K'_\rho \rangle /Ri_b^s \equiv \langle \rho '^2 \rangle$, which are more straightforward measures of scalar variance. High values of mean variance $\langle \bar {\rho }^2 \rangle \approx 0.4$ (panel d) confirm that very little mixing takes place in ${\rm L}$ and ${\rm H}$ flows beyond molecular diffusion (close to the no-mixing upper bound of $0.5$, see dotted lines). Mixing increases in ${\rm I}$ flows, where an intermediate layer of approximately uniform density achieves ‘more’ mixing than a uniformly linear stratification ($\langle \bar {\rho }^2 \rangle <1/3$, see dotted lines), while ${\rm T}$ flows are halfway between linear and full mixing ($\langle \bar {\rho }^2 \rangle \approx 1/6$). The turbulent variance $\langle \rho '^2 \rangle$ (panel e) is, surprisingly, higher in some ${\rm H}$ flows than in most ${\rm I}$ and ${\rm T}$ flows. This reflects the fact that Holmboe waves on a sharp interface can generate very large perturbations on either side of it (due to high $|\rho '|$ values), compared with a well-mixed turbulent layer (low $|\rho |$ values). This effect partially disappears when considering the relative square-root of turbulent-to-mean variance $\sqrt {\langle K'_\rho \rangle /\langle \bar {K}_\rho \rangle } = \sqrt {\langle \rho '^2 \rangle /\langle \bar {\rho }^2 \rangle }$ (panel f), typically higher in ${\rm I}$ and ${\rm T}$ flows, and reaching a maximum of 25 %, just like the kinetic energies (panel c).

We further see that the mean scalar variance $\langle \bar {\rho }^2 \rangle$ is closely correlated to the mean kinetic energy $\langle \bar {K} \rangle$ (panel g), especially in ${\rm I}$ and ${\rm T}$ flows, where they become equal (dashed line). This reflects our observation in Part 1 of self-similarity $\langle \bar {u} \rangle _y(z) \approx \langle \bar {\rho } \rangle _y(z)$ in ${\rm T}$ flows, i.e. that momentum and density become equally mixed. This general correlation in ${\rm I}$ and ${\rm T}$ flows also extends to the turbulent energies $\langle \rho '^2 \rangle$ and $\langle K' \rangle$ (panel h).

We now turn to the potential-to-kinetic energy partitions. The mean partition $\langle \bar {K}_\rho \rangle /\langle \bar {K} \rangle$ (panel i) drops from ${\approx } 1$ (equipartition, see dashed line) in L1 and H1 (where $Ri_b^s \approx 0.5\text{--} 1$) down to ${\approx }0.1$ (1/10 partition, see dotted line) in the late ${\rm I}$ and ${\rm T}$ regimes (where $Ri_b^s \approx 0.1-0.2$). The turbulent partition $\langle K'_\rho \rangle / \langle K' \rangle$ (panel j) follows a similar trend of equipartition in all ${\rm H}$ flows, and asymptotic $1/10$ partition in ${\rm T}$ flows (see the zoomed-in inset for more details).

3.1.2. Energy fluxes

In figure 2 we plot the steady-state energy fluxes in all 16 data sets. All gradients were computed using second-order accurate finite differences. Limitations in the resolution of the data will be discussed in § 3.3 and Appendix B.

Figure 2. Steady-state energy fluxes in all 16 data sets. (a) Mean kinetic energy forcing (power source) as function of $\theta Ri^s_b$. (b–c) Mean kinetic energy dissipation, and turbulent kinetic energy dissipation (power sinks) as function of $\theta Re^s$ ($\approx \theta Re^h$ identified as proxy for regime transitions in LPL19). (d–f) Test of the approximate kinetic balances (2.7a), (2.8a), (2.8b), respectively. (g–i) Test of the approximate scalar variance balances (2.7c), (2.8c), (2.8d), respectively. (j–l) Test of three further commonly used ratios: $\langle \mathcal {B} \rangle / \langle \mathcal {E} \rangle \equiv \varGamma$, $\langle \mathcal {B} \rangle / \langle \mathcal {P} \rangle \equiv R_f$, and $\langle \mathcal {B} \rangle / \langle \mathcal {P}_\rho \rangle$ (${\equiv }1$ when $\partial _z \bar {\rho }=1$), respectively. All dashed lines have slope 1 and denote expected equality between fluxes. Dotted lines are labelled explicitly.

In panels (a–c) we investigate the dependence of the kinetic energy source $\langle \mathcal {F} \rangle$ and sinks $\langle \bar {\epsilon } \rangle, \langle \mathcal {E}\rangle$ with respect to two key groups of parameters $\theta Ri_b^s$ and $\theta Re^s$, respectively. Note that $\theta$ is in radians, and recall from Part 1 (see figure 2) how these output parameters depended on input parameters: $Ri_b^s \propto \theta ^{-0.9}(Re^h)^{-0.4}$ and $Re^s \propto \theta ^{0.7}(Re^h)^{1.4}$. As expected from its definition (2.4a–c), $\langle \mathcal {F}\rangle \propto \theta Ri_b^s$, with a factor $\approx 0.5$ in ${\rm L}$ and ${\rm H}$ flows, decreasing to $\approx 0.25$ in ${\rm T}$ flows (due to a lower $\langle \bar {u}\bar {\rho }\rangle$). The mean dissipation $\langle \bar {\epsilon }\rangle$ dominates over the turbulent dissipation $\langle \mathcal {E} \rangle$ at low $\theta Re^s$ (${\rm L}$ and ${\rm H}$ flows), but decreases to become comparable or lower at higher $\theta Re^s= O(100)$ (T2 and T3). These observations in panels (a–c) are key – and almost defining – features of SID flows: hydraulic control of two-layer exchange flows sets an upper bound on the magnitude of the mean flow (set by the dimensional scale $\sqrt {g'H}$ and thus $\langle |\bar {u}| \rangle \lesssim 1/2$ or $\langle \bar {K} \rangle \lesssim 1/4$), causing a plateau in $\langle \bar {\epsilon }\rangle$ in the ${\rm I}$/${\rm T}$ regimes, and thus an increase in $\langle \mathcal {E} \rangle$, which eventually dominates to match the increased $\langle \mathcal {F} \rangle$ at higher $\theta$ (see ML14 and LPL19).

In panels (d–f) we test the approximate kinetic energy balances of (2.7a), (2.8a), (2.8b), respectively. The mean balance $\langle \mathcal {P} \rangle \approx \langle \mathcal {F} \rangle - \langle \bar {\epsilon } \rangle$ is only verified (dashed line) in a subset of flows (e.g. H1, H2, H4 I8, T2, T3). The systematic underestimation of $\langle \mathcal {P} \rangle$ is due partly to the neglected boundary flux $\langle \varPhi ^{\bar {K}} \rangle$, and partly to our limited resolution of small-scale fluctuations (which are needed to measure $\langle \mathcal {P} \rangle$ but not $\langle \mathcal {F} \rangle$ and $\langle \bar {\epsilon } \rangle$). Unfortunately, it is not possible to estimate with confidence the relative importance of either source of error, since boundary fluxes are notoriously inaccurate and we cannot know how much fluctuation energy is present below our resolution. The turbulent balance of Osborn (Reference Osborn1980) $\langle \mathcal {E} \rangle \approx \langle \mathcal {P} \rangle - \langle \mathcal {B} \rangle$ is also verified in a (different) subset of flows. The general underestimation of $\langle \mathcal {E} \rangle$, especially in ${\rm I}$ and ${\rm T}$ flows, is primarily due to the limited resolution of gradients of small-scale velocity fluctuations (needed to measure $\langle \mathcal {E} \rangle$ but not $\langle \mathcal {P} \rangle$ and $\langle \mathcal {B} \rangle$). The balance $\langle \mathcal {E} \rangle \approx \langle \mathcal {F} \rangle - \langle \mathcal {B} \rangle - \langle \bar {\epsilon } \rangle$ follows from the previous two balances, and is thus the most poorly verified overall.

In panels (g–i) we test the approximate scalar variance balances (2.7c), (2.8c), (2.8d), respectively. The balance between production of turbulent variance and advective flux of mean variance (from unmixed fluid coming into the domain) $\langle \mathcal {P}_\rho \rangle \approx \langle \varPhi ^{\bar {K}_\rho } \rangle$ (panel g) is verified in most flows (e.g. H2, H4, T3 and most ${\rm I}$ flows except I4), although the cluster near $0$ is inconclusive. Some ${\rm H}$ flows (H2 and H4) even show equality between negative values, which suggests that: (i) the net effect of Holmboe wave turbulence in the measurement volume is to increase (rather than decrease) scalar variance, by sharpening (rather than broadening) the mean density interface, consistent with the findings of Zhou et al. (Reference Zhou, Taylor, Caulfield and Linden2017), Salehipour, Caulfield & Peltier (Reference Salehipour, Caulfield and Peltier2016) and our Reynolds-averaged profiles in Part 1; and/or (ii) this sharpening must be countering the net advection of mixed fluid into the volume, which means that mixing must take place outside the length of the duct occupied by Holmboe waves, presumably near the ends of the duct where plumes discharge turbulently into the reservoirs and interact with the incoming fluid, entraining mixed fluid back into the duct. Negative values of $\langle \varPhi ^{\bar {K}_\rho } \rangle <0$ in I4, T1 and T2 are, however, surprising and likely the result of experimental noise in the computation of this mean gradient. The turbulent balance of Osborn & Cox (Reference Osborn and Cox1972) $\langle \chi \rangle \approx \langle \mathcal {P}_\rho \rangle >0$ (panel h), only valid for broadening-type (${\rm I}$ and ${\rm T}$) flows (because of the neglect of $\langle \varPhi ^{K'_\rho } \rangle$), cannot be verified even in these flows. In most ${\rm H}$ flows, this balance is fundamentally impossible since $\langle \mathcal {P}_\rho \rangle <0$. The systematic and severe underestimation of $\langle \chi \rangle$ is due to our severely limited resolution of small-scale density gradients (more severe than for $\langle \mathcal {E} \rangle$, because $\rho '$ contains energetic length scales that are approximately a factor $\sqrt {Pr} \approx 25$ smaller than $\boldsymbol {u}'$). Finally, the balance $\langle \chi \rangle \approx \langle \varPhi ^{\bar {K}_\rho } \rangle >0$ (panel i) follows from the previous two balances and is thus equally poorly verified. We explain the reasons for these limitations in § 3.4.

In panels (j–l), we test the correlation of $\langle \mathcal {B} \rangle$ with the three other turbulent fluxes $\langle \mathcal {E} \rangle, \langle \mathcal {P} \rangle, \langle \mathcal {P}_\rho \rangle$ respectively, in order to assess the relevance and numerical value of the following ratios:

(3.1a–c)\begin{equation} \frac{ \langle \mathcal{B} \rangle}{ \langle \mathcal{E} \rangle } \equiv \varGamma, \quad \frac{\langle \mathcal{B} \rangle}{ \langle \mathcal{P} \rangle} \equiv R_f, \quad \frac{\langle \mathcal{B} \rangle}{ \langle \mathcal{P}_\rho \rangle} \equiv 1\ \text{when} \ \partial_z \bar{\rho}={-}1 \ \text{everywhere}. \end{equation}

The flux parameter $\varGamma$ and the flux Richardson number $R_f$ date back to Osborn (Reference Osborn1980) and have been extensively used in the literature to parameterise the ‘taxation rate’ of stratification on turbulent dissipation (Caulfield Reference Caulfield2020). Although often assumed constant, dimensional analysis suggests that $\varGamma$ and $R_f$ are functions $(\theta,Re^s,Ri^s_b,R,Pr)$ until proven otherwise. First, our data show that $\langle \mathcal {B} \rangle \propto \langle \mathcal {E} \rangle$ only in late ${\rm I}$ flows and in all ${\rm T}$ flows (panel j), where the slope indicates an asymptotic ratio $\varGamma \approx 0.1$ (dotted line), approximately half the commonly used value of 0.2 in the literature. The slightly negative values of $\langle \mathcal {B} \rangle$ can be explained by the slight non-periodicity of exchange flows at low tilt angles $0<\theta \lesssim \arctan A^{-1} \approx 1/30 \approx 2^\circ$: the convective acceleration of each layer ($u'\partial _x u'>0$) caused by a tilting interface produces downward flow $(w'<0)$ in the dense layer ($\rho '>0$) and vice versa, resulting in a net volume-averaged $\langle \mathcal {B}\rangle = Ri_b^s \langle w'\rho '\rangle <0$ in the absence of turbulence. This effect vanishes in more turbulent flows at larger tilt angles, where we instead tend to slightly overestimate $\varGamma$ by our underestimation of its denominator $\langle \mathcal {E} \rangle$ (compared with its numerator $\langle \mathcal {B} \rangle$, due to limitations in our computation of small-scale gradients). Second, we see that $\langle \mathcal {B} \rangle \propto \langle \mathcal {P} \rangle$ in most ${\rm I}$ and ${\rm T}$ flows (panel k), where the slope indicates an asymptotic ratio $R_f \approx 0.05$ (dotted line), approximately a third of the commonly used value of 0.15 in the literature. Third, we see that $\langle \mathcal {B} \rangle \approx \langle \mathcal {P}_\rho \rangle$ (dashed line) in most ${\rm I}$ and ${\rm T}$ flows (panel l), which is consistent with the theory under linear stratification (where $\partial _z \bar {\rho }=-1$), despite such a stratification being only achieved approximately in T3 (see Part 1, figure 3p). We return to these parameters in more detail in § 5.

3.1.3. Estimations of $\langle \mathcal {E}\rangle$ and $\langle \chi \rangle$ from non-dimensional parameters

In this section we combine the steady-state energy balances of § 2.5 and the experimental results of § 3.1.2 to propose indirect estimations (or proxies) of $\langle \mathcal {E}\rangle$ and $\langle \chi \rangle$ that are insightful and more accurate than their direct computations, which rely on small-scale gradients.

From (2.8a) and (3.1), we take advantage of the fact that $\mathcal {P}$ is measured with better accuracy than $\mathcal {E}$ to propose

(3.2a)\begin{align} \langle \mathcal{E}\rangle &\approx \langle \mathcal{P}\rangle - \langle \mathcal{B}\rangle \end{align}

(3.2b)\begin{align} &\approx (1-R_f) \langle\mathcal{P}\rangle \end{align}

(3.2c)\begin{align} &\approx \frac{1}{1+\varGamma} \langle \mathcal{P}\rangle, \end{align}

which means that

(3.3)\begin{equation} R_f \approx \frac{\varGamma}{1 + \varGamma} \quad \text{or} \quad \varGamma \approx \frac{R_f}{1 - R_f}. \end{equation}

These estimations depend on the balance (2.8a) and the assumption (3.1) that the fluxes $\langle \mathcal {E}\rangle, \langle \mathcal {B}\rangle, \langle \mathcal {P}\rangle$ are proportional to one another, approximately verified in ${\rm T}$ flows. Note, however, that our measurements gave slightly incompatible values of $\varGamma \approx 0.1$ and $R_f \approx 0.05$.

To address this, we first note that $\varGamma =0.1$ is mostly likely an overestimate due to the greater underestimation of its denominator $\langle \mathcal {E}\rangle$ (which relies on gradients of fluctuations which have higher energy content at higher, unresolved wavenumbers), compared with its numerator $\langle \mathcal {B}\rangle$ (which relies on fluctuations only). Figure 2(e–f) suggests that $\langle \mathcal {E}\rangle$ may be underestimated by as much as a factor of 2, bringing the actual value of $\varGamma$ closer to $0.05$. Second, we note that $R_f\approx 0.05$ (and therefore $\varGamma \approx 0.05/0.95 \approx 0.05$) appears more trustworthy because both its numerator $\langle \mathcal {B}\rangle$ and its denominator $\langle \mathcal {P}\rangle$ rely on fluctuations only, rather than gradients. Third, we add to these main comments a weaker nuance that $R_f=0.05$ might actually be a slight underestimate due to its numerator $\langle \mathcal {B}\rangle$ involving $\rho '$, whose energy spectrum extends further in wavenumber space at $Pr=700$ than its numerator $\langle \mathcal {P}\rangle$ involving $u'$ and $v'$, and is thus comparatively more poorly resolved (see § 3.1.4 for more details). In summary, we conclude that 0.05 and 0.1 are robust lower and upper bounds, i.e. $\varGamma \approx R_f \approx 0.05\text{--} 0.1$, with a greater confidence for the lower range $0.05\text{--} 0.07$.

From (2.8b), we take advantage of the fact that $\mathcal {F}$ is measured with even better accuracy than $\mathcal {P}$ to propose a series of further approximations of $\langle \mathcal {E}\rangle$ valid in the limit of very turbulent flows ($\theta Re^s \gg 100$)

(3.4a)\begin{align} \langle \mathcal{E}\rangle &\approx \langle \mathcal{F}\rangle - \langle \mathcal{B}\rangle - \langle \bar{\epsilon} \rangle \end{align}

(3.4b)\begin{align} &\approx \langle \mathcal{F}\rangle - \langle \mathcal{B}\rangle \quad \ \ \, \,\quad \qquad \text{if} \ \langle \mathcal{E}\rangle \gg \langle \bar{\epsilon} \rangle \ \text{(hydraulic control, figure 2 b--c)} \end{align}

(3.4c)\begin{align} &\approx (1-R_f) \langle\mathcal{F}\rangle \end{align}

(3.4d)\begin{align} &\approx 0.25 \,(1-R_f) \, \theta\, Ri_b^s \quad \,\text{if} \ \langle \mathcal{F} \rangle \approx 0.25 \, \theta \, Ri_b^s \ \text{(figure 2 a)} \end{align}

(3.4e)\begin{align} &\approx 0.037 \, (1-R_f) \, \theta \qquad \,\, \text{if} \ Ri_b^s \approx 0.15 \ \text{(Part 1 figure 2 b)} \end{align}

(3.4f)\begin{align} &\approx 0.035 \, \theta \qquad \qquad \qquad \text{if} \ R_f \approx 0.05 \ \text{(figure 2 k)} , \end{align}

where we recall that $\theta$ is in radians. Note that using the upper bound corresponding to $\varGamma \approx 0.1$ in the last line (3.4f) (see figure 2j) would give an almost identical expression $\langle \mathcal {E}\rangle \approx 0.034 \theta$.

From (2.8c), we propose the corresponding approximation of $\langle \chi \rangle$, in the limit of very turbulent flows with linear stratification where $\langle \mathcal {B}\rangle \approx \langle \mathcal {P}_\rho \rangle$ (figure 2l)

(3.5)\begin{equation} \langle \chi \rangle \approx R_f \langle \mathcal{F} \rangle \approx 0.037 \, R_f \, \theta \approx 0.0019\, \theta. \end{equation}

We also note that, under all the above assumptions, our estimations (3.4f) and (3.5) yield the following ratio of scalar variance to kinetic energy dissipation:

(3.6)\begin{equation} \frac{\langle \chi \rangle}{\langle \mathcal{E}\rangle } \approx \frac{\langle \mathcal{B} \rangle}{\langle \mathcal{P}\rangle -\langle \mathcal{B}\rangle} \approx \frac{R_f}{1-R_f}\approx \varGamma , \end{equation}

which, as we have seen, gives values between 0.05 and 0.1. This expression has the merit of linking $\varGamma,R_f$ with a natural measure of the irreversible ‘tax’ levied by stratification on turbulence. The key question becomes: How do $\varGamma,R_f$ scale with the non-dimensional flow parameters? We tackle this parameterisation of mixing in § 5.

Finally, we note that LPL19 explained the transitions between flow regimes by using the simple approximation $\langle \mathcal {E} \rangle \approx \langle \mathcal {P} \rangle \approx \langle \mathcal {F} \rangle \approx (h^2 \delta u ) \, \theta /8 \approx 0.04 \theta$ (where the factor $h^2 \delta u \approx 3$ converts their hydraulic non-dimensionalisation to our shear-layer non-dimensionalisation). Their expression is in good agreement with (3.4f). They argued that regime transitions are caused by thresholds in the normalised turbulent strain rate, which we write as

(3.7)\begin{equation} \langle s'_{ij} s'_{ij} \rangle \equiv \frac{Re^s}{2} \langle \mathcal{E} \rangle \approx 0.02 \, (1-R_f) \, \theta Re^s \propto \theta Re^s, \end{equation}

assuming $R_f =$ const., highlighting the key role of the group of parameters $\theta Re^s$.

The above data on mean energy reservoirs and fluxes confirm and extend LPL19's findings that flows with a similar product $\theta Re^s$ (but different individual values of $\theta$ and $Re^s$) behave similarly. Note that LPL19's hydraulic formulation used the product $\theta Re^h$ (where $Re^h$ is defined in Part 1, (3.2a)), while our more accurate shear-layer formulation uses the product $\theta Re^s \propto \theta ^{1.7} (Re^h)^{1.4}$.

Our data are also consistent with the findings in Part 1 that quantitative turbulent fractions scale strongly with both $\theta$ and $Re^s$ (enstrophy fraction $\propto \theta ^{2.7}(Re^s)^{2.8}$, and overturn fraction $\propto \theta ^{3.2}(Re^s)^{1.8}$). Since the production of perturbation enstrophy by vortex stretching is given by $s'_{ij}\,\omega '_{i}\omega '_j$, there is in fact a direct relation between an increasingly large turbulent strain rate $s'_{ij} s'_{ij}$ (slaved to $\theta Re^s$) and increasingly extreme enstrophy events, and thus enstrophy fraction (Johnson & Meneveau Reference Johnson and Meneveau2016). The relation to density overturns is more indirect; first because vorticity can be decomposed into a rotating and a shearing part (Tian et al. Reference Tian, Gao, Dong and Liu2018) (the rotating part being associated with overturns but not the shearing part), and second because overturns feed back into the enstrophy production through the baroclinic term.

3.1.4. Kolmogorov and Batchelor length scales

The estimation of the viscous dissipation of turbulent kinetic energy $\langle \mathcal {E} \rangle$ in (3.4f) allows us in turn to give a practical volume-averaged estimate of the Kolmogorov length scale $\ell _K$, marking the end of the inertial subrange for $K'$ and $K'_\rho$. Defined dimensionally as $(\nu ^3/ \langle \mathcal {E}\rangle )^{1/4}$, its non-dimensional expression in shear-layer units is

(3.8a)\begin{align} \ell_K &\equiv \langle \mathcal{E} \rangle^{{-}1/4}(Re^s)^{{-}3/4} \end{align}

(3.8b)\begin{align} &\approx 2 \theta^{{-}1/4} \, (Re^s)^{{-}3/4}, \end{align}

assuming for simplicity that $1-R_f\approx 1$.

We also estimate the Batchelor length scale $\ell _B$, marking the end of the viscous convective sub-range for $K'_\rho$, as

(3.9)\begin{equation} \ell_B \equiv\ell_K Pr^{{-}1/2} \approx 0.1 \, \theta^{{-}1/4} \, (Re^s)^{{-}3/4} \quad \text{for} \ Pr=700. \end{equation}

These estimates give $\ell _K \approx 0.02$ and $\ell _B \approx 0.0007$ for T2 and T3. For these data sets, we thus only have suitable resolution in $x,z$ for the velocity field (since ${\textrm {d}\kern0.06em x} = \textrm {d}z \approx 1.5\ell _K \approx 40 \ell _B$ and ${\textrm {d} y} \approx 5 \ell _K \approx 130 \ell _B$, see Part 1, Appendix B).

These estimates also suggest that while the magnitude of energy reservoirs and fluxes are strong functions of $\theta$, the Kolmogorov and Batchelor scales are stronger functions of $Re^s$ than of $\theta$. In particular, we note that in flows having identical ‘$\theta Re^s$ intensity’, still have $\ell _K$, $\ell _B \propto (Re^s)^{-1/2}$, suggesting inherently different small-scale dynamics in the same flow ‘regime’. This is consistent with the different ‘flavours’ of stratified turbulence described in Part 1, § 6.4, wherein high-$\theta$, low-$Re^{s}$ flows have greater overturns, while low-$\theta$, high-$Re^{s}$ flows have more extreme enstrophy events.

3.2. Spatio-temporal profiles

In figure 3 we plot the vertical, spanwise and temporal structure of the turbulent energy reservoirs $(K',K'_\rho )(\boldsymbol {x},t)$ and the volumetric fluxes $\mathcal {E}(\boldsymbol {x},t)$ and $(\mathcal {F},\bar {\epsilon },\mathcal {P},\mathcal {P}_\rho,\mathcal {B})(y,z)$. We show $z$ profiles in the left column (averaged in $x,y,t$ or $y$), the $y$ profiles in the middle column (averaged in $x,z,t$ or $z$) and the $t$ profiles in the right column (averaged in $x,y,z$ or $y,z$). We only show six data sets whose energetics previously revealed interesting aspects representative of ${\rm H}$ flows (H1 and H4, first and second rows), ${\rm I}$ flows (I7 and I8, third and fourth row) and ${\rm T}$ flows (T1 and T3, fifth and sixth row, noting that T2 was omitted because it is similar to T3). The mean energy reservoirs $\bar {K},\bar {K}_\rho$ are omitted for clarity (but can be visualised by mentally squaring $\bar {u},\bar {\rho }$ in Part 1, figure 3). Note that $\chi$ and $\varPhi ^{\bar {K}_\rho }$ are omitted too; the former because of its severe underestimation and low values (typically below the axis limits), and the latter as a consequence of our focus on kinetic energy budgets.

Figure 3. Profiles of turbulent energy reservoirs and fluxes in the vertical direction $z$ (a,d,g,j,m,p); the spanwise direction $y$ (b,e,h,k,n,q); and time $t$ (c, f,i,l,o,r) in six data sets: (a–c) H1; (d–f) H4; (g–i) I7; (j–l) I8; (m–o) T1; (p–r) T3. Axis limits and labels are identical in all panels of the left, middle and right column, respectively. Note the semi-log scale in all panels. Data that are inferior to the lower axis limit are omitted (e.g. $\mathcal {F}$ partially ${<}0$ near $z=0$ in the left columns, and $\mathcal {P}_\rho, \mathcal {B}$ typically ${<}10^{-4}$ in the middle and right columns except in T3). Also note that $\mathcal {F},\bar {\epsilon },\mathcal {P},\mathcal {P}_\rho,\mathcal {B}$ are by definition time independent (c, f,i,l,o,r).

First, looking at the vertical profiles, $K'$ (in solid black) becomes nearly flat and symmetric over most of the shear layer as the flow becomes increasingly turbulent (panels g,k,n,q). The forcing $\mathcal {F}$ (in green) is always highest near the top and bottom edges of the shear layer (where $|\bar {u}|$ and $|\bar {\rho }|$ are highest) and vanishes in the middle (where it reaches slightly negative values, not shown on the log scale, where the $\bar {u}=0$ and $\bar {\rho }=0$ levels are offset). The turbulent dissipation $\mathcal {E}$ (in blue) closely matches the structure of $K'$ in all flows, albeit with approximately $1/10$ magnitude (giving an approximate turbulent dissipation time scale $K'/\mathcal {E} = O(10$ A.T.U.$)$). In T3 only, the turbulent dissipation exceeds the mean dissipation $\bar {\epsilon }$ (in cyan) throughout most of the shear layer (panel q). The mean dissipation $\bar {\epsilon }$ highlights the structure of the mean shear $\partial _z \bar {u}$, typically higher on either side of the layer of mixed density, which matches more closely the structure of $K'_\rho$ than of $K'$. The scalar variance $K'_\rho$ (in dotted black) has a much sharper and sometimes asymmetric peak than $K'$, as seen in ${\rm H}$ flows (symmetric Holmboe waves in panel (a), asymmetric Holmboe waves in panel (d)) and some ${\rm I}$ flows (larger variance at the lower edge of the mixed layer in panel (g)). In ${\rm I}$ and ${\rm T}$ flows, $K'_\rho$ tends to exhibit two peaks on either side of the mixed layer, due to overturning motions entraining fluid from the unmixed layers. In these flows the buoyancy flux $\mathcal {B}$ (in magenta) and production of scalar variance $\mathcal {P}_\rho$ (in dotted red) also tend to be nearly equal (as would be the case under linear stratification), and to closely match the structure of $K'_\rho$ (albeit with smaller magnitude, see panels (g,k,n,q)). Finally, in ${\rm T}$ flows, the buoyancy flux $\mathcal {B}$ (in magenta) and the production of turbulent energy $\mathcal {P}$ (in red) have very similar profiles, corresponding to a uniform flux Richardson number $R_f(z)\approx 0.05$. This may be another hallmark of the self-organising equilibrium of stratified turbulent shear layers, related to the convergence of the gradient Richardson number to an equilibrium value ${\approx }0.10\text{--} 0.15$ as shown in Part 1.

Second, looking at the spanwise profiles, $K'$ nearly always has a sharper peak than the nearly flat $K'_\rho$ (panels e,h,l,o,r), a situation exactly opposite to that of their vertical profiles. The peak in $K'$ near $y=0$ is also much sharper than that of the mean flow $\bar {u}$ (see Part 1, figure 3), suggesting a peak in the ratio of turbulent-to-mean energy $K'/\bar {K}$ near $y=0$. This dichotomy between peaked vs flat spanwise profiles also extends to the turbulent fluxes $\mathcal {P},\mathcal {P}_\rho,\mathcal {B}$ vs the mean fluxes $\mathcal {F}$ and $\bar {\epsilon }$. Moreover, we know that outside the shear layer ($|y| > L_y$, $|z| > 1$) the turbulent fluxes decay to zero whereas the mean fluxes remain high.

Third, in our interpretation of the $z$ and $y$ profiles, we recall that assuming a steady state and negligible boundary fluxes $\varPhi ^{\bar {K}},\varPhi ^{K'}$ should yield local (point-wise) equality of the following fluxes: $\langle \mathcal {F}\rangle _y \approx \langle \mathcal {P}\rangle _y + \langle \bar {\epsilon } \rangle _y$ and $\langle \mathcal {E}\rangle _y \approx \langle \mathcal {P}\rangle _y - \langle \mathcal {B} \rangle _y$ at all $z$ (and vice versa, equality of $z$ averages at all $y$, as in (2.7)). In other words, these fluxes need to approximately balance everywhere both in $z$ and $y$ for $K'(\boldsymbol {x},t)$ to be steady (in term of curves: ‘$\textrm {green} = \textrm {red} + \textrm {cyan}$’ and ‘$\textrm {blue}=\textrm {red}-\textrm {magenta}$’). As we see in the left and middle columns, this is rarely exactly the case in our data, presumably due, again, to our under-resolution of turbulent variables (and their gradients), and because of the existence of non-negligible boundary fluxes $(\varPhi ^{\bar {K}}, \overline {\varPhi ^{K'}})(y,z) \neq 0$ due to (i) the slight non-periodicity of SID flows in $x$: $\overline {\partial _{x}(uK)}$, $\overline {\partial _{x}(u'K')}$; (ii) the inevitable advective transport of $K'$ across our artificial ‘shear layer’: $\overline {\partial _{y}(v' K')+\partial _{z}(w' K')}$; (iii) the unknown work of the mean and turbulent pressures in $x,y,z$.

Fourth, looking at the temporal profiles, the amplitude of the fluctuations in $K',K'_\rho,\mathcal {E}$ is small in H1, H4, T3 (panels c, f,s), and much larger in I7, I8, T1 (panels i,m,p). This is consistent with our nomenclature of the ${\rm I}$ regime as intermittently turbulent, and with our previous finding that T1 is actually closer to ${\rm I}$ flows than to T2 and T3, whose fluctuations are steadily large. Moreover, $K'$ and $K'_\rho$ are not generally correlated in ${\rm I}$ and ${\rm T}$ flows, i.e. the intensity of velocity and density fluctuations do not temporally vary hand in hand (as might be incorrectly generalised from the time- and volume-averaged statement (3.6)). Finally, we recall that the temporal profiles of $K'$ and $K'_{\rho }$ should only reflect (with opposite correlation) the profiles of $\mathcal {E}$ and $\chi$, respectively, since all other fluxes plotted are (by definition) time independent, and boundary fluxes are neglected. The negative correlation between $K'$ and $\mathcal {E}$ is, however, not always observed in our data (in fact, they appear almost positively correlated in most panels). These last two findings are not surprising, especially in light of our findings in Part 1 (figure 7) that turbulent fractions based on enstrophy or overturning can be largely uncorrelated, due to spatial heterogeneity of turbulent patches and the non-periodicity of our measurement volume along $x$.

3.3. Spectra

We now delve deeper into the flow energetics by investigating spectra. We start with spectra of the turbulent kinetic energy and scalar variance along $x$, before focusing on individual velocity components and all variables $x,y,z,t$. The limitations in our measurements of turbulent energetics, frequently hinted at in the above sections, will be summarised in light of spectral results in § 3.4.

3.3.1. Spectra of $K',K'_\rho$ in $x$

We define the spectral densities in $x$ of the mean turbulent kinetic energy $E^x_{K'}$ and scalar variance $E^x_{K'_\rho }$ such that

(3.10a,b)\begin{equation} \int_0^{k_{x, max}} E^x_{K'} \, \,\mbox{d} k_x = \langle K' \rangle , \quad \int_0^{k_{x, max}} E^x_{K'_\rho} \, \,\mbox{d} k_x =\langle K'_\rho \rangle. \end{equation}

Their unambiguous definitions and the details of their practical computation from our discrete gridded data are given in Appendix A. In the above, $k_{x, max} \equiv {\rm \pi}/{\textrm {d}\kern0.06em x}$ is the maximum (Nyquist) wavenumber that can be resolved in $x$. The (unusual) need to integrate from $k_x=0$ rather than from the minimum wavenumber $k_{x\,min} \equiv {\rm \pi}/L_x$ comes from the fact that energy is contained in the mean ($k_x=0$), an inevitable consequence of the above definitions and of our definition of fluctuations around $x-t$ averages (more details in Appendix A.3).

In figure 4, we plot the densities $E^x_{K'}$ (black solid) and $E^x_{K'_\rho }$ (grey dashed) for all data sets. To correct for errors inherent to computing Fourier transforms of noisy and non-periodic data (over-estimating high frequencies), here we plot estimations of these densities (i.e. periodograms) using Welch's averaging method. This standard method divides each original signal along $x$ into a series of overlapping segments, applies a window function to render them periodic, and returns the average square magnitude of their discrete Fourier transform (more details in Appendix A.4).

Figure 4. Spectral density along $x$ of the turbulent kinetic energy $E^x_{K'}$ and scalar variance $E_{K'_\rho }$ for all data sets, calculated with Welch's averaging method. The range of non-zero wavenumbers shown $[k_{x, min}, k_{x, max}] =[{\rm \pi} /L_x, {\rm \pi}/{\textrm {d}\kern0.06em x}]$ varies slightly among data sets because of different domain lengths $2L_x$ and resolutions ${\textrm {d}\kern0.06em x}$ (see Part 1, Appendix B). Note the $k_x=0$ energy content (see text and Appendix A for details). Axis limits and labels are identical in all panels.

Since we do not expect any turbulent signal in our laminar data set, the L1 spectra (panel a) are plotted as a ‘control’, i.e. a baseline measure of inevitable artefacts due to the nature of our data and analysis. In panel (a), $E^x_{K'}$ exhibits a distinct hump at intermediate wavenumbers ($k_x \approx 2\text{--} 30$), correlated with a distinct hump or flattening of $E^x_{K'_\rho }$. This artefact is also found in varying degrees in most other data sets around $k_x \approx k_{x, max}/2$, affecting our most turbulent data to a lesser degree (panels l–p).

Putting this ‘hump’ artefact aside, most ${\rm H}$, ${\rm I}$ and ${\rm T}$ spectra exhibit relatively similar shapes. The kinetic energy spectrum $E^x_{K'}$ is flat in the energy-containing range $k_x \lesssim 1$ (length scales $\gtrsim 6$), decaying approximately as $k_x^{-\beta }$ in the inertial sub-range, with typical values around $\beta \approx 2.0\text{--} 3.5$ (on either side of the ‘hump’), and a slightly different power law decay near $k_{x,max}$. This decay exponent $\beta$ is thus significantly larger than the classical Kolmogorov value of $\beta =5/3\approx 1.7$ expected in the forward cascade of isotropic turbulence (Pope Reference Pope2000, § 6.5), as well as in the horizontal spectrum of strongly stratified turbulence (Lindborg Reference Lindborg2006). At low wavenumbers, these steep horizontal spectra, indicative of an incomplete cascade, may be due to the relatively weak turbulent intensity of our data sets, in which ‘the energy is mainly dissipated by vertical shearing at scales close to the forcing scales’, as explained in Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007, § 5). At high wavenumbers, this steepness may be due to the low-pass filtering effects of particle image velocimetry (PIV), which obtains velocity vectors from ‘discrete interrogation windows’ that artificially reduce high-wavenumber content, as explained in Appendix B.2.

The scalar variance spectra $E^x_{K'_\rho }$ exhibit a shape similar to $E^x_{K'}$, albeit with slightly smaller amplitude (as expected from the asymptotic 1/10 partition in figure 1j). These spectra also have a smoother inertial sub-range decay extending all the way to $k_{x, max}$, at least in the most turbulent data (panels k–p), where the scaling $k_x^{-\beta }$ is in better agreement with the expected value $\beta = 5/3$ (Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2016, § 12.11).

3.3.2. Spectra of all components in $x,y,z,t$

In figure 5 we plot these spectral densities in $x$ (top row), $y$ (second row), $z$ (third row) and $t$ (bottom row) for three representative data sets H1 (left column), I2 (middle column) and T3 (right column). To investigate the effects of the non-periodicity of our data on the energy spectra, we plot in $x$ the energy densities obtained using the standard discrete Fourier transform (DFT) periodogram (thin lines) and using Welch's estimated periodogram (thick lines, as in figure 4). Note that we only show the Welch in $t$ for conciseness, and we only show the DFT in $y$ and $z$ because the smaller numbers of data points in these directions render Welch's segmentation inappropriate (more details in Appendix A.4).

Figure 5. Spectral density of energy in individual velocity components $u'$ (blue), $v'$ (green), $w'$ (red) and of $K'_\rho$ (grey) in all directions $x$ (a–c), $y$ (d–f), $z$ (g–i) and $t$ (j–l) for three representative data sets H1 (a,d,g,j), I2 (b,e,h,k) and T3 (c, f,i,l). The mean energies $(1/2) \langle u'^2 \rangle$, $(1/2) \langle v'^2 \rangle$, $(1/2) \langle w'^2 \rangle$ and $\langle K'_\rho \rangle$ are given by one-dimensional integration of any respective density, e.g. $(1/2) \langle u'^2 \rangle = \int _0^{k_x, max} E^x_{u'} \,\mbox {d} k_x = \int _0^{k_y, max} E^y_{u'} \,\mbox {d} k_y = \int _0^{k_z, max} E^z_{u'} \,\mbox {d} k_z = \int _0^{\omega, max} E^t_{u'} \,\mbox {d} \omega$, etc. Mean values at $k_x,k_y,k_z,\omega =0$ are not shown for clarity. The spectral range depends on domain length and resolution $(k_x,k_y,k_z,\omega ) \in [{\rm \pi} /L_x, {\rm \pi}/{\textrm {d}\kern0.06em x}]\times [{\rm \pi} /L_y, {\rm \pi}/{\textrm {d} y}] \times [{\rm \pi}, {\rm \pi}/\textrm {d}z] \times [2{\rm \pi} /L_t, {\rm \pi}/\textrm {d}t]$. Note the different axis scales between (a–i) and (j–l). In $x$, we compare spectra obtained by the standard DFT (thin lines) and by Welch's method (thick lines). In $y,z$ we only show the former, and in $t$ we only show the latter (see text for more details).

We see in panels (a–c) that the standard DFT (thin lines) consistently overestimates the high-wavenumber content ($k_x\gtrsim 30$) compared with the Welch (thick lines), as expected from the fact that the latter is designed to minimise the effects of edge discontinuities in our data, incorrectly rendered as high-wavenumber energy by the standard DFT (called spectral leakage, or the Gibbs phenomenon, see Appendix B.1). Given these observations, we should remain critical in our interpretation of DFT spectra in $y,z$ (panels d–i), where spectral leakage is expected to be more significant than in $x$ due to the shorter size of the signals. We hypothesise that these spectra exhibit an inertial sub-range decay closer to $k^{-5/3}$ in $y$ and $z$ than in $x$ as a result of greater spectral leakage (overestimating high-wavenumber energy) countering the effects of PIV filtering (underestimating high-wavenumber energy).

We further see that $u'$ usually has the most energy across all wavenumbers and frequencies, and $E_{u'}>E_{v'}>E_{w'}$. The lowest energy in $w'$ is consistent with the expectation that vertical motions are partially hindered by the stable stratification at $Ri_b^s>0$. The higher energy in $u'$ than in $v'$, particularly clear at very low streamwise wavenumbers $k_x\lesssim 1$, is partly due to our definition of fluctuations around $x-t$ averages and to the fact that the flow is not perfectly parallel (i.e. $u'$ can have a slight residual large-scale variance along $x$, as explained in Appendix A.3).

The above observation that $E_{u'}>E_{v'}>E_{w'}$ has a few notable exceptions. First, we diagnose that the hump artefact in $E^x_{K'}$ observed in most panels of figure 5 appears primarily caused by $v'$ since $E^x_{v'}>E^x_{u'}$ at medium and high $k_x$ (green lines in panels a–c), independently of the method (DFT or Welch). This artificial medium-scale structure in $v'(x)$ may come from the delicate stereo PIV calculation of $v$ (the component perpendicular to the laser sheet). Second, H1 exhibits $E_{v'}\approx E_{w'}$ across most wavenumbers (panels a,d,g) and even $E_{v'}< E_{w'}$ in the frequency range $\omega \approx 0.3\text{--} 1.5$ (panel j), which is consistent with the presence of Holmboe waves, known to generate vigorous vertical motions even in the presence of strong stratification (here $Ri_b^s = 0.567$).

The signature of Holmboe waves is indeed clear in the H1 temporal spectra at $\omega \approx 0.5$ (panel j), and also detectable in the longitudinal spectra $E^x_{w'}, E^x_{K'_\rho }$ around $k_x \approx 1$ (panel a, thin red and grey lines). These peaks suggest a typical phase speed $c\approx \omega /k_x \approx 0.5$ in agreement with observations in the spatio-temporal domain (not shown here). We also note in intermittent flow I2 a similar, albeit fainter, peak in all longitudinal spectra (panel b, thin lines), suggesting the faint presence of similar waves, in agreement with observations near the laminar/turbulent transitions (not shown here). Such spectral peaks are absent in the turbulent flow T3 (right column), suggesting dynamics across a broader range of spatio-temporal scales.

4.1. Reynolds stresses and Lumley triangle

We recall that the specific turbulent kinetic energy $K' \equiv (1/2) \, \textrm {tr} \, \overline {\boldsymbol {u}' \otimes \boldsymbol {u}'}$ is the isotropic part of the Reynolds stress tensor (half the trace of the one-point, one-time velocity cross-correlation tensor). By the Cauchy–Schwartz inequality, this diagonal part (isotropic ‘pressure’) sets a bound on the magnitude of the off-diagonal part (deviatoric stresses): $K' \geqslant |\overline {u'v'}|, |\overline {u'w'}|,$ or $|\overline {v'w'}|$ (Pope Reference Pope2000, eq. (5.109)). In idealised isotropic turbulence, all deviatoric stresses are zero, thus there is no transfer between mean and turbulent kinetic energy, hence $\mathcal {P}=0$. By contrast, in shear-driven turbulence, this bound becomes more meaningful due to the crucial production of $K'$ at rate $\mathcal {P}>0$ resulting from $\overline {u'w'} \neq 0$, i.e. from the net correlation of anisotropic eddies at large (energy-containing) scales.

To quantify this anisotropy of the Reynolds stresses, we consider the widely used normalised velocity anisotropy tensor $\boldsymbol{\mathsf{b}}$, defined as the deviatoric part of the normalised Reynolds stress tensor $\overline {\boldsymbol {u}' \otimes \boldsymbol {u}'}/(2\overline {K'})$ with components

(4.1)\begin{equation} b_{ij}(y,z) \equiv \frac{\overline{u'_i u'_j}}{\overline{u'_l u'_l}}- \frac{\delta_{ij}}{3}, \end{equation}

as in Pope (Reference Pope2000) (§ 11.3.2). Since by definition $\textrm {tr} \, \boldsymbol{\mathsf{b}} = b_{ii} = 0$ (first invariant), this tensor has only two independent invariants: $\textit {II}_b \equiv \textrm {tr} \, \boldsymbol{\mathsf{b}}^2/2$ (second invariant) and $\textit {III}_b \equiv \det \boldsymbol{\mathsf{b}} = \textrm {tr} \, \boldsymbol{\mathsf{b}}^3/3$ (third invariant), which are more conveniently defined in normalised form as

(4.2a,b)\begin{equation} \eta(y,z) \equiv \sqrt{\frac{\textit{II}_b}{3}}, \quad \xi(y,z)\equiv \sqrt[3]{\frac{-\textit{III}_b}{2}}. \end{equation}

The local state of large-scale anisotropy of a turbulent flow at any point $y,z$ can therefore be described by a point in the $\xi \text{--} \eta$ plane, lying inside the so-called Lumley triangle (Lumley Reference Lumley1978), drawn with thick lines in figure 6(a). The point $\xi =\eta =0$ corresponds to isotropic turbulence; the left (respectively right) straight edge $\eta =\mp \xi /2$ correspond to oblate (respectively prolate) axisymmetric turbulence, i.e. one principal eigenvalue being smaller (respectively larger) than the other two; and the top curved edge $\eta =\sqrt {1/27+2\xi ^3}$ corresponds to two-component turbulence (one principal eigenvalue being zero). In summary, the vertical axis $\eta$ quantifies the degree of anisotropy, while the horizontal axis $\xi$ quantifies its shape (oblate $\xi <0$ vs prolate $\xi >0$).

Figure 6. Degree and shape of Reynolds stress anisotropy in all 16 data sets. (a,b) Mean values $\langle \xi \rangle _{y,z},\langle \eta \rangle _{y,z}$ (zoomed in detail in b). The Lumley triangle is highlighted by thick lines, and the limiting cases of turbulence are shown schematically with principal axis coordinates. (c–r) All $n_y n_z$ data points of $(\xi,\eta )(y,z)$, coloured with the absolute vertical coordinate $|z|$ within the shear layer.

In figure 6(a,b) we plot the mean $\langle \xi \rangle _{y,z}$ and $\langle \eta \rangle _{y,z}$ in all 16 data sets. First, we observe in panel (a) that all points are clustered in a narrow top-right region of strong prolate anisotropy, shown in greater detail in panel (b). This is consistent with our prior spectral observation that the streamwise velocity perturbations dominate over the other two: $|u'|^2>|v'|^2,|w'|^2$. Second, our ‘control’ data set L1, being non-turbulent and therefore primarily affected by unphysical artefacts, distinguishes itself by being the only data set lying outside of Lumley's realisability triangle (though all $y,z$ points are by construction inside, the $y\text{--} z$ average does not have to be since this ‘curved triangle’ is not convex). Third, asymmetric ${\rm H}$ flows (H2, H4) lie closer to the two-component (top) limit, while symmetric ${\rm H}$ flows (H1, H3) and most ${\rm I}$/${\rm T}$ flows lie closer to the prolate axisymmetric (right) limit. Fourth, almost all ${\rm I}$ flows exhibit stronger anisotropy than ${\rm H}$ and ${\rm T}$ flows, and lie closer to the one-component limit. This is a result of greater temporal variability (intermittency) in the streamwise component $u' \equiv u - \langle u \rangle _{x,t}$, defined with respect to the streamwise and temporal average. We verified that removing intermittency effects by defining perturbations with respect to the streamwise average alone ($\boldsymbol {u} - \langle \boldsymbol {u} \rangle _{x}$) did move ${\rm I}$ flows slightly away from the one-component limit, but it did not change the qualitative picture of panels (a,b). We also verified that removing streamwise variance effects ($\partial _x \bar {u} \neq 0$) by defining perturbations with respect to the temporal average alone ($\boldsymbol {u} - \langle \boldsymbol {u} \rangle _{t}$) changed the picture very little. In other words, $u'$ always dominates and anisotropy is not significantly biased by our definition of $u' \equiv u - \langle u \rangle _{x,t}$ (used throughout Part 1 and Part 2).

In figure 6(c–r) we plot the underlying $n_y n_z$ data points within the triangle for each data set (of which panels a and b show the centre of mass), to highlight the anisotropy over the full range of y and more particularly of z (in colour). In all flows, we observe large spatial variations around the mean, closely following the prolate axisymmetry limit (right edge), i.e. a state in which $v',w'$ have nearly (but not exactly) equal magnitude, while being dominated by $u'$, no matter the $y,z$ location. This general trend is nuanced by the following subtleties. First, ${\rm H}$ flows exhibit the greatest variations, and are unique in that they include pockets of oblate anisotropy ($\xi <0$) for $|z|\approx 0.3\text{--} 1$. Second, some ${\rm I}$/${\rm T}$ flows (I3, I6, T1, T2) have data points at $|z|\lesssim 0.3$ which deviate significantly away from the right edge (axisymmetry) and lie closer to the centre of the triangle. Third, the data points closest to the mid-point of the shear layer ($|z|\lesssim 0.2$, in black) tend to be the most anisotropic (largest $\eta$) in ${\rm H}$/${\rm I}$ flows but the least anisotropic in ${\rm T}$ flows (smallest $\eta$).

These findings are qualitatively consistent with the unforced direct numerical simulations (DNSs) of Smyth & Moum (Reference Smyth and Moum2000a), who observed oblate axisymmetry during the initial growth of the Kelvin–Helmholtz instability and the turbulent transition, followed by prolate axisymmetry during the turbulent and decay phases (see their figure 6).

4.2. Spatial profiles

To delve deeper into these tantalising observations, we plot in figure 7 the spatial structure of $\eta (y,z)$, $\xi (y,z)$, and the vertical structure of the six individual components of the (symmetric) tensor $\langle (b_{ij}) \rangle _y(z)$ for the four representative data sets H4, I4, I6 and T2. In the contrasting cases of H4 and T2, we also plot the underlying probability density function (p.d.f.) of the $(u',w')$ data at three distinct vertical locations.

Figure 7. Spatial structure of the anisotropy tensor $\boldsymbol{\mathsf{b}}$ of data sets H4, I4, I6, T2 (top to bottom row). Left to right columns: second invariant $\eta (y,z)$, third invariant $\xi (y,z)$ (including the averages in each direction, superposed in white); vertical structure of the diagonal components $b_{11},b_{22},b_{33}$ averaged in $y$; vertical structure of the off-diagonal components $b_{12},b_{13},b_{23}$ averaged in $y$. In (e,r) we also show the p.d.f. of the $(u',w')$ clouds (rescaled histogram, here four equidistant contours at $20,40,60,80\,\%$) for H4 and T2 at three different vertical locations flagged by asterisks in (d,q) (these are $z\in [-1, -0.9], [-0.55, -0.45], [0.05, 0.15]$ in H4, and $z\in [-0.9, -0.8], [-0.55, -0.45], [-0.05, 0.05]$ in T2, noting that in both we restrict the region to $|y|\leqslant 0.5$ to show a stronger signal).

First, starting with the $y\text{--} z$ structures (colour plots in the left two columns), we find that the region of weak oblate anisotropy in H4 (light grey in panel a and light blue in panel b) lies at the periphery of a core of strong prolate anisotropy (this subtle structure is lost in the $y$ and $z$ averages superimposed in white). We explain this oblate pocket by the particular structure of confined Holmboe waves described in Lefauve et al. (Reference Lefauve, Partridge, Zhou, Caulfield, Dalziel and Linden2018) in this same H4 flow, and in particular by the large values of $v'$ and its odd symmetry about the $y=0$ axis, responsible for the divergence and convergence of streamlines in $z=\textrm {const.}$ planes around the upward-pointing crests of the density interface (see their figure 8(k,l) and point (v) in their § 6.1.2). By contrast, ${\rm I}$/${\rm T}$ flows have a more uniform structure, ranging from strong prolate anisotropy in I4 (panels f,g) to weaker and less prolate anisotropy in I6 (panels j,k) and T2 (panels n,o) especially in the most turbulent region $|z|\lesssim 0.5$.

Second, moving on to the diagonal components (third column), we confirm our above claims that $b_{11}=\overline {u'^2}/(2\overline {K'})-1/3$ (and therefore $u'$) dominates in all flows, where it approaches its upper bound of $2/3$ in the most anisotropic $z$ locations (darkest colours in the first two columns), while the complementary $b_{22},b_{33}$ approach their lower bound of $-1/3$. Furthermore, $b_{22}>b_{33}$ in all flows (and therefore $v'$ dominates over $w'$), a natural consequence of stratification inhibiting vertical motions.

Third, moving on to the off-diagonal components (fourth column), we recall that they are all bounded above and below by $\pm 1/2$ by the Cauchy–Schwarz inequality. We find that $b_{13}\equiv \overline {u'w'}/(2\overline {K'})$ dominates over $b_{12},b_{23}$ almost everywhere in all flows, being always positive and even reaching $\langle b_{13}\rangle _y(z=0)= 0.125$ in T2 (panel q), i.e. no less than 25 % of its upper bound. This proves that while $w'$ contributes less than $v'$ to the reservoir $K'$, it contributes much more than $v'$ to the production of $K'$ because of its much greater correlation with $u'$, recalling that the logarithmic production rate $\partial _t (\ln \sqrt {K'}) = \partial _t K'/(2K')$ is $\mathcal {P}/(2K')\equiv -b_{12}\partial _y \bar {u}-b_{13}\partial _z \bar {u}$. The broad peak of $b_{13}$ in the vigorous flows I6 and T2 (panels m,q) is absent in the less vigorous flow I4 and in the Holmboe flow H4 (panels d,i), the latter having instead two narrower peaks at $z\approx -1$ and $z\approx 0$.

Fourth, these contrasting $\langle b_{13}\rangle _y(z)$ profiles in H4 and T2 can be understood by their respective p.d.f.s in the $u'-w'$ plane (fifth column). In the Holmboe flow (panel e) the $z\approx -1$ peak (denoted by $\ast$) is due to a compact but fairly tilted distribution towards the first and third quadrant ($u'w'>0$), while the $z\approx -0.5$ trough ($\ast \ast$) is due to a broader but more up–down symmetric distribution, and the $z\approx 0$ peak ($\ast \!\ast \!\ast$) is due to a yet broader, but thinner and more tilted distribution. In the turbulent flow (panel r) the increase of $\langle b_{13}\rangle _y(z)$ with decreasing $|z|$ is due to a broader distribution (compare $\ast$ and $\ast \ast$) followed by an increased tilt (compare $\ast \ast$ and $\ast \!\ast \!\ast$).

Fifth, it is possible to improve the quantification of this tilt of $(u',w')$ distributions by investigating the orientation of the principal axes of $\boldsymbol{\mathsf{b}}$, given by its eigenvectors assembled in a matrix $\boldsymbol {V}$ such that $\boldsymbol {b}\equiv (b_{ij})=\boldsymbol {V}\boldsymbol {\varLambda } \boldsymbol {V}^{-1}$ (where $\boldsymbol {\varLambda }$ is the diagonal matrix of principal eigenvalues). To avoid the intricate analysis of three Euler angles describing the three-dimensional rotation matrix $\boldsymbol {V}$, we take advantage of the fact that $|b_{12}|,|b_{23}|\ll |b_{13}|$ in T2 to simplify the analysis to a single angle $\beta$ describing the two-dimensional rotation $\tilde {\boldsymbol {V}}$ (around the $v'$ axis) of the reduced $\tilde {\boldsymbol {b}}\equiv [b_{11}\ b_{13}; b_{13}, b_{33}]$. This angle $\beta (y,z)$ locally quantifies the tilt between the $u'$ axis and the major principal axis, and we therefore expect $0<\beta \ll 90^\circ$ based on panel (r). Although we do not plot it for conciseness, the profile of $\langle \beta \rangle _y(z)$ follows almost exactly that of $\langle b_{13}\rangle _y(z)$, with a minimum value of $5^\circ$ at $|z|=1$ and a maximum value of $16^\circ$ at $|z|=0$, in excellent agreement with the qualitative insights derived from panel (r).

4.3. Velocity gradients and dissipation surrogates

We now investigate the anisotropy of velocity gradients, controlling the rate of turbulent dissipation, which is by definition the sum of 12 squared gradient terms belonging to three key groups

(4.3)\begin{align} \langle \mathcal{E} \rangle \equiv \dfrac{2}{Re^s} \langle s_{ij}'s_{ij}' \rangle &= \dfrac{2}{Re^s} \left\langle \underbrace{(\partial_x u')^2 + (\partial_y v')^2 + (\partial_z w')^2}_\text{longitudinal}\right. \nonumber\\ &\quad + \underbrace{(\partial_y u')^2 + (\partial_z u')^2 + (\partial_x v')^2 + (\partial_z v')^2 + (\partial_x w')^2 + (\partial_y w')^2}_\text{transverse}\nonumber\\ &\quad \left. + \underbrace{\partial_y u'\partial_x v' + \partial_z u'\partial_x w' + \partial_z v'\partial_y w'}_\text{asymmetric} \right\rangle. \end{align}

In idealised homogeneous isotropic turbulence, all terms belonging to the same group (longitudinal, transverse or asymmetric) are equal. Using the continuity equation $\partial _{x_i} u'_i=0$, it can further be shown that any transverse term is twice as large as any longitudinal term (e.g. $\langle (\partial _y u')^2 \rangle = 2\langle (\partial _x u')^2 \rangle$, etc) while any asymmetric term is negative and only half as large (e.g. $\langle \partial _y u' \partial _x v' \rangle = (-1/2)\langle (\partial _x u)^2 \rangle$, etc.) (Almalkie Reference Almalkie2012, pp. 22–23). Plugging in these relations into (4.3) allows us to estimate $\langle \mathcal {E} \rangle$ under the assumption of isotropy using only one term (instead of 12), as follows:

(4.4)\begin{equation} \langle \mathcal{E} \rangle \approx \begin{cases} \dfrac{15}{Re^s}\langle (\partial_{x_i} u_i')^2\rangle & \text{(longitudinal surrogate)} \\ \dfrac{15}{2Re^s}\langle (\partial_{x_i} u_j')^2\rangle, \quad i\neq j & \text{(transverse surrogate)} \\ -\dfrac{30}{Re^s}\langle \partial_{x_i} u_j'\partial_{x_j} u_i' \rangle, \ i\neq j & \text{(asymmetric surrogate)}, \end{cases} \end{equation}

where, importantly, we do not sum over repeated indices here. These simple one-dimensional and one-component surrogates have been used for decades in laboratory and field measurements due to the difficulty of measuring more than one or two terms (although there exist more sophisticated multi-component models that relax isotropy, and e.g. assume axisymmetry instead). Our data sets provide us with the complete set of twelve terms in all directions $(x,y,z,t)$ and thus allow us to test the validity of these (time- and volume-averaged) surrogates and, thus, of the underlying assumption of small-scale isotropy.

To do so, we arrange the above 12 candidates into the following ‘surrogate dissipation matrix’ $\langle \varepsilon _{ij}\rangle$, and define its relative estimation error $\langle \tilde {\varepsilon }_{ij}\rangle$ as

(4.5a,b)\begin{equation} \langle \varepsilon_{ij} \rangle \equiv \frac{15}{Re^s} \begin{bmatrix} \displaystyle\langle (\partial_x u')^2 \rangle & \displaystyle\frac{1}{2}\langle (\partial_y u')^2\rangle & \dfrac{1}{2}\langle (\partial_z u')^2\rangle \\ \displaystyle\frac{1}{2}\langle (\partial_x v')^2\rangle & \langle (\partial_y v')^2 \rangle & \dfrac{1}{2}\langle (\partial_z v')^2\rangle \\ \displaystyle\frac{1}{2}\langle (\partial_x w')^2\rangle & \dfrac{1}{2}\langle (\partial_y w')^2\rangle & \langle (\partial_z w')^2 \rangle \\ -2\langle \partial_y u'\partial_x v' \rangle & -2\langle \partial_z u'\partial_x w' \rangle & -2\langle \partial_z v'\partial_y w' \rangle \end{bmatrix}, \quad \langle \tilde{\varepsilon}_{ij} \rangle \equiv \frac{\langle \varepsilon_{ij} \rangle - \langle \mathcal{E} \rangle }{\langle \mathcal{E} \rangle } , \end{equation}

where the top $3\times 3$ block contains the three longitudinal terms (diagonal) and the six transverse terms (off diagonal), as in Portwood, de Bruyn Kops & Caulfield (Reference Portwood, de Bruyn Kops and Caulfield2019). The fourth row contains the three asymmetric terms; these are rarely used in applications because they are impractical, but we include them nonetheless to obtain a complete picture of small-scale anisotropy.

In figure 8 we visualise the relative error matrix $\langle \tilde {\varepsilon }_{ij} \rangle$ in all ${\rm H}$, ${\rm I}$ and ${\rm T}$ flows, together with their ‘intra-regime mean’ on the left. Each entry in each $4\times 3$ matrix is coloured according to its departure away from isotropy; a negative value means that the surrogate is an underestimation (in blue, bounded below by $-1$), while a positive value means that the surrogate is an overestimation (in red, not bounded above but always ${<}3$ here). We keep in mind the potential inaccuracies of these data based on velocity gradients, given the limitations highlighted in § 3.4 and Appendix B.

Figure 8. Anisotropy of the 12 density gradients in (4.3), measured by the error made by using them as surrogates for $\langle \mathcal {E} \rangle$ based on the assumption of isotropy (as in (4.4)) in all H, I and T data sets. Colours show the value of each entry of the $4 \times 3$ matrix of relative error $\langle \tilde {\varepsilon }_{ij} \rangle$ defined in (4.5a,b). All H, I, T data sets are shown, together with their mean across each regime in the left-most panels (a, f,o). Blue indicates an underestimation, red indicates an overestimation and darker shades indicate poorer estimation, i.e. stronger anisotropy. (s) Matrix norm quantifying dissipation anisotropy vs $\theta Re^s$ (isotropy corresponds to $||\langle \tilde {\varepsilon }_{ij} \rangle || = 0\,\%$.

Focusing first on the individual data sets (right part of the figure) and considering the global pattern of each matrix, we find strong similarities between all flows, and more specifically, between all flows within a same regime (${\rm H}, {\rm I}$ or ${\rm T}$), with perhaps only one exception in H3 (panel d). Importantly, terms that are clearly positive (respectively negative) are robustly so across most flows. This implies that intra-regime means of each matrix entry do not artificially cancel out values of opposite signs (which would incorrectly imply isotropy) and, therefore, that these means give a meaningful representative picture of each regime.

Focusing then on these robust means (panels a, f,o), we also find similarities between them. First, $(\partial _z u')^2$ (top right term) is consistently overwhelming, and overestimates $\langle \mathcal {E}\rangle$ by as much as $200\,\%$ (${\rm H}$ flows), $230\,\%$ (${\rm I}$ flows) or $140\,\%$ (${\rm T}$ flows). This can be attributed to the influence of the mean shear $(\partial _z u')^2$. In ${\rm I}$/${\rm T}$ flows, $(\partial _y u')^2$ and $(\partial _z v')^2$ also tend to consistently overestimate $\langle \mathcal {E}\rangle$, whereas they are reliable in ${\rm H}$ flows. Second $(\partial _x w')^2$ consistently underestimates by 80 %–90 %, while all four other terms involving $w'$ gradients (bottom right $2\times 2$ block) consistently underestimate by 20 %–70 %. This can be attributed to the stable mean stratification, hindering vertical motion. Third, all $x$ gradients (first column of each matrix) are generally weak. This can be attributed to the elongation of flow structures along $x$ by the mean shear. Fourth, the best estimates (lightest shade) vary slightly from regime to regime, but three terms stand out as consistently reliable: $(\partial _y v')^2$, $(\partial _z v')^2$, and $\partial _y u' \partial _x v'$ (having ${<}35\,\%$ relative error everywhere, sometimes much less).

Finally, we plot the Euclidian (Frobenius) norm of each matrix $||\langle \tilde {\varepsilon }_{ij} \rangle || \equiv (\sum _{i,j} \langle \tilde {\varepsilon }_{ij} \rangle ^2)^{1/2}$ (panel s) against the product of parameters $\theta Re^s$, identified in (3.7) as controlling the dissipation (norm of the strain rate tensor). The general trend is a decrease of small-scale anisotropy with increasing $\theta Re^s$ (stronger turbulence), from typical values of 200 %–300 % in ${\rm H}$-${\rm I}$ flows (except H3) to values below $200\,\%$ in ${\rm T}$ flows. This trend suggests that even stronger turbulence ($\theta Re^s \gg 100$) would continue to approach greater isotropy. This was indeed observed by Itsweire et al. (Reference Itsweire, Koseff, Briggs and Ferzinger1993), Smyth & Moum (Reference Smyth and Moum2000a) (see their figure 14), Hebert & de Bruyn Kops (Reference Hebert and de Bruyn Kops2006), Lang & Waite (Reference Lang and Waite2019) (see their figure 4b) and Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019) (see their figure 2) with increasing ‘dynamic range’, quantified by the buoyancy Reynolds number $Re_b$. We define $Re_b$ in the next section, explain its relation to $\theta Re^s$, and introduce other ratios of kinematic and dynamic scales to tackle parameterisations.

5.1. Background: measures of mixing

5.1.1. Direct measures: eddy diffusivities

Stratified turbulent mixing is usually modelled in large-scale circulation models by a single parameter, the eddy (or turbulent) diffusivity for the stratifying agent (heat or salt) $\kappa _T$, and for the momentum $\nu _T$. This turbulence closure scheme relies on the simple turbulent flux–mean gradient relations (see Pope Reference Pope2000, Chap. 10)

(5.1a,b)\begin{equation} \frac{\kappa_T}{Re^s} \equiv \frac{-\overline{w'\rho'}}{\partial_z\bar{\rho}} \equiv \frac{\mathcal{B}}{\bar{\bar{N}}^2}, \quad \frac{\nu_T}{Re^s} \equiv \frac{-\overline{u'w'}}{\partial_z\bar{u}} \approx \frac{\mathcal{P}}{\bar{\bar{S}}^2}. \end{equation}

The approximation in $\nu _T$ reflects the fact that production is dominated by the vertical shear $|\overline {u'w'}\partial _z\bar {u}| \gg |\overline {u'v'}\partial _y\bar {u}|$ in our flows. Importantly, the $Re^s$ factor comes from the fact that we choose to define both eddy diffusivities as non-dimensional ratios relative to the molecular value for momentum $\nu$, rather than relative to the (default and implicit) inertial scale $(\Delta U \delta u H h)/16$ (recall Part 1, § 3.2–3.3). Other authors legitimately choose to define $\kappa _T$ relative to the molecular value for the scalar $\kappa$, which then gives $\kappa _T/(Re^s Pr) \equiv \mathcal {B}/\bar {\bar {N}}^2$. Also recall (Part 1, § 5) the definitions of the square buoyancy frequency $\bar {\bar {N}}^2 \equiv -Ri_b^s \, \partial _z\bar {\rho }$ and square shear frequency $\bar {\bar {S}}^2 \equiv (\partial _z \bar {u})^2$ based on the mean flow (the double overline avoids confusion with $\overline {\partial _z\rho }$ and $(\overline {\partial _zu})^2$ which are different quantities not discussed here). The gradient Richardson number based on the mean flow was defined as $\overline {\overline {Ri_g}} \equiv \bar {\bar {N}}^2/\bar {\bar {S}}^2$.

Despite being used as the ‘direct’ (or ‘ultimate’) measures of mixing in most practical models, eddy diffusivities are necessarily simplistic descriptions of the process of stratified turbulent mixing. They have been criticised for their apparent inability to address the complex underlying energetics, in particular to disentangle the partition between irreversible mixing and reversible stirring in $\mathcal {B}$ (Salehipour & Peltier Reference Salehipour and Peltier2015). However, upon inspection of the budget equation (2.3), we find that under linear stratification ($\mathcal {B}\approx \mathcal {P}_\rho$) and neglecting boundary fluxes ($\varPhi ^{K'_\rho }\approx 0$), the buoyancy flux appears to be in ‘lock step’ with the irreversible dissipation of scalar variance $\mathcal {B}\approx \chi$, which, again under linear stratification, is equivalent to the dissipation of perturbation available potential energy, i.e. irreversible mixing (Caulfield Reference Caulfield2020). This led some authors to argue that defining $\kappa _T$ using $\chi /\bar {\bar {N}}^2$ was generally more appropriate than using $\mathcal {B}/\bar {\bar {N}}^2$, an approach known as the ‘Osborn–Cox method’ after Osborn & Cox (Reference Osborn and Cox1972) (see Salehipour & Peltier (Reference Salehipour and Peltier2015); Gregg et al. (Reference Gregg, D'Asaro, Riley and Kunze2018) and Taylor et al. (Reference Taylor, de Bruyn Kops, Caulfield and Linden2019) for more details). Following this line of thought, some authors define the flux coefficient $\varGamma$ as $\chi / \mathcal {E}$, which agrees with our approximation (3.6). However, unlike DNS data, our experimental data do not allow us to access $\chi$ directly with good accuracy, which is why we pursue an indirect approach, discussed next.

5.1.2. Indirect measures: flux coefficients, mixing lengths

Starting with the original definition (5.1a,b), we attempt to relate the elusive $\mathcal {B}$ to the more tangible $\mathcal {E}$ (the ‘turbulence intensity’). This approach proposes equivalent definitions for $\kappa _T,\nu _T$ using our previous definitions of $\varGamma, R_f, \overline {\overline {Ri_g}}$ and of a new turbulent Prandtl number $Pr_T$

(5.2a,b)\begin{equation} \kappa_T \equiv Re^s \frac{\varGamma\langle\mathcal{E}\rangle}{\langle \bar{\bar{N}}^2 \rangle} \equiv \varGamma Re_b, \quad Pr_T \equiv \frac{\nu_T}{\kappa_T} \equiv \frac{\langle\mathcal{P}\rangle}{\langle\mathcal{B}\rangle} \frac{\langle\bar{\bar{N}}^2 \rangle}{\langle\bar{\bar{S}}^2\rangle} \approx \frac{\langle \overline{\overline{Ri_g}} \rangle}{R_f} \approx \frac{1+\varGamma}{\varGamma}\, \langle\overline{\overline{Ri_g}}\rangle, \end{equation}

where the buoyancy Reynolds number $Re_b \equiv Re^s \langle \mathcal {E}\rangle / \langle \bar {\bar {N}}^2 \rangle$ is a measure of the ‘turbulence intensity’ that we will return to in § 5.1.3. The first approximation in $Pr_T$ comes from $\langle \bar {\bar {N}}^2 \rangle / \langle \bar {\bar {S}}^2 \rangle \approx \langle \bar {\bar {N}}^2/\bar {\bar {S}}^2 \rangle$, and the second approximation comes from the approximate link between $\varGamma$ and $R_f$ in (3.3), valid under the simplified balance (2.8a) of Osborn (Reference Osborn1980).

Eddy diffusivities can also be expressed using the Prandtl mixing length model, which posits that the turbulent fluxes depend quadratically on the mean gradients

(5.3a–c)\begin{equation} L_\rho^2 \equiv \frac{-\overline{w'\rho'}}{|\partial_z\bar{u}|\partial_z\bar{\rho}} \equiv \frac{\mathcal{B}}{\bar{\bar{S}} \bar{\bar{N}}^2}, \quad L_m^2 \equiv \frac{-\overline{u'w'}}{|\partial_z\bar{u}|\partial_z\bar{u}} \approx \frac{\mathcal{P}}{\bar{\bar{S}}^3}, \quad \frac{L_\rho^2}{L_m^2} = Pr_T, \end{equation}

and therefore that $\kappa _T = L_{\rho }^2 \bar {\bar {S}}$, $\nu _T=L_m^2 \bar {\bar {S}}$ where $L_\rho, L_m$ are the non-dimensional ‘mixing lengths’ for density and momentum, respectively. They can be interpreted as the typical distance travelled by a fluid parcel before its density or momentum becomes mixed with its surroundings (analogous to the mean free path of a molecule in the kinetic theory of gases). The stratified shear flow experiments of Odier et al. (Reference Odier, Chen, Rivera and Ecke2009), Odier, Chen & Ecke (Reference Odier, Chen and Ecke2012) and Znaien et al. (Reference Znaien, Hallez, Moisy, Magnaudet, Hulin, Salin and Hinch2009) showed that $L_\rho,L_m$ were approximately uniform in $z$ (instead of $\kappa _T,\nu _T$), i.e. that the quadratic flux-gradient relationships (5.3a–c) were better approximations than the linear flux-gradient relationships (5.1a,b).

Nevertheless, putting this aside for now and assuming the validity of the widely used eddy diffusivity model (5.1a,b), the key challenge of parameterising $\kappa _T$ (and its related $\nu _T$) using (5.2a,b) becomes equivalent to parameterising the dependence of the indirect (or ‘proximate’) parameter $\varGamma$ (or its related $R_f$) on a few key non-dimensional parameters best characterising the flow, an approach known as the ‘Osborn method’ after Osborn (Reference Osborn1980) (see Salehipour & Peltier (Reference Salehipour and Peltier2015), Gregg et al. (Reference Gregg, D'Asaro, Riley and Kunze2018) and Taylor et al. (Reference Taylor, de Bruyn Kops, Caulfield and Linden2019) for more details). To achieve this, different dynamical balances have been proposed, based on the ratios of relevant length scales or time scales which we discuss next.

5.1.3. Parameters based on length scales and time scales ratios

Further to our definitions in § 3.1.4 of the microscopic Kolmogorov length scale $\ell _K$ (see (3.8)) and Batchelor length scale $\ell _B$ (see (3.9)), we now define the Ozmidov length scale $\ell _O$ and the Corrsin length scale $\ell _C$, which represent the smallest scales at which the distorting influences of background stratification and shear, respectively, are felt (Smyth & Moum Reference Smyth and Moum2000b). Their non-dimensional expressions in shear-layer units are

(5.4a,b)\begin{equation} \ell_O \equiv \left( \frac{\langle \mathcal{E}\rangle}{\langle \bar{\bar{N}}^3\rangle}\right)^{1/2} \equiv \frac{\langle \mathcal{E} \rangle^{1/2}}{ (Ri_b^s)^{3/4} \langle |\partial_z \bar{\rho}|^{3/4}\rangle}, \qquad \ell_C \equiv \left( \frac{\langle \mathcal{E}\rangle}{\langle \bar{\bar{S}}^3\rangle}\right)^{1/2} \equiv \frac{\langle \mathcal{E} \rangle^{1/2}}{ \langle |\partial_z \bar{u}|^{3/2}\rangle}. \end{equation}

Note the subtle fact that the $y-z$ averaging (integration) is made after raising the power in the denominator, contrary to the numerator. This choice is often ambiguous in the literature, and in the following we average $\bar {\bar {N}},\bar {\bar {S}}$ sometimes before, and sometimes after raising the power, for notational convenience. However, this (common) abuse of notation is justifiable in the more strongly turbulent flows I6–T3 in which $\bar {\bar {N}}, \bar {\bar {S}} \approx \textrm {const.}$ (thus the power and integration operators commute with good accuracy). In these flows we have the following separation of scales:

(5.5)\begin{equation} \frac{\ell_O}{\ell_C} \approx \frac{ \langle |\partial_z \bar{u}|^{3/2}\rangle}{\langle |\partial_z \bar{\rho}|^{3/4}\rangle} (Ri_b^s)^{{-}3/4} \rightarrow (Ri_b^s)^{{-}3/4} \rightarrow 5, \end{equation}

using $|\partial _z\bar {u}|\approx |\partial _z\bar {\rho }| \rightarrow 1$ and $Ri_b^s \rightarrow 0.15$. In other words, there exists a moderate range of eddy sizes that are significantly more influenced by shear than by stratification, i.e. the turbulence is slightly dominated by shear.

The separation between the Ozmidov and the Kolmogorov scales is usually quantified by the buoyancy Reynolds number $Re_b$ (first mentioned in (5.2a,b))

(5.6a)\begin{align} Re_b \equiv Re^s \frac{\langle \mathcal{E} \rangle}{\langle \bar{\bar{N}}^2\rangle} = \left(\frac{\ell_O}{\ell_K}\right)^{4/3} &= Re^s (Ri_b^s)^{{-}1}\frac{\langle \mathcal{E} \rangle}{ \langle |\partial_z \bar{\rho}|\rangle} \end{align}

(5.6b)\begin{align} &\rightarrow 0.2 \theta Re^s \end{align}

(5.6c)\begin{align} &\approx 10\text{--}20 \quad \text{for} \ \theta Re^s = 50\text{--}110 . \end{align}

The turbulent estimate (5.6b) assumes: (i) $Ri_b^s\rightarrow 0.15$ (Part 1, figure 2b), (ii) $\langle \mathcal {E}\rangle \rightarrow 0.035 \theta$ in (3.4f) and (iii) $\langle |\partial _z \bar {\rho }|\rangle \rightarrow 1$, the latter being verified to better than $5\,\%$ in I6–T3. This expression is slightly different from that of LPL19 who proposed $Re_b \rightarrow 0.12 \, \theta Re^h$ (see their (6.9) and (6.10)), using the hydraulic Reynolds number (instead of the shear-layer Reynolds number), and averaging data across the whole duct cross-section (instead of the more vigorous ‘core’ shear layer only). Our previous estimate yields (5.6c) in I6–T3 (see figure 2(c)). We conclude from this separation of scales that there exists a range of eddy sizes that are too small to be significantly affected by stratification but too large to be dominated by viscous dissipation, which is a requirement for the existence of stratified turbulent dynamics (although $Re_b=30$ is usually viewed as a minimum threshold).

The indirect measure of mixing $\varGamma$ has often been assumed constant ${\approx}\,0.2$ by physical oceanographers (corresponding to the upper bound set by Osborn (Reference Osborn1980), as mentioned in § 3.1.2). The DNSs of Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) suggested that this constant value was indeed accurate in ‘transitional’ turbulence ($Re_b\approx 7\text{--} 100$), but that $\varGamma \propto Re_b^{-1/2}$ in ‘energetic’ turbulence ($Re_b> 100$); a scaling that has been much debated and reinterpreted since.

It is now widely acknowledged that the challenge of isolating the key non-dimensional parameters controlling turbulent mixing was due to the pervasive tendency for these parameters to be correlated in often-unsuspected, flow-specific and potentially misleading ways. As an example, Maffioli, Brethouwer & Lindbord (Reference Maffioli, Brethouwer and Lindbord2016) and Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) recently argued that $\varGamma$ should not be a function of the (ambiguous) parameter $Re_b$, but of a (more fundamental) turbulent Froude number $Fr$ instead. This Froude number is defined as the ratio of the turbulent kinetic energy dissipation frequency $\mathcal {E}/K'$ to the buoyancy frequency (or the ratio of the buoyancy time scale to the dissipation time scale)

(5.7a)\begin{align} Fr \equiv \frac{\langle \mathcal{E} \rangle}{ \langle K' \rangle \langle \bar{\bar{N}} \rangle} &\equiv (Ri_b^s)^{{-}1/2} \langle |\partial_z\rho|^{{-}1/2} \rangle (\langle K' \rangle/\langle \mathcal{E} \rangle )^{{-}1} \end{align}

(5.7b)\begin{align} &\approx 0.3, \end{align}

using $Ri_b^s\rightarrow 0.15$, $\langle |\partial _z \bar {\rho }|\rangle \rightarrow 1$ and the approximate energy dissipation time scale $K' / \mathcal {E} \approx 10$ observed in § 3.2. Their scaling analyses and triply periodic, spectrally forced DNSs suggest that $\varGamma \approx 0.5 \propto Fr^0$ in strongly stratified flows ($Fr\ll 1$); that $\varGamma \propto Fr^{-1}$ in moderately stratified flows ($Fr \approx 1$); and that $\varGamma \propto Fr^{-2}$ in weakly stratified flows (${Fr\gg 1}$). Note that their argument relies on a definition of $\varGamma$ using the ratio of irreversible components $\chi / \mathcal {E}$, which is only consistent with our definition $\mathcal {B} / \mathcal {E}$ under conditions of ‘lock step’ between $\mathcal {B}$ and $\chi$ explained in § 5.1.1 (asymptotically satisfied at large $\theta Re^s$).

This turbulent Froude number is connected to a further key scale, the Ellison scale

(5.8)\begin{equation} \ell_E \equiv \frac{\langle \rho'^2 \rangle^{1/2} }{\langle |\partial_z \rho|\rangle} \rightarrow 0.07\text{--}0.12, \end{equation}

using $|\partial _z \rho |\rightarrow 1$, and $\langle \rho '^2 \rangle \approx 0.005\text{--} 0.015$ for I6–T3 as observed in figure 1(e). It measures the typical vertical distance travelled by fluid parcels to achieve a stable equilibrium density profile through adiabatic sorting. It is closely related to the Thorpe scale $\ell _T$, defined directly on any instantaneous vertical density profile as the root mean square of these sorting displacements (Mater, Schaad & Venayagamoorthy Reference Mater, Schaad and Venayagamoorthy2013). Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) argued that $Fr \approx (\ell _O/\ell _E)^2$ when $Fr\ll 1$; that $Fr \approx L_O/L_E$ when $Fr \approx 1$; and that $Fr \approx (L_O/L_E)^{2/3}$ when $Fr\gg 1$. Our estimate (5.7b) suggests that I6–T3 are relatively strongly stratified, hence that

(5.9)\begin{equation} \frac{\ell_E}{\ell_O} \approx Fr^{{-}1/2} \rightarrow 2, \end{equation}

i.e. that the separation between the Ellison and Ozmidov scales is very modest.

Figure 9 summarises the above estimates by showing the relative position of these length scales expected in the asymptotic turbulent regime. In addition to the general separation factors between different scales, we also give the corresponding specific values for each scale in data set T2 (non-dimensional value in shear-layer units and dimensional value in mm). In dark grey we highlight scales that are fixed (the reference shear-layer half-height $1$) or slaved by approximately constant parameters ($\ell _E$ in (5.8), $\ell _O$ in (5.9), $\ell _C$ in (5.5) and $\ell _B$ in (3.9)). In blue we highlight scales that are specified by the apparatus and thus subject to change by the experimenter (the duct height $4/h$ and length $120/h$, as well as the PIV vector resolution, where $h,{\textrm {d}\kern0.06em x},{\textrm {d} y},\textrm {d}z,n_x,n_y,n_z$ were given in Part 1, table 1). In maroon we highlight the only scale, $\ell _K$, which is directly controlled through the two key variable flow parameters $\theta, Re^s$ ($\ell _K$ can be visualised as a ‘slider’, unconstrained by other scales). Its definition (3.8) and the use of our turbulent estimate (3.4f) provided a (likely correct) scale separation factor of $\ell _K^{-1} \approx 0.5 \theta ^{1/4} (Re^s)^{3/4}$ (with respect to the shear-layer scale $1$). However, the scaling arguments in this section yield a slightly incompatible $\ell _K^{-1} \approx Re_b^{3/4} Fr^{-1/2} \langle \rho '^{2}\rangle ^{-1/2} \approx 20 \, Re_b^{3/4} \approx 6 \theta ^{3/4} (Re^s)^{3/4}$ (combining three factors to reach the shear-layer scale $1$). This (likely incorrect) scaling in $\theta ^{3/4}$ can likely be explained by a weak, neglected dependence of $Fr$ and $\langle \rho '^{2}\rangle$ on $\theta$ (i.e. $\ell _O$ and $\ell _E$ are not exactly constant).

Figure 9. Key length scales in SID turbulence. Relative positions and scale separation factors are based on the simple volume-averaged estimates of § 5.1.3 in the asymptotic turbulent regime ($\theta Re^s \gg 100$). Note the logarithmic axis. Non-dimensional values (in shear-layer variables) and dimensional values (in mm) are for data set T2. The Kolmogorov scale $\ell _K$ (in red) has two slightly incompatible scalings (3.8) and (5.6), flagged by ‘?’ (the lighter shade highlights the less trustworthy (5.6)). Ideally the PIV resolution would approach $\ell _K$, and the LIF resolution (here assumed equal to the PIV) would approach $\ell _B$.

This figure highlights that in the SID experiment, increasing $\theta$ and $Re$ (in red) will lead to a decreasing Kolmogorov scale $\ell _K$, while (importantly) keeping the other scales above it approximately constant (in grey). As a result, ratios such as $\ell _C/\ell _K$ and $\ell _O/\ell _K$ will increase with $\theta Re$, and SID turbulence will become increasingly isotropic.

5.2. Eddy diffusivities, mixing lengths, turbulent Prandtl number

In figure 10 we test the flux-gradient relations (5.1), (5.3) with the full clouds of $n_y n_z$ data points (left three columns). Linear fits with enforced zero intercept are plotted in blue, and provide the eddy diffusivities $\nu _T, \kappa _T$, while quadratic fits with enforced zero intercept are shown in purple, and provide the mixing lengths $L_m,L_\rho$. These ‘fit’ values are then plotted in the rightmost column.

Figure 10. Eddy diffusivities and mixing lengths in data sets I6–T3 (top to bottom row). Clouds of $n_y n_z$ points of $\overline {u'w'} \equiv \mathcal {P}/\bar {\bar {S}}$ vs $-\partial _z \bar {u} \equiv \bar {\bar {S}}$ (first column); $\overline {w'\rho '} \equiv \mathcal {B}/Ri_b^s$ vs $-\partial _z \bar {\rho } \equiv \bar{\bar{N}}^2/Ri_b^s$ (second column); $\overline {w'\rho '} \equiv \mathcal {B}/Ri_b^s$ vs $\partial _z \bar {u} \, \partial _z \bar {\rho }\equiv \bar {\bar {S}} \bar {\bar {N}}^2/Ri_b^s$ (third column). Note the log–log axes, and symbol colour and size respectively indicating the $|z|$ and $|y|$ location. Linear and quadratic least-squares fits provide the eddy diffusivities (in blue, after multiplying by $Re^s$ as in (5.1a,b)) and mixing lengths (in purple). Diamonds show the volume-averaged values of the flux vs gradient (blue) or square gradient (purple). (s–x) The $\nu _T,\kappa _T,L_m,L_\rho$ values (s–t) against one another, giving the ratio $Pr_T=\kappa _T/\nu _T = L_\rho ^2/L_m^2$ ($Pr_T= 1, 3, 10$ shown); (u–x) against the input parameters $\theta Re^s$. Values obtained from the fit are indistinguishable from those obtained from the diamonds.

First, focusing on the left three columns, we find that the clouds generally have a wide spread, making the fits fairly poor (also note the log–log axes). Despite this spread, the fits capture a clear monotonic tendency, particularly visible in the upper boundary of each cloud (high flux values) which are indeed bounded by an approximately linear or quadratic flux relation. The symbol colours, indicating $|z|$, reveal that these high flux values tend to occur close to the mid-point of the shear layer ($|z|\approx 0$, dark colour), though less so in the buoyancy flux (second column). The symbol sizes, inversely proportional to $|y|$, do not reveal any clear correlation between flux-gradient behaviour and spanwise location, other than the fact that $|y|$ contributes to the spread of the clouds. Although the coefficients of determination are generally very low ($r^2<0.2$), the constant eddy diffusivity model (linear fit) does slightly better than the constant mixing length model (quadratic fit) overall. However, this is not very significant because uniform eddy diffusivities and uniform mixing lengths actually become compatible in our asymptotic case of uniform shear $\bar {\bar {S}} \approx 1$, since by definition $(\kappa _T/Re^s) \equiv L_\rho ^2 \, \bar {\bar {S}}$ and $(\nu _T/Re^s) \equiv L_m^2 \, \bar {\bar {S}}$. In other words, our range of $\bar {\bar {S}} \equiv -\partial _z \bar {u} \approx 0.5\text{--} 1.5$ (see left column) is not wide enough to convincingly argue in favour of either model.

Second, the diamond symbols show the volume average of the flux – the numerator – against the volume average of the gradient (in blue) or square gradient (in purple) – the denominator. As expected from our above comment that uniform eddy diffusivities and mixing lengths are compatible, blue and purple diamonds lie close to one another, near the horizontal values of $1 \approx |\partial_z \bar{u}| \approx |\partial_z \bar{u}|^2 \approx |\partial_z \bar{\rho}| \approx \partial_z \bar{u} \partial_z \bar{\rho} $. Moreover, most diamonds sit very close to the fits (lines) of their respective colour in the left two columns, but consistently above them in the third column. This proves that the definitions of $\nu _T,\kappa _T,L_m$ by volume averages would produce good approximations of the fit of the underlying distribution (i.e. the fit goes through the centre of mass of the cloud), whereas the definition of $L_\rho$ by volume averages would produce an overestimation.

Third, moving on to the rightmost column, we find good correlations $\kappa _T \propto \nu _T$ (panel s) and $L_\rho \propto L_m$ (panel t), corresponding to a constant turbulent Prandtl number $Pr_T \approx 3$ (dashed line), except in I8 which has $Pr_T\approx 7$. This is entirely consistent with the approximation in (5.2) and our previously quoted asymptotic values of $\overline {\overline {Ri_g}} \approx 0.15$ (Part 1, § 5) and $R_f\approx 0.05$ (§ 3.1.2) giving $Pr_T\approx 3$. This value is comfortably above 1, despite the tendency to self-similarity of the mean velocity and density profiles observed in ${\rm T}$ flows ($\langle \bar {u} \rangle _y(z) \approx \langle \bar {\rho } \rangle _y(z)$, see Part 1, figure 3). This value is, however, consistent with the DNSs of Salehipour & Peltier (Reference Salehipour and Peltier2015) (see their figure 10, at higher $Re^s$ but similar $Ri_b^s$) who found $Pr_T \approx 3$ at $Re_b\approx 5\text{--} 15$. Actual values for the diffusivities range from $\nu _T\approx 1$ in I6–T1 to $\nu _T\approx 3$ in T2–T3, a substantial but not overwhelming increase with respect to the molecular value for momentum $\nu$. The corresponding range $\kappa _T Pr \approx 700/3\text{--} 700$ is, by contrast, an overwhelming increase with respect to the molecular value for density $\kappa$, i.e. a high ‘eddy Péclet number’. Mixing lengths $L_m,L_\rho$ are of the order of the Kolmogorov length $\ell _K$ (see estimate in § 3.1.4) and of the resolution of our measurements in $x,z$ (see Part 1, table 3). Finally, all quantities typically increase monotonically with $\theta Re^s$ (panels u–x), though T1 is an outlier that appears less energetic than suggested by its $\theta Re^s$ value. We conclude that $\nu _T,\kappa _T$ appear linear or superlinear in $\theta$ and $Re^s$.

5.3. Flux coefficient, flux Richardson number and $\mathcal {B}/\mathcal {P}_\rho$

In figure 11, we test the relations in (3.1) with clouds of $n_y n_z$ data points (left three columns). The diamond coordinates are given by the numerator and denominator (volume averages) of (3.1); these were already plotted for all 16 data sets in figure 2(j–l). Linear fits with enforced zero intercept are also shown in blue, and provide the values for $\varGamma$, $R_f$, $\mathcal {B}/\mathcal {P}_\rho$ plotted in the rightmost column.

Figure 11. Flux coefficient, flux Richardson number and $\mathcal {B}/\mathcal {P}_\rho$ ratio in data sets I6–T3 (top to bottom rows). Clouds of $n_y n_z$ points of the numerator $\mathcal {B}$ vs the respective denominator: $\bar {\mathcal {E}}$ (first column); $\mathcal {P}$ (second column); $\mathcal {P}_\rho$ (third column). Log–log axes with identical vertical axis for all panels (a–r). Symbol styles and diamonds are as in figure 10. (s–v) $\varGamma, R_f,\mathcal {B}/\mathcal {P}_\rho$ vs $\theta Re^s$ (values obtained from the linear fit or from the diamonds are indistinguishable).

First, focusing on the left three columns, we find that the clouds of the leftmost column have the largest spread, followed by those of the second column, and finally those of the third column, which are tighter around the fit. As a result, though the linear fits capture a clear trend, a constant $\varGamma$ is a relatively poor model (mean $r^2=0.30$), while constant $R_f$ and $\mathcal {B}/\mathcal {P}_\rho$ are better models (mean $r^2=0.63$ and 0.66, respectively). Besides, the symbol colours or sizes do not reveal any clear pattern between this behaviour and the position $|z|,|y|$ within the shear layer. The diamonds generally lie very close to the fit, which means that our previous volume-averaged estimations of figure 2(j–l) ($\varGamma \approx 0.1, R_f\approx 0.05, \langle \mathcal {B} \rangle /\langle \mathcal {P}_\rho \rangle \approx 1$) were good approximations.

Second, the values in the rightmost column confirm indeed that $R_f \approx \varGamma /2$ (panel s, dotted line), which we recall is qualitatively sensible ($R_f<\varGamma$) but quantitatively inconsistent with (3.3). As explained in § 3.1.3, this is likely due to a systematic underestimation of $\langle \mathcal {E}\rangle$, making $\varGamma \approx R_f \approx 0.05\text{--} 0.07$ more likely at high $\theta Re^s$ (panels t–u). Moreover, we confirm that our data suggest $\mathcal {B}/\mathcal {P}_\rho \rightarrow 1$ at high $\theta Re^s$ (panel v), a necessary condition for the attractive ‘lock step’ between $\mathcal {B}$ and $\chi$.

5.4. Buoyancy Reynolds number, turbulent Froude number

In figure 12, we test the relations in (5.6a) (left two columns) and (5.7a) (right two columns) with clouds of $n_y n_z$ data points, as in the previous two figures, to explore the extent to which the definitions of $Re_b$ and $Fr$ (based on mean quantities) are satisfied pointwise. Linear fits with enforced zero intercept are also shown in blue, and provide the values for $Re_b,Fr$ plotted in the bottom row against $\theta Re^s$ to test (5.6b) and (5.7b), respectively.

Figure 12. Buoyancy Reynolds number and turbulent Froude number in data sets I6–T3. Left two columns (a–f) numerator $Re^s\bar {\mathcal {E}}$ vs denominator $\bar {\bar {N}}^2$. Right two columns (h–m) numerator $\bar {\mathcal {E}}/\overline {K'}$ vs denominator $\bar {\bar {N}}$. Log–log axes, symbol styles and diamonds are as in figures 10 and 11. Bottom row: (g) $Re_b$ vs $\theta Re^s$ to test (5.6b) (dashed line), showing with empty symbols the values obtained from the fit, and with full symbols the values obtained from the diamonds; (n) $Fr$ vs $\theta Re^s$ to test (5.7b) (fit and diamonds values are indistinguishable).

First, we find that all clouds have a large spread around the fit. The linear fit in the left two columns (panels a–f) captures a trend, especially the shape of the yellow cloud ($|z|\lesssim 1$ points at the edges of the shear layer), but a uniform $Re_b$ remains a poor model. Although there is no reason to expect a uniform $Re_b$ in general, single ‘representative’ values of $Re_b$ based on averaged quantities are known to work well to characterise the flow regime in ocean mixing (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018). The linear fit in the right two columns (panels h–m), however, even fails to capture the trend, arguing against a uniform $Fr$ model. These criticisms should be nuanced by the observation that the clouds in panels (h–m) are very compact and span a very limited range (less than a decade in the horizontal and vertical axes). This limited range reveals an asymptotic tendency to uniform linear stratification (denominator), and to uniform dissipation frequency (numerator). The corresponding turbulent kinetic energy dissipation time scale $(\bar {\mathcal {E}}/\overline {K'})^{-1} \rightarrow 10$ A.T.U. (see the vertical coordinate of the diamonds), with some scatter in I6–T1 (${\approx }3\text{--} 30$ A.T.U.), but much less scatter in T2–T3 (${\approx }5\text{--} 12$ A.T.U.), suggesting some form of turbulent self-organisation.

Second, this observation allows us to deduce the asymptotic turbulent scaling $\langle K' \rangle \rightarrow 10 \langle \mathcal {E} \rangle \rightarrow 0.35 \theta$ (using the scaling for $\langle \mathcal {E} \rangle$ in (3.4)). As a result, the parameterisation of eddy diffusivities $\nu _T,\kappa _T \propto Re^s \langle K' \rangle / \langle \bar {\bar {S}} \rangle$ proposed by van Reeuwijk, Holzner & Caulfield (Reference van Reeuwijk, Holzner and Caulfield2019) in DNSs of inclined gravity currents yields $\nu _T,\kappa _T \propto \theta Re^s$ (using $\langle \bar {\bar {S}} \rangle \approx 1$, and our non-dimensionalisation of $\nu _T,\kappa _T$ by the molecular value $\nu$). This scaling appears compatible with our data in figure 10(u,v), although prefactors do not match.

Third, returning to figure 12, we find that, despite the spread of the clouds, the linear fits (blue lines) must (by construction) approximately go through the centre of mass of each cloud (blue diamonds), giving indistinguishable values of $Fr$ in panel (n). These data confirm our estimate (5.7b) of an approximately constant $Fr\approx 0.3$. The values of $Re_b$ obtained from the fit and from the diamonds are, however, slightly distinguishable, and shown using empty and full symbols respectively in panel (g). These data suggest an approximate scaling $Re_b \approx 0.1 \theta Re^s$ (dotted line), with volume-averaged values (full symbols) being consistently higher. This scaling is consistent with our estimate $Re_b \approx 0.2 \theta Re^s$ in (5.6b) (dashed line) if we again invoke the systematic underestimation of $\langle \mathcal {E}\rangle$ by a factor of 2.

Finally, let us consider the ratio $\kappa _T/Re_b$, using our $\kappa _T$ data from figure 10(v) and our $Re_b$ data from figure 12(g), to test the parameterisation of $\kappa _T$ using $Re_b$. Although not plotted here for conciseness, we verified that this ratio is very close to our $\varGamma$ data in figure 11(t) (as expected from the definition (5.2)), with values tightly grouped between 0.093 and 0.12 (except for I8 which has a ratio of 0.043). This suggests that the parameterisation $\kappa _T\approx \varGamma Re_b$ is a good model for our data, even for $3 \lesssim Re_b\lesssim 12$, qualitatively supporting Osborn's model (although our $\varGamma$ values are a factor of 2 or 3 below the upper bound 0.2 commonly used). However, it disagrees with the parameterisation of Bouffard & Boegman (Reference Bouffard and Boegman2013) (see their tables 1 and 2), who argue that Osborn's model should only be expected in the ‘transitional’ and ‘energetic’ regimes at $Re_b \gtrsim (3 \ln \sqrt {Pr})^2 \approx 100$ in salt-stratified flows, whereas $\kappa _T \approx 0.1 \, Pr^{-1/4} Re_b^{3/2}$ should be expected in the ‘buoyancy-controlled’ regime below this threshold (a scaling that we do not observe in our data).