Experimental properties of continuously forced, shear-driven, stratified turbulence. Part 2. Energetics, anisotropy, parameterisation

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