3.1 The suite of simulations

We implemented two-dimensional direct numerical simulations (DNS) with a similar set-up as in figure 1(c) using OpenFoam v2.4 (interFoam multiphase solver), which is an open-source code developed to solve the two phase fluid (gas and/or liquid) continuity and Navier–Stokes equations. The details of the code, numerical solver, mesh grid cells, domain size, governing equations, a sketch of the simulation set-up and boundary conditions are described in appendix A. Table 2 shows the details of the direct numerical simulations performed and their initial conditions. Figure 3 depicts a typical snapshot from the DNS in the late-time period of case D1 at $t=0.38$  s. One can compare this figure with the profiles of figure 2; it is clear that both figures exhibit similar shapes, which will be compared in detail below. Note that case D1 is our basic case in which water is the gravity current and air is the ambient fluid. In the next cases, we change the length, density and viscosity of the gravity current by factors of 10 to validate our reduced-order model, which is described in § 4. We decrease the viscosity ratio in case D5 to confirm that it is not important since we focus on cases where inertial effects are dominant and the density ratio $\unicode[STIX]{x1D70C}_{c}/\unicode[STIX]{x1D70C}_{a}\gg 1$ .

Figure 3. A representative result of the direct numerical simulation of the release of a fixed volume of water at early period, case D1. $\unicode[STIX]{x1D6FC}_{water}$ denotes the fraction of liquid phase in each grid cell; see appendix A for more information.

3.2 The Froude number of the current near the edge

The Froude number at the edge is defined as $Fr\equiv u/\sqrt{gh}|_{x=L}$ . From the DNS, we can compute useful depth-averaged values such as the velocity at the edge, $x=L$ and then determine the Froude number. In particular, we have vertically averaged the horizontal velocity, $u(x=L,t)$ , at the edge and divided it by $\sqrt{gh(x=L,t)}$ to obtain the $Fr$ , shown in figure 4(a), which summarizes results from the five cases. Note that we will derive the non-dimensionalization of all the variables from the governing equations in the next section. Our DNS results indicate that the $Fr$ increases in the initial adjustment phase and then approaches an approximately constant value at late times. Moreover, $Fr$ reaches a constant value at almost the same dimensionless time except for case D2, which has a smaller value of $h_{0}/L$ so its initial adjustment process occurs faster (before the wave is reflected off the wall and exits). Nevertheless, we tested many DNS cases with different ratios of $h_{0}/L$ and found that they collapse after a time $t=c(L/\sqrt{gh_{0}})$ , where $c\approx 2{-}3$ and $Fr$ remains constant in the late-time inertially dominated regime. Furthermore, figure 4(b) exhibits the height and velocity profiles for case D1 at late times. $Fr$ increases gradually from zero at the wall (since the horizontal velocity at the wall is zero) to $\simeq 1$ at the edge, as shown in figure 4(b). We will provide the theoretical background to explain why this result for the $Fr$ , i.e. $Fr\simeq 1$ , is actually expected and will use it in the next section to characterize the right boundary condition. These results from the DNS indicate that the flow transits from a subcritical regime before the edge to a critical depth where $Fr\simeq 1$ near the edge.

Figure 4. (a) $Fr\equiv u/\sqrt{gh}|_{x=L}$ versus time from all of the DNS data, where $u$ is the vertically averaged horizontal velocity; $\unicode[STIX]{x1D70F}$ is the dimensionless time defined in equation (4.2), which is derived with the typical gravity wave speed and the horizontal length of the domain. $Fr$ is found to be a constant close to 1 after an early transition period. (b) The height and velocity profiles and the corresponding $Fr$ shown at $\unicode[STIX]{x1D70F}=4$ for case D1. The $Fr$ increases from zero at the wall to ${\approx}1$ at the edge.

4.1 Governing equations and boundary conditions

Assuming that the flow is approximately hydrostatic, the height of the current is much smaller than its length $h\ll L$ and viscous effects are negligible compared to the inertial effects, because $Re$ is high during the initial adjustment and inertial periods, the continuity and momentum equations can be vertically averaged to obtain the shallow-water equations (see e.g. Ungarish Reference Ungarish2009, ch. 2). Therefore, for the configuration of figure 1(d), the governing equations of a gravity current can be written as:

where $x$ is the horizontal coordinate, $t$ is time, $g^{\prime }\equiv g|\unicode[STIX]{x1D70C}_{c}-\unicode[STIX]{x1D70C}_{a}|/\unicode[STIX]{x1D70C}_{c}\simeq g$ is the reduced gravitational constant in our non-Boussinesq system where $\unicode[STIX]{x1D70C}_{a}\ll \unicode[STIX]{x1D70C}_{c}$ , $u(x,t)$ denotes the vertically averaged velocity and $h(x,t)$ represents the thickness of the gravity current. We seek to solve these equations with the initial conditions $u(x,0)=0$ and $h(x,0)=h_{0}$ . We define non-dimensional variables as:

where $\sqrt{g^{\prime }h_{0}}$ is the typical speed of shallow gravity waves. Now, we can rewrite equations (4.1) in terms of the dimensionless variables as:

which are to be solved with initial conditions: At the early times, the height of the fluid at the left boundary decreases slowly (the height is almost constant until influenced by the gravity wave following the release of the gate). The characteristic directions, $C_{\pm }=U\pm \sqrt{H}$ , play a significant role in the initial motion of the flow. Figure 4(a) indicates that the early period occurs when $\unicode[STIX]{x1D70F}\lesssim 2$ since it takes $t\simeq L/\sqrt{g^{\prime }h_{0}}$ for a gravity wave to travel to the end wall and about the same amount of time to return to the free end. The rarefaction wave is associated with the left-travelling characteristic direction $C_{-}\approx -\sqrt{H}$ (since at early times $U\approx 0$ ). As time evolves, $U$ increases and $H$ decreases at the edge such that the characteristic direction $C_{-}\approx 0$ prevails (choked condition, see e.g. Ungarish (Reference Ungarish2009)), which yields the $Fr\approx 1$ condition, where $Fr=U/\sqrt{H}|_{(1,\unicode[STIX]{x1D70F})}$ is the Froude number at the edge of the domain. Figure 4(a) also suggests that during the early-time transition period, the $Fr$ at the edge increases linearly as $Fr=\unicode[STIX]{x1D70F}/\unicode[STIX]{x1D70F}_{t}$ in all the numerical cases, in which $\unicode[STIX]{x1D70F}_{t}$ represents the transition time between the early and late time periods. We observed that the transition time $\unicode[STIX]{x1D70F}_{t}\approx 2$ for the DNS examples we considered (table 2), except for case D2, when the initial aspect ratio is different from the rest, which influences the transition time. The direct numerical simulations shown in figure 4 also confirm that $Fr\simeq 1$ captures the DNS results well after an initial transition period. Hence, based on the analysis of the characteristic directions that is confirmed by the DNS results, the boundary conditions of the two coupled, first-order partial differential equations (PDEs) can be written as:

Hence, we have divided our problem into two regimes: (i) an early period where the initial conditions of the flow are important and (ii) a late period, which holds when the effects of initial conditions are no longer important and the solutions describe an asymptotic state of the inertial dynamics. A solution of (4.3) subject to (4.4) and (4.5) describes the dynamics of the flow in both of these regimes and can be obtained numerically. In addition, an analytical solution for the late period is obtainable as we show later. At sufficiently long times, the film becomes so thin that viscous effects become significant; however, we do not aim to investigate this limit here since the flow regime will change into a low Reynolds number flow for which the shallow-water equations (4.1) do not hold anymore (e.g. Zheng et al. (Reference Zheng, Soh, Huppert and Stone2013) studied this viscous limit for gravity current drainage from the edge of a porous reservoir). The limit in which the shallow-water equations are still valid and an inertial phase exists can be established quantitatively. These equations are valid when the time to propagate a gravity wave a distance $L$ , $L/\sqrt{g^{\prime }h}$ , is shorter than the time at which viscous effects are felt across the film, $h^{2}/\unicode[STIX]{x1D708}$ . This means that $h/L>Fr/Re$ . Therefore, the focus of the paper is to investigate the late inertial period of this drainage problem, i.e. $\unicode[STIX]{x1D70F}\gtrsim \unicode[STIX]{x1D70F}_{t}$ .

The boundary condition for the height during the late period, equation (4.5), is given through this $Fr$ relation at the edge, which was found for all times and both regimes using our DNS data as a priori analysis in § 3.2. The $Fr$ condition has been a widely used approach to characterize high Reynolds number gravity currents in different contexts. For instance, many such studies impose a constant $Fr$ jump condition at the nose of the propagating currents (see Rottman & Simpson Reference Rottman and Simpson1983; Monaghan et al. Reference Monaghan, Meriaux, Huppert and Monaghan2009; Ungarish Reference Ungarish2009). The $Fr$ condition has also been used in other contexts such as in the hydraulic control problems for the flows over obstacles see e.g. Armi & Farmer (Reference Armi and Farmer1986), for maximal flow exchange through a contraction. The control separates a subcritical ( $Fr<1$ ) regime from a supercritical ( $Fr>1$ ) regime. In such configurations, the flow can always travel toward the control (analogous to the gate here) from the subcritical side (similar to our problem) but not from the supercritical side.

4.2 Numerical solutions of the governing PDEs

We numerically solved the shallow-water PDEs (4.3) with given initial conditions (4.4) and boundary conditions (4.5). To do so, we used a finite-volume discretization method and wrote the PDEs in a conservative form. The details of solving the PDEs and the numerical methods employed, along with the validation of the code, are provided in appendix B.

We present solutions for the time evolution of $H(X,\unicode[STIX]{x1D70F})$ and $U(X,\unicode[STIX]{x1D70F})$ . Figures 5(a) and 5(d), exhibit the height and velocity profiles respectively at early times where the gate is lifted and one can observe the left-travelling wave caused by this release. Once this wave approaches the end wall, it is reflected to the right and hence in figure 5(b,e) we observe a right-travelling wave. Finally, at late times around $\unicode[STIX]{x1D70F}\approx 5$ , the profiles shown in figure 5(c,f) reach the similar shapes, which we call the self-similar phase.

Moreover, this flow behaviour is observed in our DNS results. For example, figure 6 shows the profile shape and velocity magnitude contours in the early period, figure 6(a), and the late period, figure 6(b). The velocity magnitude increases in $X$ and peaks at the edge. The direction of the velocity vectors in the figure indicates that the flow velocity is more horizontal at late times (figure 6 a) compared to early times (figure 6 b). Hence, this figure confirms the assumption that the flow is mainly horizontal for the shallow-water PDE and it holds when $h\ll L$ . The numerical solution of the governing PDEs and the DNS results clearly demonstrate the existence of a self-similar phase. In the next subsection, we will derive the analytical solution for this phase and validate it with the shallow-water solutions. Finally, § 5 will validate the self-similar solutions with the help of experiments and DNS.

Figure 5. Numerical solution of the one-dimensional equations (4.3) in which (a,d) indicate an initial left-travelling wave caused by lifting the gate, (b,e) show the wave reflected by the wall and (c,f) exhibit the late-time profiles that reach a self-similar shape.

Figure 6. The contour of the velocity magnitude (where $U$ and $V$ respectively represent the horizontal and vertical velocities) for the DNS case D1. (a) The early-time behaviour and (b) the late-time behaviour in the self-similar phase. The black curves show the height of the flow and the arrows indicate the velocity vectors.

4.3 Self-similar solutions

Using equations (4.3), and recognizing that the domain has a fixed length, one can employ an appropriate scaling argument to infer the time evolution of the normalized height and speed of the gravity current in the self-similar phase. For instance, balancing the two terms in the continuity equation (4.3a ) yields $H/\unicode[STIX]{x1D70F}\approx -UH$ , which implies $U(X,\unicode[STIX]{x1D70F})\propto \unicode[STIX]{x1D70F}^{-1}$ . A quantitative comparison will be made below. Using this time dependence of the scaled velocity and balancing the three terms in (4.3b ) leads to $H(X,\unicode[STIX]{x1D70F})\propto \unicode[STIX]{x1D70F}^{-2}$ . Thus, we introduce

where $\widetilde{U}$ and $\widetilde{H}$ are, respectively, the similarity functions for velocity and height, and $\unicode[STIX]{x1D6FE}$ is a constant time shift that will be determined to match the shallow-water PDE solutions (e.g. the fluid volume or the fluid height at the edge). Substituting equations (4.6) into (4.3) yields two nonlinearly coupled first-order ordinary differential equations (ODEs) in which the time dependencies cancel, with the boundary conditions

We will integrate our derived ODE system in equations (4.7) analytically. The two ODEs can be rearranged as

This autonomous system can be written as an ODE for $\widetilde{U}(\widetilde{H})$ ,

Integrating this equation and using the left boundary condition $\widetilde{U}(X=0)=0$ with $\widetilde{H}(0)$ (to be determined) gives

It then follows that the $Fr$ condition at the free end, $\widetilde{U}/\sqrt{\widetilde{H}}\rightarrow 1$ as $X\rightarrow 1$ , leads to $\widetilde{H}(1)/\widetilde{H}(0)=3/4$ . We can then construct the shape of the current by substituting into equation (4.9b ) to find

Integrating this equation we obtain an analytical solution for the profile shape as

Evaluating (4.13) at $X=1$ and considering $\widetilde{H}(1)/\widetilde{H}(0)=3/4$ we find $\widetilde{H}(0)=[6/(3\sqrt{3}-\unicode[STIX]{x03C0})]^{2}\approx 8.5$ . The corresponding velocity distribution $U(X)$ follows from (4.11).

Based on the obtained analytical solutions, the constant time shift, $\unicode[STIX]{x1D6FE}$ , in (4.6b ) can be determined by matching the fluid height at the edge. Following our numerical solutions of the governing PDEs in figure 5  $H(X=1,\unicode[STIX]{x1D70F})$ can be found, and considering $\widetilde{H}(1)\approx 8.5\times 3/4$ we obtain $\unicode[STIX]{x1D6FE}\approx 2$ . The validity of these solutions compared to the DNS and experimental results will be assessed in the next section.

We note that many of the previous similarity solutions are for problems in which the volume of the gravity current does not change, e.g. the classic paper by Rottman & Simpson (Reference Rottman and Simpson1983). The dam-break problems have constant volume; however, the volume in our domain decreases and thus we have different scaling results for the time dependence and a different boundary condition, i.e. $Fr\simeq 1$ at the edge. Hence, our self-similar solutions are different from the self-similar flows in the reservoir domain of spreading gravity current problem, produced from a sudden dam-break mechanism (e.g. Ungarish Reference Ungarish2009, ch. 10). In particular, our approach in this paper was to first balance the terms in the continuity and momentum equations, considering the fixed length of the domain, to infer the time dependence in the similarity phase (4.6) and then use them to obtain the self-similar profile shapes.

We have also plotted the PDE solutions of equations (4.3) against the solutions of the ODEs in equations (4.7) at different times to investigate whether the PDE solutions will converge to the obtained self-similar solutions. As figure 7 shows, the PDE solutions first undergo an early adjustment phase and finally they approach the self-similar solution after some small oscillations. These results clearly demonstrate good agreement between the self-similar solution and the governing PDE solutions for both the velocity and height profiles.

Figure 7. Time evolution of (a) the profile shape $H(X,\unicode[STIX]{x1D70F})$ and (b) the velocity field $U(X,\unicode[STIX]{x1D70F})$ (dashed lines) from numerically solving the one-dimensional PDE (4.3) subject to (4.4) and (4.5). The self-similar solutions are also plotted (solid lines) from (4.13).

5.1 Validating the scaling arguments and self-similar solutions using experiments

In order to validate the scaling arguments in § 4, we have plotted in figure 8(a) the height time series for all experimental measurements recorded in the middle of the domain ( $x=L/2$ ). Figure 8(b) depicts the same data with dimensionless variables, in which all six curves collapse very well at late times. For comparison, the solid line indicates the $\unicode[STIX]{x1D70F}^{-2}$ scaling that our reduced-order model (4.6b ) predicts. We can observe from figure 8(b) that the height at $X=1/2$ decreases as $(\unicode[STIX]{x1D70F}+2)^{-2}$ in the similarity phase, $\unicode[STIX]{x1D70F}\gtrsim 2$ , and our scaling theory shows good agreement with the experimental results.

Figure 8. Time evolution of the height of the gravity current at the centre of the tank ( $X=0.5$ ) for all of the experiments (E1–E6): (a) dimensional data and (b) rescaled data. The theoretical prediction of the self-similar solution, i.e. $(\unicode[STIX]{x1D70F}+2)^{-2}\widetilde{H}(0.5)$ , is also plotted in (b), which agrees well with the rescaled heights after an initial transition period, $\unicode[STIX]{x1D70F}\approx 2$ .

Now we compare the experimental data for the profile shape of the current with the results of the reduced-order model. Figure 9(a) displays the dimensional values of the profile evolution while figure 9(b) shows the non-dimensional values. This figure exhibits good agreement between the experiments in the self-similar phase and the reduced-order model; the maximum relative error is less than 10 % for each measured point of case E6 that is shown. The transition period is also clear from these results where the normalized height of the flow gradually increases up to the similarity phase.

Next, we present the experimental depth profiles at late times ( $\unicode[STIX]{x1D70F}\approx 5$ ) when all data should be in the similarity phase. The results in figure 10 collapse well after non-dimensionalizing for all of the performed experiments. Note that in the self-similar phase, $Re=uh/\unicode[STIX]{x1D708}_{c}=Fr\sqrt{g^{\prime }h}h/\unicode[STIX]{x1D708}_{c}$ , which is approximately $10^{2}{-}10^{4}$ at $\unicode[STIX]{x1D70F}\approx 5$ for the different cases. This good agreement demonstrates that our reduced-order model is able to accurately predict the gravity current behaviour in the self-similar phase for different conditions.

Figure 9. Time evolution of the profile shapes in the similarity phase for experimental case E6: (a) dimensional data and (b) rescaled data. The predictions of the self-similar solutions are plotted using the analytical solution (4.13) of the ODEs (4.7).

Figure 10. The profile shapes of all the experiments in the self-similar phase: (a) dimensional data and (b) rescaled data. The prediction of the self-similar solution is plotted in (b), which is found to agree well with the rescaled experimental results after the transition period.

5.2 Validating scaling arguments and self-similar solutions using DNS

The scaling arguments implied that the heights of the profiles have to decrease as $\unicode[STIX]{x1D70F}^{-2}$ in the self-similar phase. Figure 11 shows the results of the DNS at $X=0$ , $X=1$ and an averaged height throughout the domain for case D1. The results in the figure demonstrate that the DNS data also confirm our scaling arguments.

Figure 11. Height versus time for the case D1 of the DNS results: (a) dimensional data and (b) rescaled data. The slope of the scaling argument of equation (4.6b ) is shown in (b), matching the DNS data well.

To further validate our scaling arguments, the height time series at $X=1/2$ for all the DNS results are investigated. Figure 12 shows the depth time series for all cases in terms of raw data (figure 12 a) and non-dimensionalized values (figure 12 b). All the simulations collapse well after non-dimensionalization and match with the theoretical solution obtained from (4.6).

Figure 12. Height versus time for the DNS data at $X=1/2$ for all cases: (a) dimensional data and (b) rescaled data. The prediction of the self-similar solution, i.e. $(\unicode[STIX]{x1D70F}+2)^{-2}\widetilde{H}(0.5)$ , is plotted in (b), which agrees with the rescaled heights after the transition period.

Furthermore, non-dimensional depth profiles from DNS data are in good agreement with the obtained analytical similarity solution. Figure 13 depicts our dimensional DNS results and their non-dimensional values versus the similarity solution.

Figure 13. The profile shapes of the DNS results in the self-similar phase: (a) dimensional data and (b) rescaled data. The prediction of the similarity solution is shown in (b). The oscillations for case D2 are related to the vertical resolution and numerical accuracy of the simulations.

6.1 Practical applications

One of the practical applications of this reduced-order model could be in characterizing reservoir drainage problems in any configuration that results in $Fr\approx 1$ at the boundary (free fall with a perpendicular edge as in our experiments and DNS or simply a very steep slope on the ambient fluid side of the problem). For instance, when a wall in a fluid container is removed, when a tank of fluid (e.g. a tank of water in various structures) breaks or when a spillway is damaged, the discharge rate in the similarity phase could be predicted by our reduced-order model. A potential practical application concerns a break of a spillway structure on top of a dam such as the one that recently occurred in the Oroville dam in California (http://www.cnn. com/2017/02/12/us/california-oroville-dam-failure/). Figure 14(a) displays the DNS results for the volume $V_{out}$ that is spilled from the domain. For the initial conditions we considered, the results indicate that ${\approx}50\,\%$ of the initial volume exits at early times (initial adjustment period $\unicode[STIX]{x1D70F}\lesssim 2$ ) and ${\approx}40\,\%$ of the initial volume is removed during the similarity phase by the end of the simulations (which takes approximately two times longer than the initial adjustment period). The discharge rate of this drainage is also shown in figure 14(b). The discharge gradually increases at early time periods and, after reaching a maximum, it decreases in the similarity phase. Note that at early periods, the drainage rates might be approximated by some available gravity current solutions using simplified assumptions although further research is required to establish the interface shape of the draining edge and its dynamics in the initial adjustment period. For instance, the DNS drainage value is in good agreement with the simple dam-break solutions (i.e. $Q(X=1,\unicode[STIX]{x1D70F}=2)=8/27$ see Ungarish (Reference Ungarish2009)) for 1 ${\lesssim}\unicode[STIX]{x1D70F}\lesssim 2$ . However, significant discrepancies with our problem are seen for $\unicode[STIX]{x1D70F}\lesssim 1$ and also for the predicted profile shapes of dam-break solutions (not shown), which are due to the different physics of the two problems. Alternatively, the drainage values can be obtained from the numerical solutions of the PDEs in § 4.2. The numerical solutions of the shallow-water equations in figure 14 display a good agreement with DNS results. Furthermore, figure 14 shows that our reduced-order model is in a good agreement with the DNS results during the similarity period and hence it could be used in this context for predicting the discharge rates in practical applications or engineering designs.

In fact, the geometry of the problem we have considered (figures 1(c) and 2) is a common shape for open channel flows. Therefore, the analysis presented here could also be useful for hydraulic engineering and design where edge flows occur. Furthermore, the analysis could help to understand the flow behaviour after a dam break where the released fluid falls from the broken structure or for describing the behaviour of an overflow caused by the excess inflow due to rain or other reasons. This analysis could also be extended to other cross-sectional shapes such as V-shaped or curved geometries (see e.g. Ungarish Reference Ungarish2013b ; Zemach & Ungarish Reference Zemach and Ungarish2013, for non-rectangular cross-section studies of gravity currents) and for ventilation problems where the interface shape and the discharge rates of cold (dense) air and hot (light) air systems in a building are sought.

Figure 14. An example of reservoir drainage as one of the applications of drainage from an edge. When a wall in a tank of a fluid is removed or when a reservoir such as an aquarium breaks, the dynamics of the interface shape, the discharge and the volume removal rate could be described using our results. The ratio of the removed volume, $V_{out}$ , over the initial volume, $V_{0}$ , for three cases of our DNS results are shown in (a). The prediction of the self-similar solutions for the discharge, i.e. $\widetilde{Q}(1)=(\unicode[STIX]{x1D70F}+2)^{-3}\widetilde{H}(1)^{3/2}$ , is plotted in (b), which agrees with the rescaled discharge in the similarity phase. An approximate dam-break solution for the initial slumping phase is also plotted; it works well for $1\lesssim \unicode[STIX]{x1D70F}\lesssim 2$ . The numerical shallow-water PDE solutions obtained in § 4.2, which are plotted in (b), agree well with the DNS results.

6.2 Model assumptions and limitations

In our analysis, we have employed the shallow-water equations, which are derived from the Euler equations. It is assumed that the motion in the ambient fluid has negligible effects on the flow of the dense fluid. Furthermore, the shallow-water equations are vertically averaged and hence two-dimensional effects are neglected, which could be important especially during the early-time period when the heights could be comparable with the length but not during the similarity phase that we focused on here. At sufficiently long times, viscous effects become important and eventually dominate the inertial effects; however, we do not consider such long-time behaviour. We also ignored entrainment and surface tension effects, which are typically small in the kinds of problems studied here.