3.1 Model formulation

In order to obtain the vertical profiles of mean velocity and temperature in fully developed turbulent stratified plane Couette flow, we integrate the following set of equations of $U$ and $\unicode[STIX]{x1D6E9}$ in time (using the laminar profiles as initial conditions) until reaching a steady state:

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}U}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\unicode[STIX]{x1D708}_{t}\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}\right), & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D705}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\unicode[STIX]{x1D705}_{t}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}y}\right), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{t}$ and $\unicode[STIX]{x1D705}_{t}$ are the turbulent (eddy) viscosity and diffusivity respectively.

The closure for $\unicode[STIX]{x1D708}_{t}$ in the Reynolds-averaged momentum equation (3.1) can be obtained by specifying a mixing length (see e.g. Pope Reference Pope2000):

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{t}=\ell _{m}^{\ast 2}\left|\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}\right|=\ell _{m}^{\ast }u^{\ast },\end{eqnarray}$$

where $\ell _{m}^{\ast }$ is the mixing length for momentum and the fluctuation velocity

(3.4) $$\begin{eqnarray}u^{\ast }=\ell _{m}^{\ast }\left|\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}\right|.\end{eqnarray}$$

Similarly, the turbulent flux of scalar in (3.2) can be modelled as

(3.5) $$\begin{eqnarray}-\langle v^{\prime }\unicode[STIX]{x1D703}^{\prime }\rangle =\unicode[STIX]{x1D705}_{t}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}y}=u^{\ast }\unicode[STIX]{x1D703}^{\ast }=\ell _{m}^{\ast }\left|\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}\right|\ell _{s}^{\ast }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}y},\end{eqnarray}$$

where $\ell _{s}^{\ast }$ is the scalar mixing length, and it follows that

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D705}_{t}=\ell _{s}^{\ast }\ell _{m}^{\ast }\left|\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}\right|=\ell _{s}^{\ast }u^{\ast }.\end{eqnarray}$$

It remains to specify the two mixing lengths $\ell _{m}^{\ast }$ and $\ell _{s}^{\ast }$ .

To do this, we start by considering unstratified flows, i.e. $L\rightarrow \infty$ . We define $y_{w}$ as the wall-normal (vertical) distance to the closer wall, i.e. $y_{w}\equiv \min (h-y,h+y)$ . The length $y_{w}$ can be normalized in wall units as $y^{+}\equiv y_{w}/(\unicode[STIX]{x1D708}/u_{\unicode[STIX]{x1D70F}})$ . The ‘law of the wall’ of unstratified wall-bounded flows (see e.g. Bradshaw & Huang Reference Bradshaw and Huang1995) prescribes the wall-normal gradients of $U$ and $\unicode[STIX]{x1D6E9}$ in the log-law region, i.e. $y^{+}>30$ (Pope Reference Pope2000), as

(3.7a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}=\frac{u_{\unicode[STIX]{x1D70F}}}{k_{m}y_{w}}\quad \text{and}\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}y}=\frac{\unicode[STIX]{x1D703}_{\unicode[STIX]{x1D70F}}}{k_{s}y_{w}}=\frac{\unicode[STIX]{x1D703}_{\unicode[STIX]{x1D70F}}\hat{Pr}_{t}}{k_{m}y_{w}},\end{eqnarray}$$

where $k_{m}$ and $k_{s}$ are the von Kármán constants for momentum and scalar respectively, and $\hat{Pr}_{t}=k_{m}/k_{s}$ is a turbulent Prandtl number which applies for the log-law region. With

(3.8a,b ) $$\begin{eqnarray}u_{\unicode[STIX]{x1D70F}}^{2}\cong \unicode[STIX]{x1D708}_{t}\left|\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}\right|\quad \text{and}\quad q_{w}=\unicode[STIX]{x1D703}_{\unicode[STIX]{x1D70F}}u_{\unicode[STIX]{x1D70F}}\cong \unicode[STIX]{x1D705}_{t}\left|\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}y}\right|,\end{eqnarray}$$

in the log-law region and following the model prescriptions in (3.3) and (3.6), the mixing lengths $\ell _{m}^{\ast }$ and $\ell _{s}^{\ast }$ corresponding to (3.7) read

(3.9a,b ) $$\begin{eqnarray}\ell _{m}^{\ast }=k_{m}y_{w}\quad \text{and}\quad \ell _{s}^{\ast }=k_{s}y_{w}=\ell _{m}^{\ast }\hat{Pr}_{t}^{-1}.\end{eqnarray}$$

As a result, the velocity scale $u^{\ast }$ in (3.3) and (3.6) can be specified as

(3.10) $$\begin{eqnarray}u^{\ast }=u_{\unicode[STIX]{x1D70F}}.\end{eqnarray}$$

When the fluid is stratified, Monin–Obukhov similarity theory prescribes the vertical gradients of $U$ and $\unicode[STIX]{x1D6E9}$ as

(3.11a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}y}=\frac{u_{\unicode[STIX]{x1D70F}}}{k_{m}y_{w}}\unicode[STIX]{x1D6F7}_{m}(\unicode[STIX]{x1D709})\quad \text{and}\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E9}}{\unicode[STIX]{x2202}y}=\frac{\unicode[STIX]{x1D703}_{\unicode[STIX]{x1D70F}}}{k_{s}y_{w}}\unicode[STIX]{x1D6F7}_{s}(\unicode[STIX]{x1D709}).\end{eqnarray}$$

In these expressions, $\unicode[STIX]{x1D6F7}_{m}$ and $\unicode[STIX]{x1D6F7}_{s}$ are Monin–Obukhov functions which are linear in the non-dimensional variable $\unicode[STIX]{x1D709}\equiv y_{w}/L$ for stable stratification:

(3.12a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{m}(\unicode[STIX]{x1D709})=1+\unicode[STIX]{x1D6FD}_{m}\unicode[STIX]{x1D709}\quad \text{and}\quad \unicode[STIX]{x1D6F7}_{s}(\unicode[STIX]{x1D709})=1+\unicode[STIX]{x1D6FD}_{s}\unicode[STIX]{x1D709}.\end{eqnarray}$$

Here we take $k_{m}=0.41$ and $k_{s}=0.48$ following Bradshaw & Huang (Reference Bradshaw and Huang1995). The choice of $\unicode[STIX]{x1D6FD}_{m}=4.8$ follows the recommendation of Wyngaard (Reference Wyngaard2010) and $\unicode[STIX]{x1D6FD}_{s}=5.6$ is used following the specific choice of $k_{m}$ and $k_{s}$ . These model constants are determined empirically using field observations of stable atmospheric boundary layers, and their values can exhibit some uncertainties (see Foken Reference Foken2006 for a review). The form of the similarity functions may also require additional corrections in order to match the field situations (see e.g. Tastula et al. Reference Tastula, Galperin, Sukoriansky, Luhar and Anderson2015), such as varying fluxes with height. In the idealized situation considered here, where the entire flow between the walls is a constant-flux layer by construction, we use the classical canonical forms of Monin–Obukhov functions described in (3.12) for clarity and simplicity.

The mixing length formulation corresponding to Monin–Obukhov theory becomes

(3.13a,b ) $$\begin{eqnarray}\ell _{m}^{\ast }=k_{m}y_{w}\unicode[STIX]{x1D6F7}_{m}^{-1}(\unicode[STIX]{x1D709})\quad \text{and}\quad \ell _{s}^{\ast }=k_{s}y_{w}\unicode[STIX]{x1D6F7}_{s}^{-1}(\unicode[STIX]{x1D709}).\end{eqnarray}$$

Taking $\unicode[STIX]{x1D709}\rightarrow 0$ in (3.13), one recovers the unstratified formulation (3.9). Since $\ell _{m}^{\ast }$ and $\ell _{s}^{\ast }$ are specified in (3.13) in very similar ways, the ratio $\ell _{m}^{\ast }/\ell _{s}^{\ast }$ is expected to be of order unity.

3.2 The near-wall layer

Here we focus on the viscous wall region, i.e. $y^{+}<50$ (Pope Reference Pope2000), where the (molecular) Prandtl number $Pr$ plays a critical role. The mean velocity and temperature differences relative to the closer wall can be written in wall units as

(3.14a,b ) $$\begin{eqnarray}U^{+}=\frac{\min (U+U_{w},U_{w}-U)}{u_{\unicode[STIX]{x1D70F}}}\quad \text{and}\quad \unicode[STIX]{x1D6E9}^{+}=\frac{\min (\unicode[STIX]{x1D6E9}+T_{w},T_{w}-\unicode[STIX]{x1D6E9})}{\unicode[STIX]{x1D703}_{\unicode[STIX]{x1D70F}}},\end{eqnarray}$$

where the velocity and temperature at the upper and lower walls are fixed at $\pm U_{w}$ and $\pm T_{w}$ , respectively. In the viscous/conductive sublayer near the wall (as shown in figure 2 for $\unicode[STIX]{x1D6E9}^{+}$ ),

(3.15a,b ) $$\begin{eqnarray}U^{+}=y^{+}\quad \text{and}\quad \unicode[STIX]{x1D6E9}^{+}=y^{+}Pr.\end{eqnarray}$$

As $y^{+}$ increases, the viscous/conductive sublayer transitions into the log-law region for which the mean profiles can be obtained by integrating (3.7) to yield

(3.16a,b ) $$\begin{eqnarray}U^{+}=\frac{1}{k_{m}}\ln y^{+}+C_{m}\quad \text{and}\quad \unicode[STIX]{x1D6E9}^{+}=\frac{1}{k_{s}}\ln y^{+}+C_{s}=\frac{\hat{Pr}_{t}}{k_{m}}\ln y^{+}+C_{s}.\end{eqnarray}$$

DNS of stratified plane Couette flows recover such behaviour in the near-wall region, as shown in figure 2. Unlike $C_{m}$ which is a constant (we take $C_{m}=5.0$ following Bradshaw & Huang (Reference Bradshaw and Huang1995)), $C_{s}$ is thought to be a function of $Pr$ , e.g. following Schlichting & Gersten (Reference Schlichting and Gersten2003),

(3.17) $$\begin{eqnarray}C_{s}=13.7Pr^{2/3}-7.5,\end{eqnarray}$$

or following Davidson (Reference Davidson2004),

(3.18) $$\begin{eqnarray}C_{s}=1.67(3Pr^{1/3}-1)^{2}.\end{eqnarray}$$

Figure 2 confirms such an effect of $Pr$ on the log-law layer. The empirical estimates of $C_{s}$ as a function of $Pr$ , i.e. (3.17) and (3.18), agree well with DNS, as shown in figure 3.

Figure 2.

Normalized temperature difference from the wall value, $\unicode[STIX]{x1D6E9}^{+}$ as defined in (3.14), plotted as a function of normalized wall distance $y^{+}$ . (a) $Pr=0.7,Ri=0$ (simulation 1); (b) $Pr=7,Ri=0$ (simulation 4). Circles show DNS data; the conductive law (3.15) is plotted with a solid line; the logarithmic law (3.16) in which the additive constant $C_{s}$ varies with $Pr$ , is plotted with a dashed line; and the dot-dashed line shows the location of $y^{+}=\unicode[STIX]{x0394}y^{+}$ , where $\unicode[STIX]{x0394}y^{+}$ marks the characteristic height of the conductive sublayer as defined in (3.19).

The value of $C_{s}$ effectively determines the height of the conductive sublayer which can be measured by $\unicode[STIX]{x0394}y^{+}$ (as marked with vertical dot-dashed lines in figure 2), the intersect of the conductive law (3.15) and the log law (3.16), i.e.

(3.19) $$\begin{eqnarray}\unicode[STIX]{x0394}y^{+}Pr=\frac{1}{k_{s}}\ln \unicode[STIX]{x0394}y^{+}+C_{s}(Pr).\end{eqnarray}$$

The quantity $\unicode[STIX]{x0394}y^{+}$ is observed to decrease with $Pr$ (see figures 2 and 3), and, in particular, for $Pr\gg 1$ (Davidson Reference Davidson2004),

(3.20) $$\begin{eqnarray}\unicode[STIX]{x0394}y^{+}\propto Pr^{-1/3}.\end{eqnarray}$$

With (3.15), the temperature difference across the conductive sublayer, i.e.

(3.21) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D6E9}^{+}\sim Pr\unicode[STIX]{x0394}y^{+}\propto Pr^{2/3},\end{eqnarray}$$

varies strongly with $Pr$ .

Figure 3.

Effect of $Pr$ on the near-wall layer. (a) Variation with $Pr$ of the additive constant $C_{s}$ in the log law for scalar (3.16) as determined by simulations 1 and 4 of stratified plane Couette flows and empirical relations (3.17) and (3.18). (b) The height of the conductive sublayer $\unicode[STIX]{x0394}y^{+}$ as a function of $Pr$ . $\unicode[STIX]{x0394}y^{+}$ values obtained by solving (3.19) and (3.17) are plotted with a solid line, while the dashed line shows the scaling $\unicode[STIX]{x0394}y^{+}\propto Pr^{-1/3}$ .

It is thus shown that $Pr$ has a significant effect on the near-wall structure of the mean scalar field. A thinner conductive layer is expected at higher values of $Pr$ , as suggested by (3.20). Moreover, as the temperature gradients (in wall units) are sharper at a larger $Pr$ , as quantified by (3.15), the temperature jump across the conductive sublayer increases with $Pr$ , as quantified by (3.21). This generic behaviour of the ‘law of the wall’ for varying $Pr$ has implications for the overall temperature profile across the channel gap in stratified plane Couette flows, as we discuss in detail in § 4.

3.3 Damping functions

To complete the mixing length specifications by taking into account the near-wall layer and the effect of $Pr$ mentioned above, one can apply the van Driest damping functions (van Driest Reference van Driest1956) to the mixing lengths in (3.13). This near-wall correction improves the modelling of the turbulent fluxes in terms of their dependence on $y_{w}$ in the viscous/conductive sublayer (Pope Reference Pope2000). The momentum mixing length is corrected by the damping function $\text{D}_{m}(y^{+})$ to become

(3.22) $$\begin{eqnarray}\ell _{m}^{\ast }=k_{m}y_{w}\unicode[STIX]{x1D6F7}_{m}^{-1}(\unicode[STIX]{x1D709})\text{D}_{m}(y^{+})=k_{m}y_{w}\unicode[STIX]{x1D6F7}_{m}^{-1}(\unicode[STIX]{x1D709})[1-\exp (-y^{+}/A_{m}^{+})],\end{eqnarray}$$

where the van Driest constant for momentum $A_{m}^{+}$ is set to be 26 (van Driest Reference van Driest1956; Pope Reference Pope2000).

Similarly, the scalar mixing length becomes

(3.23) $$\begin{eqnarray}\ell _{s}^{\ast }=k_{s}y_{w}\unicode[STIX]{x1D6F7}_{s}^{-1}(\unicode[STIX]{x1D709})\text{D}_{s}(y^{+})=k_{s}y_{w}\unicode[STIX]{x1D6F7}_{s}^{-1}(\unicode[STIX]{x1D709})[1-\exp (-Pr^{-1}y^{+}/A_{s}^{+})],\end{eqnarray}$$

where the constant $A_{s}^{+}$ is inherently related to the $Pr$ -dependent additive constant $C_{s}$ in (3.16) (Pope Reference Pope2000) and is thus also a function of $Pr$ .

As $y^{+}\rightarrow 0$ , the turbulent diffusivity $\unicode[STIX]{x1D705}_{t}$ in the conductive sublayer, following (3.23), scales as

(3.24) $$\begin{eqnarray}\unicode[STIX]{x1D705}_{t}=\ell _{s}^{\ast }u^{\ast }=\ell _{s}^{\ast }\ell _{m}^{\ast }\left|\frac{\text{d}U}{\text{d}y}\right|\sim k_{s}k_{m}\frac{y^{+4}}{A_{s}^{+}A_{m}^{+}}\frac{\unicode[STIX]{x1D708}}{Pr}\sim k_{s}k_{m}\frac{y^{+4}}{A_{s}^{+}A_{m}^{+}}\unicode[STIX]{x1D705}.\end{eqnarray}$$

Note that (3.24) does not yield the expected power law, i.e. $\unicode[STIX]{x1D705}_{t}\propto y^{3}$ , that describes the near-wall variation of $\unicode[STIX]{x1D705}_{t}$ , which is a shortcoming of the van Driest model (see Pope Reference Pope2000, p. 305). We use the standard van Driest model for its simplicity. More sophisticated near-wall treatments for large Prandtl (Schmidt) number can be found in e.g. van Reeuwijk & Hadžiabdić (Reference van Reeuwijk and Hadžiabdić2015). The inclusion of $Pr^{-1}$ in the scalar damping function $\text{D}_{s}(y^{+})$ in (3.23) is such that $\unicode[STIX]{x1D705}_{t}$ in the near-wall limit is proportional to the molecular diffusivity $\unicode[STIX]{x1D705}$ (rather than $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D705}Pr$ ).

The quantity $A_{s}^{+}$ is, by definition, a dimensionless wall distance below which the damping takes place. A natural choice for $A_{s}^{+}$ is to take $A_{s}^{+}\sim \unicode[STIX]{x0394}y^{+}$ where the latter is a characteristic height of the conductive sublayer as defined by (3.19). In this model, we take $A_{s}^{+}=0.65\unicode[STIX]{x0394}y^{+}$ . This results in $A_{s}^{+}$ values of $\{7.9,4.3,2.1\}$ respectively for $Pr\in \{0.7,7,70\}$ .

Figure 4.

Comparison of the model prediction of $L^{+}$ and $Re_{\unicode[STIX]{x1D70F}}$ with DNS data from the present study and Deusebio et al. (Reference Deusebio, Caulfield and Taylor2015).  $L_{model}^{+}$ and $Re_{\unicode[STIX]{x1D70F},model}$ are the results of the mixing length model as described in § 3; $L_{DNS}^{+}$ and $Re_{\unicode[STIX]{x1D70F},DNS}$ are the results of DNS which are tabulated in table 1. Varying Reynolds numbers are used in the simulations with $Pr=0.7$ (plotted with circles) and the fill colour is made darker for larger values of $Re$ .

3.4 Comparisons with DNS

This Monin–Obukhov mixing length model as described above can be implemented to produce predictions of wall fluxes $u_{\unicode[STIX]{x1D70F}}^{2}$ and $q_{w}$ and the dimensionless parameters defined in terms of the various fluxes, given the external parameters $Re$ , $Ri$ and $Pr$ . Figure 4 shows the comparisons between the model predictions of $L^{+}\equiv L/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ and $Re_{\unicode[STIX]{x1D70F}}\equiv h/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ (see definitions in (1.3) and (2.6) respectively) and DNS results at $Re=4250$ (the present study) and a (crucially) wider range of $Re$ values (Deusebio et al. Reference Deusebio, Caulfield and Taylor2015) as listed in table 1. Given the considerable range of parameters, the agreement of the model predictions with DNS data is reasonable, with the $L_{2}$ norm of percentage relative errors being 16.4 % for $L^{+}$ and 13.9 % for $Re_{\unicode[STIX]{x1D70F}}$ over all simulations tested in figure 4. We believe that the model is thus validated and can be employed to produce an estimate of the wall fluxes given the (externally set) $(Re,Ri,Pr)$ parameters.

5.1 Equilibrium Richardson number

We first revisit the midgap gradient Richardson number $Ri_{g}|_{y=0}$ for fully developed stationary (equilibrium) stratified plane Couette flows as prescribed by Monin–Obukhov scaling. The concept of just such a characteristic equilibrium $Ri_{g}$ value was discussed by Turner (Reference Turner1973) in the context of constant flux layers. There also exists a large body of literature considering the ‘stationary Richardson number’ in homogeneous sheared stratified turbulence, e.g. see Shih et al. (Reference Shih, Koseff, Ferziger and Rehmann2000), where the particular value of the gradient Richardson number is imposed by construction, and the references therein. A more recent discussion by Galperin, Sukoriansky & Anderson (Reference Galperin, Sukoriansky and Anderson2007) questioned whether such a unique ‘critical Richardson number’ exists, although the flows considered there differed in several significant ways from the flows considered here. Specifically, and most importantly, stratified plane Couette flow exhibits intermittency for the bulk Richardson number $Ri\lesssim O(1)$ . Also, as we discuss in more detail below, the turbulent Prandtl number, i.e. the ratio of eddy diffusivities of heat and momentum, behaves in a qualitatively different manner in stratified plane Couette flow from the behaviour of the ‘quasi-normal scale elimination’ (QNSE) model used in Galperin et al. (Reference Galperin, Sukoriansky and Anderson2007). Under the plausible assumption that a critical Richardson number exists at least in the flow geometry under consideration here, it may help us to assess if the turbulence would be self-sustained if the externally imposed Richardson number matches the equilibrium condition, or if the flow would self-adjust under non-equilibrium conditions (Turner Reference Turner1973). Examples of the adjustment in the latter scenario include the formation of ‘layer’ and ‘interface’ structures through the rearrangement of velocity and density profiles so that the equilibrium Richardson number is maintained everywhere in the vertical direction (see § 10, Turner Reference Turner1973).

Figure 7(a) compares the midgap equilibrium $Ri_{g}|_{y=0}$ values from DNS data (from both the present study and those by DCT crucially at a range of $Re$ ) with the model prediction (B 3) derived in appendix B, i.e.

(5.1) $$\begin{eqnarray}Ri_{g}|_{y=0}=\frac{k_{m}}{k_{s}}\frac{(h/L)^{-1}+\unicode[STIX]{x1D6FD}_{s}}{[(h/L)^{-1}+\unicode[STIX]{x1D6FD}_{m}]^{2}},\end{eqnarray}$$

which suggests that such an equilibrium $Ri_{g}$ value is determined solely by the length scale ratio $h/L$ (note that $k_{m}$ , $k_{s}$ , $\unicode[STIX]{x1D6FD}_{s}$ and $\unicode[STIX]{x1D6FD}_{m}$ are model constants defined in § 3). The data points indeed collapse in figure 7(a) for the wide range of external parameters (in particular, Prandtl number $Pr$ , but also Reynolds number $Re$ ) examined, and the DNS results compare well with the Monin–Obukhov prediction (5.1).

Two scenarios in stratified plane Couette flows arise from (5.1) when $h/L$ approaches different limits. First, when $h/L\rightarrow \infty$ , the midgap equilibrium $Ri_{g}$ saturates at

(5.2) $$\begin{eqnarray}Ri_{g}|_{y=0}=\frac{k_{m}}{k_{s}}\frac{\unicode[STIX]{x1D6FD}_{s}}{\unicode[STIX]{x1D6FD}_{m}^{2}}\simeq 0.21.\end{eqnarray}$$

This scenario is at least superficially similar to the discussion of constant-flux layers in ‘very stable’ stratification (Ellison Reference Ellison1957; Turner Reference Turner1973), although as discussed further below, the behaviour of the turbulent Prandtl number is qualitatively different from that assumed by Ellison (Reference Ellison1957). When $\unicode[STIX]{x1D709}=h/L\gg 1$ , the linear dependence of Monin–Obukhov functions $\unicode[STIX]{x1D6F7}_{m}$ and $\unicode[STIX]{x1D6F7}_{s}$ on $\unicode[STIX]{x1D709}$ dominates (see (3.11) and (3.12)). Fluid in the channel gap interior does not ‘feel’ the impact of the wall directly (but still indirectly though the wall fluxes $u_{\unicode[STIX]{x1D70F}}^{2}$ and $q_{w}$ ), because the vertical motions are strongly damped by stratification. In the channel gap interior, the distance to the wall $y_{w}$ (or the channel gap half-height $h$ ) becomes irrelevant, as shear and temperature gradients both become constant (by taking the limit of (3.12) at $\unicode[STIX]{x1D709}\rightarrow \infty$ ), which renders the turbulence close to homogeneous in the wall-normal direction. Interestingly, the maximum stationary $Ri_{g}$ reported in homogeneous sheared turbulence is also approximately 0.2 (see e.g. Shih et al. Reference Shih, Koseff, Ferziger and Rehmann2000). This reinforces the notion that Monin–Obukhov scaling may also apply to such homogeneous triply periodic flows (Chung & Matheou Reference Chung and Matheou2012). There remains some debate as to whether the standard Monin–Obukhov theory holds in the $\unicode[STIX]{x1D709}\rightarrow \infty$ limit in a stable atmospheric boundary layer, see for example the discussion on ‘ $z$ -less’ stratification by Mahrt (Reference Mahrt1999), and any such differences between the standard theory and boundary layer flow are likely linked to the variation of fluxes with height in a real boundary layer. However, the statistically stationary stratified plane Couette flows examined here, which are constant-flux layers by construction, appear to be consistent with the standard Monin–Obukhov theory.

As an aside, we note that this maximum observed Richardson number is close to $Ri_{g}=1/4$ which arises in the well-known Miles–Howard criterion for linear normal mode stability of inviscid parallel steady stratified shear flows (Howard Reference Howard1961; Miles Reference Miles1961). This closeness is apparently fortuitous, as the arguments leading to the prediction of the value in (5.2) are entirely constructed under the assumption of statistically stationary turbulent flow. Therefore, it is at least possible that observations of $Ri_{g}$ close to $1/4$ , as, for example, in the equatorial undercurrent (Smyth & Moum Reference Smyth and Moum2013), are due to turbulent balances, not ‘marginal stability’ of the flow, as argued by Thorpe & Liu (Reference Thorpe and Liu2009), although, it is also important to remember, as shown for example by Pham, Sarkar & Winters (Reference Pham, Sarkar and Winters2013), that the dynamics of the Equatorial Undercurrent is inevitably non-stationary, due to diurnal forcing.

Second, when $h/L$ is $O(1)$ or smaller, the equilibrium $Ri_{g}$ at midgap varies strongly with $h/L$ , which can be seen from figure 7(a). Under this scenario, the stabilising effects due to stratification are relatively weak. The direct influence of the walls on the interior turbulence becomes significant, and both $h$ and $L$ become relevant scales for the channel gap interior.

5.2 $L^{+}$ , $Re_{b}$ and intermittency

The parameter $L^{+}$ is a useful diagnostic quantity to predict if stratified plane Couette flows can sustain a fully turbulent state or become intermittent (as discussed by DCT). On the other hand, the buoyancy Reynolds number $Re_{b}\equiv \unicode[STIX]{x1D700}/(\unicode[STIX]{x1D708}N^{2})\sim (\ell _{O}/\unicode[STIX]{x1D702})^{4/3}$ , which describes the scale separation between the Ozmidov scale $\ell _{O}$ and the Kolmogorov scale $\unicode[STIX]{x1D702}$ , is often used to predict whether small scale turbulence can exist given the level of turbulent dissipation and stratification (see e.g. Riley & Lindborg Reference Riley, Lindborg, Davidson, Kaneda and Sreenivasan2012), typically in homogeneous simulations (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007).

A natural question to ask is then whether $L^{+}$ and $Re_{b}$ are related to each other, at least in stratified plane Couette flows. The analysis in appendix B has, through Monin–Obukhov similarity theory, predicted a linear scaling between $L^{+}$ and $Re_{b}$ as given by (B 5) shown in appendix B, i.e.

(5.3) $$\begin{eqnarray}Re_{b}\sim L^{+}k_{m}.\end{eqnarray}$$

In figure 7(b) this scaling is confirmed from DNS data (shown for simulations 1–12), and has already been noted by Scotti & White (Reference Scotti and White2016) in a more limited range of $Ri\leqslant 0.1$ , $Re\leqslant 55\,000$ and $Pr=1$ . $Re_{b}$ estimates presented here are based on $\unicode[STIX]{x1D700}$ and $N$ values that are sampled pointwise in the vertical direction $y$ . However, in open flows, there are different possible choices of averaging volumes for $\unicode[STIX]{x1D700}$ and $N$ (see e.g. Salehipour, Caulfield & Peltier Reference Salehipour, Caulfield and Peltier2016), and caution needs to be exercised when comparing specific numerical values of $Re_{b}$ between different flow geometries, or indeed between different analyses. A reanalysis of DCT’s data (simulations 13–23, not shown) suggests the same linear scaling for a wide range of $Re$ and $Ri$ . This indicates that the $L^{+}$ criterion for predicting intermittency, which is specific to wall-bounded flows, is also linked to this more general $Re_{b}$ argument. The critical (minimum) $Re_{b}$ for fully developed turbulence, as inferred from the $L^{+}>200$ criterion reported by DCT and the scaling (5.3), is approximately 80 (as $k_{m}\approx 0.4$ ) for stratified plane Couette flows. This critical $Re_{b}$ of 80 is close to the cutoff value $Re_{b}=100$ between the ‘intermediate’ and ‘energetic’ regimes of Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) which is discussed in detail in § 6, although one needs to be careful about whether the $Re_{b}$ value is a ‘bulk’ or local estimate when comparing the numerical values. Here simulation 12 is in the intermediate regime ( $Re_{b}<35$ , see figure 7 b), and in what follows, we focus instead on the other simulations ( $Re_{b}>60$ ) which are close to or within this ‘energetic’ regime in terms of the $Re_{b}$ value.

5.3 Turbulent Prandtl number

In appendix B, it is shown through scaling arguments that the flux Richardson number $Ri_{f}$ is proportional to $Ri_{g}$ . The particular scaling derived in appendix B is given by (B 7) i.e.

(5.4) $$\begin{eqnarray}Ri_{f}\sim Ri_{g},\end{eqnarray}$$

and is compared to DNS results (simulations 1–12) in figure 7(c). In general, $Ri_{f}$ is proportional to $Ri_{g}$ with a multiplicative constant of approximately unity, which is consistent with DCT. The group of points which appear to be outliers, correspond to simulation 12 ( $Pr=70,Ri=1.44$ ). As discussed previously, the atypical behaviour associated with this simulation is likely to be due to low- $Re_{b}$ , and hence inherently viscously dominated effects.

With the turbulent viscosity $\unicode[STIX]{x1D708}_{t}$ defined through the flux gradient relation

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{t}\equiv -\frac{\langle u^{\prime }v^{\prime }\rangle }{S},\end{eqnarray}$$

and turbulent diffusivity $\unicode[STIX]{x1D705}_{t}$ defined in (A 5), the turbulent Prandtl number $Pr_{t}\equiv \unicode[STIX]{x1D708}_{t}/\unicode[STIX]{x1D705}_{t}$ can be expressed as

(5.6) $$\begin{eqnarray}Pr_{t}=\frac{Ri_{g}}{Ri_{f}}.\end{eqnarray}$$

The $Ri_{f}\simeq Ri_{g}$ scaling can thus be interpreted alternatively as the turbulent Prandtl number $Pr_{t}$ being approximately unity, which is consistent with the Reynolds analogy, as noted independently by Scotti & White (Reference Scotti and White2016). This result can be derived from Monin–Obukhov theory (appendix B) and is consistent with DNS data for the present study (simulations 1–12) shown in figure 7(c), as well as from revisited DCT datasets (simulations 13–23) which exhibit the same behaviour (not shown).

$Pr_{t}$ is often parameterized as a function of $Ri_{g}$ in the literature (see, for example, Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2010). $Pr_{t}$ being unity, as we observe in stratified plane Couette flows (see figure 7 c), appears to be typical for gradient Richardson numbers $Ri_{g}<0.2$ which are sufficiently small in this context – again, one needs to be careful about the exact definition of $Ri_{g}$ when comparing across different studies, and also it is necessary to remember that this is distinct from keeping the bulk Richardson number $Ri$ (set by the boundary conditions) small. This observation is consistent with previous studies of stably stratified wall-bounded flow simulations (Armenio & Sarkar Reference Armenio and Sarkar2002; García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; García-Villalba et al. Reference García-Villalba, Azagra, Uhlmann and Nagel2011a ) and in homogeneous stratified turbulence (Rohr & Van Atta Reference Rohr and Van Atta1987; Chung & Matheou Reference Chung and Matheou2012).

The behaviour of $Pr_{t}$ becomes more complex at higher values of $Ri_{g}$ , i.e. for $Ri_{g}>0.2$ (Taylor, Sarkar & Armenio Reference Taylor, Sarkar and Armenio2005; Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2010; Karimpour & Venayagamoorthy Reference Karimpour and Venayagamoorthy2014, Reference Karimpour and Venayagamoorthy2015; Salehipour & Peltier Reference Salehipour and Peltier2015; Wilson & Venayagamoorthy Reference Wilson and Venayagamoorthy2015). However, turbulent flows with larger gradient Richardson numbers $Ri_{g}>0.2$ do not appear to be accessible in stratified plane Couette flows, for reasons that have been discussed in § 5.1. There also exist $Re_{b}$ -based parameterizations for $Pr_{t}$ in the literature. Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005) and Salehipour & Peltier (Reference Salehipour and Peltier2015) reported $Pr_{t}$ approaching order unity for intermediate to large values of $Re_{b}$ , which is consistent with our observations. Salehipour & Peltier (Reference Salehipour and Peltier2015) also observed larger than $O(1)$ values of $Pr_{t}$ when the values of $Re_{b}$ are small, i.e. $O(1)$ to $O(10)$ . This is consistent with our outlier group (simulation 12) in figure 7(c) whose $Re_{b}$ value is $O(10)$ (see figure 7 b) and the $Pr_{t}$ value is larger than unity ( $Ri_{g}\gg Ri_{f}$ ).

Crucially, all the evidence points towards $Pr_{t}\sim O(1)$ while the flow is turbulent, with the flow becoming intermittent before $Ri_{g}$ reaching large values. This is qualitatively different behaviour to that assumed by Ellison (Reference Ellison1957), who stated that ‘it seems more likely’ that turbulence could be ‘maintained’ at large values of $Ri_{g}$ with still finite $Ri_{f}<1$ , and so, from (5.6) and consequences derived from it with further turbulence modelling assumptions, Ellison (Reference Ellison1957) was led to the conclusion that $Pr_{t}$ inevitably reaches large values. Galperin et al. (Reference Galperin, Sukoriansky and Anderson2007) analogously arrived at the conclusion that $Pr_{t}$ reaches large values in strongly stratified, yet still ‘turbulent’ flows, in the relatively weak sense that the eddy diffusivities (particularly in the horizontal) remain elevated above molecular values. A potential major point of difference leading to this apparent qualitative discrepancy in the behaviour of the turbulent Prandtl number is the central role played in open flows of propagating internal waves, which is not possible in stratified plane Couette flow.

5.4 Realizability of strongly stratified regime

Finally, we test the scaling in (B 11) in appendix B, i.e.

(5.7) $$\begin{eqnarray}Fr_{h}\sim \frac{1}{\sqrt{Ri_{g}}},\end{eqnarray}$$

for the turbulent Froude number $Fr_{h}$ . Figure 7(d) shows the DNS results for which an empirical scaling of

(5.8) $$\begin{eqnarray}Fr_{h}\simeq \frac{0.95}{\sqrt{Ri_{g}}}\end{eqnarray}$$

applies, which is consistent with the Monin–Obukhov prediction in appendix B. The outlier once again corresponds to simulation 12 for which the $Re_{b}$ value may not be high enough for the inertially dominated forward cascade assumption underlying (B 8), i.e. $\unicode[STIX]{x1D700}={U^{\prime }}^{3}/\ell _{h}$ , to hold.

Given that the maximum $Ri_{g}$ in fully developed stratified plane Couette flow is approximately 0.2 (see § 5.1), the minimum $Fr_{h}$ that can be obtained in the interior of an stratified plane Couette flow (at large enough $Re_{b}$ ) is approximately 2, following (5.8). However, for the turbulence to reach the strongly stratified regime, it is typically argued that $Fr_{h}$ needs to be smaller than 0.02 (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). Therefore, the strongly stratified regime, which is characterized by layering in the density field with characteristic vertical length scale $U^{\prime }/N$ , may be fundamentally non-realizable in stratified plane Couette flows, at least under the equilibrium conditions we have been considering. Once again, it is important to emphasize that it is the midgap gradient Richardson number $Ri_{g}$ which cannot become large in quasi-steady turbulent stratified plane Couette flow, for any choice of $Re$ and $Ri$ set by the boundary conditions.

5.5 Summary

To summarize the results in § 5, we have identified certain generic characteristics of the turbulence in the interior regions of stratified plane Couette flows. We find that: the length scale ratio $h/L$ determines the midgap $Ri_{g}$ ; $Re_{b}$ scales linearly with $L^{+}\equiv L/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ ; $Pr_{t}$ is of order unity for the range of accessible $Ri_{g}$ associated with turbulence; and $Fr_{h}$ is proportional to $Ri_{g}^{-1/2}$ . The scalings, consistent with, and extending the observations of Scotti & White (Reference Scotti and White2016) into the crucially important regime where the externally set bulk $Ri>0.1$ , apply not only to the DNS performed for the present study which cover a wide range of $Pr$ (simulations 1–12), but also to those by DCT as listed in table 1 which covered a wider range of $Re$ (simulations 13–23). These characteristics of stratified plane Couette flows fundamentally relate to the fact that it is the upper and lower walls which impose momentum and heat fluxes on the fluids. These fluxes then dictate the self-similar behaviour of both the mean flow (as characterized by $Ri_{g}$ ) and the turbulence (as characterized by $Re_{b}$ , $Pr_{t}$ and $Fr_{h}$ ) in the interior. These results are expected to hold not only for stratified plane Couette flows but also for other constant-flux layers to which the Monin–Obukhov scaling applies. These Monin–Obukhov scalings are intended for regions sufficiently far from the walls. Through the wide range of Prandtl numbers examined, our DNS data suggest that the dynamics away from the walls is $Pr$ -independent for given wall fluxes.

6.1 Osborn formulation for stratified plane Couette flow

Now we turn our attention to the third main theme of interest, namely the parameterization of mixing. Here we use the framework proposed by Osborn (Reference Osborn1980) to formulate a parameterization for the turbulent diffusivity $\unicode[STIX]{x1D705}_{t}\equiv -\langle \unicode[STIX]{x1D70C}^{\prime }v^{\prime }\rangle /(\text{d}\overline{\unicode[STIX]{x1D70C}}/\text{d}y)=-B/N^{2}$ . As described in appendix C, key to this formulation is the turbulent flux coefficient, $\unicode[STIX]{x1D6E4}\equiv B/\unicode[STIX]{x1D700}\approx Ri_{f}/(1-Ri_{f})$ . With $\unicode[STIX]{x1D6E4}$ appropriately parameterized, the Osborn formulation yields an expression for $\unicode[STIX]{x1D705}_{t}$ , i.e.

(6.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D705}_{t}}{\unicode[STIX]{x1D708}}\approx \frac{Ri_{f}}{1-Ri_{f}}\frac{\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D708}N^{2}}=\unicode[STIX]{x1D6E4}Re_{b}.\end{eqnarray}$$

It is important to appreciate that key aspects of the Osborn (Reference Osborn1980) framework are based on the theoretical considerations of Ellison (Reference Ellison1957) and the experimental data of Britter (Reference Britter1974), both associated with stratified flows in the presence of boundary forcing and thus expected to have at least some similar properties to the turbulence in stratified plane Couette flows. Osborn (Reference Osborn1980), following Ellison (Reference Ellison1957), postulated that $\unicode[STIX]{x1D6E4}\leqslant 0.2$ , or equivalently $Ri_{f}\leqslant 0.15$ , although the inequality in Osborn’s original paper has often been ignored subsequently. Interestingly, the experimental data by Britter (Reference Britter1974) (see e.g. pp. 8–37 of the thesis) led to his conclusion that ‘a critical Richardson flux number (i.e. $Ri_{f}$ ) of approximately 0.2 is predicted’. This is entirely consistent with our results presented in § 5 that

(6.2) $$\begin{eqnarray}Ri_{f}\simeq Ri_{g}\lesssim 0.2\end{eqnarray}$$

in stratified plane Couette flows for turbulence to be maintained, although as already noted we observe the turbulent Prandtl number remaining of order one, unlike in the model developed by Ellison (Reference Ellison1957). Indeed, using this scaling, $\unicode[STIX]{x1D6E4}$ can be written as a function of the gradient Richardson number $Ri_{g}$ :

(6.3) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}\approx \frac{Ri_{g}}{1-Ri_{g}},\end{eqnarray}$$

remembering that $Ri_{g}$ appears to have an upper bound above which turbulence cannot be maintained, even for asymptotically large $Re$ (see figure 18 of DCT), In the literature, however, $\unicode[STIX]{x1D6E4}$ is often parameterized as a function of $Re_{b}$ (see e.g. Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005). The connection between the $Ri_{g}$ -based and $Re_{b}$ -based scalings for $\unicode[STIX]{x1D6E4}$ is discussed further in § 6.3.1. It follows from (6.1) and (6.3) that

(6.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D705}_{t}}{\unicode[STIX]{x1D708}}\approx \frac{Ri_{g}}{1-Ri_{g}}Re_{b}\end{eqnarray}$$

in the context of stratified plane Couette flows. Noting that $Pr_{t}\equiv \unicode[STIX]{x1D708}_{t}/\unicode[STIX]{x1D705}_{t}\approx 1$ in stratified plane Couette flows (as shown in § 5) and as also noted by Scotti & White (Reference Scotti and White2016), we can also approximate the turbulent viscosity $\unicode[STIX]{x1D708}_{t}$ with the same scaling for $\unicode[STIX]{x1D705}_{t}$ in (6.4), i.e.

(6.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D708}_{t}}{\unicode[STIX]{x1D708}}\approx \frac{Ri_{g}}{1-Ri_{g}}Re_{b}.\end{eqnarray}$$

Figure 8.

$\unicode[STIX]{x1D705}_{t}$ and $\unicode[STIX]{x1D708}_{t}$ , as defined by (A 5) and (5.5), and both normalized by $\unicode[STIX]{x1D708}$ , as a function of $Re_{b}$ ; $\unicode[STIX]{x1D705}_{t}$ , $\unicode[STIX]{x1D708}_{t}$ and $Re_{b}$ values are computed pointwise in $y$ in the channel gap interior with $y^{+}>50$ . Scaling laws of $\unicode[STIX]{x1D705}_{t}/\unicode[STIX]{x1D708}=0.2Re_{b}$ (Osborn Reference Osborn1980) plotted with a solid line, and $\unicode[STIX]{x1D705}_{t}/\unicode[STIX]{x1D708}=2Re_{b}^{1/2}$ (Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005) plotted with a dashed line, are also shown.

6.2 Numerical results

These $\unicode[STIX]{x1D705}_{t}$ and $\unicode[STIX]{x1D708}_{t}$ values are estimated directly using their definitions through the flux gradient relation (A 5) and (5.5) at all locations in the wall-normal direction $y$ that are at least 50 wall units ( $y^{+}>50$ ) away from the walls, where the local equilibrium (A 3) is expected to hold (García-Villalba et al. Reference García-Villalba, Azagra and Uhlmann2011b ). These results are first plotted in figure 8 to test the $Re_{b}$ -based parameterizations that are commonly seen in the literature, e.g. those reviewed by Ivey et al. (Reference Ivey, Winters and Koseff2008) and also discussed in Scotti & White (Reference Scotti and White2016). Our results are plotted in figure 9 to validate the scalings (6.4) and (6.5). Simulation 12, in which the flow is viscously controlled and exhibits spuriously small ( $O(1)$ or smaller) or negative (counter-gradient) values of $\unicode[STIX]{x1D705}_{t}/\unicode[STIX]{x1D708}$ or $\unicode[STIX]{x1D708}_{t}/\unicode[STIX]{x1D708}$ , is not included in the plots to allow the discussion to stay focused on the fully turbulent simulations.

Figure 8 compares the DNS results of $\unicode[STIX]{x1D705}_{t}/\unicode[STIX]{x1D708}$ against the classical $Re_{b}$ -based parameterizations of Osborn (Reference Osborn1980) and Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005). The DNS data points in figure 8 are sampled locally (pointwise) at various $y$ locations across the channel gap interior of stratified plane Couette flows. Within each simulation, the $Re_{b}$ value stays relatively constant, while the diffusivities span a wider range – the latter is somewhat expected because $\unicode[STIX]{x1D705}_{t}$ and $\unicode[STIX]{x1D708}_{t}$ scale linearly with the mixing lengths $\ell _{s}^{\ast }$ and $\ell _{m}^{\ast }$ respectively, both of which increase with the wall distance $y_{w}$ , as described in § 3.1. These $Re_{b}$ -based scalings are effective in describing the homogeneous flow dataset of Shih et al. (Reference Shih, Koseff, Ivey and Ferziger2005), but they do not provide a good agreement with our DNS data from stratified plane Couette flows which are inherently inhomogeneous due in particular to the presence of the wall. The data for $\unicode[STIX]{x1D708}_{t}$ , which are also plotted in figure 8, behave similarly to $\unicode[STIX]{x1D705}_{t}$ , since the turbulent Prandtl number $Pr_{t}\equiv \unicode[STIX]{x1D708}_{t}/\unicode[STIX]{x1D705}_{t}$ is approximately unity (as shown in § 5.3).

Figure 9 compares the DNS data against the scalings (6.4) and (6.5). The collapse of the DNS data improves significantly when $Ri_{g}$ is included in the parameterizations, as they capture the critical (linear) dependence of $Ri_{f}$ on $Ri_{g}$ . At sufficiently large values of $Re_{b}$ , i.e. $Re_{b}\gtrsim 60$ , the $\unicode[STIX]{x1D705}_{t}/\unicode[STIX]{x1D708}\sim \unicode[STIX]{x1D708}_{t}/\unicode[STIX]{x1D708}\sim Re_{b}Ri_{g}/(1-Ri_{g})$ scaling, based on the turbulent kinetic energy budget argument by Osborn (Reference Osborn1980) and incorporating Monin–Obukhov scaling for constant-flux layers to account for the importance of the (coupled) value of $Ri_{g}$ , provides an accurate description of the turbulent diffusivity in stratified plane Couette flows. It is certainly of interest that the Osborn scaling appears to hold, at least qualitatively, even though the underlying assumption of Ellison (Reference Ellison1957) (on which the Osborn scaling is at least partially based) that $Pr_{t}$ becomes large is violated in stratified plane Couette flow. We further discuss this scaling with respect to other previously proposed scalings in the next subsection.

Figure 9.

$\unicode[STIX]{x1D705}_{t}$ and $\unicode[STIX]{x1D708}_{t}$ , both normalized by $\unicode[STIX]{x1D708}$ , as a function of $Re_{b}Ri_{g}/(1-Ri_{g})$ ; $\unicode[STIX]{x1D705}_{t}$ , $\unicode[STIX]{x1D708}_{t}$ , $Re_{b}$ and $Ri_{g}$ values are computed pointwise in $y$ in the channel gap interior with $y^{+}>50$ . The dashed line marks equality between the abscissa and the ordinate.