2.1. Simulation set-up

We consider a negatively buoyant gravity current flowing down a slope of constant angle $\alpha$ after the passage of any transient ‘head’, such that the current exhibits slow vertical growth induced by entrainment at its upper interface, as shown in figure 1. Periodic boundary conditions are imposed for all flow variables on the lateral boundaries of a finite-sized computational domain. Consequently, the simulations are statistically homogeneous in the streamwise $(x)$ and spanwise $(y)$ directions, but evolve with time. The simulation set-up follows the framework established by van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019), with the exception of the bottom boundary condition, which in this study is specified as no-slip rather than free-slip. This set-up leads to the evolution of a temporal gravity current, resulting in significant computational savings compared to simulations of a spatially evolving gravity current, especially for shallow-angle cases involving a long evolution process. A detailed description of temporal gravity currents is provided in van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019).

If the flow is assumed to be Boussinesq, then the governing equations in the coordinate system in figure 1 may be written as

(2.1) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}&=-\rho _a^{-1}\,\boldsymbol{\nabla }p+v\,{\nabla} ^2\boldsymbol{u} +b\boldsymbol{e}, \end{align}

(2.2) \begin{align} \frac {\partial b}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }b &=\kappa\, {\nabla} ^2 b, \end{align}

(2.3) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &=0, \end{align}

where $ \boldsymbol{u}=(u,v,w)$ is the velocity vector, $p$ is the pressure, $b=(\rho _a-\rho )g/\rho _a$ is the buoyancy, and $\rho _a$ is a reference density (taken to be that of the ambient fluid). The vertical unit vector resolved in the coordinate system is $\boldsymbol{e} = (-\sin {\alpha }, 0, \cos {\alpha })$ , and $\nu$ and $\kappa$ are the kinematic viscosity and diffusivity, respectively.

We solve numerically the governing equations using an in-house DNS code SPARKLE, which employs a conservative fourth-order-accurate differencing scheme (Verstappen & Veldman Reference Verstappen and Veldman2003) for spatial discretisation, and an adaptive third-order Adams–Bashforth scheme for explicit time advancement. The code is described in detail by Craske & van Reeuwijk (Reference Craske and van Reeuwijk2015) and has been widely used in simulations of gravity currents (Krug et al. Reference Krug, Holzner, Marusic and van Reeuwijk2017; van Reeuwijk et al. Reference van Reeuwijk, Holzner and Caulfield2019; Dieu Reference Dieu2020). The grid size $\Delta x$ of the domain varies between cases to ensure $\Delta x/\eta _K \lt 3/2$ . Here, $\eta _K$ is a characteristic Kolmogorov length scale defined as $(\nu ^3/\varepsilon _T)^{1/4}$ , where $\varepsilon _T=h^{-1}\int \varepsilon\, {\textrm{d}} z$ is a characteristic dissipation rate. Here, $h$ is a length scale defined in (2.10a ). The initial conditions ( $t=0$ ) can be written as

(2.4) \begin{equation} u(x,y,z) = \begin{cases} u_0, & z \leqslant h_0, \\ 0, & z \gt h_0, \end{cases}\quad v(x,y,z)=0,\quad w(x,y,z)=0,\quad b= \begin{cases} b_0, & z \leqslant h_0, \\ 0, & z \gt h_0, \end{cases} \end{equation}

where $u_0$ and $h_0$ each maintain the same value across all simulations, and $b_0$ is varied to keep $B_0 = \int -b_0h_0\sin \alpha\, {\textrm{d}} z$ constant across the cases. The initial profiles are shown in figure 5(a,b) below. Small random perturbations are introduced to the initial velocity field to trigger turbulence. The computational domain size is $20h_0 \times 20h_0 \times 20h_0$ for all cases. In the streamwise direction, the domain spans approximately ten integral length scales $L_T$ on average, ensuring sufficient data for statistical analysis. Here, $L_T = e_T^{3/2} / \varepsilon _T$ , with $e_T$ denoting the characteristic TKE as defined in (2.11). A free-slip wall is imposed on the top boundary to maintain compatibility with the temporal framework. The vertical extent of the domain is significantly larger than the current height so that the top boundary condition has negligible influence on the current itself. Further details are provided in table 1.

2.2. Characteristic quantities

Given the statistical homogeneity in the $x$ and $y$ directions, we spatially average the governing equations, and write the Reynolds-averaged momentum and buoyancy equations as

(2.5) \begin{align} &\frac {\partial \overline {u}}{\partial t} +\frac {\partial \overline {w^{\prime}u^{\prime}}}{\partial z}=v\frac {\partial ^2\overline {u}}{\partial z^2}-\overline {b}\sin \alpha , \end{align}

(2.6) \begin{align} &\frac {\partial \overline {b}}{\partial t} + \frac {\partial \overline {w^{\prime}b^{\prime}}}{\partial z}=\kappa \frac {\partial ^2\overline {b}}{\partial z^2}, \end{align}

where $\overline {*}=\iint *\,{\textrm{d}} x\, {\textrm{d}} y/(L_xL_y)$ represents the spatial averaging operator for the quantity $\ast$ , $L_x$ and $L_y$ are the dimensions of the domain in the $x$ and $y$ directions, respectively, and a prime represents the departure from the corresponding average, i.e. $*^{\prime} = * - \overline {*}$ . Taking the dot product of $\boldsymbol{u}$ with the momentum (2.1), subtracting the mean kinetic energy ( $\overline {\boldsymbol{u}}\boldsymbol{\cdot }\overline {\boldsymbol{u}}/2$ ) and averaging over $x$ and $y$ directions, we obtain the TKE budget

(2.7) \begin{equation} \frac {\partial e}{\partial t} = - \overline {w^{\prime}u^{\prime}}\, \frac {\partial \overline {u}}{\partial z} + \overline {w^{\prime}b^{\prime}} \cos \alpha - \overline {u^{\prime}b^{\prime}} \sin \alpha - \varepsilon , \end{equation}

where $e=\overline {{u^{\prime}_{i}}^2}/2$ is the TKE, $- \overline {w^{\prime}u^{\prime}} ({\partial \overline {u}}/{\partial z})$ is the shear production of TKE, and $\varepsilon =\nu\, \overline {(\partial u_i^{\prime}/\partial x_j)^2}$ is the dissipation rate of TKE. Note that the TKE transport terms are neglected.

Table 1. Simulation details: $\textit{Ri}_0 = -b_0h_0\cos \alpha /u_0^2$ is the initial Richardson number, where $b_0, h_0, u_0$ are the initial buoyancy, layer thickness and velocity, respectively (see (2.4)); $\textit{Ri}_\infty$ represents the stabilised value of Richardson number $\textit{Ri}$ when the flow is fully developed, where $\textit{Ri}$ is defined in (2.10d ); the Reynolds number ${{Re}}_\tau =u_\tau z_{um}/\nu$ characterises the inner layer, where $u_\tau =\sqrt { \nu\, ({\partial \overline {u}}/{\partial z})|_{z=0}}$ is the friction velocity, and $z_{um}$ is the vertical coordinate of the velocity maximum; $t_{ave}$ is a time interval towards the end of the simulation over which the numerical results are averaged; and $t^* = h_0/\sqrt {B_0}$ is a typical time scale, where $B_0 = \int -b_0h_0\sin \alpha\, {\textrm{d}} z$ is the initial buoyancy forcing. The initial Reynolds number ${Re}_0 = u_0h_0/\nu$ is 3800, and the Prandtl number $Pr_T = \nu /\kappa$ is 1.

The integral volume flux $Q$ , momentum flux $M$ and integral buoyancy forcing $B$ of the gravity current are defined here as

(2.8a–c) \begin{align} Q=\int _0^\infty \overline {u}\,{\textrm{d}} z,\quad M=\int _0^\infty \overline {u}^2\,{\textrm{d}} z,\quad B=\int _0^\infty -\overline {b}\sin {\alpha }\,{\textrm{d}} z. \end{align}

We decompose these quantities into inner and outer components, denoted with subscripts $i$ and $o$ , respectively, i.e.

(2.9a–c) \begin{align} \begin{aligned} Q =&\underbrace {\int _0^{z_{um}}\overline {u}\, {\textrm{d}} z}_{Q_i} + \underbrace {\int _{z_{um}}^\infty \overline {u}\, {\textrm{d}} z}_{Q_o},\quad M = \underbrace {\int _0^{z_{um}}\overline {u}^2\,{\textrm{d}} z}_{M_i} + \underbrace {\int _{z_{um}}^\infty \overline {u}^2\,{\textrm{d}} z}_{M_o},\\ B=& \underbrace {\int _0^{z_{um}}-\overline {b}\sin {\alpha }\,{\textrm{d}} z}_{B_i} + \underbrace {\int _{z_{um}}^\infty -\overline {b}\sin {\alpha }\,{\textrm{d}} z}_{B_o}, \end{aligned} \end{align}

where $z_{um}$ is the vertical coordinate of the velocity maximum. Note that buoyancy is conserved in the flow, thus $B$ remains constant (equal to $B_0$ ). With $B_0$ prescribed identically across all cases (by adjusting $b_0$ ), $B$ is therefore invariant over the full range of slope angles. The characteristic velocity scale $u_{T*}$ , layer thickness $h_*$ , buoyancy $b_{T*}$ and bulk Richardson number $\textit{Ri}_*$ are defined as

(2.10a–d) \begin{align} h_*=\frac {Q_*^2}{M_*},\quad u_{T*}=\frac {Q_*}{h_*},\quad b_{T*}=-\frac {B_*}{h_*\sin \alpha },\quad Ri_*=\frac {-b_{T*}h_*\cos \alpha }{u_{T*}^2}=\frac {B_*}{u_{T*}^2\tan \alpha }, \end{align}

respectively, where the subscript $*$ is either omitted, $i$ or $o$ , and is used to characterise the entire current, the inner layer or the outer layer, respectively. In a similar manner, the characteristic scales for TKE $e_{T*}$ are given by

(2.11) \begin{equation} \underbrace {\int _0^\infty e\, {\textrm{d}} z}_{e_Th} = \underbrace {\int _0^{z_{um}} e\, {\textrm{d}} z}_{e_{Ti}h_i} + \underbrace {\int _{z_{um}}^\infty e\, {\textrm{d}} z}_{e_{To}h_o}. \end{equation}

Figure 2. Temporal evolution of (a) dimensionless TKE: $e/B_0$ for 1N, where the dotted line denotes the boundary between the inner and outer layers, together with the instantaneous buoyancy field $\overline {b}/b_0$ and profiles of horizontally averaged velocity and buoyancy, $\overline {u}/u_{To}$ and $\overline {b}/b_{To}$ , at (b) $t/t^* = 69.5$ , (c) $t/t^* = 83.25$ , (d) $t/t^* = 90$ and (e) $t/t^* = 101$ . Temporal evolution of normalised (f) $e_T, {e_T}_i,{e_T}_o$ , (g) $u_T, {u_T}_i, {u_T}_o$ , (h) $b_T, {b_T}_i, {b_T}_o$ (i) $h, h_i, h_o$ , (j) $B_i, B_o$ and (k) $\textit{Ri}, Ri_i, Ri_o$ . Note that the subscripts $o$ and $i$ in the legend denote the results of the outer layer (displayed with symbol $\circ$ ) and the inner layer (displayed with symbol $\times$ ), respectively.

3.1. Inner–outer flow dynamics

As we observe essentially similar evolution processes for all the slope angles considered, we use case 1N as an example to illustrate the dynamics. Figure 2(a) shows the temporal evolution of $e(z,t)$ , normalised by $B_0$ , against $t/t^*$ for case 1N, where $t^*=h/\sqrt {B_0}$ is a typical time scale. Also shown is the location of the velocity maximum $z_{um}$ normalised by $h_0$ (black dashed line). Instantaneous snapshots at different time intervals of the normalised buoyancy field in figure 2(b–e) show the development of turbulent structures in the flow. The evolution of layer-specific characteristic flow variables for case 1N are presented in figure 2(f–k); as above, an omitted subscript, $i$ or $o$ is used to denote the overall current, the inner layer or the outer layer, respectively.

An intense initial burst of turbulence associated with shear instabilities is observed for $20\lt t/t^*\lt 40$ in figure 2(a), caused by the sharp initial acceleration (figure 2 g) from the initial conditions. This initial burst leads to a noticeable plunge in velocity (see figure 2 g), as the mean flow kinetic energy is converted to TKE and potential energy (see the increase of $h_*$ in figure 2 i). Consequently, large bulk Richardson numbers arise after the initial burst (see figure 2 k) as damping of turbulence and even relaminarisation occurs in the outer layer (see figure 2(a), $50\lt t/t^*\lt 80$ ). Meanwhile, turbulence is sustained in the inner layer over the whole evolution process.

During the period of damping ( $50\lt t/t^*\lt 80$ ), the outer layer again accelerates (figure 2 g) due to the buoyancy forcing, leading to increased shear and reduced Richardson number (see figure 2 k). The outer layer eventually transitions to a turbulent state as the shear instabilities overcome the restoring stratification. The temporal sequence of instantaneous buoyancy fields (figure 2 b−e) illustrates the transition to a turbulent regime in the outer layer, which initiates with the onset of instabilities near the velocity maximum, followed by the growth of eddies and vortices. Figure 2(f) quantifies the turbulence level throughout the evolution process, showing the first burst of turbulence and subsequent damping (confined to the outer layer), followed by the eventual transition to a nearly constant turbulence level.

Restricting our attention to the time period $t/t^*\gt 100$ , it is noteworthy that the characteristic velocities ${u_T}_i$ , ${u_T}_o$ and $u_T$ attain nearly constant values (see figure 2 g). This behaviour is consistent with the so-called equilibrium state commonly assumed to exist for inclined gravity currents, in which a bulk force balance is achieved (Ellison & Turner Reference Ellison and Turner1959; Britter & Linden Reference Britter and Linden1980; Odier, Chen & Ecke Reference Odier, Chen and Ecke2014; Martin, Negretti & Hopfinger Reference Martin, Negretti and Hopfinger2019).

The buoyancy variables ${b_T}_i$ , ${b_T}_o$ and $b_T$ shown in figure 2(h) continue to reduce gradually as entrainment of ambient fluid continues to dilute the current and increase the layer thickness (figure 2 i). However, the integral buoyancy forcings $B_o$ and $B_i$ (figure 2 j) are remarkably invariant, i.e. the total buoyancy in each of the inner and the outer layers is approximately conserved. This behaviour is closely linked to the interaction between the two layers, and is discussed in § 5.1. The bulk Richardson numbers shown in figure 2(k) also attain approximately constant values in the turbulent regime. We will term this state, in which multiple flow parameters take constant characteristic values, the dynamically equilibrated regime, which implies self-similar behaviour in the outer layer and a quasi-steady state in the inner layer; see §§ 4 and 6, respectively.

Figure 3. Temporal evolution of overall quantities $(a)\, h/h_0$ and $(b)\, Ri$ , outer-layer quantities $(c)\, h_o/h_0$ and $(d)\, Ri_o$ , and inner-layer quantities $(e)\, h_i/h_0$ and $(f)\, Ri_i$ .

3.2. Slope angle dependence

Figures 3(a), 3(c) and 3(e) show the evolution of the characteristic thicknesses ( $h, h_o, h_i$ ) of the currents for different slope angles. The overall thickness $h$ (see figure 3 a) and outer-layer thickness $h_o$ (see figure 3 c) grow more rapidly on steeper slopes because of relatively vigorous entrainment. Notably, the overall current and the outer layer exhibit a similar normalised growth rate (approximately linear in time in the dynamically equilibrated regime), whilst the inner layer exhibits a much slower growth rate over time.

Figures 3(b), 3(d) and 3(f) plot the bulk Richardson numbers ( $\textit{Ri}, Ri_o, Ri_i$ ) as functions of time for different slope angles. In all cases, the various Richardson numbers become essentially constant in the dynamically equilibrated regime and are negatively correlated with slope angle, suggesting a greater role of stratification at a shallower angle.

Figure 4. Variation of averaged $(a)$ normalised $u_T$ , $(b)\, Ri,$ and $(c)$ normalised $e_T$ over $t_{ave}$ in the dynamically equilibrated regime against $\alpha$ in degrees, including the overall, outer and inner quantities for all the slope angles considered. The solid lines in $(a)$ and $(b)$ denote the theoretical predictions from (7.4) and (7.6), respectively. The solid lines in $(c)$ represent the prediction in (4.7) and (6.8).

Figure 4 shows the overall and layer-averaged normalised velocities, Richardson numbers and normalised TKE (the subscript $*$ denoting whether the ordinate in a given panel pertains to an overall or layer-specific quantity) for the currents as a function of slope angle in the dynamically equilibrated regime. The characteristic velocities $u_{T*}$ shown in figure 4 $(a)$ attain larger values at shallower angles ( $u_0$ is set to the same value across all the slope angles), reflecting the need for greater shear to overcome stronger stratification and transition the current to the turbulent state. We also observe that $u_T\approx u_{To}$ for all the slope angles.

Figure 4 $(b)$ shows that the Richardson numbers in the dynamically equilibrated state increase as the slope angle decreases; $\textit{Ri}_o$ appears to approach an asymptotic value at a small angle (note the logarithmic scale). This is consistent with the conjecture that stratified shear flows adjust to a ‘marginally stable’ state (Thorpe & Liu Reference Thorpe and Liu2009) characterised by a critical Richardson number. Nevertheless, further simulations at milder slopes are needed to confirm the asymptotic behaviour. The overall and inner Richardson numbers $\textit{Ri}$ and $\textit{Ri}_i$ , respectively, increase rapidly as the slope angle decreases (see § 6.1 and further discussion in § 7). The normalised TKE is seen in figure 4 $(c)$ to be approximately constant at the three smallest angles ( $1^\circ, 2^\circ, 5^\circ$ ) in each of the outer and inner layers.

Figure 5. Profiles of $(a)$ $\overline {u}/u_T$ , $(b)$ $\overline {b}/b_T$ , $(c)\, e/e_T$ and $(d) \overline {w^{\prime}u^{\prime}}/e_T$ against scaled height $z/h$ , and $(e)$ $\overline {u}/{u_T}_o$ , $(f)$ $\overline {b}/{b_T}_o$ , $(g)\, e/{e_T}_o$ and $(h)\, \overline {w^{\prime}u^{\prime}}/e_{To}$ against scaled distance from the velocity maximum $(z-z_{um})/h_o$ (outer-layer scaling). Profiles for each case at a series of times in the dynamically equilibrated regime are plotted. Strong self-similarity and collapse of profiles are observed in the outer layer for all the cases considered when normalisation is based on outer-layer integral quantities. The results for no-slip boundaries (from this study), and free-slip boundaries (adapted from van Reeuwijk et al. Reference van Reeuwijk, Holzner and Caulfield2019) are appended with ‘N’ and ‘F’, respectively, in the legend (e.g. 5F for a $5^\circ$ slope with free-slip boundaries). The solid lines in (e–h) represent the predictions from (4.1) and (4.5). The black dashed lines in (a,b) indicate the initial conditions used for all simulations presented in this study. The black dash-dotted lines show results from DNS of a spatially evolving current over slope $2.86 ^\circ$ , adapted from Salinas et al. (Reference Salinas, Balachandar and Cantero2021a ), which almost collapse onto the velocity profiles of case 2N in (a).

4.1. Self-similar profiles

Figure 5 shows normalised profiles of velocity $\overline {u}$ , buoyancy $\overline {b}$ , TKE $e$ and turbulent shear stress $\overline {w^{\prime}u^{\prime}}$ during the dynamically equilibrated regime for all the cases considered, including both no-slip (solid lines) and free-slip (dashed lines) boundary conditions. The results for the free-slip boundary conditions are from van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019), and the corresponding dataset is provided in van Reeuwijk (Reference van Reeuwijk2019). Each profile is scaled by the appropriate integral quantity ( $u_T$ , $b_T$ or $e_T$ ) at the time of sampling in the dynamically equilibrated regime to give figure 5(a–d). Note that $\overline {w^{\prime}u^{\prime}}$ is scaled with $e_T$ , and this is discussed in § 4.3. The normalised profiles for the flow variables largely collapse for all free-slip cases, whereas deviations become evident with no-slip boundary conditions.

The observed deviations suggest use of a local scaling based on the integral flow quantities in the outer layer (i.e. ${u_T}_o$ , ${b_T}_o$ and ${e_T}_o$ ). Figure 5(e–g) show the profiles of velocity, buoyancy, TKE and Reynolds stress rescaled with the appropriate outer-layer integral quantities. In addition, a normalised vertical (‘outer’) coordinate $(z-z_{um})/h_o$ is used to facilitate meaningful comparison between the results for free-slip and no-slip boundary conditions (i.e. the vertical coordinate has its origin at the level of the maximum velocity for both types of boundary conditions). Remarkably, all the normalised profiles nearly collapse in the outer layer whether or not an inner layer (corresponding to the region $(z-z_{um})/h_o \lt 0$ ) is present, indicating analogous self-similar dynamics dominate there. We will therefore refer to the dynamically equilibrated regime in the outer layer as the self-similar regime.

Although there are some deviations from universal forms, these look to be associated primarily with the presence of a strong stratification localised near the velocity maximum for currents on a low-angled no-slip boundary (see § 6.1 for more details). An important observation is that $ \overline {w^{\prime}u^{\prime}}$ is approximately zero at $z_{um}$ for all cases (see figure 5 h), as $z_{um}$ is defined via the velocity maximum, and the gradient diffusion hypothesis (see (4.4)) works reasonably well. This implies that $\partial \overline {u}/\partial z = 0$ and $\overline {w^{\prime}u^{\prime}}\approx 0$ hold at the bottom boundary ( $z=z_{um}$ ) of the outer layer in this study and of the currents in van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019), which underpins the consistency between them as shown in figure 5(e–h).

The DNS results of Salinas et al. (Reference Salinas, Balachandar and Cantero2021a ) for a spatially evolving current over a $2.86 ^\circ$ slope show consistent behaviour with the temporally evolving simulations presented in this study in terms of both bulk scaling in figure 5(a,b) and outer-layer scaling in figure 5(c,d). This comparison suggests that temporal currents (at long time) and spatial currents (sufficiently far from the head) share the same self-similar dynamics.

4.2. Approximate self-similar solutions

Given the observed near-collapse of the outer-layer profiles in § 4.1, we now examine the usefulness of the approximate self-similar descriptions developed by van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019) for inclined currents on a free-slip boundary. We thus propose that the outer-layer profiles are modelled as

(4.1a–c) \begin{align} \overline {u}=\underbrace {a_uB_o^{1/2}}_{u_{To}} \frac {2}{\eta _1^2}(\eta _1 - \eta _o), \quad \overline {b}=\underbrace {-a_b\frac {B_o}{h_o\sin \alpha }}_{b_{To}} \frac {2}{\eta _1^2}(\eta _1 - \eta _o), \quad e = \underbrace {a_eB_o}_{e_{To}}\frac {6 \eta _o}{\eta _1^3}(\eta _1 - \eta _o), \end{align}

where the outer-layer self-similarity variable is defined as $\eta _o = (z-z_{um})/h_o \in [0,\eta _1]$ , and the shape factor is $\eta _1 = 4/3$ . Note that (4.1a – c ) reduce to the form considered by van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019) for a current on a free-slip boundary upon setting $z_{um} = 0$ and dropping the subscript $o$ .

The coefficients $a_u$ , $a_b$ and $a_e$ depend on the dimensionless parameters of the problem, and need to be determined. We adopt the results of van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019), who observed that the eddy viscosity, eddy diffusivity, shear production $P_S$ , dissipation rate and turbulent Prandtl number could be parametrised as

(4.2a–e) \begin{align} K_m = c_m \frac {e}{S}, \quad K_\rho = c_\rho \frac {e}{S}, \quad {P_s = c_m e S},\quad \varepsilon = c_\varepsilon e S, \quad Pr_T=\frac {c_m}{c_\rho }, \end{align}

respectively, where $S=|\partial \overline {u}/\partial z|$ is the absolute strain rate of the mean flow, and $c_m=0.25\pm 5.5\times 10^{-4}$ , $c_\rho =0.31\pm 2\times 10^{-3}$ and $c_\varepsilon =0.21\pm 2.1\times 10^{-3}$ are empirical coefficients based on the DNS results. It follows from the success of these scalings in terms of the strain rate $S$ and TKE $e$ that the turbulence is in the shear-dominated regime (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014; Krug et al. Reference Krug, Holzner, Marusic and van Reeuwijk2017). Armed with this turbulence closure, van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019) integrated the equations for Reynolds-averaged momentum, buoyancy and TKE, (2.5)–(2.7), and used the Von Kármán–Pohlhausen method (Lighthill Reference Lighthill1950; Spalding Reference Spalding1954; Schlichting & Gersten Reference Schlichting and Gersten2016) to find the coefficients

(4.3a,b) \begin{align} a_u = \left (\frac {9}{8}\frac {Pr_T(c_m-c_\varepsilon )\tan \alpha }{\tan \alpha\, Pr_T + c_m} \right )^{-1/2}, \quad a_b=1. \end{align}

Note, however, that the coefficient for the TKE $a_e$ (see (4.1c )) did not follow from the analysis of van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019), and we evaluate it in the next subsection.

The theoretical solutions given by (4.1)–(4.3) are shown in figure 5(e–h), and are in good agreement with data from the simulations conducted in this study. Notably, the theory predicts that $e$ will tend to zero near $\eta _o = 0$ (i.e. $z=z_{um}$ ), but because this is not a solid boundary in these simulations, the TKE does not have to be zero there. Despite this, the shear production of TKE is zero at the velocity maximum by definition, and the magnitude of $e/e_{To}$ is indeed close to zero. Therefore, the theory still provides a reasonable approximation.

Figure 6(a–c) show the temporally averaged scaling coefficients defined in (4.2) over times sampled during the self-similar phase. We observe that there is a convincing collapse of profiles across the range of angles, with the converged values matching those from the free-slip cases. We thus conclude that the theory developed by van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019) can be effectively applied to the outer layer of an inclined current on a no-slip boundary.

Figure 6. Outer-layer scaling of averaged turbulence parameters over $t_{ave}$ : $(a)\, c_m=K_mS/e=0.25\pm 5.5\times 10^{-4}, (b)\, Pr_T=c_m/c_{\rho }=0.81\pm 5.2\times 10^{-3}$ (i.e. $c_{\rho } \approx 0.31$ ) and $(c)\, c_\varepsilon =\varepsilon /(eS)=0.21\pm 2.1\times 10^{-3}$ against scaled distance to velocity maximum $(z-z_{um})/h_o$ . The converged values are denoted with the vertical dashed lines.

4.3. Scaling of TKE

In this subsection, we explore if the TKE can be scaled with the integral buoyancy forcing $B$ as suggested by van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019), thereby allowing all the turbulence quantities in the closure to be related to macroscopic flow quantities. We first assume that the turbulent shear stress $\overline {w^{\prime}u^{\prime}}$ can be parametrised in the outer layer (where $\partial \overline {u} / \partial z \lt 0$ ) using the gradient diffusion hypothesis and (4.2):

(4.4) \begin{equation} \overline {w^{\prime}u^{\prime}}= -K_m \frac {\partial \overline {u}}{\partial z} = c_m \frac {e}{S} S = c_m e. \end{equation}

Substituting for $e$ using (4.1c ), we expect $\overline {w^{\prime}u^{\prime}}$ to take a self-similar form

(4.5) \begin{equation} \overline {w^{\prime}u^{\prime}}=c_m{e_T}_o\frac {6\eta _o}{\eta _1^3}\left (\eta _1-\eta _o\right ), \end{equation}

which is plotted in figure 5 $(h)$ and shows good agreement with the DNS data. Note that in contrast with $e$ , we observe that $\overline {w^{\prime}u^{\prime}}$ does become zero at $\eta _o=0$ , as discussed in § 4.1. Thus the quadratic profile (4.5) is more appropriate for $\overline {w^{\prime}u^{\prime}}$ than for $e$ .

In order to find $a_e (= e_{To}/B_o)$ , we first substitute the self-similar expressions (4.1a ) and (4.5) into the streamwise momentum (2.5) to give

(4.6) \begin{equation} \frac {\partial \overline {u}}{\partial t}=\underbrace {\frac {12c_m{e_T}_o-2B_o\eta _1}{{h_o}\eta _1^3}\eta _o + \frac {2B_o\eta _1-6c_m{e_T}_o}{(h_o)\eta _1^2}}_{-\partial \overline {w^{\prime}u^{\prime}}/\partial z-\overline {b}\sin \alpha }. \end{equation}

Secondly, a crucial simplification is motivated by the observation that the characteristic velocity in the outer layer, $u_{To}$ , becomes approximately constant (or only evolves over a relatively large time scale; see figure 2) in the self-similar regime, consistent with the theory described in van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019). Therefore, the maximum velocity in the outer layer $\overline {u}_m=3u_{To}/2$ (see (4.1a )) is also expected to become approximately constant in the self-similar regime, consistent with figure 7 $(a)$ .

Setting $\partial \overline {u}/\partial t \approx 0$ at $\eta _o=0$ in (4.6) gives $2B_o\eta _1-6c_m{e_T}_o=0$ , thus

(4.7) \begin{equation} a_e \equiv \frac {{e_T}_o}{B_o}= \frac {\eta _1}{3c_m} \approx 1.77 \end{equation}

and

(4.8) \begin{equation} \frac {\partial \overline {u}}{\partial t}=\frac {9B_o}{8h_o}\eta _o,\quad \eta _o \in [0, \eta _1]. \end{equation}

The prediction from (4.7) is shown in figure 4( $c$ ). Although the shear-dominated scaling mainly applies in the core region ( $\eta _o\in [0.5, 1]$ ) of the outer layer, there is fairly good agreement with the DNS data over the entire outer layer.

Equation (4.8) predicts that the acceleration ${\partial \overline {u}}/{\partial t}$ increases with height above the velocity maximum, and is consistent with the instantaneous velocity profiles in figure 7 $(c)$ . Analysis of the momentum budget shows that the individual terms in (4.6) also vary linearly with height above the velocity maximum in the outer layer (figure 7 b). Notably, the gradient of the Reynolds stress and the buoyancy terms are in approximate balance (i.e. no acceleration) at the velocity maximum and in the inner layer. This observation is discussed in detail in § 6.

Over long times, however, ${\partial \overline {u}}/{\partial t}$ in the outer layer gradually decreases as $h_o$ increases through entrainment, while $B_o$ remains constant, implying a slow approach towards zero acceleration. In reality, the current is likely disrupted by external effects such as tidal motions or rough topography, which are beyond the scope of the present study.

Figure 7. $(a)$ Temporal evolution of the maximum velocity $\overline {u}_m$ . $(b)$ Momentum budgets averaged in the dynamically equilibrated regime over $t_{ave}$ for cases 1N, 2N, 5N and 10N. $(c)$ The instantaneous along-slope velocity profiles in the dynamically equilibrated regime at $t/t^*=140, 152, 164$ for case 1N. The top legend applies to $(a)$ and $(b)$ .

5.1. Integral momentum and buoyancy budgets

The integral momentum and buoyancy equations for the inner and outer layers can be obtained by integrating (2.5) and (2.6) over the respective layer to give

(5.1a–d) \begin{align} \begin{aligned} \frac {{\textrm{d}} B_o}{{\textrm{d}} t} &= f_{um} \sin {\alpha },\quad \frac {{\textrm{d}} B_i}{{\textrm{d}} t} = -f_{um} \sin {\alpha }, \\ \frac {{\textrm{d}} Q_o}{{\textrm{d}} t} &= B_o + m_{um}, \quad \frac {{\textrm{d}} Q_i}{{\textrm{d}} t} = B_i -\tau _w -m_{um}, \end{aligned} \end{align}

where $\tau _{w} = \left .\nu \,({\partial \overline {u}/\partial z})\right \rvert _{0}$ is the shear stress at the lower boundary. The respective exchanges of buoyancy and momentum between the inner and outer layers are

(5.2a,b) \begin{align} \begin{aligned} f_{um} &= -\overline {w^{\prime}b^{\prime}}|_{z_{um}} +{\overline {b}|_{z_{um}}}\frac {{\textrm{d}} z_{um}}{{\textrm{d}} t}+\left . \kappa \frac {\partial \overline {b}}{\partial z}\right |_{z_{um}},\\ m_{um} &= \overline {w^{\prime}u^{\prime}}|_{z_{um}}-{\overline {u}|_{z_{um}}}\frac {{\textrm{d}} z_{um}}{{\textrm{d}} t}. \end{aligned} \end{align}

Here, $\overline {w^{\prime}b^{\prime}}|_{z_{um}}$ and $\overline {w^{\prime}u^{\prime}}|_{z_{um}}$ are the turbulent buoyancy and momentum fluxes at the level of the velocity maximum $z_{um}$ , respectively, $- \kappa ({\partial \overline {b}}/{\partial z})|_{z_{um}}$ is the molecular buoyancy flux at $z_{um}$ , and $\overline {b}|_{z_{um}}({{\textrm{d}} z_{um}}/{{\textrm{d}} t})$ and $\overline {u}|_{z_{um}} ({{\textrm{d}} z_{um}}/{{\textrm{d}} t})$ are the Leibniz terms (Schatzmann Reference Schatzmann1978; Davidson Reference Davidson1986; van Reeuwijk et al. Reference van Reeuwijk, Vassilicos and Craske2021) representing the effective buoyancy and momentum fluxes associated with a change in the height of the velocity maximum.

Figure 8( $a$ ) plots the terms in the integral buoyancy forcing budget of the outer layer ((5.1a ) and (5.2a ), scaled by ${B_o}/t_o$ ) as a function of slope angle, where $t_o= h_o/u_{To}$ is a typical turnover time scale of the outer layer. The magnitudes of the normalised fluxes are of order $10^{-2}$ , suggesting that the exchange of buoyancy happens over time scales much longer than $t_o$ . The flux with the largest magnitude is the Leibniz term (especially for small angles), but interestingly, the turbulent and molecular terms counteract it, creating a net buoyancy flux $f_{um}$ that is practically zero for all currents under consideration.

Similarly, figure 8( $c$ ) shows the buoyancy flux terms from the inner-layer budget (of equal magnitude and opposite sign to the outer-layer budget), but instead normalised by $B_i/t_i$ , where $t_i= h_i/u_{Ti}$ is a turnover time scale of the inner layer. Here, the scaled budget terms are also of order $10^{-2}$ , and as in the outer layer, the Leibniz term is in approximate balance with the turbulent and molecular terms. It is apparent that the integral buoyancy forcing $B_*$ can be approximated as constant in each of the inner and outer layers on time scales up to at least $t_i$ and $t_o$ , respectively.

Given that the evolution of the integral buoyancy forcing ( ${\textrm{d}} B_*/{\textrm{d}} t$ ) is not a leading-order term in the inner- and outer-layer budgets at any slope angle considered, the dynamics governing the interface can therefore be regarded as quasi-steady in the self-similar regime, and (5.1a,b ) can be approximated as

(5.3) \begin{equation} \frac {{\textrm{d}} B_o}{{\textrm{d}} t}\approx \frac {{\textrm{d}} B_i}{{\textrm{d}} t}\approx 0, \end{equation}

which is also consistent with the results in figure 2 $(j)$ .

Figure 8. Normalised terms against slope angles for the budgets of $(a)$ outer integral buoyancy forcing, $(b)$ outer volume flux , $(c)$ inner integral buoyancy forcing, and $(d)$ inner volume flux. Note that these term are averaged over $t_{ave}$ in the dynamically equilibrated regime.

Equations (5.1c , d ) and (5.2b ) present the integral momentum (volume flux) budgets, where $B_*$ acts to accelerate the flow and increase the volume flux. Figure 8( $b$ ) illustrates the individual terms in the integral momentum budget of the outer layer, normalised by $B_o$ . At all angles, it is clear that the Leibniz term makes the dominant contribution to the momentum exchange between layers $m_{um}$ (see (5.2b )), while the Reynolds stress $\overline {w^{\prime}u^{\prime}}|_{z_{um}}$ plays a negligible role. However, this momentum exchange has magnitude approximately $0.15B_o$ at all angles considered, and therefore plays a minor role in modifying the rate of volume flux increase, i.e.

(5.4) \begin{equation} \frac {{\textrm{d}} Q_o}{{\textrm{d}} t}\approx B_o. \end{equation}

Figure 8 $(d)$ shows the integral momentum budget for the inner layer, and we observe a leading-order balance between the buoyancy forcing and bottom shear stress, i.e.

(5.5) \begin{equation} B_i\approx \tau _w. \end{equation}

The rate of increase of volume flux in the inner layer is of a similar order to the Leibniz term (towards which the Reynolds stress contribution is negligible, as in the outer layer).

5.2. Layer-specific force balance

Ellison & Turner (Reference Ellison and Turner1959) showed that gravity currents reach an equilibrium state where the gravitational forces are balanced by ‘entrainment drag’ and bottom friction. It is useful to interpret this finding in the light of weak interaction between the inner and outer layers.

The starting point is the integral momentum balance for the entire layer, which can be obtained by adding (5.1c ) and (5.1d ) to give the time derivative of $Q(=u_Th)$ :

(5.6) \begin{equation} \frac {{\textrm{d}} Q}{{\textrm{d}} t}=h\frac {{\textrm{d}} u_T}{{\textrm{d}} t}+u_T\frac {{\textrm{d}} h}{{\textrm{d}} t}= (B_o + B_i)-\tau _w. \end{equation}

In terms of top-hat variables (Ellison & Turner Reference Ellison and Turner1959), this equation can be written as

(5.7) \begin{equation} h\frac {{\textrm{d}} u_T}{{\textrm{d}} t}=(B_o + B_i)-Eu_T^2-\tau _w, \end{equation}

where $E=u_T^{-1} {\textrm{d}} h/{\textrm{d}} t$ is the entrainment coefficient of a temporal gravity current (van Reeuwijk et al. Reference van Reeuwijk, Krug and Holzner2018, Reference van Reeuwijk, Holzner and Caulfield2019). Since $u_T$ is expected to be constant in the dynamically equilibrated regime (Ellison & Turner Reference Ellison and Turner1959) (see also figure 2), (5.7) simplifies to

(5.8) \begin{equation} B_o+B_i=\tau _w+Eu_T^2. \end{equation}

Given (5.5) and that $m_{um}$ is of a similar magnitude to both $(B_i - \tau _w)$ and $({\textrm{d}} Q_o/{\textrm{d}} t - B_o)$ , as observed in figures 8 $(b)$ and 8 $(d)$ , we deduce that

(5.9) \begin{equation} B_o\approx Eu_T^2. \end{equation}

These results support distinct dynamics in the inner and outer layers. Buoyancy in the inner layer primarily overcomes the bottom friction, as shown in figure 7 $(b)$ , whilst the buoyancy in the outer layer overcomes drag associated with entrainment of ambient fluid. In the absence of significant exchange of momentum, the inner and outer layers are only weakly coupled (subject to the continuity condition at the velocity maximum).

Importantly, our theoretical parametrisation and DNS results show that this weak coupling is not confined to currents on small-angle slopes (with relatively strong stabilising stratification), as conjectured by Salinas et al. (Reference Salinas, Balachandar, Shringarpure, Fedele, Hoyal, Zuñiga and Cantero2021b ), but also applies on larger-angle slopes (where the stabilising stratification is relatively weak). This observation signifies that the weak coupling is not a result of a density interface forming near the velocity maximum. Instead, it appears to be a natural behaviour of inclined gravity currents.

Although it has been reported that there is a mismatch between the levels of the velocity maximum and zero turbulent shear stress (Salinas et al. Reference Salinas, Balachandar, Shringarpure, Fedele, Hoyal, Zuñiga and Cantero2021b ; Wei et al. Reference Wei, Wang and Yang2021), our results indicate that this is not a leading-order effect, and the gradient–diffusion hypothesis remains useful, i.e. both the viscous shear stress and turbulent shear stress change sign and cross zero at the level of the velocity maximum, across which there is essentially no momentum exchange by turbulence or diffusion. This ‘decoupling’ of the two layers explains why the outer layer behaves independently of bottom friction and much like a current on a free-slip slope described by van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019).

6.1. Inner-layer profiles

Figure 9(a–d) show the inner-layer profiles of $\overline {u}, \overline {b}, e, \overline {w^{\prime}u^{\prime}}$ normalised by the characteristic scales ${u_T}_i$ , ${b_T}_i$ and ${e_T}_i$ , respectively, at a series of times in the dynamically equilibrated regime. We employ $z/h_i$ as the scaled slope-normal coordinate.

We decompose the inner shear layer into two regions: a turbulent wall region (TWR) and a viscous wall region (VWR) defined in terms of the dimensionless wall distance $z^+ = zu_\tau /\nu$ . Figure 9 $(i)$ illustrates these regions using case 2N. The VWR is defined up to $z^+ = 10$ , rather than $z^+ = 50$ as in traditional boundary layers (Pope Reference Pope2000), since the viscous contribution to total shear stress becomes negligible beyond $z^+ = 10$ for small-angled gravity currents. The TWR lies above the VWR, with TKE peaking at its lower boundary, and decreasing to a minimum at its upper boundary. The upper boundary of the TWR is defined via a density interface (see figure 9 b) in the vicinity of the velocity maximum for small-angled cases. The distinct regions for case 2N are also illustrated in figure 9(a–h) with the same colour scheme as in figure 9 $(i)$ . Note that the VWR occupies only a minor area at the bottom of the inner shear layer in figure 9(a–d), as the vertical coordinate is scaled as $z/h_i$ .

Figure 9. Profiles of spatially averaged $(a)$ $\overline {u}/{u_T}_i$ , $(b)$ $\overline {b}/{b_T}_i$ , $(c)$ $e/{e_T}_i$ and $(d)$ $\overline {w^{\prime}u^{\prime}}/{e_T}_i$ against scaled height $z/h_i$ , and $(e)$ $\overline {u}/{u_\tau }$ , $(f)$ $\overline {b}/\overline {b}_w$ , $(g)$ $e/u_\tau ^2$ and $(h)$ $\overline {w^{\prime}u^{\prime}}/u_\tau ^2$ against $z^+ = zu_\tau /\nu$ . Here, $u_\tau$ is the friction velocity, and $\overline {b}_w$ is the spatially averaged buoyancy at the wall. Note that the profiles at a series of times in the dynamically equilibrated regime are plotted for each case. Here, CCF indicates the closed channel flow data adapted from Lee & Moser (Reference Lee and Moser2015), and OCF indicates the open channel flow data adapted from Yao et al. (Reference Yao, Chen and Hussain2022), both at ${Re}_\tau =550$ . Distinct regions are highlighted with shading, which are depicted in $(i)$ using case 2N with ${Re}_\tau =620$ , including the VWR, the TWR and the density interface, represented by blue, yellow and green, respectively.

The profiles of scaled mean velocity in figure 9 $(a)$ exhibit considerable collapse over a range of times. However, this collapse is not indicative of a self-similar regime like the outer layer; instead, the inner layer reaches a quasi-steady state that evolves over very long time scales. Indeed, figure 7 $(c)$ depicts (for slope angle $1^\circ$ ) profiles of mean streamwise velocity $\overline {u}$ scaled by the initial (constant) velocity $u_0$ in the inner shear layer at different times. The profiles remain nearly unchanged with time because the buoyancy and shear stress gradient are in local balance and $\partial \overline {u}/\partial t \approx 0$ throughout the inner layer, as shown in figure 7 $(b)$ , i.e.

(6.1) \begin{equation} \frac{\partial}{\partial z} \left ( \overline {w^{\prime}u^{\prime}}-\nu \frac {\partial \overline {u}}{\partial z} \right ) \approx -\overline {b}\sin \alpha . \end{equation}

Furthermore, the profiles of buoyancy in the inner layer (figure 9 b) are practically uniform due to strong turbulent mixing, especially for small angles, implying that

(6.2) \begin{equation} \overline {w^{\prime}u^{\prime}}-\nu \frac {\partial \overline {u}}{\partial z} \approx -\overline {b}z - \tau _w \approx B_i(z/h_i-1). \end{equation}

The approximately linear profile of total momentum flux suggested by (6.2) shares a clear analogy with a canonical plane turbulent channel flow subject to constant streamwise pressure gradient, where the momentum balance is given by $\partial(\overline{w'u'} - \nu \partial \overline{u}/\partial z)/ \partial z = \partial \overline{p}/ \partial x$

In order to explore this connection, profiles corresponding to a closed channel flow (CCF, adapted from Lee & Moser Reference Lee and Moser2015) and an open channel flow (OCF, adapted from Yao et al. Reference Yao, Chen and Hussain2022) with turbulent Reynolds number ${Re}_\tau =550$ have been included in figure 9. Note that CCF is subject to no-slip conditions on both the bottom and top boundaries of the channel, with symmetry about the plane at half-height. For the purposes of comparison, the channel half-height is plotted as corresponding to $z/h_i =1$ . In contrast, OCF has a shear-free surface at the top, which is taken to correspond to $z/h_i =1$ . The values of ${Re}_\tau$ that characterise the inner shear layer in the gravity currents are a function of slope angle (see table 1 for more details); however, case 2N corresponds to ${Re}_\tau =620$ , similar to that of the selected channel flow comparisons. The scaled velocity profiles of the channel flows shown in figure 9 $(a)$ nearly overlap with those for the inner layer in the gravity currents, suggesting the strong similarity of these two flow types. Notably, this similarity is also observed in the spatially evolving DNS results of Salinas et al. (Reference Salinas, Balachandar and Cantero2021a ).

Figure 9 $(c)$ shows the TKE profiles normalised by $e_{Ti}$ . For small-slope-angle currents, these profiles nearly collapse and decrease approximately linearly with height within the TWR. Deviations become apparent in the TKE profiles between gravity currents on small- and large-angled slopes, the likely explanation for which is that the density interface near the velocity maximum strengthens as the slope angle decreases (figure 9 b), acting to suppress turbulence. This explanation is supported by the near collapse of the normalised TKE profiles for the small-angle cases (1N, 2N, 5N) with that for OCF, in which no vertical transport of turbulence is possible at the free surface. The magnitude of normalised TKE for CCF is slightly smaller than for OCF in the TWR, an effect attributed to stronger ‘very-large-scale motions’ in OCF (Kim & Adrian Reference Kim and Adrian1999; Balakumar & Adrian Reference Balakumar and Adrian2007; Yao et al. Reference Yao, Chen and Hussain2022).

Figure 9 $(d)$ shows the normalised turbulent shear stress $\overline {w^{\prime}u^{\prime}}/e_{Ti}$ , which varies linearly with height throughout most of the inner layer for small-angled gravity currents and all channel flows. This arises because viscous shear stress is negligible compared to turbulent shear stress in the TWR, allowing (6.2) to be simplified as

(6.3) \begin{equation} \overline {w^{\prime}u^{\prime}} \approx B_i(z/h_i-1). \end{equation}

The observed collapse of these profiles upon scaling with $e_{Ti}$ is discussed in § 6.2.

Figure 9(e–h) show the scaled quantities in wall units with a focus on the VWR, which is highlighted using the same colour scheme as in figure 9 $(i)$ . The velocity profiles in figure 9 $(e)$ fully collapse in the VWR ( $z^+ \lt 10$ ) across all cases, conforming to the well-known relationship $\overline {u}/u_{\tau } = z^+$ in this viscosity-dominated region. Figure 9 $(f)$ shows that the buoyancy throughout the VWR is close to the wall buoyancy $\overline {b}_w$ .

The profiles of TKE and $\overline {w^{\prime}u^{\prime}}$ normalised by $u_\tau ^2$ are presented in figures 9 $(g)$ and 9 $(h)$ , respectively, in terms of wall units. The normalised TKE profiles are seen to collapse onto a single curve within the VWR that increases rapidly with height from zero on the lower boundary to a local maximum at the top of the VWR. The normalised $\overline {w^{\prime}u^{\prime}}$ profiles for the 2N slope current and the CCF and OCF cases, which all have a similar ${Re}_\tau$ , also nearly collapse in the VWR. However, the deviations from this normalised profile for the other slope currents considered suggest a possible ${Re}_\tau$ dependence in this scaling.

6.2. Approximate steady-state solution

Upon inspection, we find that an approach analogous to that for the outer layer in §§ 4.2 and 4.3 can also be applied to the inner shear layer. On the basis of the observed collapse in figure 9, approximate solutions are proposed, especially for currents on small-angled slopes:

(6.4) \begin{align} \overline {u}={u_T}_i {f_u}_i(\zeta ), \quad \overline {b}={b_T}_i {{f_b}_i}(\zeta ),\quad e={e_T}_i {{f_e}_i}(\zeta ), \quad \zeta =z/h_i \in (0, 1). \end{align}

Note that $f_{ui}(\zeta )$ , $f_{bi}(\zeta )$ and $f_{ei}(\zeta )$ are expected to be strictly valid only within the TWR. However, as the VWR volume is negligible compared with that of the TWR for the small-angled slope currents, we assume that these functions may be applied throughout the inner shear layer. Taking $f_{bi}(\zeta )\approx 1$ (consistent with figure 9 b), we assume that ${f_u}_i$ takes a logarithmic form, as for an unstratified boundary layer adjacent to a no-slip surface:

(6.5) \begin{equation} {f_u}_i(\zeta )= c_i\ln {\zeta }+c_{um}, \end{equation}

where $c_i$ is a coefficient to be determined, and $c_{um}={f_u}_i(1)=\overline {u}_m/{u_T}_i$ . Consistency with the volume transport decomposition for the inner layer in (2.9a ) requires

(6.6) \begin{equation} Q_i=\int _0^{h_i}\overline {u}\,{\textrm{d}} z = {u_T}_i h_i\int _{0}^1{f_u}_i(\zeta )\,{\textrm{d}}\zeta ={u_T}_ih_i \quad \Rightarrow \quad \int _{0}^1{f_u}_i(\zeta )\,{\textrm{d}}\zeta =1, \end{equation}

thus $c_i=c_{um}-1$ .

As in § 4.3, we propose that $e$ and $\overline {w^{\prime}u^{\prime}}$ in the inner layer (where $\partial \overline {u}/\partial z\gt0$ now) are related by the scaling

(6.7) \begin{equation} e=K_{mi}S/c_{mi}=-\overline {w^{\prime}u^{\prime}}/c_{mi}. \end{equation}

Note that an additional subscript $i$ is used to distinguish the inner-shear-layer eddy parametrisation coefficient ( $c_{mi}$ and, later, $c_{\rho i}$ and $c_{\epsilon i}$ ) from that applicable to the outer shear layer. Combining (6.3), (6.4) and (6.7) suggests that $f_{ei}$ takes a linear form (consistent with collapse of the small-angled cases in figure 9 c), and consistency with the TKE decomposition for the inner layer (2.11) requires that $\int _{0}^{1}f_{ei}(\zeta )\,{\textrm{d}}\zeta =1$ , thus

(6.8) \begin{equation} f_{ei}=2(1-\zeta ), \quad e_{Ti}=e/f_{ei}=B_i/(2c_{mi}). \end{equation}

The approximate solutions proposed in (6.4) and (6.7) can thus be summarised as

(6.9a-c) \begin{align} \frac {\overline {u}}{{u_T}_i}= \underbrace {(c_{um}-1)\ln {\zeta }+c_{um}}_{{f_u}_i(\zeta )}, \quad \frac {\overline {b}}{b_{Ti}}=\underbrace {1}_{f_{bi}(\zeta )},\quad \frac {e}{e_{Ti}} = -c_{mi}\frac {\overline {w^{\prime}u^{\prime}}}{e_{Ti}}=\underbrace {2(1-\zeta )}_{f_{ei}(\zeta )}. \end{align}

The ratio of $u_m$ to $u_{Ti}$ (i.e. $c_{um}$ ) is found from the DNS data to be approximately 1.12. It is interesting to consider the analogy with channel flow using the well-known approximations for the inner layer $\overline {u}/u_\tau =c_\kappa ^{-1} \ln z^+ +C$ and $(\overline {u}_m - u_{Ti})/u_\tau = c_k^{-1}$ (Pope Reference Pope2000, equations (7.43) and (7.52) therein), where $c_\kappa =0.41$ is the von Kármán constant, and $C$ is a constant with value approximately 5.2. Comparison with (6.9a ) suggests that

(6.10) \begin{equation} c_{um} = \frac {c_{\kappa }C + \ln {Re}_{\tau }}{c_{\kappa }C + \ln {Re}_{\tau }-1}, \quad \frac {u_{Ti}}{u_\tau } =\frac {1}{(c_{um}-1)c_\kappa }, \end{equation}

which is in excellent agreement with the DNS data, e.g. yielding $c_{um} = 1.13$ and $u_{Ti}/u_{\tau } = 18.8$ for ${Re}_{\tau } = 620$ .

Figures 10(a), 10(b) and 10(c) show the coefficients $c_{mi}$ , $Pr_{Ti}$ and $c_{\varepsilon i}$ based on the DNS data and the form of the outer-layer scaling in (4.2). We observe that a single value for each scaling coefficient can be applied with reasonable success to the core region of the TWR (where $z/ h_i \in [0.25, 0.75]$ ) in the small-angle slope currents (1N, 2N, 5N) and the channel flows. In particular, $c_{mi}$ converges to approximately 0.27, which is similar to the value of $c_m$ (namely 0.25), indicating comparable shear-driven dynamics in the core regions of the inner and outer layers. Likewise, $c_{\varepsilon i}$ also converges to approximately 0.27, suggesting an approximate balance between turbulence production and dissipation in the quasi-steady state (noting the relation in (4.2c , d )). The converged value of $Pr_{Ti}$ is close to unity, in agreement with the Reynolds analogy. However, $Pr_{Ti}$ exhibits opposite trends in small- and large-angle cases near the velocity maximum, arising from variations in the relative magnitudes of turbulent viscosity and diffusivity. This behaviour presumably suggests the existence of a critical slope angle below which the strong density interface that inhibits mixing spontaneously forms near the velocity maximum. The mechanisms underlying this critical threshold warrant further investigation.

The functions given by (6.9) with $c_{mi} = 0.27$ and $c_{um} = 1.12$ are plotted in figure 9(a–d), and show reasonable agreement with the DNS data for currents on small-angle slopes (e.g. up to $5^{\circ }$ ) within the TWR. The corresponding prediction from (6.8) that $e_{Ti}/B_i\approx 1.9$ is plotted in figure 4 $(c)$ and agrees well with the DNS results for small-angled currents. We suggest that the likely reason for the success of this scaling is that the core region of the TWR is sufficiently far from both the wall and the density interface for the turbulence to be shear-dominated (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014), and therefore that parametrising the turbulence using $|\partial \overline {u}/\partial z|$ and $e$ is reasonable.

Figure 10. Scaling of inner-layer turbulence parameters averaged over $t_{ave}$ in the dynamically equilibrated regime: $(a)$ $ c_{mi}=K_{mi}S/e= 0.27\pm 0.018$ , $(b)$ $ Pr_{Ti}=c_{mi}/c_{\rho i}=1\pm 0.048$ and $(c)$ $ c_{\varepsilon i}=\varepsilon /(eS)=0.27\pm 0.012$ against $z/h_i$ . The converged values are denoted with vertical dashed lines.

7.1. Buoyancy partition

It is clear that the buoyancies in the inner and outer layers ( $B_i$ and $B_o$ ) are the primary forcings in each layer. Therefore, a crucial step in obtaining a closed-form description of a slope current is to predict how buoyancy is partitioned between the layers. Given that $u_T\approx u_{To}$ (figure 2 g) and $\tau _w=u_\tau ^2$ , (5.5) and (5.9) can be rewritten as

(7.1a,b) \begin{align} B_o=Eu_{To}^2,\quad B_i=u_\tau ^2. \end{align}

Combining this with a velocity-matching condition in terms of the self-similar relations $\overline {u}_m=3u_{To}/2=c_{um}u_{Ti}$ gives

(7.2a–c) \begin{align} B_o=\frac {4}{9}E\overline {u}_m^2 ,\quad B_i=c_{mk}^2\overline {u}_m^2, \quad \textrm {and thus}\quad \frac {B_i}{B_o}=\frac {9}{4}c_{mk}^2E^{-1}, \end{align}

where $c_{mk}={(c_{um}-1)c_{\kappa }}/{c_{um}}$ . Recalling from (4.1a ) that $E=B_o/u_{To}^2=a_u^{-2}$ , where $a_u$ is given by (4.3a ), the predicted dependence of $B_i/B_o$ on $\alpha$ is

(7.3) \begin{equation} \frac {B_i}{B_o}=\frac {2c_{mk}^2(Pr_T+{c_m}/\tan \alpha )}{Pr_T\,(c_m-c_\varepsilon ) }, \end{equation}

which is compared in figure 11 $(a)$ with $B_i/B_o$ calculated from the DNS data. Reasonably good agreement is found across the range of slope angles $\alpha$ considered, and the ratio $B_i/B_o$ is seen to increase as both the slope angle $\alpha$ and the associated entrainment rate decrease (7.2c ). This insight could explain the long runout of submarine gravity currents over mild slopes. As the slope angle reduces: (i) a greater proportion of the buoyancy is confined in the inner layer, where it propagates largely undiluted because of a weak interaction with the outer layer; and (ii) the outer-layer buoyancy also experiences limited dilution due to the reduced entrainment rate.

Figure 11. Plots of $(a)\, B_i/B_o$ , $(b)\, {\overline {u}_m}/u_0$ and $(c)\, h_i/h_0$ averaged over $t_{ave}$ in the dynamically equilibrated regime as a function of slope angle $\alpha$ (in degrees). The solid lines in (a,b) denote the theoretical prediction in (7.3), (7.4), respectively.

The maximum velocity $\overline {u}_m$ is determined by substituting (7.3) and $B=B_i+B_o$ into (7.2) to give

(7.4) \begin{equation} \overline {u}_m(\alpha )=\frac {\sqrt {B_i}}{c_{mk}}=\sqrt {\frac {2B( Pr_T+c_m/\tan \alpha )}{2c_{mk}^2(Pr_T+c_m/\tan \alpha )+Pr_T\,(c_m-c_\varepsilon )}}=\frac {3}{2}u_{To}=c_{um}u_{Ti}. \end{equation}

The DNS data for $\overline {u}_m$ and $u_{T*}$ shown in figures 11 $(b)$ and 4 $(a)$ are seen to be well predicted by (7.4). Notably, the characteristic thickness $h_i$ remains a free parameter. Applying buoyancy conservation yields

(7.5) \begin{equation} \frac {\int _0^{z_{um}}\overline {b}(t,z)\,{\textrm{d}} z}{\int _0^\infty \overline {b}(t=0,z)\,{\textrm{d}} z}=\frac {b_{Ti}h_i}{b_0h_0}=\frac {B_i}{B_o+B_i}\quad \Rightarrow \quad \frac {h_i}{h_0} = \frac {B_i/B_o}{1+B_i/B_o}\frac {b_0}{b_{Ti}}. \end{equation}

Here, $B_i/B_o$ is a function of $\alpha$ as given in (7.3). Equation (7.5) indicates $h_i/h_0$ ( $h_0$ set to constant) depends on $\alpha$ and the ratio of the initial buoyancy $b_0$ to the inner characteristic buoyancy $b_{Ti}$ . As shown in figure 11 $(c)$ , $h_i/h_0$ attains larger values at smaller angles, with a sharp increase observed between cases 2N and 5N. This sudden rise is likely driven by the dilution due to the initial burst (see $t/t^*\in [20, 40]$ in figure 2), which substantially increases the ratio $b_0/b_{Ti}$ . This observation somewhat suggests that the history of a flow influences its subsequent evolution, as highlighted by Caulfield (Reference Caulfield2021).

7.2. Entrainment law

Using the self-similar relations for $\overline {u}_m$ , $u_{To}$ and $u_{Ti}$ in § 7.1, we can obtain expressions for $\textit{Ri}_*$ upon combining (2.10d ), (7.3), (7.4) and $u_T \approx u_{To}$ :

(7.6a–c) \begin{eqnarray} && Ri_o=\frac {9}{8}\frac {Pr_T\,(c_m-c_\varepsilon )}{\tan \alpha\, Pr_T +c_m},\quad Ri_i=\frac {c_{mk}^2c_{um}^2}{\tan \alpha },\quad \nonumber\\&& Ri=\frac {9}{8}\frac {2c_{mk}^2(Pr_T+c_m/\tan \alpha )+Pr_T\,(c_m-c_\varepsilon )}{\tan \alpha\, Pr_T+c_m}. \end{eqnarray}

The predictions for these Richardson numbers are plotted in figure 4 $(b)$ , with good agreement apparent for small slope angles. Equation (7.6a ) indicates that the outer-layer Richardson number $\textit{Ri}_o$ increases and approaches a finite limit ( $\textit{Ri}_{om} \approx 9\,{Pr}_T\, (c_m-c_\varepsilon )/8c_m$ ) as $\alpha$ decreases ( $\tan \alpha \ll c_m/{Pr}_T$ ). This limiting value suggests that the outer layer remains marginally stable and weakly stratified at all the slope angles considered here. However, $\textit{Ri}$ and $\textit{Ri}_i$ continue to increase as the slope angle decreases. Notably, $\textit{Ri}_i$ loses physical relevance for small-angled cases, as the inner layer is nearly well-mixed in these scenarios. Richardson number $\textit{Ri}$ is directly related to the buoyancy partition, as shown by the following relationship derived from (2.10d ) with $u_T \approx u_{To}$ :

(7.7) \begin{equation} Ri/Ri_o = B/B_o. \end{equation}

Therefore, the increase in $\textit{Ri}$ with decreasing $\alpha$ essentially reflects the increasing proportion of the integral buoyancy held in the inner layer.

Figure 12. Entrainment rate $E$ against $\textit{Fr}=1/\sqrt {Ri}$ . The black solid line and triangles denote the prediction in (7.9) and the present DNS data, respectively. The shaded band represents the 95 % uncertainty interval of the present theory obtained via Monte Carlo uncertainty analysis. The data and theory from van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019), Salinas et al. (Reference Salinas, Zúñiga, Cantero, Shringarpure, Fedele, Hoyal and Balachandar2022), and the previous data and fitted functions compiled by Odier et al. (Reference Odier, Chen and Ecke2014) and Salinas et al. (Reference Salinas, Cantero, Shringarpure and Balachandar2019) are also shown – incorporating the studies by Ellison & Turner (Reference Ellison and Turner1959), Ashida & Egashira (Reference Ashida and Egashira1975), Parker et al. (Reference Parker, Garcia, Fukushima and Yu1987), Cenedese et al. (Reference Cenedese, Whitehead, Ascarelli and Ohiwa2004), Wells (Reference Wells2007), Odier et al. (Reference Odier, Chen and Ecke2014) and Salinas et al. (Reference Salinas, Cantero, Shringarpure and Balachandar2019), together with field data (filled markers) collected from the Mediterranean, Lake Ogawara and the Faroe Bank Channel). The dataset for this plot is available at Cui (Reference Cui2025).

As we have shown in this study that entrainment is associated with the outer layer dynamics in a slope current, we now adapt the entrainment law derived by van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019) for currents on free-slip boundaries ((4.24) and (4.25) therein). Using (7.7) to relate $\textit{Ri}$ and $\textit{Ri}_o$ , their theoretical expression can be rewritten as

(7.8) \begin{equation} E=\frac {c_m}{Pr_T}(Ri_{om}-Ri_o)=\frac {c_m}{Pr_T}\left (Ri_{om}-Ri\,\frac {B_o}{B}\right ),\quad Ri_{om}=0.15, \end{equation}

which, with (7.2c ), gives the entrainment law as an explicit function of $\textit{Ri}$ :

(7.9) \begin{eqnarray} && E=\frac {c_m}{Pr_T}\left (Ri_{om}-\frac {E\,Ri}{9c_{mk}^2/4+E}\right ) \quad\Rightarrow \quad \nonumber\\ && E=\frac {\sqrt {\left ({c_m\,Ri/Pr_T}+c_{R1}\right )^2+c_{R2}}-c_m\,Ri/Pr_T-c_{R1} }{2}, \end{eqnarray}

where $c_{R1}={9}c_{mk}^2/{4}-{c_m}\,Ri_{om}/Pr_T$ and $ c_{R2}={9}c_{mk}^2{c_m}\,Ri_{om}/Pr_T$ are constants.

Figure 12 shows the predicted entrainment rate $E$ as a function of densitometric Froude number $\textit{Fr}=1/\sqrt {Ri}$ . The present theory (7.9) (solid black line) is shown together with entrainment models from van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019) as a dashed line, Ellison & Turner (Reference Ellison and Turner1959) as a dash-dotted line, and Parker et al. (Reference Parker, Garcia, Fukushima and Yu1987) as a dotted line, along with a broad dataset (denoted by coloured symbols) from laboratory experiments, DNS and field observations, as compiled by Odier et al. (Reference Odier, Chen and Ecke2014) and Salinas et al. (Reference Salinas, Cantero, Shringarpure and Balachandar2019) (details provided in the caption). The present data are in good agreement with the present theory, and both show consistency with the spatially evolving DNS results of Salinas et al. (Reference Salinas, Zúñiga, Cantero, Shringarpure, Fedele, Hoyal and Balachandar2022), suggesting that similar entrainment processes operate in both temporal and spatial configurations.

The classical parametrisation based on the experimental data from Ellison & Turner (Reference Ellison and Turner1959) (proposed by Turner Reference Turner1986) aligns well with the high- $\textit{Fr}$ data, but shows a different asymptotic behaviour at low $\textit{Fr}$ (with $E$ dropping to 0 at a critical $\textit{Fr}$ value beyond the range accessible to their experiments). The parametrisation fitted by Parker et al. (Reference Parker, Garcia, Fukushima and Yu1987) has asymptotic behaviour that is more consistent with the data and offers better overall performance, but its functional form lacks a solid theoretical basis.

The theoretical model of van Reeuwijk et al. (Reference van Reeuwijk, Holzner and Caulfield2019) predicts $E$ as a function of the outer-layer Richardson number $\textit{Ri}_o$ , as described by (7.8). A comparison between their prediction $E(Ri_o)$ and the present theoretical prediction $E(Ri)$ (see (7.9)), presented in figure 12, highlights the significant role that the outer-layer dynamics is likely to play in many slope current applications. As indicated by (7.7), $\textit{Ri}\approx Ri_o$ at large $\textit{Fr}$ (corresponding to steep slopes), where the integral buoyancy forcing is primarily confined within the outer layer (i.e. $B_o/B \approx 1$ ). Under these conditions, the two predictions are in close agreement, while increasing deviations are observed as $\textit{Fr}$ decreases (note, however, that the entrainment rate remains consistent for the same slope angle). Both theories suggest a maximum entrainment rate $E_m = c_m\,Ri_{om}/Pr_T\approx 0.046$ as $\textit{Fr}$ approaches infinity, which is in good agreement with the value 0.04 proposed by Wells, Cenedese & Caulfield (Reference Wells, Cenedese and Caulfield2010).

The present theory indicates that entrainment is not completely suppressed at a critical $\textit{Ri}$ ( $\textit{Fr}$ ), but rather asymptotes towards 0 as $\textit{Ri}\rightarrow \infty\ (\alpha \rightarrow 0)$ , consistent with the level of TKE in the outer layer $e_{To}$ that scales with $B_o$ (see (4.7)). Notably, this differs conceptually from the hypothesis of ‘continued (high-Richardson-number) mixing’ associated with intermittent turbulence under strong stratification (see e.g. Wells et al. Reference Wells, Cenedese and Caulfield2010). Despite $\textit{Ri}$ approaching infinity as $\alpha$ decreases towards 0, the outer layer herein remains weakly stratified, with $\textit{Ri}_o$ asymptotically approaching $\textit{Ri}_{om}$ as shown in figure 4 $(b)$ . Crucially, the present theory shows good agreement with the field data (filled symbols in figure 12), offering a physical basis and the prospect of general applicability to flows at small $\textit{Fr}$ of geophysical relevance.