2.1 Normalization of the equations and governing parameters

To express the governing equations as a function of suitable dimensionless parameters, a characteristic time, length, buoyancy and velocity scale can be defined as

where the surface buoyancy $\hat{b}_{s}$ determines whether the flow is anabatic ( $\hat{b}_{s}>0$ ) or katabatic ( $\hat{b}_{s}<0$ ). These aforementioned parameters can now be used to introduce the following normalized variables

Relations (2.7) are derived selecting $\hat{b}_{s}$ and $\hat{N}$ as repeating parameters though this choice is by no means unique and other options are possible. Substituting the expressions (2.8) into the governing equations (2.4)–(2.6) results in

where $\mathit{Gr}\equiv \hat{b}_{s}^{4}\,\hat{\unicode[STIX]{x1D708}}^{-2}\,\hat{N}^{-6}$ is the Grashof number, defined as the ratio between buoyancy and viscous forces. There are a number of advantages to adopting $Gr$ over other possible dimensionless numbers (e.g. Reynolds number). For the particulars of sloping flow with zero slip at the boundary, $Gr$ measures the relative strength of the body (or non-contact) force to the surface (or contact) force arising from friction (or viscosity). The flow adjustment to this imbalance gives rise to non-hydrostatic pressure distribution and advective acceleration.

Table 1.

Geometry and parameters for the DNS runs. $L_{i}$ and $N_{i}$ denote the domain size and the number of collocation nodes in the three coordinate directions, respectively, $T$ denotes the characteristic oscillation period characterizing the buoyancy and velocity fields (see § 3.4), $\mathit{Gr}\equiv \hat{b}_{s}^{4}\,\hat{\unicode[STIX]{x1D708}}^{-2}\,\hat{N}^{-6}$ where $b_{s}$ is the imposed (normalized) surface buoyancy. Simulation labels indicate whether a specific run is in buoyantly [u]nstable or [s]table regime ( $b_{s}=+1$ and $b_{s}=-1$ respectively), the surface sloping angles $\unicode[STIX]{x1D6FC}$ , and which among the three considered Grashof numbers ([h]igh, [m]edium and [l]ow) is used; for instance u30h denotes an anabatic flows regime with $\unicode[STIX]{x1D6FC}=30^{\circ }$ at the highest among the considered Grashof numbers.

4.1 Time evolution and structure of the flow

The time evolution of the slope-normal integrated, space-averaged ( $x,y$ directions), normalized streamwise velocity $\langle u\rangle$ and buoyancy $\langle b\rangle$ is displayed in figure 2. The system exhibits the classical quasi-periodic, low-frequency, oscillatory behaviour (surges), superimposed to a base flow, as observed in FS09 and in experiments (e.g. Doran & Horst Reference Doran and Horst1981; Monti et al. Reference Monti, Fernando, Princevac, Chan, Kowalewski and Pardyjak2002; Princevac, Hunt & Fernando Reference Princevac, Hunt and Fernando2008). It can be shown (McNider Reference McNider1982) that the slope-normal integrated $\langle u\rangle ,\,\langle b\rangle$ variables resemble a system of coupled oscillators, which in the case of laminar flow are characterized by a period $\hat{T}\approx 2\unicode[STIX]{x03C0}(\hat{N}\sin \unicode[STIX]{x1D6FC})^{-1}$ (normalized period is $T\approx 2\unicode[STIX]{x03C0}/\sin \unicode[STIX]{x1D6FC}$ ). The oscillation periods from DNS are very close to those of corresponding laminar solutions, which for typical atmospheric values $\hat{N}=10^{-2}\,(\text{Hz})$ , $\unicode[STIX]{x1D6FC}=15^{\circ },\,30^{\circ }$ and $60^{\circ }$ , are $40,\,20$ and 10 min. A recent study (Fedorovich & Shapiro Reference Fedorovich and Shapiro2017) has shown how damping of integral oscillations by surface stress in laminar Prandtl flows decays very rapidly and independently from the flow forcing mechanism. Such fast reduction in the damping factor was also observed in the turbulent slope flows simulated in FS09. Current results further support such findings, i.e. observed integral oscillations are characterized by a progressively weaker decay rate as time advances.

Figure 2.

Time evolution of slope-normal integrated $\langle u\rangle$ (solid lines) and $\langle b\rangle$ (dashed lines) fields for simulations s90h, s60h, s30h and s15h (katabatic flow regime). The total time integration period is shown for each run.

Figure 3.

Dynamic (black lines) and energy (red lines) identities (4.1) and (4.2) for the considered simulations. Profiles have been shifted on the $y$ axis to allow for proper visualization. From top to bottom, profiles correspond to u15h, u30h, u60h, u90h and s15h, s30h, s60h, s90h, respectively. We here denote $\langle \unicode[STIX]{x1D70F}_{xz}^{tot}\rangle =\mathit{Gr}^{-1/2}(\text{d}\langle u\rangle /\text{d}z)-\langle u^{\prime }w^{\prime }\rangle$ and $\langle \unicode[STIX]{x1D70F}_{bz}^{tot}\rangle =\mathit{Gr}^{-1/2}Pr^{-1}(\text{d}\langle b\rangle /\text{d}z)-\langle b^{\prime }w^{\prime }\rangle$ (sum of molecular and turbulent kinematic fluxes of streamwise momentum and buoyancy in the slope-normal direction).

Averaging equations (2.9)–(2.11) in time and over directions of statistical homogeneity $(x,y)$ results in

where $\langle \unicode[STIX]{x1D70F}_{xz}^{tot}\rangle =\mathit{Gr}^{-1/2}(\text{d}\langle u\rangle /\text{d}z)-\langle u^{\prime }w^{\prime }\rangle$ and $\langle \unicode[STIX]{x1D70F}_{bz}^{tot}\rangle =\mathit{Gr}^{-1/2}Pr^{-1}(\text{d}\langle b\rangle /\text{d}z)-\langle b^{\prime }w^{\prime }\rangle$ are the normalized total (molecular $+$ turbulent) slope-normal kinematic momentum and buoyancy fluxes. Equations (4.1) and (4.2) can be used to test the ‘quality’ of computed statistics (steady state is guaranteed only if the two identities hold). Numerical results are displayed in figure 3, and certify that averaging over $3T$ , after a transient of at least $3T$ is sufficient to satisfy both equations (4.1) and (4.2).

Figure 4.

Pseudocolour plot of instantaneous katabatic flow streamwise velocity $u$ (a,b) and buoyancy $b$ (c,d), on the plane $y=L_{y}/2$ for simulations s60h (a,c) and s15h (b,d). The $z$ -axis denotes the slope-normal coordinate direction, whereas the $x$ -axis denotes the along-slope direction. The displayed $u(x,z)$ and $b(x,z)$ fields correspond to the crest of the last simulated oscillation for both runs. For detailed viewing, only the near-surface region of the total domain is shown.

Figure 5.

Pseudocolour plot of instantaneous anabatic flow streamwise velocity $u$ (a,b) and buoyancy $b$ (c,d), on the plane $y=L_{y}/2$ for simulations u60h (a,c) and u15h (b,d). The displayed $u(x,z)$ and $b(x,z)$ fields correspond to the crest of the last simulated oscillation for both runs. For detailed viewing, only the near-surface region of the total domain is shown.

Figures 4 and 5 display a pseudocolour plot of instantaneous streamwise normalized velocity field $(u)$ and normalized buoyancy field $(b)$ for simulations s60h, s15h and u60h, u15h, respectively. The TBL appears shallower than the DBL, as previously noted in FS09. A reversed flow characterizes the above-jet regions, resulting from the interaction between the flow and the background stably stratified environment, in qualitative agreement with the predictions of the Prandtl model (Prandtl Reference Prandtl1942). Progressive dissimilarity in DBL and TBL boundary layer thickness arises as the sloping angle decreases from the reference $\unicode[STIX]{x1D6FC}=90^{\circ }$ vertical wall set-up (recall that at $\unicode[STIX]{x1D6FC}=90^{\circ }$ the two flow regimes are equal up to a sign in $u_{i}$ ). As $\unicode[STIX]{x1D6FC}$ decreases, the LLJ regions in the anabatic flow solution experience a significant thickening induced by the convective type regime characterizing the flow at small sloping angles. Katabatic flows are instead characterized by a strong static stability at small sloping angles (see figure 4), which damps positive slope-normal velocity fluctuations, thus maintaining the LLJ relatively close to the wall and reducing the overall mixing of momentum and buoyancy in the near-wall regions. In addition, the strong stability induced by the imposed surface buoyancy in the katabatic flow regime at small $\unicode[STIX]{x1D6FC}$ results in partial laminarization of the LLJ.

4.2 Mean flow

Mean profiles of kinematic momentum $\langle u\rangle$ and buoyancy $\langle b\rangle$ are displayed in figure 6. Computed quantities qualitatively resemble those from the classic constant-K Prandtl analytic solution (Prandtl Reference Prandtl1942). The most significant features are a peak velocity $(u_{j})$ , the nose of the LLJ (occurring at $z_{j}$ ), and a return flow region capping both the DBL and TBL. As previously observed in FS09, such features are sensitive to the sloping angle $(\unicode[STIX]{x1D6FC})$ , and as $\unicode[STIX]{x1D6FC}$ decreases from the vertical wall set-up ( $\unicode[STIX]{x1D6FC}=90^{\circ }$ ), the anabatic and katabatic flow solutions progressively depart from each other.

Figure 6.

Comparison of streamwise mean velocity $\langle u\rangle$ (a) and $\langle u^{+}\rangle$ (c) and of mean buoyancy $\langle b\rangle$ (b) and $\langle b^{+}\rangle$ (d) for anabatic (dashed lines) and katabatic (solid lines) flows. Symbols: $\unicode[STIX]{x1D6FC}=90^{\circ }$ , black lines; $\unicode[STIX]{x1D6FC}=60^{\circ }$ , blue lines; $\unicode[STIX]{x1D6FC}=30^{\circ }$ , green lines; $\unicode[STIX]{x1D6FC}=15^{\circ }$ , red lines. All cases are characterized by $\mathit{Gr}=2.1\times 10^{11}$ . The $z$ -axis denotes the slope-normal coordinate direction.

The constant-K Prandtl solution provides a useful framework for the interpretation of results, and is thus here briefly summarized. Based on the proposed normalization (see Sect. 2.1), the Prandtl one-dimensional solution for imposed constant surface buoyancy reads (FS09)

where $\unicode[STIX]{x1D70E}=(\mathit{Gr}\,Pr^{-1})^{1/4}\sin (\unicode[STIX]{x1D6FC})^{1/2}$ and $b_{s}=\pm 1$ . Note that (4.3) and (4.4) correspond to a laminar flow solution of the time and space averaged (2.9)–(2.11). It predicts a velocity maximum $u_{j}=\pm 1/\sqrt{2}\exp (-\unicode[STIX]{x03C0}/4)$ that is independent of the sloping angle $(\unicode[STIX]{x1D6FC})$ , whereas the characteristic length scale of the flow $L\propto \sin (\unicode[STIX]{x1D6FC})^{-1/2}$ whereby $z_{j}\propto \sin (\unicode[STIX]{x1D6FC})^{-1/2}$ (Grisogono & Axelsen Reference Grisogono and Axelsen2012; Grisogono et al. Reference Grisogono, Jurlina, Večenaj and Güttler2014).

The location of the LLJ for the katabatic flow DNS solution is in good agreement with prediction from the Prandtl laminar flow solution, i.e.

Conversely, $u_{j}$ is significantly smaller (a reduction from ${\approx}30\,\%$ to ${\approx}50\,\%$ depending on the slope angle) when compared to the laminar solution, mainly due to additional diffusion of momentum caused by turbulent motions at $z_{j}$ (as shown in § 4.5). In addition, $u_{j}$ is not independent of $\unicode[STIX]{x1D6FC}$ as in the Prandtl model, but is characterized by a modest increase as $\unicode[STIX]{x1D6FC}$ decreases. This behaviour can be explained by considering that as the sloping angle decreases, the stable stratification induced by the imposed surface buoyancy damps turbulent motions in the LLJ region, thus reducing turbulent mixing of momentum and resulting in a higher peak velocity $u_{j}$ .

The anabatic flow solution is more sensitive to variations in the sloping angle, when compared to its katabatic counterpart. As $\unicode[STIX]{x1D6FC}$ is reduced, a simultaneous increase in the height of the LLJ and a reduction in its peak speed are observed. Variations are significant when compared to those characterizing the katabatic solution and the laminar case. This pattern is related to the strengthening of the slope-normal component of the buoyancy force as $\unicode[STIX]{x1D6FC}$ decreases ( $b\unicode[STIX]{x1D6FF}_{i3}\cos \unicode[STIX]{x1D6FC}$ in (2.9)). This results in an unstable near-wall stratification, which enhances the slope-normal flux of momentum, and leads to progressively more mixed profiles of velocity and buoyancy in the LLJ regions, in agreement with the findings of FS09.

Figure 6 also features mean velocity and buoyancy profiles scaled with inner units, i.e. $z^{+}=\hat{z}\hat{u} _{\star }/\hat{\unicode[STIX]{x1D708}},\,u^{+}=\hat{u} /\hat{u} _{\star },\,b^{+}=\hat{b}/\hat{b}_{\star }$ (where $\hat{u} _{\star }\equiv \sqrt{\langle \unicode[STIX]{x1D70F}_{xz}^{tot}\rangle |_{z=0}}$ and $\hat{b}_{\star }\equiv (\langle \unicode[STIX]{x1D70F}_{bz}^{tot}\rangle |_{z=0})/u_{\star }$ , with $\unicode[STIX]{x1D70F}_{xz}^{tot}$ and $\unicode[STIX]{x1D70F}_{bz}^{tot}$ being the total kinematic momentum and buoyancy fluxes). The nose of the LLJ is located at $z^{+}\approx 25$ in katabatic flows, separating the viscosity dominated wall regions from the turbulent outer layers. Anabatic flow solutions are characterized by a progressively more mixed LLJ as the sloping angle decreases, with peak speed at $z^{+}\approx 100$ for $\unicode[STIX]{x1D6FC}=15^{\circ }$ , again highlighting a broadening of the DBL as the slope is reduced.

Another notable difference between the katabatic and the anabatic flow solutions is the sensitivity of the DBL and TBL thickness to $\unicode[STIX]{x1D6FC}$ . The DBL and TBL thickness ( $\unicode[STIX]{x1D6FF}_{d}$ and $\unicode[STIX]{x1D6FF}_{t}$ respectively), which can be identified as the first and second zero crossing of $\langle b\rangle$ and $\langle u\rangle$ respectively are insensitive to $\unicode[STIX]{x1D6FC}$ for the katabatic flow regime, whereas they vary by a factor of three across the considered $\unicode[STIX]{x1D6FC}$ -range for the anabatic flow solution. In line with this finding, the slope-normal integrated horizontal momentum flux $(I_{u})$ in the anabatic flow regime is also strongly sensitive to variations in $\unicode[STIX]{x1D6FC}$ , whereas it is insensitive in the katabatic cases.

Figure 7.

$\unicode[STIX]{x1D6FC}$ dependence of the surface buoyancy flux $B_{w}$ for katabatic (a) and anabatic (b) flow cases. Simulations correspond to cases s90h, s60h, s30h, s15h and u90h, u60h, u30h, u15h for the katabatic and anabatic regimes, respectively. Predictions from the Prandtl laminar solution are included for comparison (black lines).

Such behaviour can be understood by integrating (4.2),

where $B_{w}$ is the surface buoyancy flux. From (4.6) is apparent how larger $B_{w}$ and shallower slopes are characterized by larger $I_{u}$ . Variations of the surface buoyancy flux as a function of $\unicode[STIX]{x1D6FC}$ for the considered anabatic and katabatic cases are displayed in figure 7. As apparent, the anabatic flow solution is characterized by modest variations in $B_{w}$ across the considered $\unicode[STIX]{x1D6FC}$ range, when compared against the katabatic flow cases. The latter experience a significant decrease in the surface buoyancy flux as the sloping angle decreases from $90^{\circ }$ to $15^{\circ }$ , slowly approaching the surface flux predicted by the laminar flow solution (at $\unicode[STIX]{x1D6FC}=15^{\circ }$ the difference between the DNS surface flux and the laminar solution surface flux is only $20\,\%$ the magnitude of the DNS surface flux itself). Such a behaviour is closely related to the strengthening of the near-surface inversion layer, which damps turbulent motions and leads to an apparent laminarization of the flow. In contrast, anabatic flows depart from the laminar solution as $\unicode[STIX]{x1D6FC}$ decreases, due to the unstable stratification and likely onset of convective motions. Based on results shown in figure 7 one could introduce the equally crude approximations

which justifies the observed variations in the vertically integrated horizontal momentum flux in the anabatic and katabatic flow regimes based on (4.6).

Note that since $I_{u}={\mathcal{U}}{\mathcal{L}}$ , where ${\mathcal{U}}$ and ${\mathcal{L}}$ are the characteristic (normalized) velocity and length scale of the flow, and given that ${\mathcal{U}}\approx const.$ for both anabatic and katabatic flows, from (4.7) and (4.8) follows that ${\mathcal{L}}\propto 1/\sin \unicode[STIX]{x1D6FC}$ for anabatic flows, and ${\mathcal{L}}\approx const.$ for katabatic flows. Such information can be of use for scaling purposes.

Given the strong $\unicode[STIX]{x1D6FC}$ dependence of $B_{w}$ in the katabatic flow solution, profiles are expected vary significantly – when compared to those proposed herein – if a constant $B_{w}$ is used to drive the flow. Conversely, $B_{w}$ in the anabatic flow solution is weakly dependent on $\unicode[STIX]{x1D6FC}$ , and profiles are therefore expected to be poorly sensitive on the specific flow driving mechanisms (i.e. constant surface buoyancy or constant surface buoyancy flux). Because of this, the proposed anabatic flow solutions share strong similarities with those reported in the FS09 study, especially when considering the $\unicode[STIX]{x1D6FC}$ dependence of $z_{j},u_{j}$ and $\unicode[STIX]{x1D6FF}_{d},\unicode[STIX]{x1D6FF}_{t}$ , whereas the katabatic flow solutions differ significantly.

When a constant surface buoyancy flux is applied, the laminar flow solution is characterized by $z_{j}=\unicode[STIX]{x03C0}/(4\unicode[STIX]{x1D70E})$ and $u_{j}\propto (\sin (\unicode[STIX]{x1D6FC}))^{-1/2}$ (Grisogono & Axelsen Reference Grisogono and Axelsen2012; Grisogono et al. Reference Grisogono, Jurlina, Večenaj and Güttler2014). Hence, $u_{j}$ increases as the sloping angle decreases and is not $\unicode[STIX]{x1D6FC}$ independent as predicted by (4.3) and (4.4).

Similar variations in $z_{j}$ and a stronger $\unicode[STIX]{x1D6FC}$ dependency of $u_{j}$ and $\unicode[STIX]{x1D6FF}_{d},\unicode[STIX]{x1D6FF}_{t}$ have been reported in FS09 for the katabatic flow solution when compared to current DNS results, qualitatively resembling predictions of the analytic laminar flow solution.

Note that when the flow is forced through a constant surface buoyancy flux, as in the FS09 study, anabatic and katabatic flow solutions are constraint to share the same slope-normal integrated horizontal momentum flux $I_{u}^{\unicode[STIX]{x1D6FC}}/I_{u}^{\unicode[STIX]{x1D6FC}=90^{\circ }}=1/\sin \unicode[STIX]{x1D6FC}$ , hence the stronger variability of $\unicode[STIX]{x1D6FF}_{d}$ and $\unicode[STIX]{x1D6FF}_{t}$ in the FS09 katabatic flow solution, when compared to the current DNS results.

4.3 TKE and buoyancy variance

Slope-normal variations of TKE and of buoyancy variance $\langle b^{\prime }b^{\prime }\rangle$ are featured in figure 8. The buoyancy variance may also be interpreted as a form of turbulent potential energy. The stronger surface stable stratification that characterizes katabatic flows as the sloping angle decreases results in a weakening of TKE in the inner flow regions (below the LLJ), whereas the observed increase of TKE in the outer flow region as $\unicode[STIX]{x1D6FC}$ decreases is likely related to the broadening of the flow length scales, which despite being modest for the katabatic flow solution, have an apparent effect on TKE. TKE profiles from the anabatic flow solution are again more sensitive to $\unicode[STIX]{x1D6FC}$ when compared to their katabatic counterparts. In the anabatic flow cases, the peak TKE location increases as $\unicode[STIX]{x1D6FC}$ decreases, but its magnitude shows non-monotonic behaviour, thus suggesting a more complex dependence on $\unicode[STIX]{x1D6FC}$ . Furthermore, the TKE in the neighbourhood of the LLJ is approximately constant throughout the considered flow regimes.

Figure 8.

Comparison of turbulent kinetic energy $(1/2)\langle u_{i}^{\prime }u_{i}^{\prime }\rangle$ (a) and buoyancy variance $(\langle b^{\prime }b^{\prime }\rangle )$ (b) for the katabatic (solid lines) and the anabatic flow (dashed lines) regimes at $\unicode[STIX]{x1D6FC}=90^{\circ }$ (black), $\unicode[STIX]{x1D6FC}=60^{\circ }$ (blue), $\unicode[STIX]{x1D6FC}=30^{\circ }$ (green), and $\unicode[STIX]{x1D6FC}=15^{\circ }$ (red). Simulations correspond to the highest $\mathit{Gr}=2.1\times 10^{11}$ cases. Recall that the $z$ -axis denotes the slope-normal coordinate direction.

Figure 9.

Normal stress components $\langle u^{\prime }u^{\prime }\rangle$ (solid lines), $\langle v^{\prime }v^{\prime }\rangle$ (dashed lines) and $\langle w^{\prime }w^{\prime }\rangle$ (dot-dashed lines) for the katabatic (a) and the anabatic (b) flow regimes at $\unicode[STIX]{x1D6FC}=90^{\circ }$ (black), $\unicode[STIX]{x1D6FC}=60^{\circ }$ (blue), $\unicode[STIX]{x1D6FC}=30^{\circ }$ (green) and $\unicode[STIX]{x1D6FC}=15^{\circ }$ (red). Simulations correspond to the highest $\mathit{Gr}=2.1\times 10^{11}$ cases (simulations s90h, s60h, s30h, and s15h and u90h, u60h, u30h and u15h, respectively).

The buoyancy variance $\langle b^{\prime }b^{\prime }\rangle$ peaks in the near-wall regions for both flow regimes where strong buoyancy gradients occur, in agreement with findings from FS09. Variations in $\langle b^{\prime }b^{\prime }\rangle$ as a function of $\unicode[STIX]{x1D6FC}$ in the below-LLJ region are significant only for the katabatic flow regime, with peak value and its location increasing and decreasing as $\unicode[STIX]{x1D6FC}$ is reduced. The above-LLJ regions of the boundary layer are characterized by a rapid decay in $\langle b^{\prime }b^{\prime }\rangle$ , most evident for the katabatic flow regime. The anisotropic nature of turbulence in slope flows is apparent from figure 9, where normal stress components $\langle u^{\prime }u^{\prime }\rangle$ , $\langle v^{\prime }v^{\prime }\rangle$ and $\langle w^{\prime }w^{\prime }\rangle$ are compared. The boundary layer character of the system is apparent with the wall providing an effective damping of the $\langle w^{\prime }w^{\prime }\rangle$ central moment, in both anabatic and katabatic flow regimes. It is to be noted the strong sensitivity of $\langle w^{\prime }w^{\prime }\rangle$ with respect to $\unicode[STIX]{x1D6FC}$ for both wind regimes. This behaviour is related to the direct effect of stratification (background $+$ perturbation) on $\langle w^{\prime }w^{\prime }\rangle$ , given that buoyancy is effective at damping/exciting slope-normal velocity fluctuations  $w^{\prime }$ . Because of this, as $\unicode[STIX]{x1D6FC}$ decreases, the turbulence characterizing katabatic flows becomes more anisotropic (the strong, effective, stable stratification damps $\langle w^{\prime }w^{\prime }\rangle$ ) in contrast to its anabatic counterpart, where the positive surface buoyancy leads to more isotropic turbulent motions. The observed trend here supports the recently proposed scaling of Shapiro & Fedorovich (Reference Shapiro and Fedorovich2014) based on the assumption of large-scale separation between slope-normal and slope-parallel motions populating katabatic flows, which might indeed hold at small sloping angles.

Figure 10.

Total (solid lines) and turbulent (dashed lines) momentum flux for the katabatic (a) and the anabatic (b) flow regimes, and total (solid lines) and turbulent (dashed lines) buoyancy slope-normal flux for the katabatic (c) and the anabatic (d) flow regimes. All cases are characterized by $\mathit{Gr}=2.1\times 10^{11}$ . $\langle \unicode[STIX]{x1D70F}_{xz}^{tot}\rangle$ denotes the total slope-normal flux of momentum, whereas $\langle \unicode[STIX]{x1D70F}_{bz}^{tot}\rangle$ denotes the total slope-normal buoyancy flux. The height of the LLJ ( $z_{j}$ ) is displayed (horizontal lines) for the different cases to provide a reference.

4.4 Momentum and buoyancy fluxes

Turbulent and total (turbulent $+$ molecular) vertical buoyancy and momentum fluxes are displayed in figure 10 for each of the considered runs. Turbulent buoyancy fluxes obey Fick’s law of diffusion, whereas momentum fluxes consistently feature a thin layer just below the LLJ where counter gradient fluxes are observed. This feature is related to the extrema of the mean $\langle u\rangle$ -profiles not being co-located with the zero crossing of the turbulent momentum fluxes. At the LLJ location for instance, a positive (negative) total slope-normal momentum flux for the katabatic (anabatic) flow regime are observed. Note that counter gradient fluxes have been previously reported in experimental studies (Smeets et al. Reference Smeets, Duynkerke and Vugts2000; Oldroyd et al. Reference Oldroyd, Pardyjak, Higgins and Parlange2016a ). The peak magnitude of both $\langle \unicode[STIX]{x1D70F}_{xz}^{tot}\rangle$ and $\langle \unicode[STIX]{x1D70F}_{bz}^{tot}\rangle$ in the above-LLJ regions is dependent on $\unicode[STIX]{x1D6FC}$ for both flow regimes (it decreases as $\unicode[STIX]{x1D6FC}$ decreases). In addition, katabatic flow solutions are characterized by a modest upward shift of flux extrema as the sloping angle decreases, whereas the upward shift is more significant for the anabatic flow cases. Such trends are directly related to the combined effects of $\unicode[STIX]{x1D6FC}$ on the scales of the flow (recall that the analytic Prandtl solution predicts $L\propto (\sin (\unicode[STIX]{x1D6FC}))^{-1/2}$ ) and of the surface induced stratification in katabatic and anabatic flows, whose strength increases as $\unicode[STIX]{x1D6FC}$ decreases, progressively damping (for the katabatic cases) or enhancing (for the anabatic cases) turbulent fluctuations. In the FS09 study, the peak magnitude of $\langle \unicode[STIX]{x1D70F}_{xz}^{tot}\rangle$ and $\langle \unicode[STIX]{x1D70F}_{bz}^{tot}\rangle$ in the katabatic flow regime was found to be approximately constant and independent of $\unicode[STIX]{x1D6FC}$ . Conversely, it was sensitive to $\unicode[STIX]{x1D6FC}$ in the anabatic flow cases, in apparent contrast with current findings. Such a mismatch is related to the different surface forcing approach that characterizes the current and the FS90 study (in FS90 a constant surface buoyancy flux is applied to drive the flow).

Figure 11.

Total (solid lines) and turbulent (dashed lines) momentum flux for the katabatic (a) and the anabatic (b) flow regimes in inner units. All cases are characterized by $\mathit{Gr}=2.1\times 10^{11}$ . $\langle \unicode[STIX]{x1D70F}_{xz}^{tot}\rangle$ denotes the total slope-normal flux of momentum. The height of the LLJ ( $z_{j}$ ) is displayed (horizontal lines) for the different cases to provide a reference.

Figure 11 features the total and the turbulent momentum fluxes for the considered katabatic and anabatic flow regimes in inner units. For katabatic flows the viscous contribution to the total flux is dominant in the below-LLJ regions, and drops to less than $10\,\%$ as $z^{+}>50$ for the $\unicode[STIX]{x1D6FC}=90^{\circ }$ case. As the sloping angle decreases the viscous contribution in the above-LLJ region increases, due to the strengthening of the stable stratification resulting from the imposed surface buoyancy, and amounts to approximately $35\,\%$ the total at $z^{+}=50$ for the shallower of the considered slopes $(\unicode[STIX]{x1D6FC}=15^{\circ })$ . In contrast, anabatic flows are characterized by a decrease of viscous contributions to the total momentum flux as the sloping angle decreases throughout the boundary layer, with viscous effects remaining dominant up to $z^{+}\approx 12$ , where viscous and turbulent fluxes are roughly of the same magnitude, and becoming negligible $({<}10\,\%)$ as $z^{+}>50$ . From figure 11 is also apparent how the surface contribution in terms of total stress is larger than the peak stress in the above-LLJ region for katabatic flows, whereas for anabatic flows as the sloping angle decreases the surface stress becomes progressively smaller when compared to the outer stress peak.

Figure 12.

$\unicode[STIX]{x1D6FC}$ dependence of the kinematic surface average stress $\unicode[STIX]{x1D70F}_{w}$ for katabatic (a) and anabatic (b) flow cases. Simulations correspond to cases s90h, s60h, s30h, s15h and u90h, u60h, u30h, u15h for the katabatic and anabatic regimes, respectively. Predictions from the Prandtl laminar solution are included for comparison (black lines).

Turbulent fluctuations contribute to the overall buoyancy flux in the below-LLJ region of katabatic flows, whereas they provide a negligible contribution to the overall momentum flux. This latter behaviour has direct effect on the overall surface stress, which is nearly equal to the laminar flow prediction for $\unicode[STIX]{x1D6FC}>90^{\circ }$ , as displayed in figure 12. Turbulent motion in the below-jet regions of katabatic flows are thus ‘inactive’ from a momentum transport perspective (in the Townsend (Reference Townsend1956) sense), but are relatively effective in vertically transporting buoyancy. To the contrary, turbulent fluctuations in the anabatic flow solution contribute to both the overall surface buoyancy and momentum fluxes that progressively depart from their Prandtl ‘laminar flow’ counterparts as $\unicode[STIX]{x1D6FC}$ is reduced. An interesting feature of the considered anabatic flow simulations is the magnitude of the time- and space-averaged surface total momentum flux ( $\unicode[STIX]{x1D70F}_{w}$ ), as displayed in figure 12. $\unicode[STIX]{x1D70F}_{w}$ is consistently smaller than its laminar flow prediction in the anabatic cases, and its magnitude decreases as the sloping angle is reduced. Note that such feature might be related to the relatively modest $Gr$ value characterizing the various runs.

4.5 The mean kinetic energy budget

The governing equation for the time- and space-averaged velocity $\langle u_{i}\rangle$ is readily derived by applying Reynolds decomposition to the instantaneous flow variables (e.g. $u_{i}=\langle u_{i}\rangle +u_{i}^{\prime }$ ) and space $+$ time averaging (2.9). The budget equation for MKE is then obtained by multiplying the equation for $\langle u_{i}\rangle$ by $\langle u_{i}\rangle$ itself. Assuming horizontal homogeneity $(\unicode[STIX]{x2202}\langle \cdot \rangle /\unicode[STIX]{x2202}x=\unicode[STIX]{x2202}\langle \cdot \rangle /\unicode[STIX]{x2202}y=0)$ and no subsidence ( $\langle w\rangle =0$ ) it reads

where the left-hand side of (4.9) is the storage term of MKE, ${\mathcal{P}}_{s}\equiv \langle u^{\prime }w^{\prime }\rangle (\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z)$ denotes shear production/destruction of MKE, ${\mathcal{P}}_{b}\equiv -\langle u\rangle \langle b\rangle \sin (\unicode[STIX]{x1D6FC})$ denotes buoyancy production/destruction of MKE, transport of MKE by turbulent motions is ${\mathcal{T}}_{t}\equiv -\unicode[STIX]{x2202}(1/2\langle u\rangle \langle u^{\prime }w^{\prime }\rangle )/\unicode[STIX]{x2202}z$ and dissipation of MKE by viscous diffusion is ${\mathcal{E}}\equiv \mathit{Gr}^{-1/2}\langle u\rangle (\unicode[STIX]{x2202}^{2}\langle u\rangle /\unicode[STIX]{x2202}z^{2})$ . When time averaging over a sufficiently long time period then $\unicode[STIX]{x2202}\langle \cdot \rangle /\unicode[STIX]{x2202}t=0$ and the storage term can be neglected.

The normalized MKE budget terms for the considered anabatic and katabatic runs are displayed in figure 13. The choice of $\hat{b_{s}}^{2}\hat{N}^{-1}$ as a normalizing factor is not critical for the interpretation of the budget, since the relative magnitude of the terms is unchanged. As expected, the overall main source of MKE is from buoyancy production $({\mathcal{P}}_{b})$ , which peaks in the below-jet regions, and is characterized by a gradual decrease throughout the boundary layer. In the outer regions of the flow, ${\mathcal{P}}_{b}$ becomes a sink of MKE in both flow regimes starting from the zero crossing of $\langle b\rangle$ and up to the start of the return flow region. Here, energy is provided by turbulent transport $({\mathcal{T}}_{t})$ , which balances shear production $({\mathcal{P}}_{s})$ and buoyant production $({\mathcal{P}}_{b})$ (both ${\mathcal{P}}_{s}$ and ${\mathcal{P}}_{b}$ are a sink term of MKE in such layer). At the wall, buoyant production is overcome by dissipation for both upslope and downslope flows, and transport from turbulent motions is responsible to close the MKE budget. ${\mathcal{T}}_{t}$ acts as a sink of MKE in the highly energetic LLJ regions, displacing it toward the wall to balance the enhanced dissipation, and also into the outer layer of the flow.

Figure 13.

Slope-normal structure of the MKE budget for the katabatic (a) and for the anabatic (b) flow regimes at $\unicode[STIX]{x1D6FC}=90^{\circ }$ (solid lines), $\unicode[STIX]{x1D6FC}=60^{\circ }$ (dashed lines), $\unicode[STIX]{x1D6FC}=30^{\circ }$ (dot-dashed lines) and $\unicode[STIX]{x1D6FC}=15^{\circ }$ (dotted lines). All cases are characterized by $\mathit{Gr}=2.1\times 10^{11}$ . The location of the LLJ is highlighted by horizontal lines for the various runs to provide a reference height (note that as $\unicode[STIX]{x1D6FC}$ decreases the LLJ height increases). All terms are normalized by $\hat{U} ^{3}\,\hat{L}^{-1}\equiv \hat{b_{s}}^{2}\,\hat{N}^{-1}$ .

In both anabatic and katabatic flow regimes, shear production of MKE $({\mathcal{P}}_{s})$ acts as a sink of MKE in the above-jet regions, draining energy from the mean flow and transferring it to turbulence through the classical energy cascade process. Interestingly, for both regimes and all the considered sloping angles, the below-jet regions are characterized by ${\mathcal{P}}_{s}>0$ , highlighting a region of global energy backscatter, i.e. energy is transferred from the turbulent eddies to the mean flow. Forward scatter is known to be mainly caused by vortex stretching by the mean strain rate, whereas backscatter indicates vortex compression by the mean strain rate, which is not commonly observed in canonical wall-bounded flows.

4.6 The turbulent kinetic energy budget

Under the assumptions leading to (4.9), the budget equation for TKE is given as

where $\unicode[STIX]{x2202}\langle (1/2)u_{i}^{\prime }u_{i}^{\prime }\rangle /\unicode[STIX]{x2202}t$ is the storage of TKE term, shear production of TKE is denoted as $P_{s}\equiv -\langle u^{\prime }w^{\prime }\rangle (\unicode[STIX]{x2202}\langle u\rangle /\unicode[STIX]{x2202}z)$ , buoyant production/destruction of TKE is composed of two terms, namely $P_{b,1}\equiv \langle b^{\prime }u^{\prime }\rangle \sin (\unicode[STIX]{x1D6FC})$ and $P_{b,3}\equiv \langle b^{\prime }w^{\prime }\rangle \cos (\unicode[STIX]{x1D6FC})$ , turbulent transport of TKE is $T_{t}\equiv -\unicode[STIX]{x2202}\langle (1/2)u_{i}^{\prime }u_{i}^{\prime }w^{\prime }\rangle /\unicode[STIX]{x2202}z$ , pressure transport $T_{p}\equiv -\unicode[STIX]{x2202}\langle \unicode[STIX]{x03C0}^{\prime }w^{\prime }\rangle /\unicode[STIX]{x2202}z$ , viscous diffusion of TKE is $T_{\unicode[STIX]{x1D708}}\equiv \mathit{Gr}^{-1/2}\unicode[STIX]{x2202}^{2}\langle (1/2)u_{i}^{\prime }u_{i}^{\prime }\rangle /\unicode[STIX]{x2202}z^{2}$ and viscous dissipation $\unicode[STIX]{x1D716}\equiv -\mathit{Gr}^{-1/2}\langle (\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j})(\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j})\rangle$ . With regard to the buoyancy production/destruction terms, $P_{b,1}$ accounts for production/destruction of TKE due to cross-correlation between along-slope velocity ( $u$ ) and buoyancy ( $b$ ), whereas $P_{b,3}$ accounts for production/destruction of TKE due to cross-correlation between normal-to-slope velocity ( $w$ ) and buoyancy ( $b$ ). The splitting of the buoyancy production term is clearly a result of the inclined reference system that is adopted to describe the evolution of the system.

Figure 14.

Comparison of TKE budged terms for katabatic (a,c) and anabatic (b,d) flow regimes at $\unicode[STIX]{x1D6FC}=90^{\circ }$ (solid lines), $\unicode[STIX]{x1D6FC}=60^{\circ }$ (dashed lines), $\unicode[STIX]{x1D6FC}=30^{\circ }$ (dot-dashed lines) and $\unicode[STIX]{x1D6FC}=15^{\circ }$ (dotted lines). All cases are characterized by $\mathit{Gr}=2.1\times 10^{11}$ . Production and destruction terms (a,b) have been separated from transport and residual terms (c,d). The $z$ -axis represents the slope-normal coordinate direction. The location of the LLJ is highlighted by horizontal lines for the various runs to facilitate interpretation (note that $\unicode[STIX]{x1D6FC}\propto z_{j}$ ). All terms are normalized by $\hat{U} ^{3}\hat{L}^{-1}\equiv \hat{b}_{s}^{2}\hat{N}^{-1}$ .

Figure 15.

Comparison of return-to-isotropy terms for katabatic (a) and anabatic (b) flow regimes at $\mathit{Gr}=2.1\times 10^{11}$ . We denote $\unicode[STIX]{x1D6F7}_{1}\equiv \langle p^{\prime }(\unicode[STIX]{x2202}u^{\prime }/\unicode[STIX]{x2202}x)\rangle ,\unicode[STIX]{x1D6F7}_{2}\equiv \langle p^{\prime }(\unicode[STIX]{x2202}v^{\prime }/\unicode[STIX]{x2202}y)\rangle$ and $\unicode[STIX]{x1D6F7}_{3}\equiv \langle p^{\prime }(\unicode[STIX]{x2202}w^{\prime }/\unicode[STIX]{x2202}z)\rangle$ . The location of the LLJ is once again highlighted with horizontal lines and the $\unicode[STIX]{x1D6FC}=90^{\circ },\,60^{\circ },\,30^{\circ }$ and $15^{\circ }$ runs (simulations s90h, s60h, s30h, s15h for the katabatic regimes; u90h, u60h, u30h, u15h for the anabatic regimes) are denoted with solid, dashed, dot-dashed and dotted lines respectively. The $z$ -axis represents the slope-normal coordinate direction and all terms are normalized by $\hat{U} ^{3}\,\hat{L}^{-1}\equiv \hat{b}_{s}^{2}\,\hat{N}^{-1}$ .

TKE budget terms for the considered runs (characterized by $\mathit{Gr}=2.1\times 10^{11}$ ) are displayed in figure 14. Shear production $(P_{s})$ appears with opposite signs in the budgets of MKE and TKE as expected. It represents the net transfer from MKE to TKE as the result of their interactions that often sustains turbulence in classical boundary layer theory on flat slopes. For both anabatic and katabatic flow regimes, $P_{s}$ is characterized by two positive peaks, one in the above jet regions and one in the very near-wall regions, and is negative in a small interval just below the LLJ, where global energy backscatter occurs. Occurrence of negative $P_{s}$ is related to the presence of local (in $z$ ) counter-gradient turbulent momentum flux, in line with findings of § 4.4. In both katabatic and anabatic flow regimes, $\unicode[STIX]{x1D716}$ peaks at the wall, is approximately constant in the near-LLJ region and decreases logarithmically in the core of the flow. The $P_{b,3}$ is a sink of TKE for the katabatic regime and a source of TKE for the anabatic regime, as expected. In the anabatic regime $P_{b,3}=0$ at $\unicode[STIX]{x1D6FC}=90^{\circ }$ , but gains considerable importance (as a TKE source term) in the overall budget as $\unicode[STIX]{x1D6FC}$ decreases. For instance, considering the $\unicode[STIX]{x1D6FC}=15^{\circ }$ run, $P_{b,3}$ alone overcomes TKE dissipation in the core of the LLJ. To the contrary, the modest magnitude of $P_{b,3}$ highlights how buoyant destruction of TKE is not the primary mechanism through which buoyancy acts to suppress turbulence in katabatic flows. Following the same reasoning of Shah & Bou-Zeid (Reference Shah and Bou-Zeid2014) (where stability effects on the Ekman layer were studied through DNS), it is argued here that negative buoyancy directly reduces $\langle w^{\prime }w^{\prime }\rangle$ , thus reducing local production of $\langle u^{\prime }w^{\prime }\rangle$ . A reduction in $\langle u^{\prime }w^{\prime }\rangle$ would ultimately result in the observed decrease in $\langle P_{s}\rangle$ and related TKE magnitude as $\unicode[STIX]{x1D6FC}$ decreases (i.e. as the static stability of the environment increases). $P_{b,1}$ is the major source of TKE at the LLJ for the katabatic flow regime at all the considered $\unicode[STIX]{x1D6FC}$ . On the other hand, in the anabatic flow regime $P_{b,3}$ overcomes $P_{b,1}$ as $\unicode[STIX]{x1D6FC}$ decreases, becoming the leading buoyant production term. Overall, the sum of production terms ( $P_{s}+P_{b,1}+P_{b,3}$ ) overcome dissipation in the above-jet regions (approximately up to $10z_{j}$ ), and transport terms are responsible for dislocating this excess in TKE down towards the wall and towards the outer regions of the flow. Turbulent transport ( $T_{t}$ ) is a more effective carrier of TKE in the outer regions of the flow, whereas pressure fluctuations $(T_{p})$ are more effective in transporting TKE down toward the wall, to balance dissipation and viscous diffusion. The viscous diffusion term $T_{\unicode[STIX]{x1D708}}$ resembles its pressure-driven boundary layer analogue, where $T_{\unicode[STIX]{x1D708}}$ is a sink of TKE in the buffer sublayer, and a source of TKE in the laminar sublayer, below $z^{+}=5$ (corresponding to $z=5\times 10^{-4}$ in current units).

Not shown here is the vertical structure of flux Richardson number $(Ri_{f}=\langle b^{\prime }w^{\prime }\rangle /(\text{d}\langle u\rangle /\text{d}z)\langle u^{\prime }w^{\prime }\rangle )$ , which is positive (negative) throughout the core of the boundary layer (above the LLJ and below the return flow region) in the katabatic (anabatic) flow regimes. In the katabatic flow regime $Ri_{f}$ is lower in magnitude than its critical value of 0.25, except in the neighbourhood of the LLJ the stable stratification and low velocity gradients result instead in $Ri_{f}\gg 1$ .

It is worth noting the significant contribution of turbulent and pressure transport terms in the neighbourhood of the LLJ and in the wall regions in both flow regimes. Pressure fluctuations are relevant in the near-wall regions when compared to corresponding values observed in neutral (Moser, Kim & Moin Reference Moser, Kim and Moin1999) and stably stratified (Iida, Kasagi & Nagano Reference Iida, Kasagi and Nagano2002) pressure-driven channel flow DNS. In the neighbourhood of the LLJ $Ri_{f}$ is relatively large in magnitude (not shown) for both flow regimes, promoting strong internal gravity wave activity, which does not transfer buoyancy, but can be effective in transporting momentum through the action of pressure force, thus justifying such relevant contribution of pressure transport to the TKE budget.

The return-to-isotropy term (also known as pressure redistribution term) contracts to zero, and so vanishes from the TKE budget equation (4.10). However, when analysed for the single TKE budget components displayed in figure 15, the role of turbulence in distributing turbulent energy becomes clear. For instance, wall damping effects on slope-normal velocity fluctuations are apparent in the below-LLJ regions, where $\unicode[STIX]{x1D6F7}_{3}<0$ indicating energy redistribution from the slope-normal component $(\langle w^{\prime }w^{\prime }\rangle )$ to the horizontal components ( $\langle u^{\prime }u^{\prime }\rangle$ and $\langle v^{\prime }v^{\prime }\rangle$ respectively). In the above-jet regions for the katabatic flow regime, a consistent energy redistribution among the TKE components are observed across the sloping angles with energy being transferred from the streamwise component $(\langle u^{\prime }u^{\prime }\rangle )$ to the spanwise and slope-normal components ( $\langle v^{\prime }v^{\prime }\rangle$ and $\langle w^{\prime }w^{\prime }\rangle$ respectively). For the anabatic flow regime, the return-to-isotropy terms in the above-jet regions highlight a transition in the system as a function of $\unicode[STIX]{x1D6FC}$ . When the two highest sloping angles are considered ( $\unicode[STIX]{x1D6FC}=60^{\circ }$ and $\unicode[STIX]{x1D6FC}=90^{\circ }$ ), energy transfer is qualitatively equivalent to that characterizing the katabatic flow regime, i.e. the streamwise variance feeds the spanwise and slope-normal variance components. For $\unicode[STIX]{x1D6FC}=15^{\circ }$ and $\unicode[STIX]{x1D6FC}=30^{\circ }$ , the return-to-isotropy term becomes a sink for $\langle w^{\prime }w^{\prime }\rangle$ and a source for $\langle u^{\prime }u^{\prime }\rangle$ and $\langle v^{\prime }v^{\prime }\rangle$ , indicative of energy transfer from the slope-normal TKE component to the streamwise and spanwise TKE components. This transition suggests that at low sloping angles, anabatic flow regimes are characterized by slope-normal elongated eddies as apparent from figure 5, which feed $\langle u^{\prime }u^{\prime }\rangle$ and $\langle v^{\prime }v^{\prime }\rangle$ from $\langle w^{\prime }w^{\prime }\rangle$ , the latter being directly sustained by the slope-normal projection of the buoyancy force $b\cos (\unicode[STIX]{x1D6FC})$ . Conversely, katabatic flow eddies are streamwise elongated and remove energy from $\langle u^{\prime }u^{\prime }\rangle$ – directly fed by the streamwise component of the buoyancy force – to transfer it to $\langle w^{\prime }w^{\prime }\rangle$ and $\langle v^{\prime }v^{\prime }\rangle$ . This energy redistribution term ensures some self-preservation of the slope-normal velocity variance in katabatic flows despite the adverse role of stability.

Overall, the proposed TKE budget analysis suggests a subdivision of the boundary layer into four distinct regions, for the considered $\unicode[STIX]{x1D6FC}$ - and $Gr$ -ranges, namely

1. (i) an outer layer, corresponding approximately to the return flow region, where turbulent transport $(T_{t})$ is the main source of TKE and balances dissipation $(\unicode[STIX]{x1D716})$ ;

2. (ii) an intermediate layer, bounded below by the LLJ and capped above by the outer layer, where the sum of shear and buoyant production ( $P_{s}+P_{b,1}+P_{b,3}$ ) overcomes dissipation $(\unicode[STIX]{x1D716})$ , and where turbulent and pressure transport terms ( $T_{t},T_{p}$ ) are a sink of TKE;

3. (iii) a buffer layer, corresponding to $5\lessapprox z^{+}\lessapprox 30$ , where TKE is provided by turbulent and pressure transport terms, to balance viscous diffusion and dissipation;

4. (iv) a laminar sublayer, corresponding to $z^{+}\lessapprox 5$ , where the influence of viscosity is significant and the flow is approximately laminar.