3.1 Growth rates and eigenfunctions

To explore the linear instability mechanism dependent on the longitudinal and transverse wavenumbers, the maximal real value of the spectrum for a range of $(k_{x},k_{y})$ is calculated at various fixed values of $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6F1}_{s}$. Unless otherwise stated, we consider a constant $Pr$ of 0.71 corresponding to air at a temperature of around $2\,^{\circ }\text{C}$. For a slope angle of $\unicode[STIX]{x1D6FC}=2^{\circ }$ and $\unicode[STIX]{x1D6F1}_{s}=1.7$, the maximal possible growth rates for wavenumber vectors $(k_{x},k_{y})$ within the interval $[0,2.0]\times [0,3.2]$ are shown in figure 2(a). Only the positive growth rates, i.e. unstable modes, are highlighted. We can see that the growth rates tends to grow with decreasing $k_{x}$ component such that the maximal instabilities occur at $k_{x}=0$, i.e. the most unstable modes are purely along the direction transverse to the main velocity component. We will designate this instability as the ‘transverse mode’ due to the fact that its only non-zero wavenumber $k_{y}$ is transverse to the base-flow direction. It should be noted that similar modes such as observed in Rayleigh–Bénard convection or Görtler flows are also commonly referred to as longitudinal vortices or rolls in other literature. This mode is almost the same as the shallow-slope instability for katabatic flows described in Xiao & Senocak (Reference Xiao and Senocak2019), but with a smaller transverse wavenumber. Thus, it also deviates from Squire’s theorem, which asserts that, in the absence of buoyancy effects, the most unstable mode of parallel base flows (e.g. the Prandtl model flow profile) must be two-dimensional and vary along the streamwise direction (Schmid & Henningson Reference Schmid and Henningson2001). This apparent disagreement is due to the presence of the buoyancy and stable stratification, which act as additional body forces both orthogonal and parallel to the main slope flow, similar to the emergence of centrifugal instabilities in Görtler vortices or Taylor–Couette flows (Taylor Reference Taylor1923; Görtler Reference Görtler1959; Schmid & Henningson Reference Schmid and Henningson2001). The imaginary part of the eigenvalues of the pure transverse instability are all zero, which indicates that the vortex modes in this case are all spatially stationary.

At a moderately steep slope angle of $\unicode[STIX]{x1D6FC}=30^{\circ }$ and stability with $\unicode[STIX]{x1D6F1}_{s}=6.1$, the maximal possible growth rates for wavenumber vectors $(k_{x},k_{y})$ within the interval $[0,0.8]\times [0,0.25]$ are displayed in figure 2(d). In contrast to the situation at $\unicode[STIX]{x1D6FC}=2^{\circ }$, figure 2(d) shows that the growth rate tends to increase with decreasing $k_{y}$ component such that the maximal instabilities occur at $k_{y}=0$, thus the most unstable modes are purely along the slope direction, parallel to the main flow. Following the earlier convention, we will designate this instability as the ‘longitudinal mode’ due to the fact that its only non-zero wavenumber $k_{x}$ is along the base-flow direction. This behaviour is in line with Squire’s theorem and agrees with the well-known behaviour in other types of parallel flow such as plane mixing layers with KH instability. Another difference from the gentle slope case is that, for $\unicode[STIX]{x1D6FC}=30^{\circ },\unicode[STIX]{x1D6F1}_{s}=6.1$, the imaginary part of the most unstable eigenvalue $\unicode[STIX]{x1D714}$ (i.e. the angular frequency of the most unstable mode) is non-zero with value approximately $\text{Im}(\unicode[STIX]{x1D714})\approx 0.27$. Thus, in contrast to the transverse mode consisting of stationary vortices, the longitudinal mode is a transverse wave propagating along the slope. This shows that, at the given slope angle $\unicode[STIX]{x1D6FC}=30^{\circ }$, the misalignment of shear and stratification can lead to spontaneous emission of streamwise-travelling waves from the equilibrium state given by Prandtl’s profile. It is also important to note that the most unstable eigenvalue is single, i.e. not as part of a conjugate pair, which means that the instability has a preferred direction of propagation due to the existence of the slope. The switch from the stationary transverse mode to longitudinal wave modes in slope flows with growing slope angle could be attributed to the stronger along-slope gravity component, which dominates the stable stratification orthogonal to the surface at steep inclination angles.

In order to determine which type of stratified flow instability the longitudinal wave mode could be identified with, the base-flow profiles displayed in figure 1(b) and given by equations (2.4)–(2.5) can be used to ascertain which kind of waves are supported by the Prandtl base flow. Firstly, from the given buoyancy profile under Prandtl’s slope-flow model, it can be verified that the buoyancy gradient attains its maximum at the surface and decays exponentially with growing height, hence Prandtl’s base flow lacks two distinct buoyancy interfaces to create resonating internal gravity waves which would be required for Taylor–Caulfield type of instability, as described in Eaves & Balmforth (Reference Eaves and Balmforth2019). For the identification of a Holmboe-type instability, which arises from resonance between a vorticity wave and an internal gravity wave as described in Carpenter, Balmforth & Lawrence (Reference Carpenter, Balmforth and Lawrence2010) and Eaves & Balmforth (Reference Eaves and Balmforth2019), we can apply the characterisation given in Salehipour, Caulfield & Peltier (Reference Salehipour, Caulfield and Peltier2016), which asserts that the base-flow profile supportive of Holmboe modes should consist of a sharp density or buoyancy interface immersed within a broader shear region. This characterisation is also backed up by plots of typical base profiles prone to the Holmboe instability shown in Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011) and Eaves & Balmforth (Reference Eaves and Balmforth2019). From figure 1(b) and equations (2.4)–(2.5), however, it can be deduced for the Prandtl base-flow profiles that both the velocity shear and buoyancy gradient attain their maximum at the surface $z_{n}=0$ and exhibit the same exponential decay. Thus, the buoyancy interface is not embedded within the shear layer of the Prandtl base profile, and accordingly this means that the Prandtl base profile is not conducive to Holmboe-type instability. Hence, the longitudinal travelling mode observed in this work could be more reasonably attributed to a shear instability.

The growth-rate contours at two intermediate angles between $2^{\circ }$ and $30^{\circ }$ are displayed in figure 2(b,c), and a gradual transition from the transverse mode to the longitudinal mode can be clearly observed. For all the slope-flow configurations with angles shown above, however, the most unstable instability mode propagates either parallel ($k_{y}=0$) or orthogonal ($k_{x}=0$) to the along-slope direction, but never in an oblique direction, which is a different behaviour from the oblique instabilities observed in the spanwise-stratified Couette flow studied by Facchini et al. (Reference Facchini, Favier, Le Gal, Wang and Le Bars2018).

Figure 3. Eigenfunctions of the neutral state at $\unicode[STIX]{x1D6F1}_{s}=2.1$ for the stationary transverse mode at $\unicode[STIX]{x1D6FC}=5^{\circ }$. Dash-dotted lines represent the real part, asterisks represent the imaginary part, and the dashed line is the magnitude. The disturbance magnitudes have been normalised with the maximal occurring buoyancy disturbance magnitude in each case. The solid line represents the normalised base Prandtl velocity profile $u_{n}$.

The computed eigenfunctions for the neutral state of both instability types are displayed in figures 3 and 4. From the shape of the eigenmodes in figure 3, it is clear that the transverse mode at $5^{\circ }$ has strong disturbances in all three velocity components, which decay rapidly with growing height. Another key observation from figure 3 is that the node of the velocity disturbance, i.e. centre of the instability rolls, is located approximately near the maximum of the base velocity profile. On the other hand, the longitudinal mode at slope angle $30^{\circ }$ as shown in figure 4 is two-dimensional with zero cross-flow velocity component; its main flow and slope-normal flow disturbances are a lot weaker compared to its buoyancy disturbance, which, however, persist at further distances away from the surface.

In contrast to the stationary transverse mode, the longitudinal mode is propagating with wave speed $c_{x}=\text{Im}(\unicode[STIX]{x1D714})/k_{x}$, which is approximately $0.4$ at angle $\unicode[STIX]{x1D6FC}=30^{\circ }$, as can be seen from figure 6(b). The first critical layer, i.e. the position $z_{c}$ within the shear layer where $c_{x}$ equals the Prandtl base-flow velocity, can be found just below the location of the near-surface jet with maximal velocity. Since Prandtl’s profile is not monotonic, there exists also one other critical point slightly above the near-surface jet where $c_{x}$ is attained. Upon closer inspection of figure 4, it can also be seen that the maximum of the longitudinal velocity disturbance $\hat{u}$ is located closely above the surface where the shear of the base profile is maximal. On the other hand, the slope-normal velocity disturbance ${\hat{w}}$ attains its maximum near the first inflection point of the base velocity profile (located at $z_{n}\approx 2.22$), which is also a local maximum of the base profile’s shear. These observations suggest that the longitudinal mode is principally a shear instability driven by resonating waves formed at vorticity interfaces.

These differences between the transverse and longitudinal instabilities strongly imply that the two modes are driven by different flow-physical mechanisms. Following our previous arguments, the transverse mode exhibits hallmarks of a centrifugal instability such as Görtler or Taylor–Couette vortices (e.g. Taylor Reference Taylor1923; Görtler Reference Görtler1959), whereas the oscillatory behaviour of the longitudinal mode could be attributed to an inviscid shear instability modulated by buoyancy, which can be regarded as a generalised KH instability (Drazin & Reid Reference Drazin and Reid2004).

Figure 4. Eigenfunctions of the neutral state at $\unicode[STIX]{x1D6F1}_{s}=5.9$ for the propagating longitudinal mode at $\unicode[STIX]{x1D6FC}=30^{\circ }$. Dash-dotted lines represent the real part, asterisks represent the imaginary part, and the dashed line is the magnitude. The disturbance magnitudes have been normalised with the maximal occurring buoyancy disturbance magnitude in each case. The solid line represents the normalised base Prandtl velocity profile $u_{n}$.

Figure 5. Growth-rate contours at $\unicode[STIX]{x1D6FC}=5^{\circ }$ for (a) transverse and (b) longitudinal modes and at $\unicode[STIX]{x1D6FC}=30^{\circ }$ for (c) transverse and (d) longitudinal modes. The neutral curves are given by $\text{Re}(\unicode[STIX]{x1D714})=0$.

3.2 Neutral curves and critical stability

The results from the previous subsection clearly show that the most unstable modes at each slope angle $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6F1}_{s}$ are propagating along either the main flow direction or the transverse direction, i.e. only one component of the wavenumber vector $(k_{x},k_{y})$ is non-zero in order to attain maximal growth rates. This observation is consistent with our previous discovery as outlined in Xiao & Senocak (Reference Xiao and Senocak2019). The critical $\unicode[STIX]{x1D6F1}_{s}$ for the onset of instability at a specific slope angle $\unicode[STIX]{x1D6FC}$ and $Pr$ number can be found by plotting the growth-rate contours over a range of $\unicode[STIX]{x1D6F1}_{s}$ separately for $k_{x}$ and $k_{y}$, assuming that the other wavenumber is zero. For $\unicode[STIX]{x1D6FC}=5^{\circ }$, the results are shown in figure 5(a,b). We can see that the minimal $\unicode[STIX]{x1D6F1}_{s}$ for the transverse mode is approximately 2.05, whereas the longitudinal mode requires a minimal $\unicode[STIX]{x1D6F1}_{s}=2.25$ to become dynamically unstable. Thus, the most dominant instability in this case is the transverse mode, in agreement with the growth-rate contour shown in the previous subsection. In figure 5(c,d), the results for a steep angle of $\unicode[STIX]{x1D6FC}=30^{\circ }$ are displayed. This time, we can see that the critical $\unicode[STIX]{x1D6F1}_{s}$ for the longitudinal mode is approximately $5.9$, whereas the transverse mode requires a minimal $\unicode[STIX]{x1D6F1}_{s}\approx 13.5$ to become unstable. This means that, in contrast to the case of shallower slope angle, the instability to be triggered first at $\unicode[STIX]{x1D6FC}=30^{\circ }$ is the longitudinal mode in the along-slope direction, which is also supported by the growth-rate contour shown the previous subsection.

As these results have shown, the separation between the minimal $\unicode[STIX]{x1D6F1}_{s}$ values required for the onset of each instability mode is a lot larger at the steeper slope angle $\unicode[STIX]{x1D6FC}=30^{\circ }$ than at $\unicode[STIX]{x1D6FC}=5^{\circ }$. The implication of this is that, while the longitudinal mode is clearly the preferred instability to be triggered first at large slope angles, for shallow slopes, the transverse mode is initiated only slightly before the onset of the longitudinal mode. As a result, we should expect that, under most configurations in which the anabatic slope flow is unstable at low angles, both types of instability modes can coexist simultaneously. Moreover, since it is known from Xiao & Senocak (Reference Xiao and Senocak2019) that the gradient Richardson number of the base profile is related to its $\unicode[STIX]{x1D6F1}_{s}$ value via the relation $Ri=Pr/\unicode[STIX]{x1D6F1}_{s}^{2}$, an interesting observation pertaining to the stability behaviour at shallow slopes less than $3^{\circ }$ is that, at those angles, the implied gradient Richardson number corresponding to the stability threshold value of $\unicode[STIX]{x1D6F1}_{s}$ is larger than the critical value of $Ri_{c}=0.25$. Thus, the anabatic base-flow profiles at low slope angles under Prandtl’s model can serve as another set of counter-examples to the celebrated stability theorem presented in Miles (Reference Miles1961).

For both the shallow-slope and steep-slope cases, we can determine from figure 5 that the unstable growth rates increase with growing $\unicode[STIX]{x1D6F1}_{s}$ within the range of shown $\unicode[STIX]{x1D6F1}_{s}$ values. A more detailed analysis of the $\unicode[STIX]{x1D6F1}_{s}$ dependence of the growth rate will be given in § 3.3.

The aforementioned switch from a stationary transverse mode at slope angle of $5^{\circ }$ to a travelling longitudinal mode at around $30^{\circ }$ can be attributed to the increasing longitudinal gravity component driving the flow. To further explore the influence of slope angle $\unicode[STIX]{x1D6FC}$ on the instability, the stability thresholds for both longitudinal and transverse instabilities (i.e. the critical $\unicode[STIX]{x1D6F1}_{s}$ value at which the growth rate is zero) as functions of slope angle over the range of $[2^{\circ },60^{\circ }]$ are determined and displayed in figure 6(a,c). We observe that the critical $\unicode[STIX]{x1D6F1}_{s}$ threshold for the transverse mode of instability is higher than its longitudinal counterpart at shallow angles less than $9^{\circ }$. However, it increases for growing angles, whereas the critical value for the longitudinal mode stays almost constant over the same range of slope angles. This implies that the transverse mode becomes increasingly stable with growing $\unicode[STIX]{x1D6FC}$, and the maximal angle at which its stability threshold is less than or equal to its longitudinal counterpart is approximately $9^{\circ }$, as shown by figure 6(a). The slope-angle dependence of the wavenumber of each instability mode at its critical stability threshold $\unicode[STIX]{x1D6F1}_{s}$ is displayed in figure 6(b), which shows that the critical transverse mode characterised by $k_{y}$ has a shorter wavelength than its longitudinal counterpart $k_{x}$ for all angles within $[2^{\circ },60^{\circ }]$. The frequency of the oscillatory critical longitudinal mode grows monotonically for all slope angles larger than $9^{\circ }$, i.e. when the longitudinal mode is the dominant instability with lower stability threshold.

Figure 6. Critical curves of both instability types dependent on slope angle at$Pr=0.71$ for anabatic slope flow: (a) critical perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$ for slopes$\unicode[STIX]{x1D6FC}\leqslant 30^{\circ }$; (b) critical wavenumber $k_{x},k_{y}$ and frequency $\text{Im}(\unicode[STIX]{x1D714})$ for slopes $\unicode[STIX]{x1D6FC}\leqslant 30^{\circ }$; (c) critical perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$ for slopes $30\leqslant \unicode[STIX]{x1D6FC}\leqslant 60^{\circ }$; and (d) critical wavenumbers $k_{x}$ and $k_{y}$ and frequency $\text{Im}(\unicode[STIX]{x1D714})$ for slopes $30\leqslant \unicode[STIX]{x1D6FC}\leqslant 60^{\circ }$. Since the unstable transverse mode is always stationary, its critical frequency (i.e. $\text{Im}(\unicode[STIX]{x1D714})$) is zero for all angles and not shown here.

3.3 Instability growth at different stratification perturbation regimes

We have learned from our preceding analysis that, for each slope angle and at sufficiently high $\unicode[STIX]{x1D6F1}_{s}$, two instability modes emerge. Thus, using the value of $\unicode[STIX]{x1D6F1}_{s}$ for a given slope angle, it is possible to delineate a stratification regime with small $\unicode[STIX]{x1D6F1}_{s}$ where all linear instabilities are suppressed from more dynamically unstable regimes at higher $\unicode[STIX]{x1D6F1}_{s}$.

Mahrt (Reference Mahrt1998) is credited for proposing the weakly stable and very stable regimes as a qualitative classification scheme to identify distinct flow conditions in stratified atmospheric boundary layer flows. A weakly stable regime refers to conditions where fully turbulent flows prevail, and a very stable regime is associated with flow intermittency and patchy turbulence. Based on our analysis, we deem it necessary to further refine Mahrt’s classification and introduce a super-stable regime to distinguish flow conditions where there are no flow instabilities and the flow maintains a steady state. The super-stable regime is equivalent to a linearly stable regime. The reason we call it ‘super’-stable is because we prefer to complement Mahrt’s existing terminology and emphasise a natural progression from weakly stable towards more stable regimes.

As we have shown in our previous work (Xiao & Senocak Reference Xiao and Senocak2019), the extent of the linearly stable flow conditions is defined based on the slope angle and the stratification perturbation number for a given fluid. Therefore, a super-stable regime can be proposed as a quantitative lower bound for the very stable regime. Next, we investigate how the properties of these two instability modes evolve over a large range of $\unicode[STIX]{x1D6F1}_{s}$ covering the super-stable, very stable and weakly stable regimes. We note that these regimes are strictly for the Prandtl model for slope flows, where there is no externally imposed flow shear.

3.3.2 Weakly stable regime: $\unicode[STIX]{x1D6F1}_{s}\rightarrow \infty$

Under weak background stratification relative to a very strong thermal forcing at the surface, the perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$ approaches infinity. In this limit, all terms involving $\unicode[STIX]{x1D6F1}_{s}$ in (3.1)–(3.5) disappear and the disturbance equations simplify to Rayleigh’s equation (see e.g. Drazin & Reid (Reference Drazin and Reid2004) and Schmid & Henningson (Reference Schmid and Henningson2001)). The buoyancy disturbance equation reduces to the following direct relation between buoyancy and the slope-normal velocity:

(3.8)$$\begin{eqnarray}\hat{b}=\frac{\text{i}b_{n}^{\prime }}{u_{n}k_{x}}{\hat{w}}.\end{eqnarray}$$

This means that neither the buoyancy $\hat{b}$ nor the slope angle $\unicode[STIX]{x1D6FC}$ appears in the momentum equations, hence they have no influence on other flow variables. This shows that, independent of the slope angle, the unstable modes at the weakly stable limit are KH modes driven purely by shear in the base-flow profile, and not by density or buoyancy gradient. The decoupling of buoyancy from the momentum equations in the weakly stable limit further supports the assertion that, in this regime, the two-dimensional longitudinal instabilities are not to be identified with Holmboe modes.

Since the Prandtl model for slope flow has infinitely many inflection points, it satisfies Rayleigh’s necessary (but not sufficient) instability condition (Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2016, p. 573). This analysis coincidentally shows that, in contrast to the configuration described by Candelier et al. (Reference Candelier, Le Dizès and Millet2011), the longitudinal mode of instability at the weakly stable regime is dependent on the shape of the base-flow profile; and that, for a base flow other than given by the Prandtl model, which does not satisfy the well-known necessary instability conditions by Rayleigh and Fjørtoff, there can be no unstable modes at the limit $\unicode[STIX]{x1D6F1}_{s}\gg 1$. Because this is a pure shear instability, we may conclude via Squire’s theorem that the most unstable mode must be two-dimensional, which we have already showed in figure 4, and, therefore, longitudinal (i.e. the transverse wavevector and velocity component are zero in that case). The nature of the longitudinal mode arising here is also different from the travelling waves discovered over heated inclined surfaces by Sparrow & Husar (Reference Sparrow and Husar1969) and on a slope-normal wall within a stable background stratification studied by Chen et al. (Reference Chen, Bai and Le Dizès2016), in which the instabilities were identified as Tollmien–Schlichting waves on a boundary layer, decaying at the inviscid limit.

In figure 7(a), the normalised growth rate of the most unstable longitudinal mode with $k_{y}=0$ is plotted for three different angles from $10^{\circ }$ to $30^{\circ }$. It is clear that, in all cases, the growth rates gradually decrease with growing $\unicode[STIX]{x1D6F1}_{s}$ to the same limit value of approximately 0.013 for $\unicode[STIX]{x1D6F1}_{s}>300$. A similar behaviour can be seen for the normalised wavevector (figure 7b) and angular frequency (figure 7c), which converge to an angle-independent value at high $\unicode[STIX]{x1D6F1}_{s}$. Both the wavenumber and angular frequency of the most unstable longitudinal mode decrease with growing $\unicode[STIX]{x1D6F1}_{s}$, i.e. the instability tends to have a larger wavelength and lower frequency.

The behaviour of the most unstable transverse mode with $k_{x}=0$ is displayed in figure 8. All such modes are stationary (i.e. $\text{Im}(\unicode[STIX]{x1D714})=0$). It shows that, for $\unicode[STIX]{x1D6F1}_{s}<1000$, the growth rate of the transverse mode is at least of the order of 0.01, and hence comparable to that of the longitudinal mode, which is $0.013$. At larger values of $\unicode[STIX]{x1D6F1}_{s}$, the transverse growth rate slowly decreases to zero, in agreement with Squire’s theorem and the previous analysis for $\unicode[STIX]{x1D6F1}_{s}\rightarrow \infty$. In contrast to the behaviour of the longitudinal modes, with growing perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$, figure 8(b) indicates that the wavenumber of the most unstable transverse mode increases, i.e. the wavelength of the instability decreases, with the result that the longitudinal vortices get closer to each other.

Figure 7. Most unstable longitudinal mode for three different slope angles: (a) normalised growth rate; (b) normalised longitudinal wavenumber; and (c) normalised angular frequency depending on perturbation $\unicode[STIX]{x1D6F1}_{s}$. Insets show the behaviour near the critical $\unicode[STIX]{x1D6F1}_{s}$ of the neutral state.

3.3.3 Very stable regime

In § 3.3.1 we have argued that the most unstable growth rate for $\unicode[STIX]{x1D6F1}_{s}\ll 1$ cannot be positive. We then infer that, at each slope angle, there has to be a specific finite value of $\unicode[STIX]{x1D6F1}_{s}$ where the growth rate of the unstable mode attains an extreme value. As we observe from figures 7(a) and 8(a), this turns out to be the maximal growth rate over all stratification perturbation numbers $\unicode[STIX]{x1D6F1}_{s}$ and holds true for both the longitudinal as well as the transverse mode.

Figure 8. Most unstable transverse mode for three different slope angles: (a) normalised growth rate; and (b) normalised transverse wavenumber depending on $\unicode[STIX]{x1D6F1}_{s}$. Insets show the behaviour near the critical $\unicode[STIX]{x1D6F1}_{s}$ value with zero growth rate.

Figure 9. Comparison of most unstable longitudinal and transverse modes depending on slope angle. (a) Stratification perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$ at which maximal growth rate is attained. The inset shows the ratio between $\unicode[STIX]{x1D6F1}_{max}$ of maximal growth rate and neutral stability threshold $\unicode[STIX]{x1D6F1}_{0}$. (b) Maximal growth rate and corresponding wavevector.

The properties of the most unstable longitudinal and transverse mode are plotted over the slope angle in figure 9(a) ($\unicode[STIX]{x1D6F1}_{s}$ with maximal growth) and in figure 9(b) (wavenumber). The inset of figure 9(a) shows that the maximal growth rates for both instability modes are attained at $\unicode[STIX]{x1D6F1}_{max}$ that are less than four times the critical stability threshold value $\unicode[STIX]{x1D6F1}_{0}$; thus the linear growth rates begin to decrease with growing $\unicode[STIX]{x1D6F1}_{s}$ from $\unicode[STIX]{x1D6F1}_{s}=\unicode[STIX]{x1D6F1}_{0}$ onwards. We see that, for angle values ranging from gentle slopes at $2^{\circ }$ to moderately steep slopes at $30^{\circ }$, the most unstable transverse mode has a larger growth rate than its longitudinal counterpart. The wavenumber plot (i.e. figure 9b) indicates that the most unstable longitudinal mode has a longer wavelength than its transverse counterpart for the range of slope angles considered here, and this divergence grows with increasing angle such that, at $\unicode[STIX]{x1D6FC}=30^{\circ }$, the wavelength of the most unstable longitudinal mode is approximately three times as large as that of its transverse counterpart. This finding is in agreement with the behaviour of the wavelength for the neutral states as shown in figure 6(a).

3.4 Prandtl-number dependence

Since thermal and momentum diffusivities are related via $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D708}/Pr$, it is obvious that, when fluid viscosity is kept constant, an increase in Prandtl number implies a diminishing effect of thermal diffusion relative to momentum diffusion, and vice versa. From the inspection of the disturbance equations (3.1)–(3.5), we observe that, for anabatic flows with positive surface buoyancy flux, the buoyancy contribution to the momentum equations is directly proportional to $Pr$ when all other parameters are held constant. We can then deduce that thermal diffusion tends to smear out or dampen the buoyancy force term. This can be investigated by varying $Pr$. The effect of the Prandtl number thus also directly measures the influence of buoyancy on the stability of the anabatic slope flow. Therefore, we investigate whether an increased buoyancy force contribution to the momentum equations tends to enhance or dampen the growth rate of different instability types.

3.4.1 Asymptotic analysis for large Prandtl number

With increasing Prandtl number, the importance of heat conduction at the surface diminishes relative to viscous diffusion, and it completely vanishes when $Pr$ reaches infinity. For large Prandtl-number values, pressure $p$ scales approximately linearly with $Pr$. Thus if we define $\bar{p}=\hat{p}/Pr$ to be the pressure disturbance normalised by the Prandtl number, then the eigenvalue problem for the flow disturbances (3.1)–(3.5) at very large $Pr$ for wavenumbers $k_{x},k_{y}\ll Pr$ and perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$, which is of a similar order of magnitude as $\sin \unicode[STIX]{x1D6FC}$, can be approximated by the following simpler set of equations:

(3.9)$$\begin{eqnarray}\displaystyle & \displaystyle 0=\text{i}k_{x}\hat{u} +\text{i}k_{y}\hat{v}+\frac{\unicode[STIX]{x2202}{\hat{w}}}{\unicode[STIX]{x2202}z}, & \displaystyle\end{eqnarray}$$

(3.10)$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\text{i}k_{x}\bar{p}+\frac{\sin \unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6F1}_{s}}\left(-(k_{x}^{2}+k_{y}^{2})\hat{u} +\frac{\unicode[STIX]{x2202}^{2}\hat{u} }{\unicode[STIX]{x2202}z^{2}}+\hat{b}\right), & \displaystyle\end{eqnarray}$$

(3.11)$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\text{i}k_{y}\bar{p}+\frac{\sin \unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6F1}_{s}}\left(-(k_{x}^{2}+k_{y}^{2})\hat{v}+\frac{\unicode[STIX]{x2202}^{2}\hat{v}}{\unicode[STIX]{x2202}z^{2}}\right), & \displaystyle\end{eqnarray}$$

(3.12)$$\begin{eqnarray}\displaystyle & \displaystyle 0=-\frac{\unicode[STIX]{x2202}\bar{p}}{\unicode[STIX]{x2202}z}+\frac{\sin \unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6F1}_{s}}\left(-(k_{x}^{2}+k_{y}^{2}){\hat{w}}+\frac{\unicode[STIX]{x2202}^{2}{\hat{w}}}{\unicode[STIX]{x2202}z^{2}}+\hat{b}\cot \unicode[STIX]{x1D6FC}\right), & \displaystyle\end{eqnarray}$$

(3.13)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}\hat{b}+\text{i}uk_{x}\hat{b}+b_{n}^{\prime }{\hat{w}}=\frac{\sin \unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6F1}_{s}}\left(-(k_{x}^{2}+k_{y}^{2})\hat{b}+\frac{\unicode[STIX]{x2202}^{2}\hat{b}}{\unicode[STIX]{x2202}z^{2}}-(\hat{u} +{\hat{w}}\cot \unicode[STIX]{x1D6FC})\right). & \displaystyle\end{eqnarray}$$

The above system of simplified equations can be solved in the same way as the complete set of disturbance equations (3.1)–(3.5) to obtain the limit growth curves for $Pr$ approaching infinity.

Figure 10. Growth rate of the most unstable transverse mode for different Prandtl numbers at (a) $30^{\circ }$ and (b) $5^{\circ }$. The insets show the behaviour near the neutral state.

Figure 11. Growth rate of the most unstable longitudinal mode for different Prandtl numbers at (a) $30^{\circ }$ and (b) $5^{\circ }$. The insets show the behaviour near the neutral state.

The influence of the Prandtl number on the most unstable mode for the transverse and longitudinal modes of instabilities at two different angles $\unicode[STIX]{x1D6FC}=5^{\circ }$ and $30^{\circ }$ are shown in figures 10 and 11, respectively. The figures display the maximal growth rate as a function of the perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$. We observe from these figures that, for both the transverse and the longitudinal modes, an increase in $Pr$ leads to a larger growth rate. In other words, despite a growing viscous contribution due to an increase in $Pr$, the buoyancy force in the momentum equations further strengthens both types of instabilities.

We also observe from figures 10 (transverse mode) and 11 (longitudinal mode) that, close to the neutral state with zero growth rate and low $\unicode[STIX]{x1D6F1}_{s}$ values, the growth curves for different Prandtl numbers approach the limit curve with $Pr=\infty$ for both types of instabilities. Therefore, at sufficiently high $Pr$ numbers, the neutral stability threshold value for $\unicode[STIX]{x1D6F1}_{s}$ become independent of $Pr$. However, away from the neutral state, Prandtl-number curves begin to diverge from each other with increasing $\unicode[STIX]{x1D6F1}_{s}$. Since this analysis does not explicitly depend on the sign of the surface buoyancy flux, a similar conclusion can be reached for the katabatic case, albeit with different values for the critical stability threshold.

For the transverse mode, it can be seen in figure 10(b) that, at the lower angle $\unicode[STIX]{x1D6FC}=5^{\circ }$, the Prandtl number has very little effect on the stability threshold value of $\unicode[STIX]{x1D6F1}_{s}$, which is the intercept of the growth-rate curve with the abscissa. However, figure 10(a) shows that, at the steeper angle $\unicode[STIX]{x1D6FC}=30^{\circ }$, a larger $Pr$ corresponds to a lower stability threshold. The situation is similar for the longitudinal mode where the value of $Pr$ at $\unicode[STIX]{x1D6FC}=5^{\circ }$ has little effect on the stability threshold, as shown by figure 11(b). However, at the larger angle $\unicode[STIX]{x1D6FC}=30^{\circ }$, according to figure 11(a), an increase of $Pr$ raises the stability threshold.

3.4.2 Effect on neutral stability

The neutral stability curves for different values of Prandtl number ranging from $Pr=0.71$ for air at night-time temperature of $2^{\circ }$ to $Pr=6.1$ for water at room temperature are displayed in figure 12, including the limit case of infinite Prandtl number. The curves shown in figure 12(a) indicate that, for low angles $\unicode[STIX]{x1D6FC}<9^{\circ }$, the Prandtl-number effect on the critical curves is negligibly small. For larger angles, however, both modes of instabilities are starting to be visibly affected by a variation of the Prandtl number ranging from 0.71 to 6.1. At $\unicode[STIX]{x1D6FC}=30^{\circ }$, the critical $\unicode[STIX]{x1D6F1}_{s}$ of the transverse mode for $Pr=6.1$ is approximately 75 % of the corresponding value at $Pr=0.71$, whereas for the longitudinal mode the critical $\unicode[STIX]{x1D6F1}_{s}$ at $Pr=6.1$ is around 20 % higher than at $Pr=0.71$. For angles $\unicode[STIX]{x1D6FC}>10^{\circ }$, we can see that the critical stability threshold as well as the critical wavenumber at the neutral state of the transverse mode decreases with growing $Pr$, whereas the opposite behaviour is observed for the same quantities of the longitudinal mode.

Figure 12. Critical curves of both instability types dependent on slope angle for anabatic slope flow at different Prandtl numbers. (a) Critical perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$. The inset shows a detailed view at low slope angles from $2^{\circ }$ to $5^{\circ }$. (b) Critical wavenumber. (c) Critical angular frequency for the longitudinal mode. Since the unstable transverse mode is always stationary, its critical frequency is zero for all angles and is hence not displayed in panel (c).

3.4.3 Effect on maximal unstable mode

The behaviours of the most unstable longitudinal and transverse modes for different values of Prandtl number ranging from $Pr=0.71$ to $Pr=6.1$ as a function of slope angle $\unicode[STIX]{x1D6FC}$ are displayed in figure 13. Figure 13(a) clearly indicates that, with growing Prandtl number, the maximal attainable growth rate at a given angle $\unicode[STIX]{x1D6FC}$ also increases. It is obvious from figure 13(c) that the wavenumber of the most unstable transverse mode is almost independent of $Pr$. On the other hand, the wavenumber as well as oscillation frequency of the most unstable longitudinal mode decrease with growing Prandtl number, as displayed in figure 13(c,d).

Figure 13. Comparison of most unstable longitudinal and transverse modes depending on slope angle for different Prandtl numbers: (a) $\unicode[STIX]{x1D6F1}_{s}$ at which maximal growth rate is attained; (b) maximal growth rate; (c) wavenumber at most unstable state; and (d) oscillation frequency of the longitudinal mode at most unstable state.

4.1 Linear stability regions

The anabatic version of the Prandtl model embodies the same distinct instability modes as observed computationally in the katabatic version, namely stationary longitudinal vortices (transverse mode) and two-dimensional transverse waves (longitudinal mode) as identified in Xiao & Senocak (Reference Xiao and Senocak2019). At any given slope angle, the mode with the smaller stability threshold, i.e. the critical $\unicode[STIX]{x1D6F1}_{s}$ value at which zero growth rate is reached, is also the more dominant mode that gets triggered. A comparison between the slope-angle-dependent stability regions as determined in the above way for both the anabatic and katabatic cases is displayed in figure 14. We can observe that, while for both cases, an increase of slope angle gradually shifts the dominant instability from the transverse to the longitudinal mode, the transition angle at which the stability thresholds of both modes are equal is much lower in the anabatic case, i.e. around $9^{\circ }$, compared to approximately $62^{\circ }$ for the katabatic case. Thus, while the longitudinal instability mode propagating along the slope could be observed already at moderate angles in anabatic flows, a very steep slope is required for that to happen in katabatic flows.

The critical wavenumbers for both modes and critical frequencies for the wave mode obtained at the stability threshold are displayed in figure 15. The comparison between anabatic and katabatic cases within the angle range of $[2^{\circ },30^{\circ }]$ clearly shows that, for both types of instabilities, the vortices and waves produced in katabatic flows tend to have wavelengths that are three times as large as in anabatic flows, whereas the angular frequency of the wave mode is multiple times higher in the anabatic case. Despite the opposite base-flow directions for katabatic and anabatic flows, the longitudinal modes for the two flow types both propagate downslope, as shown by the opposite signs of their oscillation frequencies and their base-flow velocities. This divergence between the critical wavelengths and frequencies grows with decreasing slope angle and becomes most pronounced at the shallowest slopes. From a flow-physical point of view, this phenomenon can be readily explained by the fact that the anabatic flow configuration is more unstable due to the stronger destabilising influence of the surface heating acting against the stable background stratification compared to the opposite effect of surface cooling in the katabatic case. Unsurprisingly, the discrepancy between the stability effects of equal-magnitude, opposite-sign buoyancy fluxes in the Prandtl model increases for smaller slope inclinations due to the larger buoyancy force component normal to the surface, as shown in figure 15.

Figure 14. Stability regions for slope flows under Prandtl’s model: (a) anabatic case and (b) katabatic case. The four delineated stability regions are the stable regime (S), longitudinal wave mode (L), transverse stationary vortex mode (T), and higher-order instabilities (H).

Figure 15. Comparison between anabatic and katabatic instabilities. (a) Normalised wavenumbers for neutral anabatic and katabatic modes (transverse and longitudinal) depending on slope angle $\unicode[STIX]{x1D6FC}$. (b) Normalised angular frequency for the neutral longitudinal modes of both katabatic and anabatic flows as a function of slope angle. The transverse modes are stationary for both flow types.

4.2 Direct numerical simulations

We carry out DNS to independently validate the findings of LSA and to visualise the flow structures induced by the flow instabilities. We solve the buoyancy-driven incompressible flow equations (2.1)–(2.3) using a Cartesian mesh three-dimensional Navier–Stokes solver (Jacobsen & Senocak Reference Jacobsen and Senocak2013). The code adopts a second-order-accurate Adams–Bashforth scheme for time advancement and a second-order central difference scheme for spatial derivatives. The pressure Poisson equation is solved with a geometric multigrid technique. Umphrey et al. (Reference Umphrey, DeLeon and Senocak2017) validated the current code using the Prandtl model and demonstrated globally second-order-accurate solutions. This same code was also used to study the katabatic flow instabilities in Xiao & Senocak (Reference Xiao and Senocak2019). The simulation domains are rectangular boxes of dimensional size $L_{x}\times L_{y}\times L_{z}$. Periodic boundary conditions are imposed in both the longitudinal and transverse directions, whereas no-slip conditions with a constant buoyancy flux are applied on the lower surface at $z=0$, and the top surface is subject to adiabatic free-slip conditions. The top boundary was placed no less than 50 times the characteristic slope length scale $l_{0}$ given by (2.6) to capture quiescent conditions aloft. The longitudinal and transverse size of the domain are chosen to be an integer multiple of the targeted wavelength of an instability in the particular direction. The initial conditions are Prandtl’s laminar solution (2.4) and (2.5) without any disturbances.

Figure 16. Transverse instability mode for (a) anabatic and (b) katabatic flows at$\unicode[STIX]{x1D6FC}=5^{\circ }$ and $\unicode[STIX]{x1D6F1}_{s}=2.15$ after nonlinear saturation. Colour shows the normalised slope-normal velocity component. The contours are drawn for half of the respective maximum and minimum slope-normal velocity values.

4.2.1 Transverse mode of instability

At a slope angle $\unicode[STIX]{x1D6FC}=5^{\circ }$, $\unicode[STIX]{x1D6F1}_{s}=2.15$ and $Pr=0.71$, DNS is carried out over a domain with a dimensionless size of $36\times 13.5\times 100$. The results display a purely transverse instability, in agreement with the prediction of the LSA carried out in the previous subsection. In figure 16(a), the steady along-slope velocity disturbance after initial growth and nonlinear saturation along the cross-slope direction is shown. We can see that three full periods of the instability mode have been established in the simulation, which for a total length of 13.5 in the cross-slope dimension of the simulation domain leads to a transverse wavelength of $\unicode[STIX]{x1D706}=4.5$ for the instability, corresponding to $k_{y}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}=1.4$. This wavenumber value falls within the range predicted by LSA for the transverse mode at $\unicode[STIX]{x1D6FC}=5^{\circ }$ at which the growth rate is positive, as shown in figure 5. The DNS also shows that the instability does not oscillate along the transverse direction, in agreement with LSA results predicting a stationary mode. For comparison, DNS of a saturated transverse mode for katabatic conditions at the same slope angle and $\unicode[STIX]{x1D6F1}_{s}$ value is shown in figure 16(b).

To effect a more quantitative analysis of the simulation results, we have conducted several DNS of the transverse mode at various slope angles $\unicode[STIX]{x1D6FC}$ where this mode is expected to be the dominant one, and imposed different $\unicode[STIX]{x1D6F1}_{s}$ values via surface buoyancy heat flux. The non-dimensional simulation domain size used is $3\times 43\times 70$, where the along-slope dimension is restricted to a small size in order to focus on the cross-slope extent along which the transverse mode is established. The cross-slope domain extent translates to a resolution of approximately $\unicode[STIX]{x0394}k_{y}\approx 0.15$ in the transverse wavenumber space. Thus, the total number of spatial transverse mode periods $n_{y}$ within the simulation domain has a wavenumber $k_{y}=\unicode[STIX]{x0394}k_{y}\,n_{y}$, and $n_{y}$ is such that the transverse mode with wavenumber $\unicode[STIX]{x0394}k_{y}\,n_{y}$ obtains maximal growth rate compared to all other integers. The Cartesian grid is chosen to ensure that at least two points resolve one length scale $l_{0}$ in each direction. The simulation results for the wavenumber selection, normalised growth rate and angular frequency of the transverse mode are displayed in table 1. It can be seen that the wavenumbers of the simulation agree with the values predicted by LSA. The simulated transverse modes are all stationary during the linear growth phase, agreeing with the zero imaginary part of the eigenvalue predicted by linear stability theory.

Figure 17. Longitudinal instability mode after nonlinear saturation. Normalised slope-normal velocity component is plotted. (a) Anabatic case with $\unicode[STIX]{x1D6FC}=30^{\circ }$ and $\unicode[STIX]{x1D6F1}_{s}=6.1$. (b) Katabatic case with $\unicode[STIX]{x1D6FC}=66^{\circ }$ and $\unicode[STIX]{x1D6F1}_{s}=18$. The contours are plotted at half the respective maximum and minimum slope-normal velocity values.

Table 1. Quantitative DNS results for transverse modes (LSA values in brackets).

4.2.3 Mixed-mode instability

At the shallow slope of $\unicode[STIX]{x1D6FC}=5^{\circ }$, if the flow is subjected to more unstable conditions as discussed previously, i.e. $\unicode[STIX]{x1D6F1}_{s}$ is further increased, then according to linear stability theory both the transverse and longitudinal modes become unstable. At the shallow slope angle of $5^{\circ }$, perturbation parameter $\unicode[STIX]{x1D6F1}_{s}$ is chosen to be $19.0$ and we can see from figure 14 that, in both the anabatic and katabatic cases, the maximal growth rates for transverse and longitudinal instability modes are positive. Simulation of both katabatic and anabatic flows for this set of parameters is carried out on a comparatively larger domain of non-dimensional size $100\times 100\times 50$ for $\unicode[STIX]{x1D6FC}=5^{\circ }$ to capture multiple vortex rolls along all directions. Instantaneous snapshots of the flow fields are shown in figure 18(a,b), where the isocontour of the $Q$-criterion at a positive value is used for vortex identification. It is clearly visible that both longitudinal and transverse rolls are simultaneously present in flow fields, and are entangled with each other. However, the major disparity in flow structures between anabatic and katabatic flows, despite having the same dimensionless parameters $\unicode[STIX]{x1D6FC}$, $Pr$ and $\unicode[STIX]{x1D6F1}_{s}$, as evident from the presence of much smaller eddies in the anabatic flow, hints that the near-surface convection rates of momentum and heat will be massively different in the two cases, requiring distinct surface flux parametrisation schemes for each case. At a steeper slope angle of $\unicode[STIX]{x1D6FC}=64^{\circ }$ and with $\unicode[STIX]{x1D6F1}_{s}=30$, using a non-dimensional domain size of $200\times 200\times 50$, the simulated mixed modes for both anabatic and katabatic flows are displayed in figure 18(c,d). In this case, the contrast between anabatic and katabatic flow is not as stark as in the shallow-angle configuration; however, the anabatic flow still contains smaller plume-like structures, which are arranged in a much more disorderly than in the katabatic flow.

Table 2. Quantitative DNS results for longitudinal modes (LSA values in brackets).

Figure 18. Comparison of mixed instability mode for anabatic and katabatic flows at two different slope angles: (a) anabatic flow at $\unicode[STIX]{x1D6FC}=5^{\circ }$ and $\unicode[STIX]{x1D6F1}_{s}=19$; (b) katabatic flow at $\unicode[STIX]{x1D6FC}=5^{\circ }$ and $\unicode[STIX]{x1D6F1}_{s}=19$; (c) anabatic flow at $\unicode[STIX]{x1D6FC}=64^{\circ }$ and $\unicode[STIX]{x1D6F1}_{s}=30$; and (d) katabatic flow at $\unicode[STIX]{x1D6FC}=64^{\circ }$ and $\unicode[STIX]{x1D6F1}_{s}=30$. The $Q$-criterion contour at 4 % of its maximal value is used for visualisation. Contours are coloured by the normalised slope-normal velocity magnitude.