5.1. Mean profiles with stability

Log-linear plots similar to the one in figure 2, but with positive Richardson numbers, are compared to the neutral case in figure 5(a,c), which correspond to $\mathit{Re}_{f}$ of 600 and 900, respectively. A correction to the log law in stably stratified cases was given by Monin & Obukhov (Reference Monin and Obukhov1954). The MOST, which theoretically holds for infinite Reynolds number only, postulates that the non-dimensional velocity and temperature gradients ( ${\it\phi}_{m}$ and ${\it\phi}_{h}$ , respectively) for non-neutral cases are functions of the height normalized by the Obukhov length $L$ . These functions have been experimentally shown to be linear in $z/L$ (for $0\leqslant z/L<1$ ) in stable conditions (Brutsaert Reference Brutsaert1982; Stull Reference Stull1988; see also theoretical justification in Li, Katul & Bou-Zeid (Reference Li, Katul and Bou-Zeid2012)):

(5.1) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}_{m}\biggl(\frac{z}{L}\biggr)=\frac{{\it\kappa}z}{u_{\ast }}\left[\left(\frac{\partial \langle u\rangle }{\partial z}\right)^{2}+\left(\frac{\partial \langle v\rangle }{\partial z}\right)^{2}\right]^{1/2}=1+{\it\beta}_{m}\frac{z}{L}, & \displaystyle\end{eqnarray}$$

(5.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}_{h}\biggl(\frac{z}{L}\biggr)=\frac{{\it\kappa}z}{{\it\theta}_{\ast }}\frac{\partial \langle {\it\theta}\rangle }{\partial z}={\it\alpha}_{h}+{\it\beta}_{h}\frac{z}{L}, & \displaystyle\end{eqnarray}$$

${\it\beta}_{m}$

${\it\alpha}_{h}$

${\it\beta}_{h}$

$z/L$

$L$

${\sim}L$

$L$

$H_{k,s}=0$

${\it\phi}_{m}$

The mean velocity gradients simulated in our statically stable DNS follow MOST despite the finite and relatively moderate $\mathit{Re}$ , as can be seen from figure 5(b,d). The value of ${\it\beta}_{m}$ at the Reynolds numbers simulated here ( ${\approx}9$ ), however, is higher than typical values observed in atmospheric flows ( ${\approx}5$ ) with higher Reynolds numbers for $z/L<1$ (Webb Reference Webb1970; Dyer Reference Dyer1974; Kondo, Kanechika & Yasuda Reference Kondo, Kanechika and Yasuda1978; Stull Reference Stull1988). This suggests that these MOST coefficients are Reynolds-number-dependent, i.e. viscosity has an effect on ${\it\beta}_{m}$ . This effect was in fact analysed by Chung & Matheou (Reference Chung and Matheou2012). By considering a stable regime where shear stress is mainly carried by eddies of sizes ranging between the Kolmogorov scale ( ${\it\eta}$ ) and $L$ , they got an expression of the form

(5.3) $$\begin{eqnarray}1\approx \frac{3}{4}\frac{{\it\gamma}_{1}}{(2{\rm\pi}{\it\kappa})^{4/3}}{\it\beta}_{m}({\it\beta}_{m}-1)^{1/3}(1-({\it\eta}/L)^{4/3}),\end{eqnarray}$$

where ${\it\gamma}_{1}\approx 0.15$ is the constant in the $u{-}w$ cospectra model for the inertial range (Lumley Reference Lumley1967; Wyngaard & Coté Reference Wyngaard and Coté1972; Chung & Matheou Reference Chung and Matheou2012). The expression above illustrates that the value of ${\it\beta}_{m}$ has to increase as the term $1-({\it\eta}/L)^{4/3}$ decreases. But since $L/{\it\eta}$ is the ratio of the largest scale $L$ to the smallest scale ${\it\eta}$ , it is proportional to $\mathit{Re}^{3/4}$ . As such, as this Reynolds number decreases, $1-({\it\eta}/L)^{4/3}$ also decreases and ${\it\beta}_{m}$ increases. This increase in ${\it\beta}_{m}$ leads to a reduction in the critical flux Richardson number ( $\mathit{Ri}_{f}=z/(L{\it\phi}_{m})$ ; Chung & Matheou Reference Chung and Matheou2012) at lower Reynolds numbers, which is indeed what we observed in our simulations as noted above. Similar arguments can be applied to suggest that ${\it\beta}_{h}$ also increases when this Reynolds number is reduced.

The profiles from the lowest Reynolds-number simulations carried out ( $\mathit{Re}_{f}=400$ ) are not shown in figure 5, since a distinct log region is not observed in the neutral case owing to the proximity of the viscous and outer regions of the flow (as can be seen in figure 2). For such flows, the separation between the production range and dissipation range eddies is not sufficient to allow the existence of an inertial range, which is associated with the observed shapes of the log region or MOST region mean profiles (Gioia et al. Reference Gioia, Guttenberg, Goldenfeld and Chakraborty2010; Katul, Koning & Porporato Reference Katul, Koning and Porporato2012; Li et al. Reference Li, Katul and Bou-Zeid2012).

A last point to underline here is that the plots of the non-dimensional gradients use the surface fluxes to define the Obukhov length. Similar plots (not shown here) with local values of the fluxes (Nieuwstadt Reference Nieuwstadt1984), and hence a local Obukhov length, at each height and friction velocity yield closer constants ( ${\it\beta}_{m}\approx 4$ and ${\it\beta}_{h}\approx 3.75$ ) to the ones measured and used in the ABL. This is relevant since it illustrates that the fluxes are not exactly constant in the MOST layer, and meteorologists often use fluxes measured locally using eddy-covariance instruments as proxies for analysis. Results from the present DNS (not shown here) on the variability of the total fluxes ( $\text{viscous}+\text{turbulent}$ ) with height for neutral and the most weakly stable cases reveal that momentum fluxes drop by approximately 20 % from their surface values within the first 5 % of the boundary layer (which is in fact remarkably close to the estimate of that drop associated with the pressure gradient term that Monin & Obukhov (Reference Monin and Obukhov1954) presented in their original work), while buoyancy fluxes drop by 20 % of their surface value at approximately 10 % of the boundary layer depth. The constant flux assumption in MOST comes from simplifications to the momentum and thermal energy conservation equations, which are not affected by the moderate Reynolds numbers if the viscous or molecular fluxes are included as in our analysis. As such, its violation is not related to the low Reynolds number of the simulations, and should be indicative of the drops expected in the ABL. The DNS results also confirm that the drop in fluxes with height is insensitive to the Reynolds number of the simulation for momentum flux and very weakly sensitive for buoyancy flux. However, increasing stability increases the rate of the flux decrease with height, particularly for the momentum flux, compared to neutral simulations.

Figure 6(a–c) depicts the separate variations in the mean profiles of $u$ and $v$ as stability increases for $\mathit{Re}_{f}=400$ , 600 and 900, respectively, while figure 6(d) shows the variation of these velocity components with time during one inertial period $2{\rm\pi}/f$ ( $=12\ \text{h}$ ) for the case $\mathit{Re}_{f}=600$ at $\mathit{Ri}_{0,0}=0.002$ . Increasing stability changes the Ekman rotation angle at a given height. It also reduces flow velocities very close to the wall (see inset in figure 6 c), while increasing them away from the wall (see figure 5). This is consistent with the reduced mixing and momentum transport from the upper levels of the boundary layer, where the flow becomes faster, to the near-wall region, where the flow becomes slower, as stability increases. The stable profiles are hence more ‘laminar-like’.

Figure 6.

(a–c) Reynolds-averaged non-dimensional velocity profiles at Reynolds number ( $\mathit{Re}_{f}$ ) of (a) 400, (b) 600 and (c) 900 at various stabilities ( $\mathit{Ri}_{0,0}=0.001$ for $\mathit{Re}_{f}=400$ ; 0.001, 0.002 for $\mathit{Re}_{f}=600$ ; and 0.001, 0.002 and 0.005 for $\mathit{Re}_{f}=900$ ) compared with the neutral case (black solid line in each panel). Corresponding bulk Richardson numbers $\mathit{Ri}_{b}$ can be found in table 1 and are shown in the legend of each panel. The inset in panel (c) shows a zoomed-in plot of normalized $\langle u\rangle$ near the wall. (d) The time variation of the velocity profiles, averaged over horizontal directions, at $\mathit{Re}_{f}=600$ and $\mathit{Ri}_{0,0}=0.002$ equivalent to $\mathit{Ri}_{b}=0.2$ . Here  $tf\approx 0$ (solid black lines) shows neutral-case profiles.

The neutral-case turbulent Ekman layer depth can be defined as ${\it\delta}_{t}=u_{\ast }/f$ (Caldwell et al. Reference Caldwell, Van Atta and Helland1972; Coleman et al. Reference Coleman, Ferziger and Spalart1990; Zilitinkevich & Mironov Reference Zilitinkevich and Mironov1996). As illustrated in table 1, the dimensional friction velocity $u_{\ast }$ and, correspondingly, the Ekman layer depth increase with increasing $\mathit{Re}$ , and decrease with increasing $\mathit{Ri}$ . As the Obukhov length and vertical mixing decrease for the more stable simulations, a shallower boundary layer is produced – a common feature of stable ABLs (Kaimal & Finnigan Reference Kaimal and Finnigan1994). The results obtained here are consistent with this trend. A good indicator of the boundary layer depth is the height at which there is a maximum in the velocity profile; this depth decreases with increasing stability as shown in figure 6. The bands of super-geostrophic velocities around the peaks observed in figures 5 and 6, where the flow speed in the boundary layer can be higher than the geostrophic value far from the wall, are also a common feature seen in stable ABLs (Stull Reference Stull1988).

5.2. Evolution of turbulence and stability

Figures 7(a), 8(a) and 9(a) show the time evolution of the surface Richardson number ( $\mathit{Ri}_{0,t}=g{\it\Gamma}_{0,t}D^{2}/{\it\theta}_{\infty }G^{2}$ ) at various bulk Richardson numbers $\mathit{Ri}_{b}=g{\it\delta}_{t}{\rm\Delta}{\it\theta}/{\it\theta}_{0}G^{2}$ , where ${\rm\Delta}{\it\theta}={\it\theta}_{\infty }-{\it\theta}_{0}$ . The bulk Richardson number $\mathit{Ri}_{b}$ remains almost constant throughout a simulation and we thus use it to characterize the bulk stability. From the plots, it can be seen that the gradient at the surface initially increases, since the initial temperature profiles are not in equilibrium with the turbulent velocity fields; therefore, the gradients at the surface need time to adjust. Subsequently, $\mathit{Ri}$ falls back and stabilizes. The cross marks on the plots show the period over which the statistics were averaged (the ones overlaid on the $x$ axis are for the neutral cases). Coleman et al. (Reference Coleman, Ferziger and Spalart1992) argue that the first $tf=1$ is enough for the development of temperature fluctuations in the field, which would be in equilibrium with the velocity field. However, in our analysis, the time before the first cross mark at $tf={\rm\pi}$ is not used for averaging, since during this period the changes in turbulence statistics due to stability seem to remain significant in our simulations. The time interval after $tf={\rm\pi}$ has a more or less constant surface Richardson number, confirming that the flow is in quasi-equilibrium and justifying the averaging we perform over this half inertial time period. For the stronger stability simulations at each value of $\mathit{Re}$ , intermittency is observed whereby turbulence is damped and then regenerated as discussed previously. This behaviour has also been observed by García-Villalba & Del Álamo (Reference García-Villalba and Del Álamo2011) in their channel flow simulations, and might be indicative of insufficiently large computational domains for the stronger stratification cases. While these intermittent simulations warrant more in-depth analysis to validate their behaviour and potentially analyse it, in this paper, only data from the continuously turbulent simulations will henceforth be analysed.

Figure 7.

Results at $\mathit{Re}_{f}=400$ . (a) Surface Richardson number ( $\mathit{Ri}_{0,t}$ ) versus non-dimensional time ( $tf$ ) at various stabilities indicated by the legend. Cross marks indicate the start and end times of the averaging period used for calculating means and the averages of budgets for neutral and weakly stable cases. (b) Non-dimensional volume-integrated TKE per unit area ( $E$ ). (c) Friction velocity ( $u_{\ast }$ ) and its variation with time and stability. (d) Variation of the angle ( ${\it\beta}$ ) between the geostrophic wind direction and the surface shear stress direction with time and stability.

Figure 8.

Results at $\mathit{Re}_{f}=600$ . (a) Surface Richardson number ( $\mathit{Ri}_{0,t}$ ) versus non-dimensional time ( $tf$ ) at various stabilities indicated by the legend. Cross marks indicate the start and end times of the averaging period used for calculating means and the averages of budgets for neutral and weakly stable cases. (b) Non-dimensional volume-integrated TKE per unit area ( $E$ ). (c) Friction velocity ( $u_{\ast }$ ) and its variation with time and stability. (d) Variation of the angle ( ${\it\beta}$ ) between the geostrophic wind direction and the surface shear stress direction with time and stability.

Figure 9.

Results at $\mathit{Re}_{f}=900$ . (a) Surface Richardson number ( $\mathit{Ri}_{0,t}$ ) versus non-dimensional time ( $tf$ ) at various stabilities indicated by the legend. Cross marks indicate the start and end times of the averaging period used for calculating means and the averages of budgets for neutral and weakly stable cases. For the highest stable stratification imposed, the period denoted by cross marks in panel (a) for $\mathit{Ri}_{b}=1.936$ shows the intermittency in that time period and hence has not been used for analysis. (b) Non-dimensional volume-integrated TKE per unit area ( $E$ ). (c) Friction velocity ( $u_{\ast }$ ) and its variation with time and stability. (d) Variation of the angle ( ${\it\beta}$ ) between the geostrophic wind direction and the surface shear stress direction with time and stability.

In figures 7(b), 8(b) and 9(b), we plot the variation of volume-integrated TKE per unit area ( $E$ ), non-dimensionalized by the laminar Ekman layer depth ( $D$ ) and the square of the geostrophic wind velocity ( $G^{2}$ ), versus time, for various stabilities. The volume-integrated TKE (per unit area) averaged over a horizontal plane is defined as

(5.4) $$\begin{eqnarray}E=\frac{1}{2}\int _{0}^{L_{z}}\langle u_{i}^{\prime }u_{i}^{\prime }\rangle _{x,y}\text{d}z.\end{eqnarray}$$

Here $\langle \;\rangle _{x,y}$ indicates averaging over $xy$ planes, the prime (  $^{\prime }$  ) indicates the fluctuating component and $L_{z}$ is the top of the domain. The TKE drops down with increasing stability (the mechanisms behind this drop are discussed later in the paper), with most of the drop occurring when the flow transitions from neutral to weakly stable. For Reynolds numbers of 600 and 900, the total TKE in the domain for the least stable cases is about half of that for the neutral simulation. For the highest Richardson numbers imposed at every value of the Reynolds number, the TKE in the domain drops almost to zero, but with strong intermittent bursts. The bulk Richardson number at which this intermittency sets in varies with $\mathit{Re}_{f}$ (recall the discussion of the variation of the critical Richardson number with $\mathit{Re}$ above). We also observe that the stronger the imposed stability, the faster the rate at which TKE drops at initial times. The decay rates of the TKE near the wall are consistent with a scaling with $L/u_{\ast }$ , as reported by Flores & Riley (Reference Flores and Riley2011). Note that all the simulations have the same initial state (fully turbulent neutral Ekman layer), and hence this comparison is meaningful.

Figures 7(c), 8(c) and 9(c) show the evolution of the friction velocity (non-dimensionalized by the geostrophic wind velocity $G$ ; $u_{\ast }/G$ is referred to in the meteorological literature as the geostrophic drag coefficient) versus time. Static stability reduces the friction velocity (and thus the surface momentum flux). The initial rate at which the friction velocity drops from the neutral-case initial value also increases with stability. At large stabilities, there are bursts in the friction velocity too, associated with the turbulent bursts, which increase $u_{\ast }$ to values comparable to those seen in the neutral cases. One should note here that, for the neutral and weakly stable simulations, $u_{\ast }/G\approx 0.05$ , which is the value used by Monin & Obukhov (Reference Monin and Obukhov1954) to estimate the rate of decrease in fluxes with height discussed above.

Figures 7(d), 8(d) and 9(d) show the variation of the angle between the geostrophic wind velocity vector and the shear stress vector at the surface versus time, for the various stabilities. Note that the angle for the turbulent neutral case is denoted by the black solid lines, while the value of the angle at the surface for the laminar limit is 45° (Stull Reference Stull1988). This limit can be obtained by solving the Reynolds-averaged equations of motion with the simplifying assumptions of steadiness, homogeneity in the horizontal directions, neutral stratification and a constant viscosity. With increasing stability, as the TKE drops and the flow becomes less turbulent, the angle at the surface increases and approaches the limit of 45° for the laminarizing intermittent cases.

The average angle, $\langle {\it\gamma}\rangle$ , between the surface shear stress (which is about aligned with the wind very close to the surface) and the wind direction with height is depicted in figure 10. In many stable boundary layer theories, it is assumed that the turning due to the Coriolis force is negligibly small near the surface (lowest 10 % of the boundary layer). From the plots, it is clear that, even if the direct effect of the Coriolis force is small on turbulence (as compared to other terms in the budgets analysed later), the angle between surface shear stress and wind direction can be significant, particularly at higher stabilities. However, $\langle {\it\gamma}\rangle$ decreases with increasing Reynolds number, and thus higher-Reynolds-number simulations will be needed to verify whether $\langle {\it\gamma}\rangle$ is indeed negligible in the lower ABL.

Figure 10.

Average angle between the shear stress at the surface and the wind velocity, $\langle {\it\gamma}\rangle =\tan ^{-1}(v/u)-{\it\beta}$ , as a function of height in the surface layer: (a) $\mathit{Re}_{f}=400$ , (b)  $\mathit{Re}_{f}=600$ and (c)  $\mathit{Re}_{f}=900$ at various stabilities indicated in the legend.

5.3. Profiles of the gradient Richardson number, TKE and stress

The vertical profiles of the gradient Richardson number are depicted in figure 11. The gradient Richardson number is defined as

(5.5) $$\begin{eqnarray}\mathit{Ri}_{gr}=\frac{g}{{\it\theta}_{\infty }}\frac{\text{d}\langle {\it\theta}\rangle /\text{d}z}{(\text{d}\langle u\rangle /\text{d}z)^{2}+(\text{d}\langle v\rangle /\text{d}z)^{2}}.\end{eqnarray}$$

As the Richardson number increases, buoyancy effects will result in complete laminarization of the flow. The exact value at which this happens, called the critical Richardson number ( $\mathit{Ri}_{c}$ , the value of $\mathit{Ri}_{gr}$ below which fluid flow becomes dynamically unstable), is determined by the state of the flow. It has been argued by Miles (Reference Miles1961) and Howard (Reference Howard1961), using linear stability analysis of two-dimensional steady stable shear flows, that a sufficient condition for dynamic stability is that the gradient Richardson number ( $\mathit{Ri}_{gr}$ ) exceed 0.25 everywhere in the flow (any small perturbation will be damped under such conditions). This value of 0.25 is hence often used as an estimate for $\mathit{Ri}_{c}$ . However, the limitations of the linear stability analysis and the possibility of sustaining turbulence if the local $\mathit{Ri}_{gr}$ is below 0.25 somewhere in the flow makes this estimate not very robust. For example, some studies (see Canuto Reference Canuto2002, and references therein) indicate that $\mathit{Ri}_{c}$ may be as large as one or beyond, and some argue that a critical limit does not exist at all (Monin & Yaglom Reference Monin and Yaglom1975; Yamamoto Reference Yamamoto1975; Lettau Reference Lettau1979; Galperin, Sukoriansky & Anderson Reference Galperin, Sukoriansky and Anderson2007). Other arguments suggest that laminarization would happen when the characteristic length of the large eddies (estimated as the Obukhov scale in Flores & Riley (Reference Flores and Riley2011) and the Ozmidov scale in Katul et al. (Reference Katul, Porporato, Shah and Bou-Zeid2014)), which is reduced by buoyancy, approaches some multiple of the Kolmogorov scale, the ratio of the two scales effectively becoming proportional to $\mathit{Re}^{3/4}$ . Chung & Matheou (Reference Chung and Matheou2012) also argue in their analysis of stationary homogeneous stratified turbulence that the Richardson number is not sufficient to determine the state of the flow; a measure of the scale separation (given by the Reynolds number) is also needed. This point of view is confirmed by our previous analysis in this paper where we illustrated that the Richardson number at which intermittency starts being observed depends on the Reynolds number.

Figure 11.

Gradient Richardson number as a function of height with varying stabilities at a Reynolds number $\mathit{Re}_{f}$ of (a) 400 and (b) 600. The normalization depth ${\it\delta}_{t}$ for all cases is the depth of the neutral boundary layer at the given Reynolds number.

From the profiles in figure 11, one can see that buoyancy starts influencing the flow significantly (large $\mathit{Ri}_{gr}$ exceeding some arbitrary threshold, e.g. 0.25) only at larger heights, but this critical height at which buoyancy becomes important decreases drastically with increasing $\mathit{Ri}_{b}$ owing to the decrease of the boundary layer depth. A very similar trend is observed in the plot with $\mathit{Re}_{f}=900$ (not shown here), indicating a decrease in the height at which the ‘critical’ Richardson number is reached. Note that the normalization of height ( $y$ axis here) is carried out using the depth (turbulent length scale) of the neutral simulations. Other scales can be used for normalization (Mironov & Fedorovich (Reference Mironov and Fedorovich2010), for example, pointed to several depth scales for stable boundary layers that are consistent with TKE and momentum budgets), but, to emphasize the effect of static stability on turbulence at a given dimensional distance from the wall, we use the fixed neutral-case turbulent Ekman layer length scale for normalization.

Figure 12 presents the vertical variation of $2q^{2}=\langle u_{i}^{\prime }u_{i}^{\prime }\rangle$ (twice the TKE) at various stabilities and at different Reynolds numbers. Increasing $\mathit{Re}$ brings the peak in $2q^{2}$ closer to the wall as expected, while increasing $\mathit{Ri}$ generally reduces the peak values of $2q^{2}$ ; for the highest stabilities, the peaks move slightly away from the wall.

Figure 12.

Profiles of the time-averaged non-dimensional sum of the variances of the three velocity components (i.e. $2q^{2}$ , where $q^{2}$ is TKE) at Reynolds number ( $\mathit{Re}_{f}$ ) of (a) 400 (b) 600 and (c) 900 at various stabilities, compared with the corresponding neutral cases (black solid line in each panel). The bulk Richardson numbers $\mathit{Ri}_{b}$ given in the legend correspond, from lowest to highest, to the following initial surface Richardson numbers: $\mathit{Ri}_{0,0}=0.001$ , 0.002, 0.005.

Figure 13 depicts the vertical profiles of the magnitude of the kinematic stress, for various Reynolds numbers and stabilities. The magnitude of this kinematic turbulent stress $|{\it\tau}/{\it\rho}|$ is defined as follows:

(5.6) $$\begin{eqnarray}|\frac{{\it\tau}}{{\it\rho}}|=(\langle u^{\prime }w^{\prime }\rangle ^{2}+\langle v^{\prime }w^{\prime }\rangle ^{2})^{1/2}.\end{eqnarray}$$

More than for the TKE, the effect of stability on the stress profiles is very pronounced. The stresses (and hence the momentum fluxes) are monotonically damped with increasing stability. This damping of fluxes by buoyancy will be a critical factor in reducing TKE in stable boundary layers, as will be illustrated in the next subsection.

Figure 13.

Time-averaged non-dimensional Reynolds stress ( $|{\it\tau}/{\it\rho}|$ ) profiles at $\mathit{Re}_{f}$ of (a) 400, (b) 600 and (c) 900 at various stabilities compared with neutral cases. The bulk Richardson numbers $\mathit{Ri}_{b}$ given in the legend correspond, from lowest to highest, to the following initial surface Richardson numbers: $\mathit{Ri}_{0,0}=0.001$ , 0.002, 0.005.

5.4. TKE and momentum flux budgets

We have seen from vertical profiles of TKE (figure 12) that buoyancy affects turbulence at all heights, but more so away from the wall where the velocity shear is weaker. It has also been observed in the evolution of TKE with time that, as static stability is increased, the total TKE in the domain decreases (figures 7 b, 8 b and 9 b). To investigate the mechanisms causing this reduction in TKE, we study the TKE budgets for stable cases in comparison with neutral ones at the same Reynolds number. Figure 14 shows TKE budgets for all the Reynolds numbers and with varying stabilities.

Figure 14.

(a) Turbulence kinetic energy ( $q^{2}=\langle u_{i}^{\prime }u_{i}^{\prime }\rangle /2$ ) budget neutral-case comparison for $\mathit{Re}_{f}=600$ (dotted-dashed lines) and $\mathit{Re}_{f}=900$ (solid lines); and (b–d) comparison of the neutral case with stable cases at Reynolds numbers ( $\mathit{Re}_{f}$ ) of (b) 400, (c) 600 and (d) 900. Different line styles indicate the various stabilities and different greyscales (colours online) denote the various terms of the budget: $\partial q/\partial t+U_{k}(\partial q/\partial x_{k})=(P_{ii}+T_{ii}+{\it\Pi}_{ii}+B_{ii}+C_{ii}+D_{ii}+{\it\epsilon}_{ii})/2$ , where turbulent transport $T_{ii}=-\partial \langle u_{i}^{\prime }u_{i}^{\prime }w^{\prime }\rangle /\partial z$ , shear production $P_{ii}=-\langle u_{i}^{\prime }w^{\prime }\rangle \partial U_{i}/\partial z$ , viscous diffusion $D_{ii}={\it\nu}\partial ^{2}\langle u_{i}^{\prime }u_{i}^{\prime }\rangle /\partial ^{2}z$ and viscous dissipation rate ${\it\epsilon}_{ii}=-2{\it\nu}\langle (\partial u_{i}^{\prime }/\partial x_{k})(\partial u_{i}^{\prime }/\partial x_{k})\rangle$ . Note that the computation of the terms is simplified to take advantage of the horizontal homogeneity of the flow and of the zero mean vertical velocity. The terms that are not plotted or shown in the legend have been computed and found to be negligibly small (including the Coriolis term $C_{ii}$ , which theoretically is zero), but the full expressions of these terms can be found in appendix C. In panel (a) we plot the residual $R$ , which should be zero for validation. We also checked the residuals in non-neutral cases, and they are sufficiently small (but non-zero since the simulations are not exactly in steady state) to validate our computations of the different terms. All terms are normalized by $u_{\ast ,neutral}^{4}/{\it\nu}$ . The vertical grey line at the centre of each panel indicates the zero gain/loss line.

Figure 14(a) shows the trends of all the terms in the TKE budget equation for the neutral simulations at Reynolds numbers of $\mathit{Re}_{f}$ of 600 and 900. Note that the peaks of the shear production $P_{ii}$ , turbulent transport $T_{ii}$ and viscous diffusion $D_{ii}$ have almost the same value near the wall, after normalization by the relevant velocity scale (friction velocity $u_{\ast }$ in this case; see caption of the figure) at both Reynolds numbers. However, the height at which those peaks occur is reduced with increasing Reynolds number, as expected. The vertical trends of all terms are in very good qualitative agreements with the budget computed by Ansorge & Mellado (Reference Ansorge and Mellado2014) (exact quantitative comparison is not possible owing to small differences in the Reynolds number). Figure 14(b–d) compares these terms in neutral and (continuously turbulent) stable cases at $\mathit{Re}_{f}=400,\ 600$ and 900, respectively. When stability is imposed, though the overall shapes remain more or less the same, the magnitudes of the individual terms drop. The friction velocity used for non-dimensionalization is the neutral-case value $u_{\ast ,neutral}$ at the corresponding Reynolds number, so that the reduction in the budget terms with stability is captured. Normalization based on respective friction velocities for stable cases has been carried out (not shown here), and it was observed that the profiles of the TKE budget terms for $\mathit{Re}_{f}=900$ (but not for lower $\mathit{Re}_{f}$ ) collapse well on to the neutral-case lines, as is observed in the ABL (see Stull Reference Stull1988).

The TKE buoyancy destruction terms are not shown in the TKE budget plots, but are plotted separately in figure 15. The term is very small in comparison to the other terms in the TKE budget, but it is not zero. In fact, one can clearly observe the increase in its magnitude with increasing static stability. This implies that the main direct impact of stability on the flow is not through TKE destruction, but rather through another mechanism: from the plots in figure 14 one can see that a drastic drop in the shear production of TKE is induced by stability. Hence, it is the inhibition of shear production that seems to be the most prominent impact of stability on the TKE budget. This is also in agreement with the conclusions of Ansorge & Mellado (Reference Ansorge and Mellado2014). Since, in a steady turbulence state with small direct buoyant destruction, domain-averaged TKE dissipation has to be approximately equal to TKE production, dissipation is also significantly reduced. Viscous diffusion terms near the wall and turbulent transport, which have smaller contributions in the TKE budget, are similarly decreased.

Figure 15.

Profiles of buoyant destruction $B_{ii}/2=g\langle w^{\prime }{\it\theta}^{\prime }\rangle /{\it\theta}_{\infty }$ normalized by the viscosity and the neutral-case friction velocity, the same non-dimensionalization used in the TKE budget, at (a)  $\mathit{Re}_{f}=600$ and (b)  $\mathit{Re}_{f}=900$ .

To elucidate the mechanism and causes of the drop in shear TKE production, we investigate the quantities involved in this term, namely the gradients and fluxes. The fluxes ( $\langle u^{\prime }w^{\prime }\rangle ,\langle v^{\prime }w^{\prime }\rangle$ ) and the gradients ( $\partial \langle u\rangle /\partial z,\partial \langle v\rangle /\partial z$ ) that enter in the two dominant terms of the TKE shear production,

(5.7a,b ) $$\begin{eqnarray}-\langle u^{\prime }w^{\prime }\rangle \frac{\partial \langle u\rangle }{\partial z}\quad \text{and}\quad -\langle v^{\prime }w^{\prime }\rangle \frac{\partial \langle v\rangle }{\partial z},\end{eqnarray}$$

are plotted in figure 16 for $\mathit{Re}_{f}=900$ . The effect of stability on velocity gradients near the wall is significant, while away from the wall, this effect is reduced. In general, gradients are reduced near the wall by a factor of approximately 2 and increased away from the wall, with higher gradients for the more stable conditions at heights where peaks in the TKE production are seen. Changes in gradients thus cannot explain the sharp drops in TKE production. The effect of stability on the Reynolds stresses ( $=-\text{vertical momentum fluxes}$ ), however, is much more significant. The fluxes are reduced due to static stability throughout the boundary layer compared to the neutral case. This indicates that the reduction due to stability in the vertical transport of momentum from the higher elevations in the boundary layer into the near-wall region, where shear production can occur, is the main mechanism for the reduction in this TKE production term in stable wall-bounded flows.

Figure 16.

(a) Magnitude of turbulent stress ${\it\tau}=(\langle u^{\prime }w^{\prime }\rangle ^{2}+\langle v^{\prime }w^{\prime }\rangle ^{2})^{1/2}$ and (b) magnitude of the viscous stress ${\it\nu}|\langle \text{d}\boldsymbol{u}/\text{d}z\rangle |$ , which represents the gradients (where $|\langle \text{d}\boldsymbol{u}/\text{d}z\rangle |=(\langle \text{d}u/\text{d}z\rangle ^{2}+\langle \text{d}v/\text{d}z\rangle ^{2})^{1/2}$ ), at $\mathit{Re}_{f}=900$ .

The normalization used in figure 16 also makes the gradients plotted in panel (b) equivalent to viscous stresses, normalized by $G^{2}$ , while panel (a) depicts the turbulent stresses also normalized by $G^{2}$ . Hence from this figure one can compare the viscous and turbulent stresses and observe the increase in the relative importance of the turbulent part as the distance to the wall increases, and as stability decreases.

To further understand the reduction of the momentum fluxes by buoyancy, we study the momentum flux budgets for $\langle u^{\prime }w^{\prime }\rangle$ (the $\langle v^{\prime }w^{\prime }\rangle$ budget was computed but is not shown here since the behaviour is very similar to the $\langle u^{\prime }w^{\prime }\rangle$ budget). Figure 17(a) shows all the terms of the flux budget for the neutral case with $\mathit{Re}_{f}=900$ . Details of the equations for the flux budgets and all the terms in it can be found in appendix C. For clarity, only the significant terms have been plotted in figure 17(b–d), where the effects of stability are illustrated. Note that, since the flux is negative, sources would be negative and sinks positive. As seen in the TKE budget plots, the peaks in vertical profiles shift towards the wall with increasing Reynolds numbers. The magnitudes of the peak values also increase with Reynolds number, even though each budget is normalized with the neutral friction velocity at its corresponding Reynolds number. The figure illustrates that the pressure correlation and return to isotropy terms, which have opposite signs near the surface, are both reduced in magnitude as stability increases; the changes in these terms near the surface for stronger stabilities partially offset each other. The flux production (which like the TKE shear production depends on the fluxes themselves) and the turbulent transport terms are also significantly damped by buoyancy. Given that the highest gradients are in the vertical direction, the dominant components of these terms will be $-\partial \langle u^{\prime }w^{\prime }w^{\prime }\rangle /\partial z$ in the transport and $-\langle w^{\prime }w^{\prime }\rangle \partial U/\partial z$ in the production, both of which are strongly dependent on $w^{\prime }$ . Therefore, damping of turbulent transport of momentum, which as we showed leads to a reduction in TKE, seems to be a result of the damping of $w^{\prime 2}$ , which is a component of the TKE. This circular dependence thus needs further elucidation. Note that the buoyancy term $B_{13}$ here remains much smaller than the other components.

Figure 17.

(a) Momentum flux ( $\langle u^{\prime }w^{\prime }\rangle$ ) budget for the neutral case at $\mathit{Re}_{f}=900$ , and (b–d) comparison of neutral and stable cases at Reynolds numbers ( $\mathit{Re}_{f}$ ) of (b) 400, (c) 600 and (d) 900. Different line styles indicate the various stabilities and different greyscales (colours online) denote the various terms of the budget: $\partial \langle u^{\prime }w^{\prime }\rangle /\partial t+U_{k}(\partial \langle u^{\prime }w^{\prime }\rangle /\partial x_{k})=P_{13}+T_{13}+C_{13}+D_{13}+{\it\Pi}_{13}+{\it\phi}_{13}+{\it\epsilon}_{13}+B_{13}$ , where turbulent transport $T_{13}=-\partial \langle u^{\prime }w^{\prime }w^{\prime }\rangle /\partial z$ , turbulent production $P_{13}=-\langle w^{\prime }w^{\prime }\rangle \partial U/\partial z$ , pressure transport ${\it\Pi}_{13}=-(1/{\it\rho}_{\infty })(\partial \langle p^{\prime }u^{\prime }\rangle /\partial z+\partial \langle p^{\prime }w^{\prime }\rangle /\partial x)$ , pressure–strain redistribution ${\it\phi}=(2/{\it\rho}_{\infty })\langle p^{\prime }S_{13}^{\prime }\rangle$ with strain-rate tensor $\mathsf{S}_{ij}=[\partial u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i}]/2$ , molecular diffusion $D_{13}={\it\nu}\partial ^{2}\langle u^{\prime }w^{\prime }\rangle /\partial ^{2}z$ , viscous dissipation ${\it\epsilon}_{13}=-2{\it\nu}\langle (\partial u^{\prime }/\partial x_{k})(\partial w^{\prime }/\partial x_{k})\rangle$ , and buoyancy $B_{13}=g\langle u^{\prime }{\it\theta}^{\prime }\rangle /{\it\theta}_{\infty }$ . All terms are normalized by $u_{\ast ,neutral}^{4}/{\it\nu}$ . The vertical grey line at the centre of each panel indicates the zero gain/loss line, but note that here sources would be negative since the momentum flux is negative. The computation of the terms is simplified to take advantage of the horizontal homogeneity of the flow and of the zero mean vertical velocity. Full budget equations can be found in appendix C.

To clarify the sequence of effects that stability induces in a turbulent boundary layer, we now plot the budget of the vertical velocity variance (figure 18). The plot in figure 18(a) compares the budget of vertical velocity variance in the neutral case to the case with stability at $\mathit{Ri}_{0,0}=0.002$ for Reynolds number $\mathit{Re}_{f}=600$ , while figure 18(b) gives the same comparison at $\mathit{Re}_{f}=900$ and $\mathit{Ri}_{0,0}=0.005$ . It is observed from the plots that the change in the buoyancy destruction here, while still small, is of the same order as the changes in the other terms. For example, the drop in the peak of the pressure redistribution term, which is the main source in this budget, from the neutral to the most stable simulation at $\mathit{Re}_{f}=900$ is only 2.5 times the increase in buoyancy destruction. Notice that we are comparing the change in the terms from neutral to stable simulations, not their magnitude (for the neutral case this buoyancy destruction is zero). This is thus the only budget where the direct effect of buoyancy, via a destruction term, is significant. This implies that buoyancy first affects motions in the vertical direction, damping $w^{\prime }$ and $\langle w^{\prime }w^{\prime }\rangle$ , which then causes a reduction in the production, and hence in the magnitude, of the momentum fluxes ( $\langle u^{\prime }w^{\prime }\rangle$ and $\langle v^{\prime }w^{\prime }\rangle$ ), which then leads to lower TKE production, particularly of the $\langle u^{\prime 2}\rangle$ and $\langle v^{\prime 2}\rangle$ components of the TKE (the shear production of $\langle w^{\prime 2}\rangle$ is zero since it involves gradients of the mean vertical velocity). The reduction in $\langle u^{\prime 2}\rangle$ and $\langle v^{\prime 2}\rangle$ could then lead to the reduction in the pressure–strain redistribution term illustrated in figure 18, further reducing $\langle w^{\prime 2}\rangle$ , but more in-depth analysis of the impact of buoyancy on this redistribution term is needed to ascertain that. We also analysed the turbulent transport efficiencies estimated via $R_{uw}=\langle u^{\prime }w^{\prime }\rangle /{\it\sigma}_{u}{\it\sigma}_{w}$ and $R_{vw}=\langle v^{\prime }w^{\prime }\rangle /{\it\sigma}_{v}{\it\sigma}_{w}$ , and, while they decrease with increasing stability away from the wall, near the wall they remain about constant. This confirms that the reduction in the momentum fluxes is related to a drop in the amplitude of the vertical velocity, rather than to a drop in the turbulent transport efficiency.

Figure 18.

Comparison of budgets of the variance of vertical velocity ( $\langle w^{\prime }w^{\prime }\rangle$ ) for neutral case with stable cases: (a)  $\mathit{Re}_{f}=600$ , where solid line corresponds to neutral case and dashed line corresponds to stable case with $\mathit{Ri}_{b}=0.200$ ; (b)  $\mathit{Re}_{f}=900$ , where solid line corresponds to neutral case and dashed line corresponds to stable case with $\mathit{Ri}_{b}=0.968$ . Different line styles indicate the various stabilities and different greyscales (colours online) denote the various terms of the budget: $\partial \langle w^{\prime }w^{\prime }\rangle /\partial t+U_{k}(\partial \langle w^{\prime }w^{\prime }\rangle /\partial x_{k})=P_{33}+T_{33}+C_{33}+D_{33}+{\it\Pi}_{33}+{\it\phi}_{33}+{\it\epsilon}_{33}+B_{33}$ , where turbulent transport $T_{33}=-\partial \langle w^{\prime }w^{\prime }w^{\prime }\rangle /\partial z$ , pressure transport ${\it\Pi}_{33}=(-2/{\it\rho}_{\infty })\ \partial \langle p^{\prime }w^{\prime }\rangle /\partial z$ , pressure–strain redistribution ${\it\Phi}_{33}=(2/{\it\rho}_{\infty })\langle p^{\prime }S_{33}^{\prime }\rangle$ with strain-rate tensor $\mathsf{S}_{ij}=[\partial u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i}]/2$ , viscous diffusion $D_{33}={\it\nu}\partial ^{2}\langle w^{\prime }w^{\prime }\rangle /\partial ^{2}z$ , viscous dissipation rate ${\it\epsilon}_{33}=-2{\it\nu}\langle (\partial w^{\prime }/\partial x_{k})(\partial w^{\prime }/\partial x_{k})\rangle$ and buoyancy destruction $B_{33}=2(g/{\it\theta}_{\infty })\langle w^{\prime }{\it\theta}^{\prime }\rangle$ . All terms are normalized by $u_{\ast ,neutral}^{4}/{\it\nu}$ ; minor terms including production, are not shown. The vertical grey line at the centre of each panel indicates the zero gain/loss line. The computation of the terms is simplified to take advantage of the horizontal homogeneity of the flow and of the zero mean vertical velocity. Full budget equations can be found in appendix C.

Figure 19 shows the profiles of the horizontal velocity variance (a) and the vertical velocity variance (b), illustrating that the reduction in $\langle w^{\prime 2}\rangle$ relative to the neutral case, which is associated with the damping of the vertical motions, is larger than the relative reduction in $\langle u^{\prime 2}\rangle$ and $\langle v^{\prime 2}\rangle$ , as can be anticipated. A particularly important implication of this finding is that even a small Richardson number associated with a small direct buoyant TKE destruction, compared to shear production, can have a significant impact on turbulence in wall-bounded flows: this small TKE destruction will primarily reduce the vertical variance, which then drastically reduces production compared to a neutral flow. This might explain the very rapid drop in TKE levels shortly after the evening transition in the ABL (Nadeau et al. Reference Nadeau, Pardyjak, Higgins, Fernando and Parlange2011). It is also worth underlining that the reduction in momentum flux under stable conditions results in a reduction in the surface drag. In a pipe or channel flow, this might not be possible if the wall drag needs to balance a fixed imposed pressure gradient. But in an Ekman boundary layer, this is always possible since it just entails a rebalancing of the equilibrium between the pressure gradient, Coriolis and surface drag forces in such a way that the Coriolis force increases and compensates for the reduced surface drag in this equilibrium.

Figure 19.

Reduction in the sum of the horizontal velocities variances (a) compared to the reduction in the vertical velocity variance (b) with increasing static stability at $\mathit{Re}_{f}=900$ . The black solid line corresponds to the neutral case; and other lines correspond to stable cases: dashed line (blue online) $\mathit{Ri}_{b}=0.193$ and dash-dotted line (red online) $\mathit{Ri}_{b}=0.387$ .

5.5. Buoyancy flux and temperature variance budgets

Heat flux budgets for all the stable cases are shown in figure 20. The shear production of heat flux given by $P^{2}=-\langle {\it\theta}^{\prime }u_{j}^{\prime }\rangle \,\partial \langle w\rangle /\partial x_{j}$ is identically zero for the flows simulated here. The main source of covariance or flux is the production term with temperature gradient $P^{1}=-\langle w^{\prime }u_{j}^{\prime }\rangle \,\partial \langle {\it\theta}\rangle /\partial x_{j}$ (note that production and other sources here are negative since the heat flux is negative or downwards; in a statically unstable boundary layer it would be positive). The two significant terms are the pressure terms, while other terms (namely the Coriolis term $C$ , dissipation ${\it\epsilon}_{23}$ , viscous diffusion $D_{23}$ and turbulent transport $T_{23}$ ) are quite small for neutral simulations. For simulations with stable stratification, the reduction in the temperature gradient production term drives the buoyancy effect on this budget. However, the buoyancy destruction term becomes important at stronger stabilities, especially for the $\mathit{Re}_{f}=900$ simulation. Note here that, for increasing stabilities, the vertical temperature gradient is larger in magnitude, and yet the gradient production term is smaller. This results from the significant drop in the magnitude of the variance of the $w$ component of velocity with stability, illustrated in the previous subsection, which appears in this production term. Expressions for all the terms involved in the heat flux budget can be found in appendix C.

Figure 20.

(a) Buoyancy flux ( $\langle w^{\prime }{\it\theta}^{\prime }\rangle$ ) budget for the stable case at $\mathit{Re}_{f}=900$ , $\mathit{Ri}_{0,0}=0.001$ and (b–d) comparison of all significant terms from the stable cases at Reynolds numbers ( $\mathit{Re}_{f}$ ) of (b) 400, (c) 600 and (d) 900. Different line styles indicate the various stabilities and different greyscales (colours online) denote the various terms of the budget: $\partial \langle w^{\prime }{\it\theta}^{\prime }\rangle /\partial t+U_{j}(\partial \langle w^{\prime }{\it\theta}^{\prime }\rangle /\partial x_{j})=P^{1}+P^{2}+T+{\it\Pi}+{\it\phi}+D^{1}+D^{2}+{\it\chi}+B$ , where turbulent production by mean temperature gradient $P^{1}=-\langle w^{\prime }w^{\prime }\rangle \partial \langle {\it\theta}\rangle /\partial z$ , buoyancy destruction $B=g\langle {\it\theta}^{\prime }{\it\theta}^{\prime }\rangle /{\it\theta}_{\infty }$ , turbulent transport $T=-\partial \langle w^{\prime }w^{\prime }{\it\theta}^{\prime }\rangle /\partial z$ , transport by pressure correlation ${\it\Pi}=-(1/{\it\rho}_{\infty })\partial \langle p^{\prime }{\it\theta}^{\prime }\rangle /\partial z$ , redistribution by return-to-isotropy term ${\it\phi}=(1/{\it\rho}_{\infty })\langle p^{\prime }\partial {\it\theta}^{\prime }/\partial z\rangle$ , diffusion by velocity gradients $D^{1}=(\partial /\partial z)({\it\nu}\partial \langle {\it\theta}^{\prime }\partial w^{\prime }/\partial z\rangle )$ , diffusion by temperature gradients $D^{2}=(\partial /\partial z)({\it\nu}_{{\it\theta}}\partial \langle w^{\prime }\partial {\it\theta}^{\prime }/\partial z\rangle )$ , and molecular dissipation ${\it\chi}=-({\it\nu}+{\it\nu}_{{\it\theta}})\langle (\partial w^{\prime }/\partial x_{j})(\partial {\it\theta}^{\prime }/\partial x_{j})\rangle$ . All terms are normalized by $u_{\ast ,neutral}^{3}({\it\theta}_{\infty }-{\it\theta}_{0})/{\it\nu}$ . The vertical grey line at the centre of each panel indicates the zero gain/loss line, but note that here sources would be negative since the buoyancy flux is negative in a stable boundary layer. The computation of the terms is simplified to take advantage of the horizontal homogeneity of the flow and of the zero mean vertical velocity. Full budget equations can be found in appendix C.

The last budget we examine is the budget of the perturbation temperature variance ( $\langle {\it\theta}^{\prime }{\it\theta}^{\prime }\rangle$ ) plotted in figure 21. Similar to TKE budget plots, from these plots one can infer that the main terms being modified by stability are the variance production terms. Variance dissipation falls with decreasing production. Turbulent transport and molecular diffusion away from the wall are small, while near the wall they are very weakly affected by buoyancy.

Figure 21.

(a) Temperature variance ( $\langle {\it\theta}^{\prime }{\it\theta}^{\prime }\rangle$ ) budget at $\mathit{Re}_{f}=900$ and $\mathit{Ri}_{0,0}=0.001$ and (b–d) comparison of all the significant terms from the stable cases at Reynolds numbers ( $\mathit{Re}_{f}$ ) of (b) 400, (c) 600 and (d) 900 and at various Richardson numbers. Different line styles indicate the various stabilities and different greyscales (colours online) denote the various terms of the budget: $\partial \langle {\it\theta}^{\prime }{\it\theta}^{\prime }\rangle /\partial t+U_{j}(\partial \langle {\it\theta}^{\prime }{\it\theta}^{\prime }\rangle /\partial x_{j})=P_{{\it\theta}}+T_{{\it\theta}}+D_{{\it\theta}}+{\it\chi}_{{\it\theta}}$ , where production $P_{{\it\theta}}=-2\langle w^{\prime }{\it\theta}^{\prime }\rangle \partial \langle {\it\theta}\rangle /\partial z$ , turbulent transport $T_{{\it\theta}}=-\partial \langle w^{\prime }{\it\theta}^{\prime }{\it\theta}^{\prime }\rangle /\partial z$ , molecular diffusion $D_{{\it\theta}}={\it\nu}_{{\it\theta}}\partial ^{2}\langle {\it\theta}^{\prime }{\it\theta}^{\prime }\rangle /\partial ^{2}z$ , and molecular dissipation ${\it\epsilon}_{{\it\theta}}=-2{\it\nu}_{{\it\theta}}\langle (\partial {\it\theta}^{\prime }/\partial x_{j})(\partial {\it\theta}^{\prime }/\partial x_{j})\rangle$ . All terms are normalized by $u_{\ast ,neutral}^{2}({\it\theta}_{\infty }-{\it\theta}_{0})^{2}/({\it\nu}/\mathit{Pr})$ . The vertical grey line at the centre of each panel indicates the zero gain/loss line. The computation of the terms is simplified to take advantage of the horizontal homogeneity of the flow and of the zero mean vertical velocity. Full budget equations can be found in appendix C.

5.6. Turbulent diffusivities under stable conditions

Most stable flows of practical engineering or geophysical relevance are at much higher Reynolds numbers than the ones we can attain here, and many can only be modelled using the Reynolds-averaged Navier–Stokes (RANS) approach (e.g. stable ABL in weather and climate models). Thus, an important follow-up on the investigation of turbulence dynamics and budgets that we performed in the previous sections is to relate the findings to the development of better turbulence closure models for stable conditions and the comparison to existing models for very high- $\mathit{Re}$ flows. Readers interested in higher-order closures will find the detailed budget analyses very useful for supporting simplifications of the budget equations that are explicitly solved in models and for guiding the development of such closures. A large fraction of weather and climate models, however, continue to rely on simple first-order or 1.5-order closures under stable conditions; therefore, we focus in this section on the eddy diffusivities and viscosities needed for these closures. We do not aim to develop new closures, but rather to test the assumptions underlying existing ones.

Figure 22 shows the vertical profiles of eddy viscosity ( ${\it\nu}_{T}$ ) and eddy diffusivity ( ${\it\alpha}_{T}$ ) for $\mathit{Re}_{f}$ of 400, 600 and 900, non-dimensionalized by molecular viscosity and diffusivity. Following Coleman et al. (Reference Coleman, Ferziger and Spalart1992), eddy viscosity and diffusivity are computed using the following relations:

(5.8) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\nu}_{T}=\frac{-\langle u^{\prime }w^{\prime }\rangle \,\text{d}\langle u\rangle /\text{d}z-\langle v^{\prime }w^{\prime }\rangle \,\text{d}\langle v\rangle /\text{d}z}{(\text{d}\langle u\rangle /\text{d}z)^{2}+(\text{d}\langle v\rangle /\text{d}z)^{2}}, & \displaystyle\end{eqnarray}$$

(5.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\alpha}_{T}=\frac{-\langle w^{\prime }{\it\theta}^{\prime }\rangle }{\text{d}\langle {\it\theta}\rangle /\text{d}z}. & \displaystyle\end{eqnarray}$$

Figure 22.

Profiles of (a,c,e) eddy viscosity and (b,d, f) eddy diffusivity non-dimensionalized by their molecular counterparts under varying stabilities for $\mathit{Re}_{f}$ of (a,b) 400, (c,d) 600 and (e, f) 900.

In figure 23, we plot the vertical eddy viscosities evaluated separately for the $x$ and $y$ momentum components ( ${\it\nu}_{Tx}=-\langle u^{\prime }w^{\prime }\rangle /(\text{d}\langle u\rangle /\text{d}z)$ and ${\it\nu}_{Ty}=-\langle v^{\prime }w^{\prime }\rangle /(\text{d}\langle v\rangle /\text{d}z)$ ) to test the assumption that the eddy viscosity is isotropic (hence a scalar, see Smagorinsky (Reference Smagorinsky1963), for example), which follows from the Boussinesq analogy. Note that the results are plotted only for values of $z/L_{z}$ ( $L_{z}$ is the domain height) less than 0.2 and where $\text{d}\langle v\rangle /\text{d}z$ is not close to zero. From the plots, it seems that this assumption is not exactly applicable: the eddy viscosity evaluated in the $y$ direction, ${\it\nu}_{Ty}$ , is consistently smaller than that evaluated in the $x$ direction, ${\it\nu}_{Tx}$ . However, one must note that the discrepancy between the two viscosities seems to be reduced at higher Reynolds numbers and also at higher Richardson numbers. These trends are in agreement with the LES results of Huang, Bou-Zeid & Golaz (Reference Huang, Bou-Zeid and Golaz2013) (at much higher Reynolds numbers), who showed that the Boussinesq analogy was reasonable and its applicability in fact improved with increasing stability. Assuming equality of the two viscosities is equivalent to hypothesizing that the fluxes and gradients or strains are aligned. Plots not shown here of the difference between the angles of the stress ( $\tan ^{-1}(\langle v^{\prime }w^{\prime }\rangle /\langle u^{\prime }w^{\prime }\rangle )$ ) and the gradient ( $\tan ^{-1}(\partial \langle v\rangle /\partial z)/(\partial \langle u\rangle /\partial z)$ ) also show improved alignment at higher Reynolds numbers (the angle differences are less than 16° for $\mathit{Re}_{f}=600$ and less than 10° for $\mathit{Re}_{f}=900$ ). These angles are much larger than the ones observed in the wall-modelled LES of Huang et al. (Reference Huang, Bou-Zeid and Golaz2013), again suggesting a Reynolds-number dependence despite the uncertainty in defining the effective Reynolds number in wall-modelled LES.

Figure 23.

Profiles of the vertical eddy viscosity for the $x$ and $y$ momentum components calculated separately ( ${\it\nu}_{Tx}=-\langle u^{\prime }w^{\prime }\rangle /(\text{d}\langle u\rangle /\text{d}z)$ and ${\it\nu}_{Ty}=-\langle v^{\prime }w^{\prime }\rangle /(\text{d}\langle v\rangle /\text{d}z)$ ) non-dimensionalized by their molecular counterparts under varying stabilities for $\mathit{Re}_{f}$ of (a) 600 and (b) 900. Missing parts in the vertical profiles for the $y$ momentum component correspond to layers where the gradients of $v$ have near-zero values, resulting in very large viscosities, and hence not shown.

The ratio of ${\it\nu}_{T}$ over ${\it\alpha}_{T}$ gives the turbulent Prandtl number $\mathit{Pr}_{T}={\it\nu}_{T}/{\it\alpha}_{T}$ , which is plotted in figure 24 for all the Reynolds numbers simulated here, at Richardson numbers of $\mathit{Ri}_{0,0}=0.001$ and 0.002 (only 0.001 for the lower $\mathit{Re}_{f}$ ) versus height (panels (a,c,e)) and versus the gradient Richardson number $\mathit{Ri}_{gr}$ (panels (b,d, f)). The turbulent Prandtl number increases, but very slightly, with increasing stability ( $\mathit{Ri}_{0,0}$ ), approaching unity. This indicates that buoyancy affects turbulent transfer of heat and momentum in the same way in statically stable flows. This has also been observed by Coleman et al. (Reference Coleman, Ferziger and Spalart1992). However, for higher Reynolds number, several studies suggest that the turbulent Prandtl number increases more significantly with increasing stability (Bou-Zeid et al. Reference Bou-Zeid, Higgins, Huwald, Meneveau and Parlange2010; Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2010; Huang & Bou-Zeid Reference Huang and Bou-Zeid2013). Despite this divergence on the effect of stability on the similarity of momentum and scalar turbulent transport, a robust conclusion from all those studies is that the Reynolds analogy holds much better in stable wall-bounded flows than in unstable (convective) ones (Li et al. Reference Li, Katul and Bou-Zeid2012; Katul et al. Reference Katul, Li, Chamecki and Bou-Zeid2013). Therefore, a first-order closure seems to be more acceptable under stable conditions.

Figure 24.

Turbulent Prandtl number ( $\mathit{Pr}_{t}$ ) variation with height and with gradient Richardson number ( $\mathit{Ri}_{gr}$ ) at Reynolds numbers (a,b) 400, (c,d) 600 and (e, f) 900.

Figure 25(a,c,e) depicts the ratio of the magnitude of the horizontal Reynolds stress to the square of the vertical velocity variance (which is defined as $b_{1}$ in the caption of the figure), with varying stabilities and Reynolds numbers, as a function of height (we only plot heights above the near-wall region where viscous stresses dominate since this analysis is not meaningful there). Figure 25(b,d, f) depicts the ratio of the same stress component to twice the TKE (defined as $a_{1}$ in the caption of the figure), again with varying Richardson and Reynolds numbers, as a function of height. These quantities represent measures of turbulent transport efficiencies in the stable boundary layer. One can observe that these fluxes, when non-dimensionalized by the vertical velocity variance, collapse very well onto a single curve in the surface layer (log or MOST layer, roughly where $z/{\it\delta}_{t}<0.2$ ), except for some discrepancy for the most stable simulation at the highest $\mathit{Re}_{f}$ . The ratio $b_{1}$ is not very sensitive to stability in the surface layer since both the numerator and denominator decrease proportionally. A reduction in $b_{1}$ can be noted at larger heights, indicating that turbulent transport becomes less efficient with increasing stability. However, when non-dimensionalized with TKE, the lines do not collapse anywhere and this measure of efficiency drops down significantly with increasing stability. Horizontal variances are not reduced significantly by stability (thus the preponderance of pancake-like structures in the stable ABL (Beckers et al. Reference Beckers, Verzicco, Clercx and Van Heijst2001; Riley & de Bruyn Kops Reference Riley and de Bruyn Kops2003; Davidson Reference Davidson2004)) and hence the TKE decreases with stability at a lower rate than vertical fluxes. These results reconfirm that the variance of the vertical velocity is the primary component of the TKE that accounts for vertical fluxes and suggest that 1.5-order closure models should solve the budget equation of the vertical velocity variance, rather than TKE, or use that variance as the basis for formulating vertical eddy viscosities and diffusivities in stable wall-bounded flows. However, solving the equation for the vertical variance can pose some challenges, particularly ones related to the closure of the pressure redistribution term, and this might explain why currently most formulations rely on the full TKE as the basis for eddy viscosity and diffusivity models.

Figure 25.

Ratio of the Reynolds stress magnitude to vertical velocity variance, $b_{1}$ , at $\mathit{Re}_{f}$ of (a) 400, (c) 600 and (e) 900, respectively, and ratio of Reynolds stress magnitude to twice the TKE, $a_{1}$ , at $\mathit{Re}_{f}$ of (b) 400, (d) 600 and ( f) 900, respectively, as a function of height and bulk Richardson number, where $b_{1}=(\langle u^{\prime }w^{\prime }\rangle ^{2}+\langle v^{\prime }w^{\prime }\rangle ^{2})^{1/2}/{\it\sigma}_{w}^{2}$ and $a_{1}=(\langle u^{\prime }w^{\prime }\rangle ^{2}+\langle v^{\prime }w^{\prime }\rangle ^{2})^{1/2}/2q^{2}$ with $2q^{2}=\langle u_{i}^{\prime }u_{i}^{\prime }\rangle$ and ${\it\sigma}_{w}^{2}=\langle w^{\prime }w^{\prime }\rangle$ . We only plot for $z/{\it\delta}_{t}>0.05$ , where ${\it\delta}_{t}$ is the turbulent boundary layer depth of each simulation, to focus on the regions where the total fluxes are not dominated by the viscous contributions.