Implications of inertial subrange scaling for stably stratified mixing

Abstract

We investigate the effects of the turbulent dynamic range on active scalar mixing in stably stratified turbulence by adapting the theoretical passive scalar modelling arguments of Beguier, Dekeyser & Launder (1978) (Phys. Fluids, vol. 21 (3), pp. 307–310) and demonstrating their usefulness through consideration of the results of direct numerical simulations of statistically stationary homogeneous stratified and sheared turbulence. By analysis of inertial and inertial–convective subrange scalings, we show that the relationship between the active scalar and turbulence time scales is predicted by the ratio of the Kolmogorov and Obukhov–Corrsin constants, provided mean flow parameters permit the two subrange scalings to be appropriate approximations. We use the resulting relationship between time scales to parameterise an appropriate turbulent mixing coefficient $\varGamma \equiv \chi /\epsilon$, defined here as the ratio of available potential energy ($E_p$) and turbulent kinetic energy ($E_k$) dissipation rates. With the analysis presented here, we show that $\varGamma$ can be estimated by $E_p,E_k$ and a universal constant provided an appropriate Reynolds number is sufficiently high. This large Reynolds number regime appears here to occur at $ {{Re_b}} \equiv \epsilon / \nu N^{2} \gtrapprox 300$ where $\nu$ is the kinematic viscosity and $N$ is the characteristic buoyancy frequency. We propose a model framework for irreversible diapycnal mixing with robust theoretical parametrisation and asymptotic behaviour in this high-$ {{Re_b}}$ limit.

1. Introduction

Stably stratified turbulence may be thought of as a model flow that is potentially useful for understanding aspects of geophysical and engineering processes in many regions of the oceans and atmosphere. In particular, stratified turbulence describes the relatively small-scale dynamics, in length and time, at which turbulence and irreversible mixing occurs. Interpreting and modelling such idealised flows is a primary avenue for the development and calibration of global circulation and basin-scale models, for which mixing plays a leading-order role in global energy budgets (e.g. Ferrari & Wunsch 2009; Jayne 2009; Gregg et al. 2018). However, even in idealised stably stratified flow configurations, parameterised modelling of mixing has emerged as a challenge due to the potential for dependence on a wide range of non-independent parameters (Gregg et al. 2018; Ivey, Bluteau & Jones 2018; Caulfield 2021).

Parallel to the growing recognition of the importance of modelling such small-scale processes in stratified turbulence for broader geophysical applications, stratified turbulence has been increasingly observed to exhibit certain quantitatively similar small-scale dynamics to isotropic turbulence in some parameter regimes, specifically when the buoyancy Reynolds number, $ {{Re_b}}$, is large enough (Gargett et al. 1981; Lindborg 2006; de Bruyn Kops 2015). Of course, this does not mean that stratified turbulence is equivalent in all respects, not least because stratification inevitably introduces anisotropy. We consider $ {{Re_b}}$ in more detail in the next section, but in simple terms it quantifies the dynamic range of turbulent length scales neither directly affected by large-scale stratification nor by molecular viscosity. Indeed, it is widely acknowledged that the details of the large-scale, or outer, mean-scale, dynamics do not strongly affect the small-scale dynamics provided there is sufficient scale separation between such large and small scales. While the small-scale dynamics is certainly not independent of large scales (Corrsin 1958; Durbin & Speziale 1991), the assumption that sufficient scale separation induces statistically independent small scales has proven to be a valuable tool in the modelling of anisotropic turbulent flows because it allows for the application of theoretical models based on statistical symmetries, e.g. isotropy and homogeneity, to dynamic modelling.

To take one example, Kolmogorov–Obukhov scaling depends on the assumptions of local isotropy and homogeneity (Kolmogorov 1941; Oboukhov 1941) yet it has been observed to be a good approximation for anisotropic flows provided the scale separation is sufficiently large between anisotropic turbulence scales and the dissipative viscous–diffusive scales (Champagne, Harris & Corrsin 1970; Gargett, Osborn & Nasmyth 1984; Saddoughi & Veeravalli 1994; Shen & Warhaft 2002). However, not all stably stratified flows exhibit sufficient scale separation for the foregoing assumption to be made (Jackson & Rehmann 2014; de Bruyn Kops 2015).

This phenomenology is useful because it motivates the adaptation of existing models for isotropic turbulence to stably stratified turbulence. For instance, by analysis of Obukhov–Corrsin and Kolmogorov scalings, Beguier, Dekeyser & Launder (1978) observed that the passive scalar time scale $\tau _\phi$ and turbulence time scale $\tau$ couple at sufficiently high Reynolds number such that they can be related by universal constants, i.e.

(1.1)\begin{equation} \frac{\tau_\phi}{\tau} \equiv \frac{E_\phi \epsilon}{\epsilon_\phi E_k} = \frac{\beta}{C}, \end{equation}

where $E_\phi$ is the turbulent scalar variance, $E_k$ is the turbulent kinetic energy, $\epsilon _\phi$ and $\epsilon$ are their respective dissipation rates and $\beta,C$ are the Obukhov–Corrsin and Kolmogorov constants, respectively. The assumption that the time-scale ratio is a constant is commonly applied to model scalar dissipation in Reynolds-averaged Navier–Stokes models (e.g. Newman, Launder & Lumley 1981; Ristorcelli 2006).

In modelling mixing in stratified flows, the applicability of such a universal relation is interesting due to its implications for modelling the irreversible diapycnal flux, $j_b$. In such models, equilibrium assumptions are typically prescribed such that the irreversible diapycnal flux is equal to a particular definition of the available potential energy dissipation rate (Osborn & Cox 1972; Peltier & Caulfield 2003; Caulfield 2021) such that coupling of the turbulent and scalar dynamics, as suggested by Beguier et al. (1978), can be an insightful relation for modelling the scalar dynamics.

The importance of accurately modelling such irreversible scalar fluxes in the ocean has motivated significant research on various idealised flows using an array of modelling frameworks. Perhaps most notable, estimating the irreversible flux from the scalar gradient via a turbulent diffusivity, $\kappa _b$, is the de facto standard approach (e.g. Ivey & Imberger 1991; Barry et al. 2001; Salehipour & Peltier 2015; Maffioli & Davidson 2016), with

(1.2)\begin{equation} j_b = \kappa_b N^{2}, \end{equation}

where $N$ is an appropriately defined large-scale buoyancy frequency. A common parametrisation of $\kappa _b$ has been suggested by Osborn (1980), which asserts that $\kappa _b$ be related to the irreversible rate of dissipation of kinetic energy via a turbulent flux coefficient $\varGamma _b$ such that

(1.3)\begin{equation} \kappa_b = \varGamma_b \epsilon/N^{2} \rightarrow \varGamma_b \equiv \frac{j_b}{\epsilon} .\end{equation}

By analysis of the energy equations for a simple stratified flow model at stationary conditions, Osborn (1980) suggested that $\varGamma _b$ is a constant $\leq 0.2$. More recently, the prescription of the parameter dependence for $\varGamma _b$ has emerged as a necessary improvement to turbulent diffusivity models due to increasing evidence that a constant coefficient for $\varGamma _b$ appears not to be at all suitable. Parametrisations of $\varGamma _b$ in this framework have been determined to be difficult due to a potential dependence on multiple parameters, which themselves may well be correlated (Caulfield 2021), which difficulties lead to often contradictory subclosures (Monismith, Koseff & White 2018). While alternative models exist, such as mixing length models (Odier et al. 2009; Ivey et al. 2018) and flux transport models, they have not as yet gained widespread usage in large-scale simulations.

Therefore, as the first principal objective of this research, we consider formally adapting the Beguier et al. (1978) relationship (1.1) to stably stratified turbulence in order to model irreversible diapycnal mixing dynamically. To evaluate this hypothesis, we consider simulations of stationary homogeneous stratified and sheared turbulence (S-HSST). S-HSST is an idealised model flow configuration that enables easy adjustment of the range of length scales (the dynamic range) available for a locally isotropic subrange to form without changing other dimensionless flow parameters such as the characteristic Froude or Richardson numbers, as defined precisely in the next section. Such a fundamental flow configuration is well suited to the evaluation of the effects of locally isotropic scaling theories (Shih et al. 2000; Chung & Matheou 2012). To test the hypothesis that the Beguier et al. (1978) relation can be applied to stably stratified turbulence, we perform numerical experiments of S-HSST in parameter space extremes previously inaccessible to computation. Therefore, a second principal objective of this work is to report on the characteristics of S-HSST at very large Reynolds numbers and time scales, quantitatively described in subsequent sections.

To achieve these objectives, the rest of the paper is organised as follows. In § 2, we discuss a length-scale-based framework for stratified and sheared turbulence and use it to derive a mixing model hypothesis based on the arguments of Beguier et al. (1978). In § 3, we then describe numerical simulation experiments of homogeneous stratified and sheared turbulence (HSST). In a Reynolds number parameter space, we then present the ensemble-averaged dynamics in § 4, evaluate the mixing model in § 5 and finally draw our conclusions in § 6.

2. Theoretical background

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 1958) 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 1994).

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 (2014) 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 2014). 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 1997), 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. (2017) 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. 2000; Holt, Koseff & Ferziger 1992; Portwood, de Bruyn Kops & Caulfield 2019).

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. 1997; Shih et al. 2000)

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 2021). 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. (2019) observed this parameter to approach unity in the high $Re_b$ limit for flows where $Pr=1$, and Mater & Venayagamoorthy (2014) 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 1975). 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 1941).

This is a prerequisite of classical inertial subrange scaling arguments (Kolmogorov 1941; Oboukhov 1941; Corrsin 1951). 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 1957; Corrsin 1958). Kolmogorov scaling (Kolmogorov 1941,  1962) in the presence of anisotropic outer scales has been revealed to coincide at scales consistent with local isotropy (Champagne et al. 1970; Saddoughi & Veeravalli 1994). The equivalent local isotropy justification of the Kolmogorov subrange may be applied to Obukhov–Corrsin passive scalar scaling (Oboukhov 1949; Corrsin 1951) such that it is fit to describe an active scalar below scales affected by stratification (Monin & Yaglom 1975, 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 (1941), 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. (1978), 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 (1980) 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 (1982) with the turbulent diffusivity model of Osborn (1980), 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 (2010) 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 2006). 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. (2009), 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. Direct numerical simulations

In order to test the underlying hypotheses described in § 2.2 leading to the irreversible mixing model embodied in (2.18), we must consider a stratified model flow wherein the behaviour of $\varPi$ may be evaluated as a function of the dynamic range associated with isotropy. The consideration of S-HSST is justified by the emergence of a Reynolds number as the only independent non-dimensional descriptive flow parameter due to two distinct phenomenologies: (i) the coupling of the outer turbulence scales to that of the shear production mechanism such that $ {{S_\ast }} \approx 5$ and (ii) the gradient Richardson number assuming a balanced apparently critical value that maintains a stationary balance of energy production and dissipation mechanisms. Both phenomenologies were discovered in direct numerical simulations for this flow configuration by the investigations of Shih et al. (2000,  2005) and subsequently verified more recently in a broader Reynolds number space in simulations where the Richardson number is controlled to maintain constant turbulent kinetic energy (Portwood et al. 2019). We therefore utilise the flows presented in Portwood et al. (2019) to verify the modelling framework presented in § 2.2 and relevant associated consequences highlighted later in this manuscript.

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 (2016). The rest of the solver is derived from the method used in Hebert & de Bruyn Kops (2006), 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 2012; Jacobitz & Moreau 2016).

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 1998; Rao & de Bruyn Kops 2011), 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 (2012) 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 2011), and then by applying a first-order approximation to the temporal evolution of kinetic energy (cf. Jacobitz et al. 1997)

(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. (1997), the characteristic frequency $\omega$ is determined by the large-eddy time scale which itself is proportional to $1/S$ (Shih et al. 2000) 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. (2019), 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 (2012) and Jacobitz & Moreau (2016).

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. (1997) and Shih et al. (2000). 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 1998; Rao & de Bruyn Kops 2011) 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. (1992), Shih et al. (2000), although the emergent critical values of the Richardson numbers are broadly consistent with those reported by Shih et al. (2000). 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. 1997; Shih et al. 2000).

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 1980; Shih et al. 2005; Portwood et al. 2016).

For further comparison with the flows considered by Shih et al. (2000) and Shih et al. (2005), 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. (2000,  2005). 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. (2012) and Gregg et al. (2018), 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 2011; Kunze et al. 2012; Gregg et al. 2018).

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. (2000) and Shih et al. (2005). 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. (2005) 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. Emergent phenomena

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 1972; Schumann & Gerz 1995; Pouquet et al. 2018). Furthermore, in ‘strongly’ stratified turbulent flow, Billant & Chomaz (2001) suggests that there should be approximate equipartition between potential and kinetic energy, i.e. $R_{PK}\approx 1$, an assumption also used by Lindborg (2006). 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. (2007), 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 (2012). 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. 2007). 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 1972; Crawford 1982). Gregg et al. (2018) 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. (2019), 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. (2005) 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 (1980). 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. (2009) 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. (2009) 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. Implications of inertial subrange scaling

Before turning to model verification, an intermediate phenomenological test of the universality hypothesis outlined in § 3.1 may be performed by analysis of the one-dimensional correlations. The one-dimensional spectrum of the potential energy is expected to be determined by the energetic dissipation rates and a wavenumber, i.e. for the cross-stream spectra

(5.1)\begin{equation} E_{p}^{OC}(k_y) = \beta_1 \chi^{1} \epsilon^{{-}1/3} k_y^{{-}5/3} f_L^{OC}(k_y L^{O}_{yy})f_v^{OC}(k_y L_B) . \end{equation}

In this expression, $\beta _1$ is the one-dimensional Obukhov–Corrsin constant. Furthermore, the functions denoted by $f^{OC}_v$, $f^{OC}_L$ are the viscous and outer-scale correction functions which are necessary if no implicit assumptions of local isotropy are made. The notation $L^{O}_{ij}$ indicates the outer scale of locally isotropic turbulence associated with the velocity component $i$ with respect to the direction $j$. This is ostensibly the Corrsin length scale but is inevitably anisotropic, as discussed in Kaimal (1973) for the specific case of the stably stratified boundary layer.

Equivalently, according to Kolmogorov (1941) for an isotropic inertial subregime, the one-dimensional cross-stream longitudinal and streamwise transverse energy spectra are

(5.2)\begin{gather} E_{y}^{K}(k_y) = C_1 \epsilon^{2/3} k_y^{{-}5/3} f_L^{K}(k_y L^{O}_{yy})f_v^{K}(k_y L_K), \end{gather}

(5.3)\begin{gather}E_{x}^{K}(k_y) = C_1' \epsilon^{2/3} k_y^{{-}5/3} f_L^{K}(k_y L^{O}_{xy})f_v^{K}(k_y L_K), \end{gather}

where $C_1$ and $C_1'$ are universal constants, and $f^{K}_L$ and $f^{K}_v$ are the correction functions. The ratio of (5.1) and (5.2) yields

(5.4)\begin{equation} \frac{E_{\rho}^{OC}(k_y)}{E_{y}^{K}(k_y)} = \frac{\beta_1}{C_1} \frac{\chi}{\epsilon} \frac{f_L^{OC}}{f_L^{K}} \frac{f_v^{OC}}{f_v^{K}} .\end{equation}

From empirical studies of turbulence, even in the absence of the effects of stratification and shear, we expect the correction functions $f_L$, $f_v$ to be non-trivial (e.g. Muschinski & de Bruyn Kops 2015). However, it is at least plausible that a subregime exists wherein

(5.5a,b)\begin{equation} \frac{f_L^{OC}}{f_L^{K}} \approx 1 \quad \text{and}\quad \frac{f_v^{OC}}{f_v^{K}} \approx 1 ,\end{equation}

either due to sufficiently high scale separation of isotropic scales, i.e. $ {{Re_s}} \gg 1$, or due to shared functional forms such that $f_L^{OC}=f_L^{K}$. Subject to either condition, it would be expected that (5.4) would reduce to

(5.6)\begin{equation} \frac{E_{\rho}^{OC}(k_y)}{E_{y}^{K}(k_y)} = \frac{\beta_1}{C_1} \frac{\chi}{\epsilon} . \end{equation}

Indeed, similar arguments would motivate the $4/3$ relation between transverse and longitudinal velocity spectra in the presence of anisotropic outer scales such that $F_L^{K}(k_yL^{O}_{xy})\neq F_L^{K}(k_yL^{O}_{yy})$ in order to obtain

(5.7)\begin{equation} \frac{E_{x}^{K}(k_y)}{E_{y}^{K}(k_y)} = \frac{C_1'}{C_1}=\frac{4}{3} .\end{equation}

The validity of the relation (5.6) is tested with appropriately scaled one-dimensional spectra in figure 7(a), where we would expect the relation to hold in the locally isotropic wavenumber regime where $L_C k_y \sim O(1)$ (Saddoughi & Veeravalli 1994). We observe a tendency for the local minima of the spectra to decrease with increasing Reynolds number. The minima for cases R8, R9 and R10 are all close to, and bounded below by 0.72, within the range of values reported for the universal constants (Sreenivasan 1995; Sreenivasan & Kailasnath 1996). Furthermore, in the stratified boundary layer, Wyngaard & Coté (1971) report measurements of each constant independently, in a flow configuration very similar to the one studied here, which imply $0.76 \pm 0.11$, also in agreement with the results observed here.

Figure 7. (a) The compensated ratio of one-dimensional energy spectra in the cross-stream direction. At high wavenumbers, the spectra lie in order of Reynolds number with R1 on the left and R10 on the right, as labelled. The coexistence of Kolmogorov and Obukhov–Corrsin scalings, even when subject to non-trivial correction functions, suggests a subregime below anisotropic scales which features a plateau corresponding to the ratio $\beta _1/C_1$ of the Obukhov–Corrsin constant to the Kolmogorov constant as predicted by (5.6), which is expected to occur at $L_Ck_y\approx 1$. The dashed line at 0.76 indicates measurements made in the stratified boundary layer by Wyngaard & Coté (1971) where the surrounding shaded region represents its uncertainty from reported standard deviations and the dotted line indicates an estimate of the asymptotic value of $\beta _1/C_1$ at 0.72. (b) The ratio of one-dimensional streamwise–transverse and cross-stream–longitudinal spectra, which have a predicted ratio given in (5.7) at high $ {{Re_s}}$, expected to occur at wavenumbers a factor of $ {{Ri}}^{-3/4}\approx 4$ greater than the locally isotropic regime shown in panel (a).

The wavenumbers associated with this regime also tend toward smaller scales as $ {{Re_s}}$ increases. In cases R8, R9 and R10, the lower bound for this regime appears to remain fixed at approximately $k_yL_C=0.4$ as the right bound begins to increase such that the flat region of the spectra broadens. We note that this apparently asymptotic behaviour in cases R8, R9 and R10 occurs while $ {{Re_s}}$ increases by a factor of two. The ratio of transverse to longitudinal one-dimensional spectra, as expressed in (5.7), is shown in figure 7(b). We observe a similar downward trend of the locally isotropic wavenumber regime with $ {{Re_s}}$, leading eventually to good agreement with the $4/3$ law for case R10.

Similarly to the data shown in figure 7(a), the locally isotropic regime at higher wavenumbers becomes more consistent with the Corrsin length scale at approximately $k_yL_C\approx 2$. The misalignment of the wavenumber regimes consistent with (5.6) might be seen to be a consequence of the anisotropy of the outer scales, as accounted for in our definitions of anisotropic length scales in $f_L^{K}$ and also in the observation that the outer scale of the scalar largely coincides with the Ozmidov length. Therefore, it is expected that $L_{OC}/L_C \sim {{Ri}}^{-3/4}$ and hence the wavenumber regime associated with the four-thirds relation is expected to be approximately a factor of four larger than (5.6), essentially as observed here.

We compare measured three-dimensional potential and kinetic energy spectra with the model spectra, (2.16) and (2.15), in figure 8(a,b). In both panels, a $|\boldsymbol {k}|^{-5/3}$ curve is plotted for reference. Our intention is to verify that the proposed spectra are approximated by our analytic model spectra at high $Re_s$, in the sense that we do not observe substantial systematic deviations which would otherwise indicate alternative two-point parameterisations would be more appropriate. Here, we observe a trend wherein the range of scales associated with three-dimensional turbulence, $|\boldsymbol {k}|L_C>1$, is approximated by the $-5/3$ power law as $ {{Re_s}}$ increases. Whereas deviations from the power-law behaviour are observed at all $ {{Re_s}}$, the quantitative sensitivity of the parameter $\varPi$ to deviations from the model spectra is unclear at this stage. Therefore, we next turn to evaluating the impact of Reynolds number on our key parameter $\varPi$.

Figure 8. Three-dimensional kinetic and potential energy spectra shown in panels (a) and (b), respectively. The dashed black line indications a $-5/3$ slope, as assumed by the model spectra (2.16) and (2.15).

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 (1991) 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é 1971; Sreenivasan 1991; Sreenivasan & Kailasnath 1996; Yeung 2002; de Bruyn Kops 2015) 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 (2006).

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é (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}

6. Concluding remarks

Inspired by increasing evidence that high dynamic-range stratified turbulence exhibits scale similarity as proposed by Kolmogorov (1941), Oboukhov (1949) and Corrsin (1951), we have considered the formal adaptation of a passive scalar mixing model (Beguier et al. 1978) to stratified scalar mixing. By deriving the model from scaling theories, we demonstrate how model parameters have explicit functional dependence and are Reynolds number independent, at sufficiently high Reynolds number. In order to verify our hypothesis, we have employed a model flow (i.e. S-HSST) which equilibrates such that the effects of decreasing dissipation scales, here parameterised by increasing shear Reynolds number $ {{Re_s}}$ (as defined in (2.7a–c)), may be investigated independently of other flow parameters.

By direct numerical simulations of S-HSST, we have shown that the outer length scales of turbulence and the active scalar become parametrically coupled such that the relationship between energetics and dissipation rates is predicted by the scaling theories of Kolmogorov and Obukhov–Corrsin similarity. Therefore, the fundamental parameter associated with this relationship, $\varPi$, as defined in (2.17), is the ratio (5.8) of the Obukhov–Corrsin constant $\beta$ and the Kolmogorov constant $C$, at sufficiently high $ {{Re_s}}$, and should be considered universal to the extent that its constitutive constants are. As detailed in the context of various mixing models in § 2.2, robust characterisations of $\varPi$ have strong implications for the calibration of a number of turbulence and mixing models.

This universality has profound, but still not fully explained, implications for the turbulent flux coefficient $\varGamma$ of such steady stratified and sheared turbulence. From the expression (2.21), $\varGamma$ tending to a universal asymptotic value can now be understood as being due to $\varPi$ tending to a constant, which follows from the scaling theory and similarity arguments presented here, that in turn rely upon classical, and relatively well-established, turbulence modelling approaches. Therefore, to ‘understand’ physically and theoretically why the empirical modelling of Osborn & Cox (1972) and Osborn (1980) works so well (with $\varGamma \sim 0.2$), the remaining interesting open question to address is why the ratio of energies $R_{PK}=E_p/E_k \simeq 0.08$ for $Re_s \gtrapprox 50$, as shown in figure 4.

The parameter $\varPi$ may be stated exactly as a function of $N_\ast, {{S_\ast }} , {{Ri}}$. However, we demonstrate that it may be effectively reduced to a function of $ {{Fr}}$ at high Reynolds number, or a function of $ {{Fr}}$ and $ {{Re_s}}$ at low Reynolds numbers for the flows considered here. The transition of the high and low Reynolds number regimes is observed approximately at $ {{Re_s}}\approx 40$, or $ {{Re_b}}\approx 300$ at stationary Richardson number, for unity Prandtl number. We demonstrate that large-scale characteristics of the active scalar, as parameterised by $N_\ast$, exhibit significant dependence on the Reynolds number for $ {{Re_s}} \lessapprox 40$. The correlation between the transition point, in Reynolds number space, of the parameter $\varPi$ and $N_\ast$ supports the underlying assumptions of the proposed model.

A stated and significant limitation of the present study is the restriction of the model to $ {{Fr}} \sim O(1)$ and $Pr=1$. An objective of future work is to perform a similar procedure to § 2.2, except using scaling laws which account for expanded parameter regimes (for instance Batchelor 1959; Lindborg 2006; Kunze 2019).