2.1. Parametric framework

HSST is assumed to satisfy the Navier–Stokes equations subject to the non-hydrostatic Boussinesq approximation. The dimensional equations for the fluctuations relative to the planar means are

(2.1a)\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u} ={-} \frac{1}{\rho_0} \boldsymbol{\nabla} p - z \frac{\partial \boldsymbol{u}}{\partial x} \frac{{{\rm d}} \bar{u}_x}{{{\rm d}}z} - u_z \frac{{{\rm d}} \bar{u}_x}{{{\rm d}}z} \hat{\boldsymbol{x}}- \frac{g_z}{\rho_0}\rho\hat{\boldsymbol{z}} + \nu \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}, \end{gather}

(2.1b)\begin{gather}\frac{\partial \rho}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \rho ={-} z \frac{\partial \boldsymbol{\rho}}{\partial x} \frac{{{\rm d}} \bar{u}_x}{{{\rm d}}z} - u_z\frac{{{\rm d}} \bar{\rho}}{{{\rm d}}z} + D \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \rho , \end{gather}

(2.1c)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

where $\boldsymbol {u} = (u_x, u_y, u_z)$ is the velocity vector in the coordinate system $(x,y,z)$ and functional dependencies have been dropped for convenience, and $\rho$ is the fluctuating density field so that the total density is

(2.2)\begin{equation} \rho_t = \rho_0 + \bar{\rho}(z) + \rho(x,y,z,t), \end{equation}

with $\bar {\rho }(z) = z \, {{\rm d}} \bar {\rho }/ {{\rm d}}z$ being the ambient density that is constant in time and has a uniform gradient antiparallel to gravitational acceleration $g_z$. Similarly, $\bar {u}_z= z \, {{\rm d}} \bar {u}_x/ {{\rm d}}z=z S$ is the mean velocity and $S\equiv {{\rm d}} \bar {u}_x/ {{\rm d}}z$ is the mean shear. The pressure decomposition is analogous to that of density and velocity so that $p$ is the fluctuating mechanical pressure relative to the temporally constant planar mean. Internal energy and scalar transport affecting density have been combined into a single evolution equation for $\rho$ with $D$ the molecular diffusivity, $\hat {\boldsymbol {x}}$ and $\hat {\boldsymbol {z}}$ are the unit vectors in the $x$- and $z$-directions respectively and $\nu$ is the kinematic viscosity.

2.1.1. Parametrisation

Dimensional analysis suggests that (2.1) can be described in terms of at least four non-dimensional parameters. As we wish to interpret the dynamics in terms of the dynamic range available for various aspects of the flow, we consider parametrisation in terms of the ratios of length scales.

Here, we use the convention of defining an approximate outer scale of turbulence, the large-eddy length scale, and the (viscous) Kolmogorov length scale as

(2.3a,b)\begin{equation} L_{LE} \equiv E_k^{3/2}/\epsilon \quad \text{and}\quad L_K \equiv (\nu^{3}/\epsilon)^{1/4} , \end{equation}

where $E_k \equiv \langle \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {u} \rangle /2$, $\epsilon \equiv \nu \langle \boldsymbol {\nabla } \boldsymbol {u} : \boldsymbol {\nabla } \boldsymbol {u} \rangle$, the operator ‘$:$’ denotes the double inner product and the notation $\langle \boldsymbol {\cdot } \rangle$ indicates an ensemble average.

The (buoyancy) Ozmidov length scale

(2.4)\begin{equation} L_O\equiv(\epsilon / N^{3})^{1/2},\end{equation}

defines the lower limit of length scales significantly affected by buoyancy forces, where $N$ is the (background) buoyancy frequency, defined here as

(2.5)\begin{equation} N^{2} \equiv{-}\frac{g_z}{\rho_0} \frac{{{\rm d}} \bar{\rho}}{{{\rm d}}z} . \end{equation}

With an explicit shear scale applied in HSST, the Corrsin length scale (Corrsin Reference Corrsin1958) is thought to characterise the lower limit of scales affected by shear, and is defined as

(2.6)\begin{equation} L_C \equiv (\epsilon / S^{3})^{1/2}.\end{equation}

Scales smaller than $L_C$ have been suggested to define a locally isotropic regime of length scales in the absence of additional smaller-scale affects of the mean flow (Saddoughi & Veeravalli Reference Saddoughi and Veeravalli1994).

Given the scales defined in (2.3a,b), (2.4) and (2.6), we choose to describe turbulence by the following three parameters:

(2.7a–c)\begin{align} {{Re_s}} \equiv \frac{\epsilon}{\nu S^{2}}= \left( \frac{L_C}{L_K}\right)^{4/3}, \quad {{Ri}} \equiv \frac{N^{2}}{S^{2}} =\left( \frac{L_C}{L_O}\right)^{4/3}\quad \text{and}\quad {{Fr}} \equiv \frac{\epsilon}{N E_k}= \left( \frac{L_O}{L_{LE}}\right)^{2/3},\end{align}

in a fluid with a fixed Prandtl number $Pr\equiv \nu /D$ as the required fourth parameter (see Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014) for more details on parameterisations). We choose the shear Reynolds number, $Re_s$, since it describes the range of length scales associated with isotropy for stationary flows with Richardson number, $Ri \leq 1$ (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014). Alternative Reynolds numbers can also be useful to consider for comparison with other flow regimes, where the turbulent Reynolds number and buoyancy Reynolds numbers are defined, respectively as

(2.8a,b)\begin{equation} {{Re}} \equiv \frac{E_k^{2}}{\nu \epsilon} =\left(\frac{L_{LE}}{L_K}\right)^{4/3} , \quad {{Re_b}} \equiv \frac{\epsilon}{\nu N^{2}} =\left(\frac{L_O}{L_K}\right)^{4/3} . \end{equation}

In flows driven by mean shear, the energetic stationarity of the flow is thought to be characterised by $Ri$ (Jacobitz, Sarkar & Van Atta Reference Jacobitz, Sarkar and Van Atta1997), measuring the scale separation between the driving shear scales and stabilising buoyancy scales. It is thought that turbulent flow configurations persist where $Ri<1$, where we distinguish between sustained turbulent flows and results derived from or interpreted in terms of linear hydrodynamic instability theory (see Zhou et al. (Reference Zhou, Taylor, Caulfield and Linden2017) for more discussion). Indeed, stratified turbulent simulations with prescribed shear scales have revealed that the characteristic gradient Richardson number associated with stationarity has been observed to take a Reynolds number insensitive value of

(2.9)\begin{equation} {{Ri}} \approx 0.15 - 0.2 ,\end{equation}

as observed in stratified homogeneous shear flow (Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000; Holt, Koseff & Ferziger Reference Holt, Koseff and Ferziger1992; Portwood, de Bruyn Kops & Caulfield Reference Portwood, de Bruyn Kops and Caulfield2019).

Furthermore, in shear-dominated flows, energy production due to shear naturally induces a strong coupling between shear and outer length scales such that the shear parameter $ {{S_\ast }}$, here defined in terms of $E_k$ (as opposed to in terms of $q\equiv \langle \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {u} \rangle$) as

(2.10)\begin{equation} {{S_\ast}} \equiv \frac{E_k S}{\epsilon}=\left( \frac{L_{LE}}{L_C}\right)^{2/3}, \end{equation}

tends to a constant between 5 and 6 at sufficiently high Reynolds numbers (Jacobitz et al. Reference Jacobitz, Sarkar and Van Atta1997; Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000)

However, the foregoing turbulent parameters are likely insufficient to characterise active scalar dynamics, a specific point we wish to investigate in this paper. The smallest scales of the scalar are thought to be characterised by the Batchelor length scale $L_B$,

(2.11)\begin{equation} L_B \equiv (\nu D^{2} / \epsilon)^{1/4}, \end{equation}

which coincides with $L_K$, at unity Prandtl number as considered in this work.

Obukhov–Corrsin similarity, which is either applicable at high Froude number or at scales smaller than those substantially affected by buoyancy, i.e. for scales smaller than $L_O$, implies an outer-scale scalar length scale

(2.12)\begin{equation} L_{OC} \equiv E_p^{3/2} \epsilon^{1/2} / \chi^{3/2},\end{equation}

where $E_p$ is the available potential energy and $\chi$ is its irreversible dissipation rate, defined in the linearly stratified limit to be

(2.13a,b)\begin{equation} E_p \equiv \left \langle \frac{g^{2}/\rho_0^{2}}{2N^{2}}\rho^{2} \right \rangle \quad \text{and}\quad \chi \equiv D \left \langle \frac{g^{2}/\rho_0^{2}}{N^{2}} | \boldsymbol{\nabla} \rho | ^{2} \right \rangle, \end{equation}

respectively, demonstrating that $\chi$ in this context may also be interpreted as the scaled destruction rate of buoyancy variance (Caulfield Reference Caulfield2021). Therefore, the range of anisotropic scales of the scalar are characterised by

(2.14)\begin{equation} N_\ast{\equiv} \left( \frac{L_{OC}}{L_O} \right)^{2/3}= \frac{N E_p}{\chi}. \end{equation}

Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019) observed this parameter to approach unity in the high $Re_b$ limit for flows where $Pr=1$, and Mater & Venayagamoorthy (Reference Mater and Venayagamoorthy2014) suggest $N_\ast$ more descriptively parameterises turbulent mixing.

2.2. A mixing model constructed from analysis of subrange scaling

Local isotropy is a state wherein statistical symmetries of multi-point statistics, such as homogeneity, isotropy and stationarity, are present in a spatio-temporally localised region (Monin & Yaglom Reference Monin and Yaglom1975). In the presence of anisotropic integral scales, the conditions of local isotropy are defined, in the spatial sense, by a quasi-equilibrium range of scales that are sufficiently smaller than integral scales (Kolmogorov Reference Kolmogorov1941).

This is a prerequisite of classical inertial subrange scaling arguments (Kolmogorov Reference Kolmogorov1941; Oboukhov Reference Oboukhov1941; Corrsin Reference Corrsin1951). These state that scales within the quasi-equilibrium regime exist wherein turbulence is independent of the effects of viscosity. In shear flows, conditions for the application of local isotropy are thought to be valid at scales smaller than $L_C$ (Uberoi Reference Uberoi1957; Corrsin Reference Corrsin1958). Kolmogorov scaling (Kolmogorov Reference Kolmogorov1941,  Reference Kolmogorov1962) in the presence of anisotropic outer scales has been revealed to coincide at scales consistent with local isotropy (Champagne et al. Reference Champagne, Harris and Corrsin1970; Saddoughi & Veeravalli Reference Saddoughi and Veeravalli1994). The equivalent local isotropy justification of the Kolmogorov subrange may be applied to Obukhov–Corrsin passive scalar scaling (Oboukhov Reference Oboukhov1949; Corrsin Reference Corrsin1951) such that it is fit to describe an active scalar below scales affected by stratification (Monin & Yaglom Reference Monin and Yaglom1975, p. 391). Therefore, the spectral density of potential energy is expected to be determined by the energy dissipation rates and a wavenumber, i.e.

(2.15)\begin{equation} \hat{E}_{p}^{OC}(|\boldsymbol{k}|) = \beta \chi \epsilon^{{-}1/3} |\boldsymbol{k}|^{{-}5/3} f_L^{OC}(|\boldsymbol{k}| L_{C})f_v^{OC}(|\boldsymbol{k}| L_B), \end{equation}

where $\boldsymbol {k}$ is the wavenumber vector, $\beta$ is the Obukhov–Corrsin constant, the functions denoted by $f^{OC}_v$, $f^{OC}_L$ are the viscous and outer-scale correction functions that are necessary without applying assumptions of local isotropy. Furthermore, according to Oboukhov (Reference Oboukhov1941), the kinetic energy spectra

(2.16)\begin{equation} \hat{E}_{k}^{K}(|\boldsymbol{k}|) = C \epsilon^{2/3} |\boldsymbol{k}|^{{-}5/3} f_L^{K}(|\boldsymbol{k}| L_{C})f_v^{K}(|\boldsymbol{k}| L_K), \end{equation}

where $C$ is the Kolmogorov constant, and $f^{K}_L$ and $f^{K}_v$ are the correction functions. Following the analysis of Beguier et al. (Reference Beguier, Dekeyser and Launder1978), by integration of (2.15) and (2.16) when $L_C \gg |\boldsymbol {k}|^{-1} \gg L_K, L_B$, it is possible to obtain

(2.17)\begin{equation} \varPi \equiv \frac{E_p \epsilon}{E_k \chi}= \frac{\beta}{C} .\end{equation}

Therefore, the explicit relation for modelling the irreversible flux from (2.17) is

(2.18)\begin{equation} j_b \approx \chi \approx \frac{E_p \epsilon}{E_k \varPi},\end{equation}

with the limiting behaviour, at high Reynolds number, that $\varPi \rightarrow \beta /C$. A critical simplifying assumption in the development of (2.18) is that $ {{Fr}} \geq O(1)$, although the imposition of more complex scaling relations, which are thought to develop when $ {{Fr}} \ll O(1)$, may alternatively be imposed. Indeed, more generally,

(2.19)\begin{equation} \varPi = \left( \frac{L_{OC}}{L_{LE}} \right)^{2/3} = \frac{N_\ast}{{{S_\ast}} \sqrt{{{Ri}}}} = N_\ast {{Fr}} . \end{equation}

We stress that the assumed model spectra (2.15) and (2.16) are approximations which are difficult to observe precisely in carefully controlled experiments or simulations. Whereas these two-point scaling models are observed in flows with anisotropic large scales, they are not claimed to be generally applicable for flows in arbitrary parameter regimes. However, the sensitivity of $\varPi$ with respect to measured deviations from the model spectra is unclear a priori, a point we discuss in § 4.

The relation (2.18) has profound implications for the calibration of a number of turbulent stratified mixing models. Generally, in gradient diffusion models,

(2.20)\begin{equation} \kappa_b \approx \frac{E_p }{\varPi N} {{Fr}}, \end{equation}

while gradient diffusion subject to the model proposed by Osborn (Reference Osborn1980) leads to

(2.21)\begin{equation} \varGamma_b \approx \varGamma \equiv \frac{\chi}{\epsilon} = \frac{E_p}{E_k \varPi} . \end{equation}

Combining the turbulent viscosity model of Crawford (Reference Crawford1982) with the turbulent diffusivity model of Osborn (Reference Osborn1980), the turbulent Prandtl number may then be expressed as

(2.22)\begin{equation} {{Pr_t}} = \frac{1+\varGamma_b}{\varGamma_b} {{Ri}} = \frac{E_k \varPi+ E_p} {E_p}{{Ri}} .\end{equation}

Alternatively, the relevance of the parameter $\varPi$ to the turbulent Prandtl number model of Venayagamoorthy & Stretch (Reference Venayagamoorthy and Stretch2010) implies that their key mixing parameter $\gamma$ may be expressed as

(2.23)\begin{equation} \gamma = \frac{\varPi}{2}, \end{equation}

where $\gamma$ too has been empirically observed to assume relatively insensitive values at moderate Reynolds numbers in HSST (Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2006). Finally, if the objective is to construct a mixing length scale for applications in Prandtl mixing length models, such as the framework of Odier et al. (Reference Odier, Chen, Rivera and Ecke2009), where $\kappa _b = l_m^{2} S$, (2.18) may be equivalently expressed as

(2.24)\begin{equation} l_m = \frac{E_p }{\varPi N^{2}} {{S_\ast}} . \end{equation}

3.1. Numerical method

The homogeneous stratified shear system defined by (2.1a) and (2.1b) is solved in a triply periodic domain with a Fourier pseudo-spectral method. Timestepping is performed with a fractional step method where nonlinear terms are advanced with a third-order Adams–Bashforth scheme. Aliasing is treated by a sharp-spectral filter at $15/16 k_{max}$, which proved to be sufficient in a posteriori analysis. Temporal integration of the inhomogeneous shear term is handled by the integrating factor approach of Sekimoto, Dong & Jiménez (Reference Sekimoto, Dong and Jiménez2016). The rest of the solver is derived from the method used in Hebert & de Bruyn Kops (Reference Hebert and de Bruyn Kops2006), with further details discussed therein.

The suite of simulations are designed to resolve sufficiently both the smallest and largest spatio-temporal turbulent flow scales. The streamwise extent, $L_x$, is chosen to be 40 times larger than the large-eddy length scale. Because of the anisotropy of integral length scales that was observed in a limited parameter space study, the cross-stream domain length $L_y$ and the vertical domain length $L_z$ are chosen such that $L_x/L_y=2$ and $L_x/L_z=4$. A visualisation of the computational domain is featured in figure 1.

Figure 1. Illustration of simulation domain and contour plots of density. The three left-most panels are slices of total density on an $x-y$ slice, normalised by three times the variance of the fluctuating density in case R1; panels are shown for cases R1, R5 and R9. Note the that fluctuations decrease with Reynolds number while small-scale structure increases. Upper-right panel indicates total density normalised by minimum and maximum values for case R9, illustrating inclined large-scale structure (cf.Chung & Matheou Reference Chung and Matheou2012; Jacobitz & Moreau Reference Jacobitz and Moreau2016).

The smallest spatial scales are resolved by enforcing $k_{max}L_{K} \approx 2$, where $k_\text {max}$ indicates the maximum wavenumber supported by the domain discretisation when fully resolved.

Finally, we have adopted a mass–spring–damper-like system (as similarly adopted in Overholt & Pope Reference Overholt and Pope1998; Rao & de Bruyn Kops Reference Rao and de Bruyn Kops2011), to tune the Richardson number to its stationary value by setting a constant kinetic energy target $E_t$. We note that this is similar to the constant power criteria as used in stationary HSST by Chung & Matheou (Reference Chung and Matheou2012) but here we include damping. The system is controlled by time-dependent variation of $Ri$ being required to satisfy the equation

(3.1)\begin{equation} c_0 S {{Ri}}'(t)+2\alpha \omega \tilde{E}_k'(t) + \omega^{2} (\tilde{E}_k(t)-1)=0, \end{equation}

where the prime notation denotes a temporal derivative, $\tilde {E}_k(t)\equiv E_k/E_t$ is the normalised turbulent kinetic energy, $\omega$ is the characteristic frequency of oscillation, $\alpha$ is a dimensionless damping factor and $c_0$ is a dimensionless parameter. This control system has been derived by assuming that the kinetic energy follows a second-order linear system (e.g. Rao & de Bruyn Kops Reference Rao and de Bruyn Kops2011), and then by applying a first-order approximation to the temporal evolution of kinetic energy (cf. Jacobitz et al. Reference Jacobitz, Sarkar and Van Atta1997)

(3.2)\begin{equation} \tilde{E}_k'(t)\approx c_0 S ({{Ri}}(t)-{{Ri}}_c) \hbox{ such that } \tilde{E}_k''(t) \approx c_0 S {{Ri}}'(t). \end{equation}

The choice for the parameter $c_0\approx -1$ is supported by Jacobitz et al. (Reference Jacobitz, Sarkar and Van Atta1997), the characteristic frequency $\omega$ is determined by the large-eddy time scale which itself is proportional to $1/S$ (Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000) and finally a damping coefficient $\alpha =1.5$ was found to work well.

3.2. Summary of experiments

The stationarity constraint and the emergent empirical observation that $S_\ast \approx 5$ lead to the natural consideration of sets of simulations which can test the effects of variation of the Reynolds number, essentially independently of the other parameters. The simulations presented here are largely equivalent to those reported in Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019), although some cases have been run for longer times to ensure convergence of statistics. The various parameters of the simulations are summarised in table 1, and some flow visualisations of the density field are shown in figure 1. Two key aspects are apparent, as noted in the caption. First, in the horizontal plane visualisations, as the Reynolds number increases, the magnitude of the density fluctuations tends to decrease, while, unsurprisingly the dynamic range of scales increases, with enhanced smaller-scale structures. Second, in the vertical plane visualisation, inclined large-scale structures are apparent, consistently with the results of Chung & Matheou (Reference Chung and Matheou2012) and Jacobitz & Moreau (Reference Jacobitz and Moreau2016).

Table 1. Simulation parameters. Other parameters may be inferred from definitions in the previous section. Here, $N_x$ is the number of streamwise grid points in the equispaced domain discretisation. All cases have $ {{Pr}} \equiv \nu /D = 1$.

Time series are shown in figure 2 to illustrate convergence. Richardson numbers tend to fluctuate within approximately 10 % of their means, at a confidence interval at least as precise as critical Richardson numbers reported by Jacobitz et al. (Reference Jacobitz, Sarkar and Van Atta1997) and Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000). The fluctuations in the Richardson number occur at large integral time scales such that effects of $Ri'(t)$ are small. Energetic fluctuations are more substantial than higher-dimensional forcing schemes (i.e. Overholt & Pope Reference Overholt and Pope1998; Rao & de Bruyn Kops Reference Rao and de Bruyn Kops2011) that dynamically control the flow with more granularity by accessing a broader range of turbulent length scales compared with a single mean-scale parameter, as is done here. Nonetheless, energy remains approximately stationary over the large time scales of the simulations.

Figure 2. (a) Sample time histories of the relative kinetic energy fluctuation about the target energy and the Richardson number. Smaller fluctuations with increasing Reynolds number were observed, but generally, the kinetic energy remains within 15 % of its target. (b) Longitudinal streamwise velocity spectra for case R5 sampled throughout the run period illustrating small spectral fluctuations, where the most significant fluctuations are about the smallest wavenumbers.

Recalling that the critical Richardson number is an emergent quantity from the control system defined by (3.1), and is not determined a priori, an increase of the critical Richardson number is not observed in these simulations as suggested by Holt et al. (Reference Holt, Koseff and Ferziger1992), Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000), although the emergent critical values of the Richardson numbers are broadly consistent with those reported by Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000). Specifically, we observe the critical gradient Richardson number which induces statistical stationarity in these flows to be approximately 0.16. Similarly, we observe $ {{Fr}}\approx 0.5$ in all simulations, which implies $ {{S_\ast }} \approx 5$, as observed in other studies (Jacobitz et al. Reference Jacobitz, Sarkar and Van Atta1997; Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000).

Curiously, the lowest Reynolds number to maintain stationarity robustly according to our stringent criteria corresponds approximately to case R1. Case R1 corresponds to $ {{Re_b}} \approx 36$, which is consistent with critical values observed and estimated for sustained three-dimensional turbulence previously (Gibson Reference Gibson1980; Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005; Portwood et al. Reference Portwood, de Bruyn Kops, Taylor, Salehipour and Caulfield2016).

For further comparison with the flows considered by Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000) and Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005), we plot their relevant cases against the solutions presented for this research in a $Fr-Re_b$ parameter space, as shown in figure 3. We remark that the configurations reported by Shih and co-authors in those publications feature transitional and, crucially, variable $Fr$ when $Re_b$ is $O(100)$. As is apparent in the figure, $Fr$ and $Re_b$ appear to be coupled, with an approximate scaling relationship $Fr \propto Re_b^{1/2}$, as shown with a dashed line. Therefore, the parameterisation of mixing by $Re_b$ cannot be made independent of $Fr$ with the data presented in Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000,  Reference Shih, Koseff, Ivey and Ferziger2005). We note that this explanation does not rely on any argument that the simulations are inadequately resolved, either at small or large scales. Our perspective is qualitatively different from the perspective presented in Kunze et al. (Reference Kunze, MacKay, McPhee-Shaw, Morrice, Girton and Terker2012) and Gregg et al. (Reference Gregg, D'Asaro, Riley and Kunze2018), which suggests the $Re_b$-sensitive regime with $Re_b>100$ is due to inadequate small-scale spatial resolution in the simulations. However, it is possible that the solutions analysed in Shih et al. suffer from inadequate large-scale resolution in the domain as suggested by other reports (see Kunze Reference Kunze2011; Kunze et al. Reference Kunze, MacKay, McPhee-Shaw, Morrice, Girton and Terker2012; Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018).

Figure 3. Flows in a $Fr{-}Re_b$ parameter space. Points indication solutions obtained in the present work. For reference, the curves indicate trajectories in parameter space of solutions for cases fz, fc, fb and fd presented in Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000) and Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005). Darker segments of those lines indicate the reporting period $St>2$ used in that research.

Therefore, disentangling the behaviour of their simulations as $Re_b$ varies from the effects of variation in $Fr$ is challenging, particularly when $Re_b$ (or $Re_s$) is large. It is important to remember that it was precisely for $Re_b \gtrsim 100$ that Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) argued that the turbulent flux coefficient $\varGamma \propto Re_b^{-1/2}$, whereas it is entirely plausible that variable $Fr$ is playing a significant role. This figure highlights the novelty of our simulations in this particular flow configuration, which allow the isolated study of the effects of variations in Reynolds number from the influence of other non-dimensional turbulent parameters, and in particular at constant $Fr$.

4.1. Energetics

Just as non-dimensional parameters are emergent quantities in our simulations, so too are the partitionings of potential energy and the various components of kinetic energy. Such partitionings are significant, not least because the ratio of potential energy to kinetic energy,

(4.1)\begin{equation} R_{PK} \equiv E_p / E_k , \end{equation}

is a critical component of mixing models wherein the Reynolds number, or $ {{Re_b}}$, dependence is often omitted and the mixing is assumed to be a function of $ {{Ri}}$ (Osborn & Cox Reference Osborn and Cox1972; Schumann & Gerz Reference Schumann and Gerz1995; Pouquet et al. Reference Pouquet, Rosenberg, Marino and Herbert2018). Furthermore, in ‘strongly’ stratified turbulent flow, Billant & Chomaz (Reference Billant and Chomaz2001) suggests that there should be approximate equipartition between potential and kinetic energy, i.e. $R_{PK}\approx 1$, an assumption also used by Lindborg (Reference Lindborg2006). It is important to remember that these models do not account for Reynolds number effects. Therefore, parameterisation of such energetic partitioning is an important a priori validation necessary for model evaluation.

Recalling that the kinetic energy is the same for all cases, by construction from our definition of stationarity, we observe a significant decrease in potential energy with increasing dynamic range, as shown in figure 4(b). In comparison, the results of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), which span up to $ {{Re_b}} \sim O(10)$ at much lower Froude number in homogeneous stratified turbulence, report $R_{PK} \approx 0.15$ at $ {{Re_b}} = {{Ri}} {{Re_s}} \approx 16$, noting an asymptotic trend to this value from $R_{PK}\approx 0.05$ at $ {{Re_b}}\approx 0.1$. Although we observe similar values at $ {{Re_s}}\approx 5$ (corresponding to $ {{Re_b}} \approx 30$), higher values of $ {{Re_b}}$ reveal a subsequent decline in $R_{PK}$, as suggested by Remmler & Hickel (Reference Remmler and Hickel2012). For Reynolds numbers larger than $ {{Re_s}} \approx 40$ ($ {{Re_b}} \approx 300$), $R_{PK}$ assumes an apparently asymptotic value of $0.08$ in these simulations. This relative decrease of potential energy as a function of $ {{Re_s}}$ can be observed in the domain slice visualisations of $\rho$ featured in figure 1, as there is a clear reduction in the variance of the density fluctuations. In terms of the important question concerning the modelling of mixing, and in particular parameterising the turbulent flux coefficient $\varGamma$, (2.21) shows that, if $R_{PK}$ approaches an asymptotic value for large $Re_s$, $\varGamma$ can exhibit parameter dependence only if the key parameter $\varPi$ (defined in (2.17)) does.

Figure 4. (a) The partitioning of kinetic energy into each velocity component where a value of $1/3$ would be expected for a statistically isotropic velocity vector. (b) The ratio of potential energy to kinetic energy, approaching asymptotic values near 0.08 for $ {{Re_s}} \gtrapprox 50$.

Furthermore, a distinguishing characteristic of sheared stratified turbulence, relative to other stratified flows, is the absence of axisymmetry normal to gravity with the result that anisotropy must be parameterised in each dimension. Inertially relevant anisotropy can be characterised by the variance of individual velocity vector components as they contribute to the mean turbulent kinetic energy. These velocity variances resulting from an anisotropic dynamics are shown in figure 4(a). We observe turbulent kinetic energy is dominated by streamwise velocity fluctuations, with the smallest contributions coming from vertical velocity fluctuations. Notably distinct from axisymmetric flows, the horizontal components of velocity variance represent approximately 80 per cent of the total contributions to kinetic energy. The cross-stream velocity variance accounts for only 30 per cent of the total energy. The vertical variance, subject to exchanges to potential energy, is typically suggested to vary as a function of Froude number (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). Here, the vertical variance accounts for approximately 20 per cent of the total energy and we observe that it remains largely constant as a function of $ {{Re_b}}$.

Perhaps surprisingly, anisotropy of velocity fluctuations increases with increasing dynamic range. That is, the streamwise velocity component becomes increasingly dominant with increasing Reynolds number, apparently at the expense of the cross-stream velocity variance. One possible interpretation is that reducing viscous effects at inertial scales as the dynamic range increases allows the flow to evolve towards an inertial equilibrium state.

4.2. Dynamics

It might be assumed that the significant energetic transitions at sufficiently high Reynolds number would be explained by accompanying clear transitions in the dynamics. Here, we analyse the stationary behaviour of the dynamics relevant to the kinetic and potential energy balances. The ensemble-averaged energy equations, as derived from (2.1), obey

(4.2)\begin{gather} \left\langle \frac{\partial {E_p}}{\partial t} \right\rangle = B - \chi \quad \text{and} \end{gather}

(4.3)\begin{gather}\left\langle \frac{\partial {E_k}}{\partial t} \right\rangle = P - B - \epsilon , \end{gather}

where the turbulent production from mean shear $P=-\langle u_x u_z \, {{\rm d}}\bar {u}_x/{{{\rm d}}z}\rangle$ and the ‘buoyancy’ flux $B=\langle ({g}/{\rho _0}) {u_z} \rho \rangle$ have anisotropic effects on the evolution of the streamwise and vertical components of the kinetic energy $\langle u_x^{2}\rangle$ and $\langle u_z^{2}\rangle$, respectively.

Note that the left-hand sides of (4.3), (4.2) reduce to approximately zero due to the stationarity constraint. Therefore, the dynamics should, at least in principle, be describable by the turbulent flux coefficient $\varGamma$ defined in terms of $\chi$ in (2.21), and a flux Richardson number, defined as

(4.4)\begin{equation} {{Rf}} \equiv B/P \simeq \frac{\chi}{\chi + \epsilon} = \frac{\varGamma}{1+\varGamma}, \end{equation}

since, if the flow actually enters an emergent stationary state the flux $B \simeq \chi$, and remembering the underlying original definition (1.3) $\varGamma _b \simeq \varGamma$. Furthermore, in this circumstance, the turbulent viscosity $\nu _b$ can be naturally defined as

(4.5)\begin{equation} \nu_b \equiv \frac{P}{S^{2}} \rightarrow {{Pr_t}} \equiv \frac{\nu_b}{\kappa_b} \simeq \frac{{{Ri}}}{{{Rf}}} \simeq \left(\frac{1+\varGamma}{\varGamma} \right) {{Ri}}, \end{equation}

and so we can obtain an expression consistent with (2.22). Henceforth, we continue to describe the evolution of the energetics by $\varGamma$ defined as in (2.21) in recognition of the long-standing application of different definitions of the turbulent flux coefficient to estimating turbulent diffusivity and turbulent viscosity (Osborn & Cox Reference Osborn and Cox1972; Crawford Reference Crawford1982). Gregg et al. (Reference Gregg, D'Asaro, Riley and Kunze2018) discuss in detail the subtleties and potential for confusion concerning different definitions for turbulent flux coefficients, and so it must always be remembered that here we concentrate on using the definition in (2.21), i.e. $\varGamma \equiv \chi /\epsilon$.

The ratio of terms in the energetic balance is shown in figure 5. We observe a slight decline of $\varGamma$ from $0.2$ to approximately $0.175$, as noted in Portwood et al. (Reference Portwood, de Bruyn Kops and Caulfield2019), until $ {{Re_s}} \approx 40$. However, $\varGamma$ remains approximately constant as $Re_s$ increases further, and in particular we do not observe the proposed $ {{Re_b}}^{-1/2}$ scaling reported in a similar flow by Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) for their energetic regime of $ {{Re_b}} > 100$. Remembering the data presented in figure 3, this is perhaps unsurprising due to the apparent correlated variation of $Fr$ and $Re_b$ in those simulations. Due to the condition of stationarity being imposed on the kinetic energy in our simulations, accompanying induced stationarity of the potential energy is also guaranteed, although still plotted for reference in figure 5, demonstrating that $\chi \approx B$. It is also apparent that $ {{Rf}} \equiv B/P$ is less than $\varGamma \equiv \chi /\epsilon$, as expected.

Figure 5. Relationships between various components of the potential and turbulent kinetic energy evolution equations as a function of Reynolds number. The dashed line at unity corresponds to the unimposed, but emergent, condition that $\chi \approx B$. The dashed line at 0.2 indicates the upper bound for $\varGamma$ postulated by Osborn (Reference Osborn1980). Note that all dynamic ratios appear approximately constant as a function of Reynolds number.

4.3. Length scales

Since we have observed $ {{Fr}}$ and $ {{Ri}}$ to be essentially independent of the Reynolds number, the relatively decreasing potential energy as $ {{Re_s}}$ increases for $Re_s \lessapprox 50$ implies that scales characteristic of the scalar should also exhibit a similar asymptotic trend. The scalar outer scale associated with Obukhov–Corrsin similarity, $L_{OC}$ is defined in (2.12). In figure 6(a), we plot the ratio of this outer scale to the Ozmidov length scale $L_O$, as defined in (2.4). As $Re_s$ increases, this ratio decreases for $Re_s \lessapprox 40.$ For larger $Re_s$, there is some evidence that the ratio remains constant near unity, suggesting of course that $L_{OC}$ anomalously scales with $L_O$. This implies that the scaling $\chi \approx E_p N$ is valid in this high Reynolds number regime and that the outer scales of the scalar adjust to the regime associated with local isotropy, which is indeed a condition of a Obukhov–Corrsin subrange in anisotropic flows.

Figure 6. (a) The variation with $Re_s$ of the ratios of: the Obukhov–Corrsin scalar outer scale $L_{OC}$, defined in (2.12), to the Ozmidov length $L_O$, defined in (2.4); the scalar mixing length $L_\rho$, defined in (4.6a,b), to the Corrsin length scale $L_C$, defined in (2.6); and the scalar mixing length $L_\rho$ to the momentum mixing length $L_m$, also defined in (4.6a,b). (b) The variation with $Re_s$ of: the scaling of total dynamic range, using ${\rm \Delta} \ell _t$ and ${\rm \Delta} \ell _\rho$, as defined in (4.8) and (4.9), respectively; and (ten times) their ratio. The expected $ {{Re_s}}^{4/3}$ scaling is plotted with a dashed line, making apparent the anomalous scaling in the range of potentially isotropic scales associated with the scalar.

Length scales associated with the mixing length models suggested by Odier et al. (Reference Odier, Chen, Rivera and Ecke2009) can also be shown to remain constant in our large $Re_s$ regime, and indeed even for smaller $Re_s$. The scalar mixing length, $L_\rho$, and the momentum mixing length $L_m$, both defined by Odier et al. (Reference Odier, Chen, Rivera and Ecke2009) as

(4.6a,b)\begin{equation} L_m \equiv \left( \frac{P}{S^{3}} \right)^{1/2}, \quad L_\rho \equiv \left( \frac{B}{N^{2} S} \right)^{1/2},\end{equation}

are coupled via the turbulent Prandtl number $ {{Pr_t}}$, which here is very close to one, since

(4.7)\begin{equation} L_\rho = L_m \left( \frac{B S^{2}}{P N^{2}} \right)^{1/2} = L_m \left( \frac{{{Rf}}}{Ri} \right)^{1/2} \simeq {{Pr_t}}^{{-}1/2} L_m .\end{equation}

Furthermore, they coincide with the Corrsin length scale by the stationarity constraint from (4.3). The expected scalings $L_\rho \approx L_C$ and $L_\rho \approx L_m$ are both clearly apparent in figure 6(a) across a wide range of Reynolds numbers.

The decrease of $L_{OC}$ implies that the total dynamic range associated with the scalar anomalously scales with the shear Reynolds number $ {{Re_s}}$. We define the total range of dynamic scales for turbulence and the scalar as

(4.8)\begin{gather} {\rm \Delta} \ell_t= (L_{L} - L_K) \quad \text{and} \end{gather}

(4.9)\begin{gather}{\rm \Delta} \ell_\rho = (L_{OC} - L_B) , \end{gather}

respectively, and show them as a function of Reynolds number in figure 6(b). With respect to the turbulence scales, after an initial small $ {{Re_s}}$ transient, the range of total dynamic scales appears to increase with $ {{Re_s}}^{4/3}$, as predicted by our definitions in § 2.1.1. For the scalar, on the other hand, there is clearly a flatter-than-predicted scaling regime at high Reynolds number which is characteristic of the transient outer scales observed in figure 6(a). Furthermore, the range of dynamic scales associated with the scalar is significantly smaller than that of the turbulence, suggesting that, if a dynamic-range threshold exists, it will not occur simultaneously for the scalar and the turbulence.

5.1. Model verification

In the previous section, we have presented evidence that the ratio of the potential energy spectrum to a longitudinal energy spectrum appears to be consistent with the underlying scaling arguments central to coexistent Kolmogorov and Obukhov–Corrsin regimes. The value of the key parameter $\varPi$ can be expressed in terms of the one-dimensional universal constants, at high Reynolds number, by

(5.8)\begin{equation} \varPi = \frac{\beta}{C} = \frac{6}{11} \frac{\beta_1}{C_1} .\end{equation}

For these flows, the parameter $\varPi$ should be assumed to be universal to the extent that $C$ and $\beta$ are. For reasons outlined in Sreenivasan (Reference Sreenivasan1991) and also observed here, measuring individual constants from a sheared flow is fundamentally problematic and, furthermore, measuring $\beta$ requires an accurate measurement of $\chi$, which is inherently difficult in experimental flows.

Nonetheless, estimates of the (three-dimensional) Kolmogorov and Obukhov–Corrsin constants in the literature would imply that $\varPi \approx 0.42$ (Wyngaard & Coté Reference Wyngaard and Coté1971; Sreenivasan Reference Sreenivasan1991; Sreenivasan & Kailasnath Reference Sreenivasan and Kailasnath1996; Yeung Reference Yeung2002; de Bruyn Kops Reference de Bruyn Kops2015) with estimates of the individual constants typically being reported within approximately 15 % of each other in the literature cited. Robust estimates of $\varPi$ from the relation (5.8) and empirical observations in figure 7 indicate $\varPi \approx 0.39$, which is certainly within the range of values reported in the literature. We demonstrate that the important expression (5.8) appears to be valid within the range of uncertainty for reported values of $\beta /C$ in figure 9. We also remark that the range of $\varPi$ is small with respect to values of $\beta /C$ for the Reynolds number space considered. Such an apparently limited range of $\varPi$ has also been reported by Venayagamoorthy & Stretch (Reference Venayagamoorthy and Stretch2006).

Figure 9. Verifications of relations (5.8) and (5.9) for sufficiently large $ {{Re_s}}$. Error bars correspond to standard deviations associated with estimates of $C$ and $\beta$ from Wyngaard & Coté (Reference Wyngaard and Coté1971). Recall that $\varPi /Fr=N_\ast$ from (2.19).

Furthermore, given the empirical observation that $L_{OC}=L_O$ (figure 6a) at sufficiently high Reynolds number,

(5.9)\begin{equation} \varPi \sim {{Fr}} .\end{equation}

Therefore, the stationary Froude number can be thought of as a simple consequence of the two universal constants at high Reynolds numbers, as validated in figure 9, which further implies $N_\ast \approx 1$, remembering (2.19). This finally suggests that the behaviour at lower Reynolds number may be parameterised by a function $f_\varPi$, i.e. at unity Prandtl number

(5.10)\begin{equation} \varPi = \frac{N_\ast}{{{S_\ast}} \sqrt{{{Ri}}}} \approx {{Fr}} f_\varPi({{Re_s}}). \end{equation}