2 Reduced model development

The budget equations for the variances of individual turbulent velocity components reveal the exchange of kinetic energy through the ‘pressure redistribution’ terms. These terms, as their alternative ‘return to isotropy’ name indicates, nudge turbulence towards an isotropic state where its kinetic energy would be equipartitioned among its three velocity components (the redistribution terms sum exactly to zero in incompressible flows and thus do not appear in the equation for the full TKE budget). However, when TKE production in one or two of the components dominate, the flow cannot reach such an isotropic state except at scales commensurate with the Kolmogorov microscale $\unicode[STIX]{x1D702}_{K}$ (provided $\unicode[STIX]{x1D702}_{K}\ll$ than the Dougherty–Ozmidov scale (Dougherty Reference Dougherty1961) so that the anisotropic influence of buoyancy at these small scales is negligible). The degree of anisotropy therefore depends on the interplay between the anisotropic production by shear or buoyancy and the effectiveness of the redistribution process that acts differently across scales of motion (Brugger et al. Reference Brugger, Katul, De Roo, Kröniger, Rotenberg, Rohatyn and Mauder2018). It is therefore reasonable to hypothesize that these return-to-isotropy terms encode significant information about the statistical and structural properties of the turbulent velocity, and how such properties are modulated by the shear-to-buoyancy balance encoded in $Ri$ . Buoyancy damping or excitation of vertical velocity must be communicated by these terms to the horizontal directions.

A logical starting point is an idealized high Reynolds number flow that remains three-dimensional with finite TKE in all components. A statistically steady state is assumed, which in geophysical settings must be interpreted as steadiness over limited periods, e.g. 15–60 min in the lower atmosphere (Stull Reference Stull1988; Salesky, Chamecki & Dias Reference Salesky, Chamecki and Dias2012). Over a horizontally homogeneous wall, the flow will also be statistically homogeneous over the streamwise and cross-stream directions, which for an incompressible flow also results in negligible mean vertical velocity. We will neglect the terms in the variance budgets that result from rotation (the Coriolis force) since they are of second-order importance (Stull Reference Stull1988); the effect of the Coriolis force mainly manifests as a rotation of the mean flow. The vertical TKE-transport terms are significant in the viscous and buffer layers, and also potentially in the outer layer (Pope Reference Pope2000), but if one is interested in the inertial sublayer (the logarithmic region), the transport terms are approximately one order of magnitude smaller than the production and dissipation terms and can be neglected in a reduced model (André et al. Reference André, De Moor, Lacarrère and du Vachat1978; Pope Reference Pope2000; Shah & Bou-Zeid Reference Shah and Bou-Zeid2014a ; Momen & Bou-Zeid Reference Momen and Bou-Zeid2017). This assumption might break down under strong wall buoyancy flux or near the density inversion that typically delineates the top of the atmospheric boundary layer (Stull Reference Stull1988; Garcia & Mellado Reference Garcia and Mellado2014), but will be adopted here since it is reasonable in the logarithmic layer under a wide range of stabilities (see appendix A for a model derivation where these terms are retained). An important side note to highlight is that, unlike for kinetic energy, the transport of density (or equivalently potential temperature) fluctuations is usually not negligible, and as a result the production–dissipation balance does not apply to scalar variance (Karimpour & Venayagamoorthy Reference Karimpour and Venayagamoorthy2015) and is not invoked here.

In summary, a high- $Re$ planar–homogeneous steady flow with an equilibrium between TKE production (by shear and buoyancy) and dissipation (by viscosity and buoyancy) is postulated. Such an idealization is in fact quite routinely used in practice when interpreting measurements in the inertial turbulent sublayer in laboratory experiments or flux measurements in the all-important atmospheric surface layer, or when collapsing field experiments in the context of Monin–Obukhov similarity theory (MOST) (Monin & Obukhov Reference Monin and Obukhov1954; Kader & Yaglom Reference Kader and Yaglom1990; Foken Reference Foken2006; Bou-Zeid et al. Reference Bou-Zeid, Higgins, Huwald, Meneveau and Parlange2010; Ghannam et al. Reference Ghannam, Katul, Bou-Zeid, Gerken and Chamecki2018; Li, Katul & Liu Reference Li, Katul and Liu2018). With these assumptions, the velocity half-variance budgets reduce to (Zilitinkevich Reference Zilitinkevich1973; Stull Reference Stull1988):

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle -\overline{u^{\prime }w^{\prime }}\frac{\unicode[STIX]{x2202}\overline{u}}{\unicode[STIX]{x2202}z}+\overline{\frac{p^{\prime }}{\unicode[STIX]{x1D70C}_{r}}\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x}}-\unicode[STIX]{x1D708}\overline{\left(\frac{\unicode[STIX]{x2202}u^{\prime }}{\unicode[STIX]{x2202}x_{j}}\right)^{2}}=S+R_{u}-\unicode[STIX]{x1D700}_{u}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{\frac{p^{\prime }}{\unicode[STIX]{x1D70C}_{r}}\frac{\unicode[STIX]{x2202}v^{\prime }}{\unicode[STIX]{x2202}y}}-\unicode[STIX]{x1D708}\overline{\left(\frac{\unicode[STIX]{x2202}v^{\prime }}{\unicode[STIX]{x2202}x_{j}}\right)^{2}}=R_{v}-\unicode[STIX]{x1D700}_{v}=0, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle -\frac{g}{\unicode[STIX]{x1D70C}_{r}}\overline{w^{\prime }\unicode[STIX]{x1D70C}^{\prime }}+\overline{\frac{p^{\prime }}{\unicode[STIX]{x1D70C}_{r}}\frac{\unicode[STIX]{x2202}w^{\prime }}{\unicode[STIX]{x2202}z}}-\unicode[STIX]{x1D708}\overline{\left(\frac{\unicode[STIX]{x2202}w^{\prime }}{\unicode[STIX]{x2202}x_{j}}\right)^{2}}=B+R_{w}-\unicode[STIX]{x1D700}_{w}=0. & \displaystyle\end{eqnarray}$$

Throughout, $u,v,w$ indicate instantaneous velocity components in the streamwise ( $x$ ), cross-stream ( $y$ ) and vertical ( $z$ ) directions, respectively; $p$ is instantaneous pressure; $\unicode[STIX]{x1D70C}$ is density and $\unicode[STIX]{x1D70C}_{r}$ its reference value; and $\unicode[STIX]{x1D708}$ is kinematic viscosity. Instantaneous quantities are decomposed into their turbulent fluctuations, denoted by a prime, and their averaged states, denoted by an overbar, where averaging is performed over coordinates of statistical homogeneity. The symbols in the middle equations refer, in the same order, to the terms on the left-hand side. The first term on the left-hand side of (2.1) is shear production $(S>0)$ , while in (2.3) the first term is buoyancy production ( $B>0$ ) or destruction ( $B<0$ ). The terms with $p$ are pressure redistribution terms (denoted by $R_{u}$ , $R_{v}$ and $R_{w}$ ), and the last terms on the left-hand side are viscous dissipation rates ( $\unicode[STIX]{x1D700}_{u}$ , $\unicode[STIX]{x1D700}_{v}$ and $\unicode[STIX]{x1D700}_{w}$ , all positive).

The dissipation is performed by viscosity acting on the smallest eddies commensurate in size to $\unicode[STIX]{x1D702}_{K}$ . If $Re$ is sufficiently high, the energy cascade dissipates the anisotropy injected at the largest scales, resulting in an approximately isotropic dissipation with

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D700}_{u}=\unicode[STIX]{x1D700}_{v}=\unicode[STIX]{x1D700}_{w}\quad \Rightarrow \quad \unicode[STIX]{x1D700}_{w}+\unicode[STIX]{x1D700}_{v}=2\unicode[STIX]{x1D700}_{u}.\end{eqnarray}$$

Viscous-range isotropy might not be exact at moderate $Re$ where the limited turnovers between the production and dissipation do not completely remove the signature of the large-scale anisotropy. An anisotropy factor could then be introduced to generalize the model. However, under sufficiently high $Re$ and ‘far’ from any wall effect (i.e. outside the buffer layer), the difference between the three components of dissipation and their mean is less than 10% (Spalart Reference Spalart1988; Pope Reference Pope2000) and is ignored here. Another potential source of viscous-range anisotropy is strong static stability at high $Ri$ where the Dougherty–Ozmidov length scale decreases and could become ${\sim}\unicode[STIX]{x1D702}_{K}$ (Grachev et al. Reference Grachev, Andreas, Fairall, Guest and Persson2015; Li, Salesky & Banerjee Reference Li, Salesky and Banerjee2016); this aspect will be analysed in future studies using direct numerical simulation data.

Equations (2.1)–(2.3) are combined with (2.4) to yield

(2.5) $$\begin{eqnarray}B+R_{w}+R_{v}=2S+2R_{u}.\end{eqnarray}$$

In addition, since the redistribution terms must sum to zero in an incompressible flow,

(2.6) $$\begin{eqnarray}R_{u}+R_{v}+R_{w}=0\quad \Rightarrow \quad R_{v}+R_{w}=-R_{u}.\end{eqnarray}$$

Replacing (2.6) into (2.5), dividing by $S$ and noting that $-B/S=Ri_{f}$ yields

(2.7) $$\begin{eqnarray}\frac{R_{u}}{S}=-\frac{2}{3}-\frac{1}{3}Ri_{f}.\end{eqnarray}$$

In a similar fashion, equations (2.4)–(2.6) can be modified (writing $\unicode[STIX]{x1D700}_{u}+\unicode[STIX]{x1D700}_{w}=2\unicode[STIX]{x1D700}_{v}$ in (2.4) and $R_{u}+R_{w}=-R_{v}$ in (2.6)) to yield

(2.8) $$\begin{eqnarray}\frac{R_{v}}{S}=\frac{1}{3}-\frac{1}{3}Ri_{f}.\end{eqnarray}$$

Then, using $R_{w}=-R_{u}-R_{v}$ results in

(2.9) $$\begin{eqnarray}\frac{R_{w}}{S}=\frac{1}{3}+\frac{2}{3}Ri_{f}.\end{eqnarray}$$

Finally, the total TKE dissipation $\unicode[STIX]{x1D700}$ (the sum of the half-variance dissipations in all three components that is used to define $\unicode[STIX]{x1D702}_{K}=(\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D700})^{1/4}$ ) can be related to $Ri_{f}$ by summing (2.1) to (2.3), and using (2.6), to obtain

(2.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D700}_{u}+\unicode[STIX]{x1D700}_{v}+\unicode[STIX]{x1D700}_{w}}{S}=\frac{\unicode[STIX]{x1D700}}{S}=1-Ri_{f}.\end{eqnarray}$$

3 Validation and implications

Equations (2.7)–(2.10) constitute a basic ‘reduced model’ that describes how energy redistribution and dissipation change with $Ri_{f}$ . Under neutral conditions where $Ri_{f}=0$ , the model assumptions and formulation are very similar to algebraic Reynolds stress models such as those proposed by Launder, Reece & Rodi (Reference Launder, Reece and Rodi1975) (see also Pope (Reference Pope2000), § 11.9). In this neutral limit, the model predicts $R_{u}=-2S/3$ and $R_{v}=R_{w}=S/3$ . This yields the expected result that two thirds of shear production in the $u$ component is redistributed equally into the $v$ and $w$ components where they are dissipated via $\unicode[STIX]{x1D700}_{v}=\unicode[STIX]{x1D700}_{w}$ , and that one third of shear production gets dissipated directly via $\unicode[STIX]{x1D700}_{u}$ . This outcome results in equal energy in the cross-stream and vertical components, which is not observed in wall-bounded flows (Stull Reference Stull1988; Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2005) due to the influence of the mean shear and to the wall pressure-blockage effect (McColl et al. Reference McColl, Katul, Gentine and Entekhabi2016) that is absent here. The effect of this wall blockage is to increase pressure ( $p^{\prime }>0$ ) for descending parcels ( $\unicode[STIX]{x2202}w^{\prime }/\unicode[STIX]{x2202}z<0$ since $w=0$ at the wall), while decreasing pressure ( $p^{\prime }<0$ ) for ascending parcels ( $\unicode[STIX]{x2202}w^{\prime }\unicode[STIX]{x2202}z>0$ ). Wall blockage thus consistently produces negative ( $p^{\prime }\unicode[STIX]{x2202}w^{\prime }/\unicode[STIX]{x2202}z$ ) contributions, or reduces the magnitude of the positive contributions, and as a result reduces the magnitude of the average redistribution of energy into the $w$ component, $\overline{p^{\prime }\unicode[STIX]{x2202}w^{\prime }/\unicode[STIX]{x2202}z}>0$ , below model predictions. The comparative influence of wall blockage relative to the influence of increasing shear in the vicinity of the wall is discussed in Launder et al. (Reference Launder, Reece and Rodi1975).

To assess the impact of neglecting wall blockage and other assumptions, the redistribution terms normalized by shear production (2.7)–(2.9) are compared to results from large eddy simulation (LES) under unstable conditions ( $Ri_{f}<0$ ) and direct numerical simulation (DNS) under stable conditions ( $Ri_{f}>0$ ). The DNS and LES details, with further references, are respectively provided in appendices B and C. It is noted that (i) the pressure redistribution terms are difficult to measure experimentally (no such measurements have been reported in the literature) and thus simulations are the only avenue for estimating them, and (ii) due to subgrid-scale modelling in the LES, these computations are deemed less accurate than in DNS. For the LES and DNS, the vertical average of the redistribution terms and of $Ri_{f}$ over the inertial sublayer (logarithmic or modified logarithmic if non-neutral) are taken for comparisons since this is the layer where our model assumptions hold and that is of most interest in high- $Re$ flows. The three redistribution terms computed from the reduced model and from the simulations are compared in figure 1. Figure 1(a) is for strongly unstable conditions ( $Ri_{f}<-1$ ) and shows that the model captures LES trends but with some discrepancies. Figure 1(b) is for moderately unstable to stable conditions ( $-2<Ri_{f}<1$ ), where the model again reasonably reproduces LES trends, but the relative errors are larger than for the strongly unstable conditions. The agreement with DNS is significantly better. These DNS were of rotating (Ekman) flows where the two horizontal-variance components are difficult to disentangle due to shear production in the cross-stream direction. Here, they are lumped and compared to $R_{w}=-R_{u}-R_{v}$ , which are in better agreement relative to the match with LES results. It is not possible to ascertain whether the discrepancies between the model and LES are due to model simplifications and/or unavoidable errors in LES computations of the redistribution. Nonetheless, the model captures the leading-order physics of variance redistribution without resorting to tuneable coefficients, and can be used to further probe their impact on the structural characteristics of eddies and on turbulence regimes.

Figure 1. Variation of the variance budgets terms with stability: (a) for strongly unstable conditions, (b) for moderately unstable and stable conditions. Lines are from the reduced model; markers are determined from LES or DNS.

Focusing on the buoyancy influence on the budget terms in figure 1, the model and simulations both indicate apparent transitions in turbulence regimes at approximate $Ri_{f}$ values, although with some discrepancy between the two as to the exact $Ri_{f}$ at which a transition occurs. These regimes will henceforth be discussed in terms of the transition $Ri_{f}$ values indicated by the model, but the reader must bear in mind the potential uncertainty in these model-derived values. Moreover, in other flow set-ups with shear and buoyancy such as a jet effluent into a stable environment or flow over inclined walls, the same turbulence regimes might emerge but the transition $Ri_{f}$ values will likely be different from horizontal wall-bounded flows.

For the current configuration under unstable conditions, as $Ri_{f}$ decreases from its zero neutral value where energy is redistributed from the $u$ component into the $v$ and $w$ components, the following can be noted:

Figure 2. Standard deviation of the perturbation cross-stream vorticity, $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}_{y}^{\prime }}$ , heat transport efficiency, $\unicode[STIX]{x1D702}_{h}$ , and momentum transport efficiency, $\unicode[STIX]{x1D702}_{m}$ , versus $Ri_{f}$ . Vertical lines are at $Ri_{f}=-0.5$ (dashed), $-1$ (solid) and $-2$ (dot-dashed). Experiments and analysis methods are described in Li & Bou-Zeid (Reference Li and Bou-Zeid2011). Purple rectangles depict the directionality of energy transfer.

(i) At $Ri_{f}<-0.5$ , $R_{w}$ becomes negative, indicating a redistribution out of the vertical component $\overline{w^{\prime 2}}$ . The streamwise component, however, remains a source and the cross-stream is now the only net recipient of energy. Furthermore, in the $\overline{u^{\prime 2}}$ budget, dissipation overtakes redistribution as the main sink. The regime transition indicated by the model at $Ri_{f}=-0.5$ can in fact be observed in very high Reynolds number atmospheric flows. One illustration is provided in figure 2, where the standard deviation of the cross-stream perturbation vorticity $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}_{y}^{\prime }}$ (where $\unicode[STIX]{x1D714}_{y}^{\prime }=\unicode[STIX]{x2202}u^{\prime }/\unicode[STIX]{x2202}z-\unicode[STIX]{x2202}w^{\prime }/\unicode[STIX]{x2202}x$ ) and the heat $\unicode[STIX]{x1D702}_{h}$ and momentum $\unicode[STIX]{x1D702}_{m}$ transfer efficiencies measured in the atmosphere above a lake are featured. The efficiencies are defined as the net flux $F$ of heat or momentum divided by the total downgradient flux from ejections and sweeps $\unicode[STIX]{x1D702}=F_{total}/F_{downgradient}$ , as detailed in Li & Bou-Zeid (Reference Li and Bou-Zeid2011). At $Ri_{f}\approx -0.5$ , the momentum transfer efficiency begins a rapid decline as the now buoyantly dominated vertical velocity starts to decorrelate from the horizontal component and vertical motions are no longer simply responding to horizontal ones. Furthermore, the decline in $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}_{y}^{\prime }}$ , which encodes information about the structural characteristics of turbulent eddies, and the increase in $\unicode[STIX]{x1D702}_{h}$ slow down. Thus, not only the turbulence energy, but also its structure is modified by buoyancy below this $Ri_{f}$ as the vertical component is now dominated by buoyancy generation.

(ii) At $Ri_{f}<-1$ , $|R_{w}|>|R_{u}|$ and the $w$ component becomes the main supplier of energy in the redistribution interplay but $R_{u}$ remains a source as well. The directionality of redistribution (indicated by arrows in figure 2) does not change at this $Ri_{f}$ . However, turbulence eddy structure and statistics are altered as can be deduced from $\unicode[STIX]{x1D702}_{h}$ , which plateaus beyond this value. This plateau suggests that buoyancy fully dominates the generation of vertical velocity perturbations and their correlation with temperature perturbations, with no effect on heat transport from shear beyond this point.

(iii) At $Ri_{f}<-2$ , there is a transition to $R_{u}>0$ implying that $\overline{u^{\prime 2}}$ now becomes a net recipient of energy originating from $\overline{w^{\prime 2}}$ . The flow can be considered buoyancy-dominated beyond this point since the streamwise fluctuations are now modulated by buoyantly generated energy. At this value of $Ri_{f}=-2$ , one can note the onset of a new turbulence regime in figure 2 as $\unicode[STIX]{x1D702}_{m}$ , $\unicode[STIX]{x1D702}_{h}$ and $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D714}_{y}^{\prime }}$ all reach a plateau.

On the stable side, as $Ri_{f}$ increases from zero:

(iv) At $Ri_{f}>0.25$ , buoyancy destruction exceeds viscous dissipation in the budget of $\overline{w^{\prime 2}}$ ( $\unicode[STIX]{x1D700}_{w}/S<-B/S=Ri_{f}$ ), and beyond this point it dominates the damping and the structure of vertical fluctuations.

(v) At $Ri_{f}>0.5$ , buoyancy destruction exceeds the total TKE viscous dissipation in all components $\unicode[STIX]{x1D700}$ .

(vi) The $Ri_{f}=1$ limit is forbidden. The dissipation and the redistribution to the $v$ -component become zero requiring complete laminarization of the flow, while the model predicts $R_{u}=R_{w}$ , requiring finite turbulent energy. This point is further examined below to show that a lower limit of $Ri_{f}\approx 0.21$ exists and that turbulence can be sustained at high stabilities ( $Ri_{g}$ values) provided $Re$ remains high.

4 Extended model with linear Rotta closure

To proceed further in variance modelling, a linear Rotta-type closure for the redistribution term is introduced and is given by

(4.1a-c ) $$\begin{eqnarray}R_{u}=-\frac{c\unicode[STIX]{x1D700}}{k}\left(\overline{u^{\prime 2}}-\frac{2}{3}k\right)\quad R_{v}=-\frac{c\unicode[STIX]{x1D700}}{k}\left(\overline{v^{\prime 2}}-\frac{2}{3}k\right)\quad R_{w}=-\frac{c\unicode[STIX]{x1D700}}{k}\left(\overline{w^{\prime 2}}-\frac{2}{3}k\right),\end{eqnarray}$$

where $k=0.5(\overline{u^{\prime 2}}+\overline{v^{\prime 2}}+\overline{w^{\prime 2}})$ is the TKE and $c$ is a model constant. The ratio $k/\unicode[STIX]{x1D700}$ represents a characteristic decay time of the turbulence, which here encodes the time required to redistribute variance between components. This model captures the ‘slow’ part of the redistribution and represents the simplest closure one can propose for the return-to-isotropy terms. More detailed models that include the effect of shear and wall blockage, buoyancy and potentially other physical processes, on the pressure–strain correlations have been proposed and used (Canuto et al. Reference Canuto, Howard, Cheng and Dubovikov2001; Heinze, Mironov & Raasch Reference Heinze, Mironov and Raasch2016). However, to assess the limits of the simplest plausible formulation (and to avoid including other budget equations in the model), we will focus on the linear Rotta closure in this study.

4.1 Model prediction for velocity variances

Inserting the expressions of (4.1) into (2.7)–(2.9), and using $S/\unicode[STIX]{x1D700}=(1-Rif)^{-1}$ , yields

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\overline{u^{\prime 2}}}{k}=\frac{S}{c\unicode[STIX]{x1D700}}\left(\frac{2}{3}+\frac{1}{3}Ri_{f}\right)+\frac{2}{3}=\frac{1}{3c}\left(\frac{2+Ri_{f}}{1-Ri_{f}}\right)+\frac{2}{3}, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\overline{v^{\prime 2}}}{k}=\frac{S}{c\unicode[STIX]{x1D700}}\left(-\frac{1}{3}+\frac{1}{3}Ri_{f}\right)+\frac{2}{3}=\frac{1}{3c}\left(\frac{-1+Ri_{f}}{1-Ri_{f}}\right)+\frac{2}{3}=-\frac{1}{3c}+\frac{2}{3}, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\overline{w^{\prime 2}}}{k}=\frac{S}{c\unicode[STIX]{x1D700}}\left(-\frac{1}{3}-\frac{2}{3}Ri_{f}\right)+\frac{2}{3}=\frac{1}{3c}\left(\frac{-1-2Ri_{f}}{1-Ri_{f}}\right)+\frac{2}{3}. & \displaystyle\end{eqnarray}$$

These model variances are shown in figure 3 using $c=0.9$ (solid lines), which has been reported as optimal in closure modelling (Pope Reference Pope2000). The $c=0.9$ (or 1.8 in the full-variance budget) has some theoretical justification (Katul et al. Reference Katul, Porporato, Manes and Meneveau2013) for neutral flows as well. A notable prediction of the model is that the fraction $\overline{v^{\prime 2}}/k$ remains independent of stability at a value of approximately 0.3, and this is indeed supported by the LES results (not plotted here) where $\overline{v^{\prime 2}}/k\approx 0.25{-}0.35$ at heights $z>\unicode[STIX]{x1D6FF}/4$ , where $\unicode[STIX]{x1D6FF}$ is the total flow depth. The ratio is higher closer to the wall, but the dependence on stability in the LES runs is weak and non-monotonic.

Figure 3. Variance to TKE ratios predicted by the model as a function of $Ri_{f}$ : dot–dash lines with $c=0.7$ , solid lines with $c=0.9$ and dashed line with $c=1.1$ .

Under neutral conditions where $Ri_{f}=0$ and with $c=0.9$ , equations (4.2)–(4.4) yield $\overline{u^{\prime 2}}/k=1.41$ and $\overline{v^{\prime 2}}/k=\overline{w^{\prime 2}}/k=0.3$ , which are reasonably close to values observed in measurements and numerical simulations in the logarithmic region away from the wall (e.g. Bou-Zeid et al. (Reference Bou-Zeid, Meneveau and Parlange2005)), although as discussed before $\overline{v^{\prime 2}}/k>\overline{w^{\prime 2}}/k$ . On the other hand, under very strongly unstable conditions when $-Ri_{f}\gg 1$ and production is dominated by buoyancy, $\overline{u^{\prime 2}}$ and $\overline{w^{\prime 2}}$ reverse roles, and $\overline{w^{\prime 2}}/k=1.41$ while $\overline{u^{\prime 2}}/k=\overline{v^{\prime 2}}/k=0.3$ as can be seen in figure 3(a) at the free convection limit of $Ri_{f}\ll 0$ . The anisotropy of the variance ratios asymptotes to $\overline{w^{\prime 2}}/\overline{v^{\prime 2}}\approx \overline{w^{\prime 2}}/\overline{u^{\prime 2}}\approx 4.6$ . These values are close to the ${\approx}3.2$ reported for DNS of free convection above a heated plate by Mellado (Reference Mellado2012), and to the ${\approx}3.4$ of Salesky et al. (Reference Salesky, Chamecki and Bou-Zeid2017) at the lowest $Ri_{f}$ they studied using LES. The better predictions of the model in the free convection limit, compared to the neutral limit, are not surprising since the energy is being injected in the vertical component and the redistributions to the horizontal components must become equal. This is indicated by the reduced model and the LES in figure 1 and is expected since physically these directions are exactly similar in that limit: the wall effect does not hinder energy redistribution into either horizontal receiving components.

Under stable conditions, a maximum attainable flux Richardson number ( $Ri_{f,max}$ ) emerges beyond which $\overline{w^{\prime 2}}<0$ (figure 3 b). This $Ri_{f,max}$ ( $=$  0.21 for $c=0.9$ ) cannot be exceeded because when it is approached, $\overline{w^{\prime 2}}$ approaches zero and buoyancy destruction ( ${\sim}\overline{w^{\prime }\unicode[STIX]{x1D703}^{\prime }}$ ) will be throttled by the rate of redistribution $R_{w}$ , maintaining $Ri<Ri_{f,max}$ . That is, the constraint of $\overline{w^{\prime 2}}>0$ limits the maximum attainable $R_{w}$ according to (4.1) (it will be the $R_{w}$ attained when $\overline{w^{\prime 2}}\rightarrow 0$ ). This maximum $R_{w}$ then limits the highest attainable energy flux into the vertical component, which must simultaneously limit the buoyant energy destruction in that component (via budget (2.3)). Thus, pressure redistribution of variance from the streamwise to the vertical velocity component becomes the primary regulator of turbulence regime at this asymptotic $Ri_{f,max}$ . This limiting value of ${\approx}\,0.21$ is quite consistent with previous estimates (ranging from 0.19 to 0.21) based on second-order moments models that make comparable assumptions to the ones we make here (Yamada Reference Yamada1975; Mellor & Yamada Reference Mellor and Yamada1982), as well as with the asymptotic limit of $Ri_{f}$ that can be obtained from its relation to the empirically determined Monin–Obukhov gradient stability correction functions (Wyngaard Reference Wyngaard2010, p. 281). The exact limiting value however might not be 0.21 if buoyancy or other factors (waves, wall-blockage, laminarization, intermittency, non-stationarity, etc.) alter the redistribution efficiency and consequently the value of $c$ . To account for this possibility and to investigate the change in variances with uncertainties in $c$ , figure 3(b) presents the profiles with $c$ values of 0.7 and 1.1. The attainable $Ri_{f,max}$ limit is reduced to approximately 0.12 for $c=0.7$ and increased to about 0.28 when $c=1.1$ . From (4.4), one can derive (again by constraining $\overline{w^{\prime 2}}$ to remain ${>}0$ ) the following expression for the dependence on $c$ of this asymptotic flux Richardson number as $Ri_{f,max}=(c-0.5)/(c+1)$ , which increases with increasing $c$ . That is, the more efficient the redistribution (higher $c$ ), the higher the attainable $Ri_{f,max}$ since a high redistribution can sustain a higher buoyancy destruction rate. The limit of $Ri_{f,max}\approx 0.25$ , consistent with the classic estimate from linear stability analysis (Miles & Howard Reference Miles and Howard1964) requires $c=1$ , although we must point out that their result was derived for a bulk $Ri$ .

Figure 4. Flux versus gradient Richardson numbers from laboratory and atmospheric field experiments, data from Katul et al. (Reference Katul, Porporato, Shah and Bou-Zeid2014).

Atmospheric boundary layer field experiments corroborate the model findings: values of $Ri_{f}>\approx 0.21$ are atypical. Gradient and flux Richardson numbers are about equal and vary together between $0<Ri_{f}\approx Ri_{g}<0.21$ ; however, beyond this point, $Ri_{f}$ reaches a plateau at $Ri_{f,max}\approx 0.21$ while $Ri_{g}$ continues to increase (Katul et al. Reference Katul, Porporato, Shah and Bou-Zeid2014; van Hooijdonk et al. Reference van Hooijdonk, Clercx, Ansorge, Moene and van de Wiel2018). This is illustrated in figure 4 where the demarcation between the buoyancy-destruction limited and the redistribution-limited regimes is indicated. Thus, even at higher bulk stabilities, the redistribution throttle mitigates the continuous increase in buoyant destruction and the collapse of turbulence (until the Dougherty–Ozmidov length scale approaches $\unicode[STIX]{x1D702}_{K}$ (Katul et al. Reference Katul, Porporato, Shah and Bou-Zeid2014)). A further evidence of the physical significance of this $Ri_{f}$ limit is that the spectrum of measured vertical velocity above an extensive ice sheet no longer exhibits an inertial subrange when $Ri_{f}$ exceeds a threshold estimated at ${\approx}0.25$ (Grachev et al. Reference Grachev, Andreas, Fairall, Guest and Persson2013).

Another conclusion that can be drawn from the analysis of the relation between $Ri_{f,max}$ and $c$ is that the value of $c$ must be constrained to $c>0.5$ (otherwise $Ri_{f,max}$ in stable flows cannot be ${>}0$ ) and $c<2.5$ (otherwise $Ri_{f,max}>1$ ). As discussed above, $Ri_{f,max}>1$ cannot be attained when a high- $Re$ flow is stationary, planar–homogeneous and lacking a mean vertical velocity. For the aforementioned conditions, $Ri_{f,max}>1$ results in a negative TKE dissipation rate (i.e. viscous production).

4.2 Model prediction for velocity anisotropy, and the implied Rotta constant

As detailed in the introduction, return to isotropy is a trend, but the turbulence remains anisotropic if production is dominated by one or two components. One important question that emerges is to what degree can return to isotropy succeed in equi-partitioning the energy, and how does buoyancy influence this equi-partitioning. The flow anisotropy is best characterized by the velocity anisotropy tensor

(4.5) $$\begin{eqnarray}\unicode[STIX]{x1D623}_{ij}=\frac{{u^{\prime }}_{i}{u^{\prime }}_{j}}{2k}-\frac{\unicode[STIX]{x1D6FF}_{ij}}{3}.\end{eqnarray}$$

Based on (4.2)–(4.4), the diagonal components of this tensor can be modelled as

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D623}_{11}=\frac{1}{6c}\left(\frac{2+Ri_{f}}{1-Ri_{f}}\right), & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D623}_{22}=-\frac{1}{6c}, & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D623}_{33}=\frac{1}{6c}\left(\frac{-1-2Ri_{f}}{1-Ri_{f}}\right), & \displaystyle\end{eqnarray}$$

where the index 1 denotes the streamwise direction, 2 the cross-stream and 3 the vertical. Under neutral conditions and with $c=0.9$ , this results in $\unicode[STIX]{x1D623}_{11}=(3c)^{-1}=0.37$ and $\unicode[STIX]{x1D623}_{22}=\unicode[STIX]{x1D623}_{33}=(-6c)^{-1}=-0.185$ (they sum to zero). These parameters are reported widely in the literature on neutral homogeneous shear flows; they depend on the Reynolds number as well as the shear rate (Pumir Reference Pumir1996; Briard et al. Reference Briard, Gomez, Mons and Sagaut2016). The values obtained here are within the reported ranges; specifically, at the highest simulated shear rate, Isaza & Collins (Reference Isaza and Collins2009) found $\unicode[STIX]{x1D623}_{11}\approx 0.38$ , $\unicode[STIX]{x1D623}_{22}\approx -0.12$ and $\unicode[STIX]{x1D623}_{33}\approx -0.26$ (notice the 22 and 33 notations in their paper are the inverse of ours). It is to be noted that their simulations were of homogeneous neutral free shear and hence there were no wall or buoyancy; therefore, the asymmetry between the cross-steam and vertical components is due to the effect of shear on the vertical component. The effect of shear is in fact broadly similar to the effect of a wall: it would reduce redistribution to the shear-normal vertical component and result in $\unicode[STIX]{x1D623}_{33}<\unicode[STIX]{x1D623}_{22}$ (both negative) or equivalently $\overline{w^{\prime 2}}<\overline{v^{\prime 2}}$ . That is, the presence of a wall or the influence of shear alter the redistribution and require a modification of the simplest linear Rotta model, and this modification due to shear is in fact often included as the second term in redistribution closure schemes (e.g. Canuto et al. Reference Canuto, Howard, Cheng and Dubovikov2001).

Some studies (Gerz, Schumann & Elghobashi Reference Gerz, Schumann and Elghobashi1989; Holt, Koseff & Ferziger Reference Holt, Koseff and Ferziger1992) also examined $\unicode[STIX]{x1D623}_{11}$ , $\unicode[STIX]{x1D623}_{22}$ and $\unicode[STIX]{x1D623}_{33}$ under non-neutral conditions. The direct measurement from simulations can be compared to predictions of the present model but the influence of buoyancy on redistribution and on the value of the Rotta constant needs to be accounted for example, as in Canuto et al. (Reference Canuto, Howard, Cheng and Dubovikov2001). An alternative approach is in fact to use the present model of the components of the velocity anisotropy tensor (or the variances or variance differences that can also be derived from (4.2) to (4.4)) to estimate the Rotta model constant and how it varies with stability. Such estimates of $c$ , while possible in principle, proved difficult. They are highly sensitive to the accuracy of the computed or measured variances near $Ri_{f}$ values that would cause the numerators or denominators of (4.6)–(4.8) to approach 0. Computations from DNS for example suggest a reasonable value of $c$ that varies between 0.8 and 1.3, but with no clear trend with increasing stability. As such, while this paper establishes the influence of redistribution of energy on buoyancy modulation of turbulent flows, the influence of buoyancy on redistribution is a topic that requires further investigation.