3.1 Dependence on boundary conditions

Before considering all the cases, we first focus on the mixing behaviour for $\unicode[STIX]{x1D6FC}=10^{\circ }$ and examine the differences between no-slip boundary conditions (simulation S10N) and free-slip boundary conditions (S10). In particular, we will demonstrate: (i) that the outer layer behaviours of S10 and S10N are practically indistinguishable with respect to their mixing behaviour; and (ii) that the mixing characteristics are practically uniform in the outer layer of the flow – something which seems to be a unique feature of temporal gravity currents.

It is clear that free-slip boundaries do not occur in nature, but if the entrainment between the current and the ambient fluid is dominated by outer (effectively shear-driven) processes, the influence of velocity boundary conditions should be small. We do not address here the important question of the interplay between turbulence generated by bottom roughness and turbulence generated by shear-driven processes at the interface between the current and ambient (see Wells et al. (Reference Wells, Cenedese and Caulfield2010) for further discussion), or indeed driving mechanisms of a current other than the buoyancy force associated with a finite (albeit potentially exceptionally small) slope angle.

Figure 2.

Time variation of (a) layer thickness $h$ ; and (b) bulk Richardson number $\mathit{Ri}$ , as defined in (2.3). Scaled wall-normal variation of scaled and spatially averaged (c) $\overline{u}$ ; (d) $\overline{b}$ ; (e) $\overline{u^{\prime }w^{\prime }}$ ; and (f) $\overline{w^{\prime }b^{\prime }}$ at time intervals $t/t_{\ast }=0.7$ . Note that the multiple profiles for $\overline{u}$ and $\overline{b}$ show strong self-similarity, are approximately linear in particular away from the wall, and extend beyond $z/h=1$ . Results for simulations S10N and S10 are shown with solid blue and dashed red lines respectively.

The evolution of $h$ and $\mathit{Ri}$ are shown in figure 2(a,b). The layer thickness $h$ grows a little faster for S10 than for S10N, indicating, perhaps unsurprisingly, slightly higher entrainment for flows with no-slip boundary conditions. The evolution of $\mathit{Ri}$ in time shows an initial growth to a maximum, after which $\mathit{Ri}$ reduces monotonically and approaches a steady state for $t/t^{\ast }>30$ ; the final value of $\mathit{Ri}$ is once again slightly higher for simulation S10N than for simulation S10.

The approach to a constant value of $\mathit{Ri}$ is qualitatively similar to the behaviour observed for turbulent plumes, where plumes attain a constant value of the appropriately defined $\mathit{Ri}$ far away from the source, i.e. they attain ‘pure plume balance’. If the flow has an excess of momentum at the source the plume is referred to as being forced (Morton & Middleton Reference Morton and Middleton1973), while if it has a deficiency of momentum flux it is referred to as being lazy (Hunt & Kaye Reference Hunt and Kaye2001), although it has also been interpreted in a geophysical context as arising from a too large or ‘distributed’ source (Caulfield & Woods Reference Caulfield and Woods1995). In either case, the pure plume balance solution is a global attractor in an unstratified ambient (Caulfield & Woods Reference Caulfield and Woods1998), and the flow inevitably adjusts itself due to the work done by gravity until it attains pure plume balance. Restricting attention to $t/t_{\ast }>7$ (when the flow has transitioned to turbulence), the gravity current is ‘lazy’, approaching ‘pure’ behaviour for $t/t^{\ast }$ greater than (approximately) 30.

The (strong) self-similarity of the spatially averaged velocity $\overline{u}$ and buoyancy $\overline{b}$ profiles are shown in figure 2(c,d), although, unsurprisingly, the solutions are very different near the wall in the ‘inner layer’. In particular, the ‘toe’ observed in the buoyancy profile for S10N is caused by the fact that the turbulence shear production $P_{S}$ defined as

is zero at the velocity maximum and that as a consequence turbulence levels are low. This implies that buoyancy is ‘trapped’ near the (smooth) wall for flows with no-slip boundary conditions, an effect reported by Ellison & Turner (Reference Ellison and Turner1959). Importantly though, the flow profiles sufficiently far from the wall are very similar for the two simulations. Figure 2(e,f) shows the wall-normal transport of the streamwise momentum and buoyancy, respectively. It is apparent that there is a large negative momentum flux for simulation S10N, associated with the particular near-wall dynamics in the inner layer, and in general, the free-slip fluxes are slightly larger than the no-slip fluxes, while there is not such a strong collapse as for the $\overline{u}$ and $\overline{b}$ profiles.

The fundamental budget of interest for mixing is the budget of turbulence kinetic energy $e$ , which for the temporal gravity current is given by

where $P_{S}$ is the shear production and $P_{B}$ is the buoyancy production of turbulent kinetic energy with dissipation rate $\unicode[STIX]{x1D700}=\unicode[STIX]{x1D708}\overline{(\unicode[STIX]{x2202}u_{i}^{\prime }/\unicode[STIX]{x2202}x_{j})^{2}}$ . Since we are considering a flow upon a slope, there is a choice to be made as to the definition of the important mixing parameters known as the flux Richardson number $\mathit{Ri}_{f}$ , the turbulent flux coefficient $\unicode[STIX]{x1D6E4}$ (Osborn Reference Osborn1980) and the gradient Richardson number $\mathit{Ri}_{g}$ . We define these quantities as follows:

For the gradient Richardson number $\mathit{Ri}_{g}$ we have chosen to use the component of the density gradient in the direction of gravity in the numerator, which is a self-consistent generalisation of the conventional interpretation of a gradient Richardson number as a quantification of the relative importance of destabilising shear to stabilising density gradients.

For the two (potentially related) flux quantities $\mathit{Ri}_{f}$ and $\unicode[STIX]{x1D6E4}$ we have chosen to use the full buoyancy production of turbulent kinetic energy, which for non-zero $\unicode[STIX]{x1D6FC}$ (as always considered here) has in general a non-zero along-slope component. This choice means that the energy-based interpretation of these terms must be done with care, particularly for larger values of $\unicode[STIX]{x1D6FC}$ . Indeed, for sufficiently steep angles, it is to be expected that the buoyancy production term will change sign and become negative as the dense flow down the slope will actually lead to an increase in the turbulent kinetic energy, rather than a decrease (ultimately due to an exchange into the potential energy reservoir). It is important to appreciate that the quantities in (3.3) are in general functions of $z$ and $t$ , involving as they do only horizontal averages, and so there is not necessarily a simple way to relate $\unicode[STIX]{x1D6E4}$ and $\mathit{Ri}_{f}$ , or indeed to compare with other studies of related flows, such as the experimental flows with the same slope angle discussed in Odier et al. (Reference Odier, Chen and Ecke2014). (See Wells et al. (Reference Wells, Cenedese and Caulfield2010) for a further discussion of this subtle, and often not appreciated point.)

For the temporal gravity current, both the mean and turbulence kinetic energy are expected to scale with $B_{0}$ and are therefore independent of time in the fully self-similar regime. The turbulence kinetic energy budget will then be in a statistically steady state, and thus the two production terms when appropriately integrated in the wall-normal direction must balance the wall-normal-integrated dissipation rate

Here, the calligraphic scripts denote quantities integrated over the wall-normal direction, i.e.  ${\mathcal{P}}_{S}=\int _{0}^{\infty }P_{S}\,\text{d}z$ , ${\mathcal{P}}_{B}=\int _{0}^{\infty }P_{B}\,\text{d}z$ and ${\mathcal{E}}=\int _{0}^{\infty }\unicode[STIX]{x1D700}\,\text{d}z$ . Note that this implies that

where the angle brackets denote quantities calculated from the integral properties. The quantity $\langle \mathit{Ri}_{f}\rangle$ quantifies the proportion of the energy injected into the flow which has led to increases in the potential energy through irreversible mixing processes, and ${\mathcal{P}}_{B}$ is expected to be a sink of kinetic energy ultimately leading to an increase in the flow potential energy, although particularly for flows where the Boussinesq approximation does not apply, the energy pathways are more convoluted (see Tailleux Reference Tailleux2013; Scotti Reference Scotti2015, for detailed discussions). Understanding and parameterising such partitioning is a central challenge in oceanography (see Ivey, Winters & Koseff Reference Ivey, Winters and Koseff2008; Ferrari & Wunsch Reference Ferrari and Wunsch2009, for reviews).

Figure 3.

Scaled wall-normal variation (over time interval documented in table 1) of (a) gradient Richardson number $\mathit{Ri}_{g}$ ; (b) flux Richardson number $\mathit{Ri}_{f}$ ; and (c) turbulent flux coefficient $\unicode[STIX]{x1D6E4}$ as defined in (3.3) for simulation S10N (red lines) and simulation S10 (blue lines).

It is apparent from figure 3 that the values of these three quantities are quite similar for the two simulations with different wall boundary conditions yet small angle $\unicode[STIX]{x1D6FC}=10^{\circ }$ , except in the near-wall inner layer $z/h\lesssim 0.3$ , and indeed are quite close to constant over $0.3<z/h<1.0$ . The evidently self-similar (and essentially linear) distributions of $\overline{u}$ and $\overline{b}$ naturally lead to constant values of the gradient Richardson number $Ri_{g}$ over much of the (self-similar) depth of the current. The perhaps more surprising observation is that the flux Richardson number $\mathit{Ri}_{f}$ and the (apparently related) turbulent flux coefficient $\unicode[STIX]{x1D6E4}$ are close to constant across the interior of the current for $0.3<z/h<1.0$ . Indeed, although as already noted there is no necessity for the simple relationship in (3.5) to hold, $\unicode[STIX]{x1D6E4}$ is indeed slightly more than $Ri_{f}$ .

There is no immediately obvious reason why $\mathit{Ri}_{f}(z,t)$ (and $\unicode[STIX]{x1D6E4}(z,t)$ ) should be approximately constant and similar to $\mathit{Ri}_{g}(z,t)$ , but it is entirely consistent with the behaviour of the flows in stratified plane Couette flow, where Deusebio et al. (Reference Deusebio, Caulfield and Taylor2015) found a close to linear relationship between appropriate definitions of a gradient Richardson number and flux Richardson number, which applies at all heights within the flow. Furthermore, this is also evidence (as discussed further below) that the turbulent Prandtl number for this flow is an $O(1)$ constant. The turbulent diffusivities for momentum and buoyancy (once again in general functions of $z$ and $t$ ) and the associated turbulent Prandtl number are defined as

and so it is straightforward to establish that

Since the angle $\unicode[STIX]{x1D6FC}$ is relatively small, the term in the square bracket is close to unity and so the observation that the flux Richardson number $\mathit{Ri}_{f}$ is close to constant and similar in magnitude to the constant gradient Richardson number $\mathit{Ri}_{g}$ implies that $Pr_{T}\sim O(1)$ and close to constant (with $z$ and $t$ ). We observe that for all heights in the current sufficiently far away from the wall, the flow ‘chooses’ a balanced coupled configuration for which $\mathit{Ri}_{f}$ and $\unicode[STIX]{x1D6E4}$ are constant across the outer layer. This inherently emergent dynamical behaviour of spatially constant $\mathit{Ri}_{f}$ and $\unicode[STIX]{x1D6E4}$ occurs for all simulations, although $\mathit{Ri}_{f}$ and $\unicode[STIX]{x1D6E4}$ do have a non-trivial dependency on $\unicode[STIX]{x1D6FC}$ as discussed further below, not least because of the increasing influence of the along-slope term in $P_{B}$ as the slope angle increases.

3.2 Dependence on the slope angle

Having established that the boundary conditions have relatively limited influence on turbulent entrainment, at least for such purely buoyancy-driven down-slope gravity currents, we restrict attention to free-slip velocity boundary conditions so that by construction we can focus on outer layer dynamics. We consider simulations with a wide variety of angles $\unicode[STIX]{x1D6FC}$ between 2 $^{\circ }$ and 90 $^{\circ }$ . It is important to appreciate that the simulation with $\unicode[STIX]{x1D6FC}=90^{\circ }$ closely resembles the flow of a temporal plume, and is in a very real sense a qualitatively different flow than a temporal gravity current.

Figure 4.

Properties of the results of free-slip simulations for a range of slope angles, as marked with different line types. Time variation of (a) $h$ ; and (b) $\mathit{Ri}$ . Scaled wall-normal variation of spatially averaged (c) $\overline{u}$ ; (d) $\overline{b}$ ; (e) $\overline{u^{\prime }w^{\prime }}$ ; and (f) $\overline{w^{\prime }b^{\prime }}$ over the time interval documented in table 1.

Figure 4(a) shows the time evolution of the current depth $h$ for simulations with a range of angles, showing that the depth grows approximately linearly with time as expected. It is also clear that $h$ grows more rapidly as the angle $\unicode[STIX]{x1D6FC}$ is increased. As shown in figure 4(b), the bulk Richardson number $\mathit{Ri}$ approaches a constant value once the flow is self-similar, with the asymptotic value decreasing with increasing angle.

The extents of the self-similarity of various spatially averaged characteristics of the flows are shown in figure 4(c–f). In contrast to figure 2, $B_{0}$ and $t$ are used to scale both the wall-normal coordinate and the various flow characteristics instead of $u_{T}$ and $b_{T}$ , so that it is possible to identify the differences in absolute amplitudes between simulations with relatively small angles (i.e. $\unicode[STIX]{x1D6FC}=2^{\circ },5^{\circ },10^{\circ }$ ) and relatively larger angles (i.e. $\unicode[STIX]{x1D6FC}=25^{\circ },45^{\circ },90^{\circ }$ ). Clearly, the currents from simulations with relatively low angles are shallower and have higher amplitudes of $\overline{u}$ and $\overline{b}$ compared to the currents along relatively steeper slopes. Flows along slopes with small $\unicode[STIX]{x1D6FC}$ require more work to be done for an increase in potential energy, and so the flow will accelerate until this can be achieved, and the self-similarity is at least approximately achieved. As is apparent in figure 4(e,f), the turbulent fluxes are also close to self-similar, and it is striking that the maximum amplitude is similar for all values of $\unicode[STIX]{x1D6FC}$ , although the location of the maximum does depend on $\unicode[STIX]{x1D6FC}$ , with the maximum being closer to the wall for flows with smaller slope angle $\unicode[STIX]{x1D6FC}$ .

Figure 5.

Variation with slope angle $\unicode[STIX]{x1D6FC}$ of the numerical data (marked with circles) and theoretical predictions (marked with dashed lines) for (a) $\mathit{Ri}$ , modelled by (4.19); (b) $\widehat{\mathit{Ri}_{g}}$ , modelled by (4.19); (c) $\hat{\mathit{Ri}_{f}}$ , modelled by (4.20); and (d) $\hat{\unicode[STIX]{x1D6E4}}$ , modelled by (4.20). A hat denotes an average over the outer layer interval $z/h=[1/2,1]$ .

In figure 5, the observed dependence of the mixing parameters on the slope angle $\unicode[STIX]{x1D6FC}$ is presented by circles. Here, the hatted quantities denote average values over the outer layer interval $z/h=[1/2,1]$ (see figure 3). The error bars show the difference between the largest and smallest value within the time window $t_{stat}$ over which statistics were gathered (the statistics were gathered at intervals $t/t_{\ast }=0.7$ ). Shown with the dashed line are results from a theoretical model which will be introduced in § 4.2. Clearly, all indicators display similar behaviour in the sense that their values increase as the slope angle reduces. Importantly, though, the values remain bounded and well below the critical linear stability threshold of $\mathit{Ri}=1/4$ . $\widehat{Ri}_{f}$ and $\hat{\unicode[STIX]{x1D6E4}}$ are positive only for the three low-angle (but still finite slope) cases. For $\unicode[STIX]{x1D6FC}>25^{\circ }$ , the along-slope turbulence kinetic energy contribution of $P_{B}$ becomes larger than the perpendicular contribution, and thus both $\widehat{Ri}_{f}$ and $\hat{\unicode[STIX]{x1D6E4}}$ will become negative. As already discussed, in such a regime, the flux Richardson number $\mathit{Ri}_{f}$ and the turbulent flux coefficient $\unicode[STIX]{x1D6E4}$ are not appropriately defined.

3.3 Turbulent entrainment

Turbulent entrainment, as quantified by the entrainment parameter (van Reeuwijk et al. Reference van Reeuwijk, Krug and Holzner2018)

is naturally closely related to turbulent mixing. In figure 6(a), we plot the temporal variation of $E$ for all six simulations. Clearly, $E$ increases with $\unicode[STIX]{x1D6FC}$ . Since $\unicode[STIX]{x1D6FC}$ is directly related to the bulk Richardson number $\mathit{Ri}$ through its definition (2.3), this relationship suggests, analogously to the classic parameterisation of Ellison & Turner discussed in § 1, that $E$ has a functional dependence on $\mathit{Ri}$ . Despite keeping $B_{0}$ the same for all simulations, quantifying $E$ in the simulations at the smallest (yet crucially still finite, as noted above) angles is challenging since there is a lower turbulence intensity (not shown) as a result of the (relatively) stronger stratification. This is particularly evident in the data extracted from simulation S2, in which the flow almost completely relaminarises after the initial burst associated with Kelvin–Helmholtz-like shear instabilities as the flow accelerates from its initial condition. However, the very low, but inherently transient, level of turbulence (evident for $t/t_{\ast }\simeq 10$ ) implies that the fluid layer will accelerate due to the low level of dissipation until a second transition occurs after which ‘pure’ behaviour emerges for $t/t_{\ast }>40$ . This second transition appears to lead to a quasi-stationary turbulent state. Therefore, it is not appropriate to interpret this flow as exhibiting ‘intermittent’ turbulence, but rather that our chosen initial conditions are in some sense imbalanced, and so there is a relatively long initial transient evolution, exhibiting non-monotonic turbulence intensity, until the flow self-organises into its eventual quasi-stationary turbulent state.

Figure 6.

(a) Variation with time of the entrainment parameter $E$ as defined in (3.8) for various simulations, plotted with different line types. (b) Variation with $\mathit{Ri}$ of the averaged entrainment parameter $E$ as defined in (3.8) for various simulations, plotted with different symbols. The classical empirical parameterisation ET59 as defined in (1.1) is plotted with a solid line, while the new proposed parameterisation (3.9) is plotted with a dashed line.

The entrainment law $E(\mathit{Ri})$ is shown in the bottom plot of figure 6, together with the Ellison and Turner entrainment law defined in (1.1). Here, $\mathit{Ri}$ is the value associated with the ‘pure’ gravity current after initial transients have decayed, and the velocity and buoyancy distributions have evolved into the clearly coupled (and close to linear) distributions as shown in figure 4(c). Precisely, the value of $E$ was calculated by averaging over the interval $t_{stat}$ (table 1), and the horizontal and vertical error bars denote the maximum and minimum observed values of $\mathit{Ri}$ and $E$ over this time interval.

Interestingly, $\mathit{Ri}$ remains constrained between $0<\mathit{Ri}<0.2$ for all angles. Crucially, when defined in this way, we observe that it is impossible to reach the higher values of $\mathit{Ri}$ observed in the Ellison & Turner (Reference Ellison and Turner1959) experiments: the coupled velocity and buoyancy distributions inevitably force the Richardson number to be bounded above by $0.2$ when there is any observable turbulence at all, consistently with the recent observations of Deusebio et al. (Reference Deusebio, Caulfield and Taylor2015) and Zhou et al. (Reference Zhou, Taylor and Caulfield2017) in plane Couette flow. Furthermore, the values of $E$ as quantified from our numerical simulations are lower than those reported experimentally by Ellison & Turner (Reference Ellison and Turner1959) and Odier et al. (Reference Odier, Chen and Ecke2014) although direct quantitative comparison is challenging, not least due to the different forcing mechanisms of the different flows.

Of course, the spatially evolving experimental flow geometry is qualitatively different from the inherently temporal flow geometry considered here. It is also at least plausible that the gravity current in their experiments was not yet ‘pure’ in a statistically steady state, which would go some way to explain their higher observed values for $\mathit{Ri}$ and $E$ . Furthermore, as noted in § 1, it is also conceivable that their experimental data were contaminated by end-effect mixing at a ‘head’ being swept back along the current, as such mixing is expected to be substantially stronger, as investigated in detail by Sher & Woods (Reference Sher and Woods2015). Whatever the reason for the mismatch, it is apparent that there is a marked decrease in $E$ with $\mathit{Ri}$ , and there is a clear ‘switch-off’ of entrainment (quantified in this way) at a maximum value of  $\mathit{Ri}$ .

Focusing on the simulations S2, S5 and S10 with relatively small slope angles, (where the buoyancy production $P_{b}<0$ , and so the expected sink of turbulent kinetic energy due to mixing occurs) the entrainment law can be reasonably well approximated by a linear relation of the form

as marked in figure 6(b) by the dashed line. We present a detailed theoretical justification for this model in § 4.2, but we once again stress that for these simulations, the mean velocity and buoyancy profile are both close to linear and appear to be coupled in a way which strongly suggests that there is a maximum value of $\mathit{Ri}$ beyond which turbulence cannot be sustained, at least in the flows considered here, purely driven by buoyancy forces down slopes of small (but still finite) angle.

4.1 Turbulence parameterisation

As demonstrated in Krug et al. (Reference Krug, Holzner, Marusic and van Reeuwijk2017), the vertical length scale is proportional to $\hat{e}^{1/2}/{\hat{S}}$ , which implies that the turbulence is ‘shear dominated’ (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014). It is thus natural to investigate whether a turbulence parameterisation can be constructed using the strain rate $S=|\unicode[STIX]{x2202}\overline{u}/\unicode[STIX]{x2202}z|$ and the turbulence kinetic energy $e$ to scale the eddy diffusivity calculated from three simulations with relatively small slope angles, i.e. the simulations S2, S5 and S10. For the turbulent diffusivities defined in (3.6), this suggests that

In figure 7(a), we plot $K_{m}S/e$ for all three small-slope simulations as a function of $z/h$ at a range of times. The vertical dashed line shows $c_{m}=0.25$ , demonstrating that this parameterisation is an excellent fit for a wide range of times for all three simulations, particularly sufficiently far away from the wall, i.e.  $0.3<z/h<1$ . Conversely, in figure 7(d) by plotting $K_{m}/(u_{T}h)$ , we demonstrate that there is no collapse when scaled with the similarity variables. Similarly, for the turbulent scalar diffusion $K_{\unicode[STIX]{x1D70C}}$ , we observe a collapse upon plotting $K_{\unicode[STIX]{x1D70C}}S/e$ , as shown in figure 7(b), with an approximately constant value sufficiently far away from the wall of $c_{\unicode[STIX]{x1D70C}}=0.31$ (as plotted with a vertical dashed line). This implies that for three lowest-angle cases, the turbulent Prandtl number $\mathit{Pr}_{T}=K_{m}/K_{\unicode[STIX]{x1D70C}}=c_{m}/c_{\unicode[STIX]{x1D70C}}=0.8$ independent of the flow angle, thus clearly suggesting that $Pr_{T}$ remains finite, and indeed is of $O(1)$ . Once more there is no collapse upon scaling $K_{m}$ by $u_{T}h$ , as shown in figure 7(e).

Figure 7.

Variation with scaled height $z/h$ for the three simulations S2 (blue lines), S5 (red lines) and S10 (yellow lines) at several time instants (see table 1 for time interval) for (a) $K_{m}S/e$ ; (b) $K_{\unicode[STIX]{x1D70C}}S/e$ ; (c) $\unicode[STIX]{x1D700}/(eS)$ ; (d) $K_{m}/(u_{T}h)$ ; (e) $K_{\unicode[STIX]{x1D70C}}/(u_{T}h)$ ; (f)  $\unicode[STIX]{x1D700}h/u_{T}^{3}$ . In (a–c) there is clear collapse, and the relevant empirical constants of proportionality $c_{m}=0.25$ , $c_{\unicode[STIX]{x1D70C}}=0.31$ and $c_{\unicode[STIX]{x1D700}}=0.21$ are marked with vertical dashed lines.

Similarly, we predict that the dissipation rate $\unicode[STIX]{x1D700}$ can also be scaled appropriately with $e$ and the strain rate $S$ , which on dimensional grounds implies that

for some constant $c_{\unicode[STIX]{x1D700}}$ . We plot $\unicode[STIX]{x1D700}/(eS)$ in figure 7(c), and once again there is a clear collapse, with an approximate value of $c_{\unicode[STIX]{x1D700}}=0.21$ away from the immediate vicinity of the wall, as marked on the figure with a dashed line. Conversely, as shown in figure 7(f), there is no collapse of $\unicode[STIX]{x1D700}$ when scaling with $u_{T}^{3}/h$ .

It is also interesting to note, using the definitions of $\mathit{Ri}_{f}$ , $\unicode[STIX]{x1D6E4}$ , $P_{S}$ (3.3) and $K_{m}$ (3.6), as well as the (verified) parameterisation (4.1), that

Comparing this relationship to the empirical observation that $c_{\unicode[STIX]{x1D700}}=0.21$ and $c_{m}=0.25$ , this suggests $\widehat{\mathit{Ri}}_{f}/\hat{\unicode[STIX]{x1D6E4}}=0.84$ and thus that $\widehat{\mathit{Ri}}_{f}<\hat{\unicode[STIX]{x1D6E4}}$ , consistently with (3.5).

It is apparent that for the key properties of the turbulence in this flow, the appropriate scalings are not based around the conventional integral velocity scale $u_{T}$ and $h$ . This is clearly demonstrated in figure 8(a), where $u_{T}$ and $e_{T}^{1/2}$ are displayed as a function of the slope angle $\unicode[STIX]{x1D6FC}$ . As mentioned earlier, these quantities are both expected (and observed) to scale according to $B_{0}^{1/2}$ in the self-similar regime. It is striking that $e_{T}$ is virtually independent of $\unicode[STIX]{x1D6FC}$ , but $u_{T}$ increases dramatically as the slope angle reduces. Indeed, the structure of the $u_{T}$ variation is consistent with $u_{T}$ diverging as $\unicode[STIX]{x1D6FC}\rightarrow 0$ . As is apparent from the definition (2.2) of the integral forcing $B_{0}$ , the forcing buoyancy $b_{0}h_{0}$ must diverge as $\unicode[STIX]{x1D6FC}\rightarrow 0$ to ensure the distinguished limit that $B_{0}$ remains constant, and so it is perhaps unsurprising that $u_{T}$ becomes very large as $\unicode[STIX]{x1D6FC}$ approaches small values. The physical picture that emerges is that the flows with higher Richardson numbers for lower slope angles $\unicode[STIX]{x1D6FC}$ require higher shear (and thus higher $u_{T}$ ) to sustain a similar level of turbulence. This will be explained in detail by the conceptual model we develop in the next section. What will not be revealed by the model is why the turbulence level $e_{T}/B_{0}^{1/2}$ is virtually independent of $\unicode[STIX]{x1D6FC}$ . In figure 8(b), we plot the variation of $(3/2)(\overline{w^{\prime }w^{\prime }}/e)$ with slope angle $\unicode[STIX]{x1D6FC}$ to quantify the extent of anisotropy within the flow. Unsurprisingly, due to the inherent greater magnitude of the streamwise velocity $u$ , this quantity is always less than one, and exhibits only a slight dependence on angle when scaled in this way, decreasing slightly as $\unicode[STIX]{x1D6FC}$ tends to small values, due naturally to the stabilising effect of stratification becoming somewhat more significant. However, there are no signs of a transition to strongly stratified turbulence, and this cannot be an explanation for the observation of $E\rightarrow 0$ as $\unicode[STIX]{x1D6FC}$ tends towards zero.

In figure 8(c,d), we plot the variation with slope angle $\unicode[STIX]{x1D6FC}$ of the three key terms ${\mathcal{P}}_{S}$ , $4{\mathcal{P}}_{B}$ (the prefactor is included for visual purposes only) and ${\mathcal{E}}$ in the integral turbulent kinetic energy equation (3.4) and averaged over the time interval $t_{stat}$ (table 1). In figure 8(c), the various terms are scaled using $(B_{0}^{3})^{-1/2}$ , which reveals the appropriate and expected dynamics for small angles, showing that the actual production of turbulence increases. This can be understood from the fact that the velocity gradient is much larger (see figure 4) as the flow will keep accelerating until an essentially equilibrium state is attained. As mentioned earlier, the buoyancy production term ${\mathcal{P}}_{B}$ is negative for the smallest three angles and changes sign around $\unicode[STIX]{x1D6FC}=20^{\circ }$ , which can be understood in terms of the relative size of the two terms in ${\mathcal{P}}_{B}$ (3.2), with the contribution of the first (generally positive) term dominating for large $\unicode[STIX]{x1D6FC}$ and the second (generally negative) term for small $\unicode[STIX]{x1D6FC}$ . Clearly, not least because of the opposing effects of these two terms, the main overall balance is between shear production ${\mathcal{P}}_{S}$ and dissipation ${\mathcal{E}}$ , and ${\mathcal{P}}_{B}$ only makes a relatively small contribution to the overall budget.

Although it is apparent from figure 8(c) that the magnitude of the various key terms in the integral turbulent kinetic energy equation (3.4) increase as $\unicode[STIX]{x1D6FC}$ decreases, their rate of increase is clearly significantly smaller than that exhibited by $u_{T}$ as shown in figure 8(a). Figure 8(d) shows once more the integral turbulence kinetic energy budget but this time scaled by $u_{T}^{3}$ . Scaled in this way, as $\unicode[STIX]{x1D6FC}\rightarrow 0$ the production of turbulent kinetic energy appears to go to zero, although it is more appropriate to interpret this figure as implying that the production of turbulent kinetic energy becomes effectively insignificant as $\unicode[STIX]{x1D6FC}\rightarrow 0$ . Significantly, even though $P_{B}$ increases in terms of $B_{0}$ , the relative production tends to zero as $\unicode[STIX]{x1D6FC}\rightarrow 0$ . This has profound implications for the entrainment as there will be no energy available for turbulent mixing.

Figure 8.

Variation with slope angle of (a) the characteristic velocities $u_{T}$ and turbulence $e_{T}^{1/2}$ , (b) the anisotropy of the turbulence and (c,d) the vertically integrated terms in the turbulent kinetic energy equation: shear production ${\mathcal{P}}_{S}$ (marked with blue triangles); dissipation rate ${\mathcal{E}}$ (red inverted triangles); buoyancy production $4{\mathcal{P}}_{B}$ (yellow circles). In (c) these quantities are scaled with $(B_{0}^{3})^{1/2}$ , while in (d) they are scaled with $u_{T}^{3}$ .

4.2 An approximate integral model for small $\unicode[STIX]{x1D6FC}$

Since we observe apparently coupled and essentially linear profiles of velocity and buoyancy, we now construct a model for entrainment based around this central observation, in an attempt to predict the apparent ‘switch-off’ of the entrainment, particularly for the simulations associated with a low slope angle. The conservation equations for mean momentum, buoyancy and turbulent kinetic energy are given by, respectively

Here, we have neglected the term $\overline{u^{\prime }b^{\prime }}\sin \unicode[STIX]{x1D6FC}$ describing the along-slope turbulent production by buoyancy in the turbulent kinetic energy equation. This limits the validity of the approximation to small values of $\unicode[STIX]{x1D6FC}$ , as we are interested in the dynamical behaviour of the system where the dominant effect of the buoyancy is to extract kinetic energy due to vertical mixing associated with entrainment. The system of equations is closed using our parameterisations of the eddy diffusivities of momentum and buoyancy, and the scaling of the dissipation rate in terms of the turbulent kinetic energy and the mean shear, i.e. (4.1) and (4.2).

Furthermore, since we observe that the profiles of $u$ and $b$ appear to be self-similar, we assume a solution of the form

where the scaled wall-normal distance $\unicode[STIX]{x1D702}=z/h$ is the natural similarity variable and $a_{u}$ , $a_{b}$ and $a_{e}$ are coefficients that will need to be determined. Given the relative complexity of the equations we will not attempt to obtain closed-form solutions. Instead we employ the von Kármán–Pohlhausen method (Lighthill Reference Lighthill1950; Spalding Reference Spalding1954; Schlichting & Gersten Reference Schlichting and Gersten2000) to construct an approximate solution to the system of coupled partial differential equations using an ansatz

on the interval $\unicode[STIX]{x1D702}\in [0,\unicode[STIX]{x1D702}_{1}]$ . These profiles are constructed such that $\int _{0}^{\infty }f_{u}\,\text{d}\unicode[STIX]{x1D702}=\int _{0}^{\infty }f_{b}\,\text{d}\unicode[STIX]{x1D702}=\int _{0}^{\infty }f_{e}\,\text{d}\unicode[STIX]{x1D702}=1$ in order to enforce consistency with $u_{T}$ , $b_{T}$ and $e_{T}$ .

Using (4.7) and (4.8), the key integrals of the volume, momentum and buoyancy in the current take the forms

From the fact that the buoyancy integral is conserved and equal to $-B_{0}$ , it follows directly that $a_{b}=1$ . From the definition of $h$ (2.1), it follows that $\unicode[STIX]{x1D702}_{1}=4/3$ . Making the further scaling assumption that $h=aB_{0}^{1/2}t$ , for a further empirical constant $a$ , we obtain the following expressions for the main quantities of interest:

Here, the hat emphasises that these quantities represent the expected average value in the outer layer.

Equations for the various empirical constants can be determined by integrating (4.4) and (4.6) from $\unicode[STIX]{x1D702}=0$ to $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{1}$ , leading ultimately to

Analogous integration of the buoyancy equation (4.5) does not provide any further information, confirming that the chosen profiles are a priori consistent with (4.5). Solving for $a_{u}$ and $a^{2}$ results in

Interestingly, the empirical constant $a_{e}$ remains free, implying that these results are not directly sensitive to the actual magnitude of the turbulent kinetic energy fluctuations.

For the mixing parameters, the model thus predicts that

Note that these expressions are fully consistent with relations (3.7), (4.3) derived above in §§ 3 and 4.1 respectively. Figure 5 shows how these theoretical predictions compare to the simulation data.

Both $\mathit{Ri}$ and $\widehat{\mathit{Ri}}_{g}$ compare well across the entire range of simulations, even – somewhat surprisingly – for the large angle cases. This suggests that the along-slope buoyancy production is not significant for these two parameters, although it also could just be a straightforward consequence of the fact that both $\mathit{Ri}$ and $\hat{\mathit{Ri}}_{f}$ tend to zero as $\unicode[STIX]{x1D6FC}\rightarrow 90^{\circ }$ . Conversely, $\widehat{\mathit{Ri}}_{f}$ and $\hat{\unicode[STIX]{x1D6E4}}$ only really compare well to the numerical simulation data for the three relatively low-angle cases. For $\unicode[STIX]{x1D6FC}>25^{\circ }$ , the theoretical and numerically simulated values start deviating strongly, due presumably to the omission of the along-slope turbulent kinetic energy production due to buoyancy, as is clear from the fact that both quantities change sign for the numerical simulation data. As already discussed, in such a regime, the flux Richardson number $\mathit{Ri}_{f}$ and the turbulent flux coefficient $\unicode[STIX]{x1D6E4}$ are not appropriately defined.

The model predicts that the maximal attainable values for the mixing parameters occur as $\unicode[STIX]{x1D6FC}\rightarrow 0$ and are given by

The entrainment parameter $E=a^{2}$ , and the angle $\unicode[STIX]{x1D6FC}$ will need to be eliminated in order to obtain the entrainment law. This can be achieved by using the equation for $\mathit{Ri}$ in (4.19), which can be rearranged to yield

Substituting this relationship into the expression for $a^{2}$ (4.18), it follows that the entrainment parameter $E$ is given by

Using the empirically determined values $c_{\unicode[STIX]{x1D700}}=0.21$ , $c_{m}=0.25$ and $\mathit{Pr}_{T}=0.8$ , we find that

Interestingly, this entrainment prediction is in very good agreement with the simulation results for the low-angle cases as is clear from figure 6(b), strongly suggesting that this model is appropriate, and that there is indeed a (particular kind of) ‘switch-off’ of entrainment in these flows for sufficiently large $\mathit{Ri}$ . Crucially, this particular ‘switch-off’ is not associated with a suppression of the turbulence: for all cases, the turbulence levels $e_{T}$ were similar, and finite. Rather it is the scaled intensity of the turbulence $e_{T}^{1/2}/u_{T}$ which tends to zero as $\mathit{Ri}\rightarrow \mathit{Ri}_{max}$ . This key phenomenology is captured by the theoretical model, which predicts that for $\unicode[STIX]{x1D6FC}\ll 1$ , $e_{T}\sim B_{0}$ and $u_{T}\sim B_{0}^{1/2}\unicode[STIX]{x1D6FC}^{-1/2}$ , from which it follows that $e_{T}^{1/2}/u_{T}\rightarrow 0$ in the limit of $\unicode[STIX]{x1D6FC}\rightarrow 0$ . This in itself can be understood by the fact that the flow will keep accelerating such that the characteristic integral velocity scale $u_{T}$ defined in (2.1) increases until the shear is sufficiently large to produce a sufficient amount of turbulence to be steady. It is apparent that there is a critical Richardson number beyond which this steady balance cannot be attained. Once again, it is important to remember that this ‘switch off’ is not associated with a complete suppression of turbulence, but rather that the turbulence becomes proportionally very small compared to the characteristic scale of velocity of the current. Therefore, though analyses of such quantities are beyond the scope of this paper, other quantities associated with turbulence, such as the buoyancy flux and turbulent drag also remain finite as this ‘switch off’ is approached. Critically, however, the central point remains that quasi-steady turbulence cannot apparently be sustained for flows with characteristic Richardson number larger than the critical value of $\mathit{Ri}_{max}=0.15$ .