2 Shallow-layer formulation

We analyse the motion of a sustained current of relatively dense fluid through a less dense environment within a horizontal channel of height $H$ . The motion is driven from a source that delivers fluid of density ${\it\rho}_{1}$ at a volume flux per unit width $Q$ , while the environmental fluid is of density ${\it\rho}_{2}$ . The thickness of the flowing dense layer is denoted by $h_{1}$ and its depth-averaged velocity is given by $u_{1}$ . The thickness of the overlying layer of less dense fluid, flowing above the dense layer is $h_{2}$ and its depth-averaged velocity is denoted by $u_{2}$ (see figure 1). Conservation of mass in each layer is then given by

where the coordinate axes are aligned so that $x$ is horizontal and $z$ is vertical. In this formulation, we have neglected mixing processes between the two layers which would generate a zone of intermediate density and which could become dynamically important sufficiently far downstream from the source when the mixed layer has grown sufficiently thick. The motion is assumed to evolve so that horizontal length scales far exceed vertical length scales and thus the pressure is hydrostatic to leading order. Denoting the gravitational acceleration by $g$ , we find that

where $p_{T}(x,t)$ is the as yet unknown pressure at $z=H$ . Then following Long (Reference Long1955), Rottman & Simpson (Reference Rottman and Simpson1983) and others, we may write the inviscid momentum balance in each layer as

Here we have neglected drag at the channel boundaries, noting that it could become a dynamically significant process sufficiently far from the source (Hogg & Woods Reference Hogg and Woods2001). The sustained source of dense fluid leads to the boundary conditions

Figure 1.

Configuration of the two-layer gravity current driven by a sustained source of dense fluid in a horizontal channel. In this sketch, the layer of dense fluid includes a bore connecting two regions within one of which the interface is of constant height and within the other it is thinning towards the front.

Since $h_{1}+h_{2}=H$ , we deduce from (2.1) and (2.2) that

We may then eliminate the unknown pressure $p_{T}(x)$ by subtracting (2.4) and (2.5) and further use (2.7) to eliminate $u_{2}$ . Writing ${\it\rho}_{1}={\it\rho}_{2}(1+S)$ with $S>0$ , we find that momentum balance in the lower layer is expressed by (cf. Rotunno et al. Reference Rotunno, Klemp, Bryan and Muraki2011, Ungarish Reference Ungarish2011b )

At this point it is convenient to introduce dimensionless variables. We non-dimensionalise using the length scale $H$ and the velocity scale $(gSH)^{1/2}$ and replace $x$ and $t$ by their dimensionless counterparts, while we define the dimensionless velocity and height of the lower flowing layer by

This choice of scales introduces a dimensionless parameter, $\hat{Q}$ , which measures the magnitude of the flux per unit width in the channel, given by

Henceforth, unless stated otherwise, we deal with dimensionless variables.

We now treat dimensionless versions of (2.1) and (2.8) as the coupled governing equations for the system, from which the other dependent variables may be deduced. Potentially they form a hyperbolic system with characteristics $\text{d}x/\text{d}t={\it\lambda}_{\pm }$ , where

and where it is convenient to define

Real-valued characteristics exist and the system is strictly hyperbolic provided ${\it\alpha}^{2}+h(1-{\it\beta})>0$ , which corresponds to

As observed by Boonkasame & Milewski (Reference Boonkasame and Milewski2011), this condition of hyperbolicity may be compactly written in terms of dimensional variables by

The loss of hyperbolicity in this two-layer system will be shown to introduce an important constraint on the possible types of gravity current motion. Loss of hyperbolicity may be interpreted as an instability of the two-layer motion in the form of Kelvin–Helmholtz instability at a sufficiently long wavelength so that it is not filtered out by the assumption of hydrostatic balance (Boonkasame & Milewski Reference Boonkasame and Milewski2011).

Riemann invariants may be evaluated along each family of characteristics (Long Reference Long1955). On the characteristic $\text{d}x/\text{d}t={\it\lambda}_{\pm }$ , they satisfy

In the Boussinesq regime $(S\ll 1)$ , this differential equation may be integrated analytically to yield

and the system evolves such that $R_{\pm }$ is constant on the characteristic $\text{d}x/\text{d}t={\it\lambda}_{\pm }$ (Long Reference Long1955).

A second source condition may be applied to the system provided both characteristics propagate into the domain $x>0$ , a condition given by ${\it\lambda}_{\pm }>0$ . We note that for two-layer flow, the condition of criticality $({\it\lambda}_{-}=0)$ corresponds in general to

Thus for the source condition analysed in this study, we write

and this corresponds to a supercritical source if $F_{0}>(1+S)^{-1/2}$ . It is noteworthy that for the boundary condition studied here, this condition is identical to the criterion of supercriticality in a single-layer shallow-water model. If the source is subcritical $(F_{0}<(1+S)^{-1/2})$ then only one boundary condition (2.6) may be imposed at the source. Finally we need a dynamic condition at the front of the gravity current, which encompasses the non-hydrostatic behaviour there. Applying the model of Benjamin (Reference Benjamin1968) as a jump condition in a frame of reference moving with the front, we find

where $h_{N}=h(x_{N},t)$ and $u_{N}=u(x_{N},t)$ . This dynamic condition at the front (2.19) may only be applied provided the speed of propagation that it demands is less than the largest, local characteristic speed $(u_{N}<{\it\lambda}_{+}(u_{N},h_{N}))$ . If this criterion is not satisfied, then the flow is ‘choked’ and the motion cannot deliver fluid rapidly enough to match the required front speed (Ungarish Reference Ungarish2009, Reference Ungarish2011a ). The onset of ‘choking’ corresponds to the time at which the fastest speed of a characteristic first equals the dynamic front speed given by (2.19), $u_{N}={\it\lambda}_{\pm }(u_{N},h_{N})$ , which implies that

The onset of ‘choked’ flow therefore occurs at

which has solution $h_{N}=h_{c}\equiv 2\sin ({\rm\pi}/18)=0.3473$ and importantly from (2.20), $u_{N}=\hat{Q}+0.5273$ . The consequence of choking is that the fluid motion satisfies (2.19) when the height of the front is sufficiently shallow $(h_{N}<h_{c})$ and that the height of the front cannot exceed $h_{c}$ . We may therefore combine these conditions to write the condition at the front as

It is important to calculate the energy fluxes transported across the front to ensure that there is net dissipation. This is most conveniently evaluated in a frame moving with the front. Then the difference between the dimensionless energy fluxes transported away from and towards the front is given by

Substituting for the combined dynamic condition at the front, (2.19), we see that there is always dissipation $(\dot{\mathscr{D}}>0)$ . We emphasise that the results (2.20)–(2.23) are valid for both Boussinesq and non-Boussinesq regimes.

In what follows we examine the unsteady solutions that develop for different values of the three dimensionless parameters that characterise the system, namely the dimensionless volume flux per unit width, $\hat{Q}$ , the Froude number at the source, $F_{0}$ and the dimensionless excess density, $S$ . The height and velocity of the lower layer at the source is then given by $u=u_{0}\equiv (\hat{Q}F_{0}^{2})^{1/3}$ and $h=h_{0}\equiv (\hat{Q}/F_{0})^{2/3}$ (provided the source is supercritical). Furthermore the condition that the source must correspond to a strictly hyperbolic system of equations demands that (2.13) is satisfied. This in turn implies that for hyperbolicity,

In what follows we first analyse Boussinesq flows $(S\ll 1)$ in detail (§ 3), drawing out the different flow types that may occur for different values of the dimensionless parameters $(\hat{Q},F_{0})$ . This analysis is straightforward due to the existence of the Riemann invariants (2.16) and so many of the results are available as simple analytical expressions. We then analyse non-Boussinesq flows $(S=O(1))$ and show that the possible types of solutions are identical to the Boussinesq case, but that the boundaries between the different regimes change (§ 5). Our analysis employs Benjamin’s formula for the dynamic condition at the front of the gravity current (2.19); we note that for Boussinesq currents Borden & Meiburg (Reference Borden and Meiburg2013a ) have developed an alternative condition, while Huppert & Simpson (Reference Huppert and Simpson1980) propose an empirical front condition. The consequences of adopting their condition for the sustained flows studied in the paper are examined in § A.1.

3 Boussinesq currents $(S\ll 1)$

4 Simulations

We now proceed to analyse the flow by means of two-dimensional direct numerical simulation (DNS) in the Boussinesq regime, based on our code TURBINS (Nasr-Azadani & Meiburg Reference Nasr-Azadani and Meiburg2011), which employs a finite-difference discretization combined with a fractional projection method and TVD-RK3 time integration. Numerical details and validation results of the code are provided in Nasr-Azadani & Meiburg (Reference Nasr-Azadani and Meiburg2011) and Nasr-Azadani, Hall & Meiburg (Reference Nasr-Azadani, Hall and Meiburg2013), so that we provide only a brief review here.

The fluid motion is modelled using the dimensionless incompressible Navier–Stokes equations in the Boussinesq regime, given by

where $\boldsymbol{u}$ , $p$ and $\hat{{\it\rho}}$ are the dimensionless, two-dimensional velocity field, pressure and excess density $(0\leqslant \hat{{\it\rho}}\leqslant 1)$ , respectively, and where $\hat{\boldsymbol{z}}$ and $Re$ denote a unit vector aligned with the direction of gravity and the Reynolds number, respectively. For the purposes of the numerical simulations, it is convenient to non-dimensionalise the variables differently from the shallow-layer analysis of § 2: in the governing equations (4.1) and (4.2), we scale all lengths with respect to the height of the inflow, $h_{0}^{\ast }\equiv (Q^{2}/(gSF_{0}^{2}))^{1/3}$ , and time with respect to $t_{0}^{\ast }\equiv (h_{0}^{\ast }/(gS))^{1/2}$ . The Reynolds number thus is obtained as

and the dimensionless height of the channel by ${\hat{H}}=H/h_{0}^{\ast }$ . Under these scalings, and keeping the height of the source constant, the rescaled inflow velocity is $F_{0}$ and variations in the dimensionless flux parameter, $\hat{Q}$ , are obtained by altering the height of the computational channel, because $\hat{Q}=F_{0}{\hat{H}}^{-3/2}$ . Although the computations are performed with these new dimensionless variables, we analyse the motion quantitatively in terms of the variables defined in § 2, to enable immediate comparison between the shallow-layer theory and the computations.

To model the evolution of the dimensionless excess density field $\hat{{\it\rho}}(\boldsymbol{x},t)$ , we employ a continuum description and evolve the density field in an Eulerian manner via

Here, $Sc$ denotes the Schmidt number associated with the diffusion of the density field $\hat{{\it\rho}}$ . In most applications $Sc\gg 1$ , although test simulations indicate that the precise value of $Sc$ influences the results only weakly as long as $Sc\geqslant O(1)$ (Härtel et al. Reference Härtel, Meiburg and Necker2000). With this in mind, we employ $Sc=14$ in all of our simulations. We remark that in most simulations $Re$ is kept at 350 to ensure a stable interface between the dense current and the ambient fluid.

The computational set-up employs a channel of dimensionless size $L_{x}\times L_{z}=50\times {\hat{H}}$ , with in/outflow boundaries in the horizontal direction, cf. figure 1. The Cartesian grid in the $x$ - and $z$ -directions is uniformly spaced with ${\rm\Delta}x=0.1$ and ${\rm\Delta}z=0.016$ . We apply free-slip conditions at the top and bottom walls. At the outflow boundary, all flow variables, $q$ , are convected out of the domain via the outflow boundary condition

where $\bar{U}$ represents the maximum $u$ velocity value in the domain. For the density field, we impose no-flux conditions at the top and bottom walls. At the inlet and for $z<1$ , the density and horizontal velocity are set to unity and $F_{0}$ , respectively. These values decay smoothly over a few grid points to a zero value for $z>1$ . The computational domain is initially filled with ambient fluid at rest.

In figure 7 we plot contours of the concentration field at dimensionless times $t=20{\hat{H}}^{1/2}$ for varying values of the flux parameters $\hat{Q}$ , while keeping the source Froude number constant $(F_{0}=1.2)$ . To reiterate, here, and in the subsequent discussion, the variables are non-dimensionalised as described in § 2. It is noteworthy that the currents maintain a relatively sharp interface between the two fluids in the channel; turbulence and other diffusive processes are insufficient to cause significant mixing on these length and time scales. The other noticeable feature is that for the higher values of the dimensionless flux, $\hat{Q}$ , the dense current has completely filled the channel depth in an expanding region close to the source, while further downstream the current re-establishes the two-layer structure.

Figure 7.

Contours of the dimensionless excess density field from the DNS at $t=20{\hat{H}}^{1/2}$ for source Froude number, $F_{0}=1.2$ , and various values of the dimensionless flux, $\hat{Q}=0.15,0.28,0.42,0.70,0.78,0.94$ . The ten contours in each panel correspond to values of $\hat{{\it\rho}}$ from 0.1 to 1 with increments of 0.1.

We calculate the rate at which the current progresses downstream by defining a scaled current height

based on the vertical integral of the density field. The front of the current, $x_{N}(t)$ is then determined by finding the maximum value of $x$ for which ${\it\Phi}>{\it\epsilon}$ , where ${\it\epsilon}$ is a prescribed threshold, typically ${\it\epsilon}=0.01$ . We plot $x_{N}(t)$ as a function of time for each computational run in figure 8(a). We observe that the currents move at constant speeds dependent on $\hat{Q}$ . By fitting a linear trend line to each set of data we measure this speed and compare with the prediction of the shallow-layer model in figure 8(b). The agreement between the model and the simulations is quite good. The Navier–Stokes currents are seen to move slightly more slowly than predicted by the shallow-layer model, which is most likely due to their relatively low Reynolds number of 350, implying the weak, but non-vanishing, action of viscous processes that retard the motion.

We examined the effects of varying the Reynolds number in two ways. First for $F_{0}=1.2$ and $\hat{Q}=0.42$ , we performed simulations with Reynolds numbers ranging from 100 to 1500 and we measured the speed of the flow. Computing at higher Reynolds numbers required greater spatial resolution to achieve converged results and the computational grids were adjusted appropriately. The results are given in table 1, where we observe that the measured speed systematically increases with increasing Reynolds number, approaching the shallow-layer value. This supports the view that the slight mismatch between the simulations and the inviscid shallow-layer prediction is due to the weak, but non-vanishing, effects of viscosity.

Figure 8.

(a) The position of the front of the current, $x_{N}(t)$ , as a function of time for each of the simulations with $F_{0}=1.2$ ; and (b) the measured velocity of the current, $u_{N}$ , as a function of the dimensionless flux, $\hat{Q}$ at $Re=350$ (○) and $Re=1500$ ( $\times$ ). Also plotted are the predictions from the shallow-layer analysis (solid line) and the asymptotic result when $\hat{Q}\ll 1$ (dotted line) and $\hat{Q}=0.2938$ at which value the flow becomes choked.

Table 1.

The dimensionless speed of the front of the current, $u_{N}$ , measured from simulations conducted at different Reynolds numbers, when the source Froude number, $F_{0}=1.2$ and the dimensionless flux per unit width, $\hat{Q}=0.42$ . The shallow-layer prediction for these conditions is $u_{N}=0.952$ .

Contours of the concentration field from each run at different Reynolds number are plotted in figure 9 at the same instant of time. The increased size and intensity of vortices on the interface between the two fluids becomes more evident as the Reynolds number is increased and these vortices lead to zones within which the density takes intermediate values between the intruding and ambient fluids. However over the time and length scales of these simulations, the emergence of this mixed zone does not significantly modify the speed of the current, presumably because the vertical integral of the excess density field over the depth of the current, $(\int _{0}^{h}\hat{{\it\rho}}\,\text{d}z)$ , remains approximately constant under mixing and this determines the speed of the density-driven spreading. Johnson & Hogg (Reference Johnson and Hogg2013) have demonstrated that entrainment alone modifies the predictions of shallow-layer models, but only at time and length scales far beyond those computed in these simulations.

We also examine the effects of Reynolds number by performing simulations at $Re=1500$ for source Froude number $F_{0}=1.2$ and a range of values of the dimensionless flux, $\hat{Q}$ . In each case, as for the simulations run with $Re=350$ , we find the front of the current moves with constant speed and that this speed compares favourably with the shallow-layer predictions (see figure 8 b). For each case, the simulations at higher Reynolds number led to slightly higher front speeds, evidencing the reduced effects of viscosity (table 1).

The shallow-layer predictions always comprise currents with rarefactions, with the transition to choked flows occurring at $\hat{Q}=0.2938$ and the loss of hyperbolicity, with the associated transition to channel-filling flows, at $\hat{Q}=0.6944$ . This latter value closely agrees with the value observed in the Navier–Stokes simulations. Also plotted is the asymptotic result for a relatively deep ambient $(\hat{Q}\ll 1)$ . In this regime, to leading order, the Riemann invariants are given by $R_{\pm }=u\pm 2\sqrt{h}$ on characteristics $\text{d}x/\text{d}t=u\pm \sqrt{h}$ . For this source condition there is an expansion fan within which $R_{+}=\hat{Q}^{1/3}(F_{0}^{2/3}+2F_{0}^{-1/3})$ . The leading edge of the expansion fan corresponds to $u_{N}=\sqrt{2h_{N}}$ and thus

a result which is quite accurate for $\hat{Q}\leqslant 0.3$ when $F_{0}=1.2$ (see figure 8). This implies that the gravity current motion is relatively unaffected by the motion of the overlying less dense fluid until the flow becomes choked at $\hat{Q}=0.2938$ . Notably when $\hat{Q}\ll 1$ , the dimensionless rate of propagation is proportional to $\hat{Q}^{1/3}$ , whereas for larger values, the speed varies linearly with $\hat{Q}$ and from (2.20) is given by $u_{N}=0.5473+\hat{Q}$ .

Figure 9.

Contours of the dimensionless excess density field from the DNS at $t=20{\hat{H}}^{1/2}$ for source Froude number, $F_{0}=1.2$ , and dimensionless flux, $\hat{Q}=0.42$ , at various Reynolds numbers. The ten contours in each panel correspond to values of $\hat{{\it\rho}}$ from 0.1 to 1 with increments of 0.1.

We may also investigate the motion of the interface between the dense and less dense fluids at the upper boundary when the flow is channel filling. In this case, we define the position $x_{r}(t)$ as the maximum value of $x$ for which ${\it\Phi}>1-{\it\epsilon}_{1}$ , where ${\it\epsilon}_{1}$ is another prescribed threshold with chosen value ${\it\epsilon}_{1}=0.05$ . We note from figure 7 that for those flows which undergo the transition to a channel-filling motion, although the interface at the top boundary is relatively sharp (i.e. the region of intermediate density is quite narrow), the $\hat{{\it\rho}}=1$ contour is somewhat displaced from the upper boundary. This means that the threshold ${\it\epsilon}_{1}$ cannot be chosen to be arbitrarily small or else the criterion for determining the position $x_{r}(t)$ would fail. However, we find that $x_{r}(t)$ depends only very weakly on ${\it\epsilon}_{1}$ for $0.04<{\it\epsilon}_{1}<0.1$ . In figure 10 we plot $x_{r}(t)$ as a function of time for those runs for which the flow becomes channel filling and we overlay the theoretical prediction. It is evident that there is some temporal offset before the channel-filling layer becomes established, presumably due to the development of the interfacial instability. Thereafter, the theoretical model captures the measured speed quite well until the basal front of the two-layer current leaves the domain. Beyond that time, the existence of the outflow boundary prevents the flow in the domain from continuing to be influenced by the motion of the current front. Longer computational domains are expected to produce agreement between the two sets of data for longer times. We examine the motion of this interface at the top of the channel for simulations at higher Reynolds number ( $Re=1500$ , see figure 10). We observe that at higher Reynolds number, the channel-filling part of the flow takes longer to be fully established and that for the lowest flux considered, $\hat{Q}=0.70$ , it is not completely initiated before the front of the flow leaves the computational domain. For the other cases in which this channel-filling transition occurs, however, the measured positions are broadly in accord with the theoretical predictions after some initial period during which the flow structure is developed.

Figure 10.

The position of the interface between the dense and less fluids at the upper boundary of the channel, $x_{r}(t)$ as a function of time for cases in which the gravity current becomes channel filling close to the source from the DNS results at $Re=350$ (solid lines), at $Re=1500$ (dotted lines) and the shallow-layer modelling (dashed lines). The DNS results are curtailed when the basal front leaves the fluid domain.

Gravity currents due to the sustained emission of relatively dense fluid within a channel have recently been studied by Shringapure et al. (Reference Shringapure, Lee, Ungarish and Balachandar2013). These authors conducted numerical simulations of the Navier–Stokes equations and compared their results with a simple analytical model of the motion in which the currents are spatially uniform. Here we quantitatively investigate their results in terms of the shallow-layer formulation developed above (§ 3). We note that the source conditions they treated were subcritical $({\it\lambda}_{-}(u_{0},h_{0})<0)$ . Thus the flows are subject to downstream influences and adjust close to source to critical conditions $({\it\lambda}_{-}(u_{0},h_{0})=0)$ , which implies $F_{0}=1$ . We may then compare the numerically determined speed of propagation and average height over a small frontal region with the shallow-layer model developed here. We plot in figure 11 the dimensionless front position and height at the front as a function of the dimensionless flux, noting that there is very good agreement between the two. This provides further confidence in the shallow-layer model since it is capable of accurately reproducing features in an independently developed numerical simulation.

Figure 11.

The dimensionless speed, $u_{N}$ , and height, $h_{N}$ , of the front as a function of the dimensionless flux per unit width, $\hat{Q}$ calculated from the shallow layer with $F_{0}=1$ (solid lines) and from the numerical simulations of Shringapure et al. (Reference Shringapure, Lee, Ungarish and Balachandar2013) ( $\times$ ). Also marked is the value of the flux, $\hat{Q}_{c}$ at which the motion first becomes choked (dashed line).

5 Non-Boussinesq flows $(S=O(1))$

For non-Boussinesq currents, the same types of flow structures arise, namely, uniform and choked flows with shocks and rarefactions. However the boundaries between the different states in the $(\hat{Q},F_{0})$ plane vary with $S$ and our method of constructing them differs from the Boussinesq regime. Instead of analytical expressions for the Riemann invariants, we must now largely rely on numerically computed solutions to the differential equation along each family of characteristics (2.11).

Non-Boussinesq currents with no spatial variations between the source and the front may be determined through the solutions of coupled algebraic equations; they satisfy (3.1) for $0\leqslant \hat{Q}\leqslant Q_{c}$ . This produces a current identical to the Boussinesq case (as plotted in figure 3). For larger values of the dimensionless flux $(\hat{Q}>\hat{Q}_{c})$ , the flow becomes choked at the front and there is a rarefaction of ${\it\lambda}_{+}$ characteristics as the dense fluid accelerates to match the front condition. The curve in the $(\hat{Q},F_{0})$ plane is mapped out by integrating

subject to $u=u_{c}$ when $h=h_{c}$ and $u=u_{0}$ when $h=h_{0}$ , thus providing an implicit relationship between the parameters $\hat{Q}$ and $F_{0}$ . The form of this curve is dependent upon the dimensionless excess density, $S$ , as illustrated in figure 12 for $S=1$ and $S=5$ . However it is evident that there always exist some parameter values for which there is ‘choked’ flow because $\hat{Q}_{c}<F_{0c}(F_{0c}^{2}+S/(1+S))^{-3/2}$ , the condition for hyperbolicity (2.24) evaluated at $\hat{Q}=\hat{Q}_{c}$ and $F_{0}=F_{0c}$ . As in § 3.1, we denote the curves on which there is uniform and choked-uniform flow by $F_{U}$ and $F_{Uf}$ , respectively.

Figure 12.

Regimes of solutions, classified as function of the dimensionless flux, $\hat{Q}$ and the source Froude number, $F_{0}$ for non-Boussinesq flows (a) $S=1$ and (b) $S=5$ . As in figure 2, uniform flows correspond to $F_{0}=F_{U}\cup F_{Uf}$ . Smaller, but still supercritical, values of $F_{0}$ yield rarefactions, while larger values lead to shocks. The curve corresponding to the loss of hyperbolicity is denoted $\mathscr{C}_{h}$ , while the transition between choked and unchoked flows of rarefactions and shocks are denoted $\mathscr{C}_{R}$ and $\mathscr{C}_{S}$ , respectively. Critical conditions are given by $F_{0}=1/\sqrt{2}$ for $S=1$ and $F_{0}=1/\sqrt{6}$ for $S=5$ ; these values are also plotted.

For $F_{0}<F_{U}\cup F_{Uf}$ , there are rarefactions close to, or attached to, the source and potentially close to the front if the motion is ‘choked’. We may compute the boundary between ‘choked’ and unchoked rarefaction by integration across the rarefaction of ${\it\lambda}_{-}$ characteristics close to source. Thus we integrate

from $u=u_{0}$ and $h=h_{0}$ to $u=u_{c}$ and $h=h_{c}$ and this provides an implicit relationship between the dimensionless source flux and source Froude number. The form of this curve, $\mathscr{C}_{R}$ , is dependent on $S$ as shown in figure 12.

When $F_{0}>F_{U}\cup F_{Uf}$ , the flow involves shocks across which the height and velocity vary discontinuously. As with the Boussinesq currents (§ 3), the jump conditions feature the conservation of mass and momentum fluxes, together with the assumption that there is no dissipation in the upper, less dense fluid. Then the currents may either comprise two uniform segments, joined by a shock, or may be choked at the front, in which case the current thins and accelerates downstream of the shock. The boundary between these two possibilities is similar to that deduced above (§ 3.3), although the momentum equations on either side of the shock are modified to take account of the potentially non-negligible difference in densities between the two fluids (see § 2). We may then compute which type of current exists for different values of the dimensionless flux, $\hat{Q}$ , and the source Froude number, $F_{0}$ (see figure 12). We note that for the extent of values for which choked flow exists diminishes as the density ratio increases.

In figure 13(a,b), we plot the profiles of height and velocity for non-Boussinesq currents for different ratios of the densities of the fluids. In particular, we evaluate the form of the currents when the source is critical, $F_{0}=(1+S)^{-1/2}$ and the flux $\hat{Q}={\it\mu}(1+S)^{-1/2}$ with $0\leqslant {\it\mu}\leqslant 1$ . (In this range of values the source is always hyperbolic.) We observe that for all of the flows the current immediately transitions from their source conditions through a rarefaction to a uniform state. Furthermore, for the cases plotted, the flow is only choked for relatively large dimensionless fluxes with $S=1$ . As the density ratio becomes large $(S\gg 1)$ , the motion of the less dense fluid in the upper layer plays only a negligible role and the motion is modelled by the single-layer shallow-water equation to leading order. In this regime, we must adopt different variables, namely $\overline{t}=(1+S)^{-1/2}t$ , $\overline{u}=(1=S)^{1/2}u$ and $\overline{Q}=(1+S)^{1/2}\hat{Q}$ . The front condition is then given by $h=0$ to leading order and the velocity and height fields satisfy

In figure 13(c) we plot the height and velocity fields for $\hat{Q}=F_{0}=(1+S)^{-1/2}$ for a range of values of $S$ in terms of these rescaled variables and we note that as $S$ increases, the single-layer results are progressively attained.

Figure 13.

The height, $h$ , and velocity, $u$ as function of $x/t$ for non-Boussinesq, supercritical source conditions $F_{0}=(1+S)^{-1/2}$ ; (a) $S=1$ and (b) $S=5$ with $\hat{Q}=n(1+S)^{-1/2}/5$ $(n=1,2,3,4,5)$ . Also plotted in (c) are the height and rescaled velocity fields as a functions of $(S+1)^{1/2}x/t$ for $F_{0}=\hat{Q}=(1+S)^{-1/2}$ .