2.1 Hodograph variables

The shallow-water equations (2.1a ) and (2.1b ) are hyperbolic and may be written in terms of quantities $\unicode[STIX]{x1D6FC}=u+2h^{1/2}$ and $\unicode[STIX]{x1D6FD}=u-2h^{1/2}$ , which are invariant along characteristics (Whitham Reference Whitham1974). Thus we find that

(2.5a,b ) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6FC}}{\text{d}t}=0\quad \text{on }\frac{\text{d}x}{\text{d}t}=u+h^{1/2}\quad \text{and}\quad \frac{\text{d}\unicode[STIX]{x1D6FD}}{\text{d}t}=0\quad \text{on }\frac{\text{d}x}{\text{d}t}=u-h^{1/2}.\end{eqnarray}$$

It is also possible to interchange the dependent and independent variables to treat $x$ and $t$ as functions of the characteristic invariants, $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ . This hodograph transformation converts the nonlinear governing equations into linear equations. In particular as shown by Hogg (Reference Hogg2006), the characteristic equations (2.5) are now given by

(2.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}x}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}=\frac{(3\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD})}{4}\frac{\unicode[STIX]{x2202}t}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}\quad \text{on }\unicode[STIX]{x1D6FC}=\text{const}., & \displaystyle\end{eqnarray}$$

(2.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}x}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}=\frac{(\unicode[STIX]{x1D6FC}+3\unicode[STIX]{x1D6FD})}{4}\frac{\unicode[STIX]{x2202}t}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}\quad \text{on }\unicode[STIX]{x1D6FD}=\text{const}. & \displaystyle\end{eqnarray}$$

$J=(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD})(\unicode[STIX]{x2202}t/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC})(\unicode[STIX]{x2202}t/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD})$

It is possible to find a fundamental solution to the adjoint governing equations and this Riemann solution can then aide the construction of the solution from conditions along the boundaries in the hodograph plane. This technique was used by Hogg and coworkers to solve some gravity current and hydraulic problems using quasi-analytical techniques (Antuono & Hogg Reference Antuono and Hogg2009; Antuono, Hogg & Brocchini Reference Antuono, Hogg and Brocchini2009; Goater & Hogg Reference Goater and Hogg2011; Hogg et al. Reference Hogg, Baldock and Pritchard2011). These analyses deploy the Riemann function, which is given by (Garabedian Reference Garabedian1986)

(2.7) $$\begin{eqnarray}B(a,b;\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=\frac{(a-b)^{3}}{(a-\unicode[STIX]{x1D6FD})^{3/2}(\unicode[STIX]{x1D6FC}-b)^{3/2}}{\mathcal{F}}\left[\frac{3}{2},\frac{3}{2};1;\frac{(a-\unicode[STIX]{x1D6FC})(\unicode[STIX]{x1D6FD}-b)}{(a-\unicode[STIX]{x1D6FD})(\unicode[STIX]{x1D6FC}-b)}\right],\end{eqnarray}$$

where ${\mathcal{F}}$ denotes a hypergeometric function. Integration around a closed curve, $\unicode[STIX]{x2202}{\mathcal{D}}$ , in the hodograph plane yields

(2.8) $$\begin{eqnarray}\int _{\unicode[STIX]{x2202}{\mathcal{D}}}\boldsymbol{f}\boldsymbol{\cdot }\text{d}\boldsymbol{x}=0,\end{eqnarray}$$

where $\boldsymbol{f}(a,b;\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=-V\hat{\boldsymbol{a}}+U\hat{\boldsymbol{b}}$ , $\text{d}\boldsymbol{x}=\text{d}a\hat{\boldsymbol{a}}+\text{d}b\hat{\boldsymbol{b}}$ and

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle U=-\frac{3}{2(a-b)}tB+\frac{B}{2}\frac{\unicode[STIX]{x2202}t}{\unicode[STIX]{x2202}b}-\frac{t}{2}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}b}, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & \displaystyle V=\frac{3}{2(a-b)}tB+\frac{B}{2}\frac{\unicode[STIX]{x2202}t}{\unicode[STIX]{x2202}a}-\frac{t}{2}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}a}. & \displaystyle\end{eqnarray}$$

Thus $a$ and $b$ are dummy variables for the evaluation of integral (2.8) around the closed path $\unicode[STIX]{x2202}{\mathcal{D}}$ in the hodograph plane.

The method of characteristics and the hodograph transformation will be used to construct the solution for this free-draining flow and will be shown to bring considerable insight to the ensuing dynamics, elucidating the nonlinear wave reflections from the back wall $(x=0)$ and the barrier or constriction $(x=1)$ . Throughout this paper we will use the following terminology: complex wave region refers to a domain within which both characteristic variables are varying; simple wave region refers to a domain within which one of the characteristic variables is constant and the other varying; and uniform region refers to a domain within which both characteristic variables are constant.

2.2 Numerical method

We also integrate the governing equations (2.1a ) and (2.1b ) numerically. This is done utilising a central-upwind scheme (Kurganov et al. Reference Kurganov, Noelle and Petrova2001) with a minmod piecewise linear reconstruction (see Kurganov et al. Reference Kurganov, Noelle and Petrova2001 equation (3.16) with $\unicode[STIX]{x1D703}=1$ ). The resultant semi-discrete system is then time stepped using a third-order Runge–Kutta scheme composed of multiple Euler time steps (Shu & Osher Reference Shu and Osher1988). This scheme is positivity preserving, i.e. we always have $h\geqslant 0$ , due to the Courant–Friedrichs–Lewy (CFL) condition and lack of source terms.

To impose boundary conditions we use three ghost cells at each domain end, for which we must find values consistent with both the governing equations and the boundary conditions. We discuss below the process of finding values at the domain end itself, but to impose values in the ghost cells some extrapolation utilising the domain endpoint values and interior values is required. Zeroth-order extrapolation is inconsistent with the first-order reconstruction in the cells, and first-order extrapolation can produce negative depths. To resolve this issue we consider the ghost cells to have size $\unicode[STIX]{x1D716}\unicode[STIX]{x0394}x$ (where $\unicode[STIX]{x0394}x$ is the cell size in the interior) and limit $\unicode[STIX]{x1D716}\rightarrow 0$ . Under this limit all consistent extrapolation procedures converge to enforce the values of $h$ and $uh$ at the domain end across all ghost cells.

To deduce how the domain endpoint values should be evolved we consider the general form of a homogeneous system of $m$ hyperbolic partial differential equations of a variable $Q(x,t)\in \mathbb{R}^{m}$

(2.11) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}Q}{\unicode[STIX]{x2202}t}+A(Q,x,t)\frac{\unicode[STIX]{x2202}Q}{\unicode[STIX]{x2202}x}=\unicode[STIX]{x1D6F9}(Q,x,t).\end{eqnarray}$$

For the purposes of imposing boundary conditions alone, (2.11) is adjusted by including a small amount of diffusion $0<{\mathcal{D}}\ll 1$ , and then multiplying by the $p$ th left eigenvector of $A(Q,x,t)$ , $l^{p}(Q,x,t)$ corresponding to eigenvalue $\unicode[STIX]{x1D706}^{p}(Q,x,t)$ , which results in the characteristic form

(2.12) $$\begin{eqnarray}l^{p}\frac{\unicode[STIX]{x2202}Q}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D706}^{p}l^{p}\frac{\unicode[STIX]{x2202}Q}{\unicode[STIX]{x2202}x}={\mathcal{D}}l^{p}\frac{\unicode[STIX]{x2202}^{2}Q}{\unicode[STIX]{x2202}x^{2}}+l^{p}\unicode[STIX]{x1D6F9}.\end{eqnarray}$$

Diffusion is required for cases when a weak solution to the system contains a discontinuity between the bulk and the boundary, that is when the value at the boundary is not the limit of the interior values. Such a situation is not permitted in the vanishing viscosity solution, and to ensure this solution is found we include a small amount of diffusion. The boundary conditions are of the form $f(Q,t)=0$ , where the function $f:\mathbb{R}^{m}\times \mathbb{R}\rightarrow \mathbb{R}^{k}$ , and specify the values advected by characteristics entering the domain. We focus our discussion on the right-hand end of the domain, that is we are at $x=x_{e}$ such that the domain is $x<x_{e}$ (and in this study $x_{e}=1$ ). This means that we must enforce (2.12) for those $p$ such that $\unicode[STIX]{x1D706}^{p}\geqslant 0$ as well as the boundary conditions to determine $Q_{e}(t)=Q(x_{e},t)$ . We discretise the $x$ derivatives using finite differences over the grid with cell size, $\unicode[STIX]{x0394}x$ , and take ${\mathcal{D}}=\max _{p}|\unicode[STIX]{x1D706}^{p}|\unicode[STIX]{x0394}x$ which is of the same form as the effective diffusion from the method in the bulk and is stable under the same CFL condition. This yields a system of differential algebraic equations (DAE) (see Ascher & Petzold Reference Ascher and Petzold1998; Kunkel & Mehrmann Reference Kunkel and Mehrmann2006) of the form

(2.13a ) $$\begin{eqnarray}\displaystyle & \displaystyle L(Q_{e},t)\frac{\text{d}Q_{e}}{\text{d}t}=C(Q_{e},t), & \displaystyle\end{eqnarray}$$

(2.13b ) $$\begin{eqnarray}\displaystyle & \displaystyle f(Q_{e},t)=0, & \displaystyle\end{eqnarray}$$

$L$

$m-k$

$\unicode[STIX]{x1D706}^{p}\geqslant 0$

$t$

$C$

$m=2$

$k=1$

$(x=0)$

$(x=1)$

$t^{n}$

$Q_{e}^{n}=Q_{e}(t^{n})$

$\unicode[STIX]{x0394}t$

$\unicode[STIX]{x0394}Q_{e}^{n}=Q_{e}^{n+1}-Q_{e}^{n}$

$f=0$

$t^{n+1}$

$\unicode[STIX]{x0394}Q_{e}^{n,0}$

(2.14) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}L(Q_{e}^{n},t^{n})\\ F(Q_{e}^{n}+\unicode[STIX]{x0394}Q_{e}^{n,0},t^{n+1})\end{array}\right]\unicode[STIX]{x0394}Q_{e}^{n}=\left[\begin{array}{@{}c@{}}C(Q_{e}^{n},t^{n})\unicode[STIX]{x0394}t\\ F(Q_{e}^{n}+\unicode[STIX]{x0394}Q_{e}^{n,0},t^{n+1})\unicode[STIX]{x0394}Q_{e}^{n,0}-f(Q_{e}^{n}+\unicode[STIX]{x0394}Q_{e}^{n,0},t^{n+1})\end{array}\right]\!,\qquad\end{eqnarray}$$

where $F=\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}Q$ and is non-vanishing, which is similar to the poststabilisation method with one iteration of projection (Ascher & Petzold Reference Ascher and Petzold1998, § 10.2).

As we are solving an ODE on a manifold (2.13) using Euler’s method, we can show the problem under consideration has sufficient smoothness conditions so that the method has a truncation error of $O(\unicode[STIX]{x0394}t)$ . Similarly, if $\overline{\unicode[STIX]{x0394}Q_{e}}^{n}$ is the limit of iterating (2.14) with the output $\unicode[STIX]{x0394}Q^{n}$ being taken as the next input $\unicode[STIX]{x0394}Q_{e}^{n,0}$ (the implicit version of (2.14)), then we can show that the Newton–Raphson portion of our method will provide $f(Q^{n+1},t^{n+1})=O(\Vert \overline{\unicode[STIX]{x0394}Q_{e}}^{n}-\unicode[STIX]{x0394}Q_{e}^{n,0}\Vert ^{2})$ and then if $\Vert \overline{\unicode[STIX]{x0394}Q_{e}}^{n}-\unicode[STIX]{x0394}Q_{e}^{n,0}\Vert =O(\unicode[STIX]{x0394}t)$ , we deduce that the convergence to the manifold is second order. Finally we can use the fact that, as time evolves, we are repeatedly performing Newton–Raphson steps on very similar problems to improve the scheme further. By taking our guess to be the value of a recent time step, for example $\unicode[STIX]{x0394}Q_{e}^{n,0}=\unicode[STIX]{x0394}Q_{e}^{n-1}$ , then $f(Q^{n+1},t^{n+1})=O(\unicode[STIX]{x0394}t^{4})$ and the scheme provides fourth order convergence to the manifold.

To use this method with the Runge–Kutta time steps from Shu & Osher (Reference Shu and Osher1988) some modifications must be made. These time steps are composed of multiple Euler time steps, which we shall denote by the operator $E$ , and the intermediate results by $Q_{e}^{(i)}$ where $Q_{e}^{(0)}=Q_{e}^{n}$ . The intermediate results are from computations of the form

(2.15) $$\begin{eqnarray}Q_{e}^{(i+1)}=\unicode[STIX]{x1D702}_{i}Q_{e}^{(0)}+(1-\unicode[STIX]{x1D702}_{i})E(Q_{e}^{(i)}),\end{eqnarray}$$

where $0\leqslant \unicode[STIX]{x1D702}_{i}\leqslant 1$ . Based on the results presented in Kunkel & Mehrmann (Reference Kunkel and Mehrmann2006, § 5.2) we argue that it is appropriate to approximately enforce $f(Q_{e}^{(i+1)},t^{(i+1)})=0$ , where $t^{(i)}$ is the time at which $Q_{e}^{(i)}$ is an approximation, rather than $f(E(Q_{e}^{(i)}),t^{(i)}+\unicode[STIX]{x0394}t)=0$ as in (2.13). We therefore use a modified version of our scheme suitable for the application of (2.15), specifically

(2.16a ) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}L(Q_{e}^{(i)},t^{(i)})\\ (1-\unicode[STIX]{x1D702}_{i})F(Q_{e}^{(i)\star },t^{(i+1)})\end{array}\right]\unicode[STIX]{x0394}Q_{e}^{(i)}=\left[\begin{array}{@{}c@{}}C(Q_{e}^{(i)},t^{(i)})\unicode[STIX]{x0394}t\\ (1-\unicode[STIX]{x1D702}_{i})F(Q_{e}^{(i)\star },t^{(i+1)})\unicode[STIX]{x0394}Q_{e}^{(i),0}-f(Q_{e}^{(i)\star },t^{(i+1)})\end{array}\right]\!,\qquad\end{eqnarray}$$

where

(2.16b ) $$\begin{eqnarray}Q_{e}^{(i)\star }=\unicode[STIX]{x1D702}_{i}Q_{e}^{(0)}+(1-\unicode[STIX]{x1D702}_{i})(Q_{e}^{(i)}+\unicode[STIX]{x0394}Q_{e}^{(i),0}),\end{eqnarray}$$

and $E(Q_{e}^{(i)})=Q_{e}^{(i)}+\unicode[STIX]{x0394}Q^{(i)}$ . The initial guesses are taken as $\unicode[STIX]{x0394}Q_{e}^{(i),0}=\unicode[STIX]{x0394}Q_{e}^{(0)}$ , except for when $i=0$ in which case $\unicode[STIX]{x0394}Q_{e}^{(0),0}=\unicode[STIX]{x0394}Q_{e}^{n-1}$ .

The numerical scheme has been validated against a variety of known solutions, but for the purposes of this paper we simply observe the excellent correspondence between the simulations and analytic solutions at early times (see §§ 3.2 and 3.3), and the convergence to the asymptotic wave structure at late times, for the particular problem considered here.

Figure 2. The characteristic plane for the unsteady draining of fluid from a partially confined reservoir with barrier height, $D=0.05$ and width $w=1$ . In this figure we plot and label the uniform, simple and complex wave regions following the initiation of the flow. We also plot the characteristics bounding the wave regions; the suffices $a$ and $b$ denote the $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -characteristics and are plotted with dashed and solid lines, respectively. (These characteristic curves are defined in § 3.) Additionally we label key points in the characteristic plane, which occur at the boundaries of the wave regions.

Figure 3. The hodograph plane for the unsteady draining of fluid from a partially confined reservoir, with barrier height, $D=0.05$ and $w=1$ . In this figure, the solid line corresponds to the free-draining boundary at $x=1$ and the dotted line to the impermeable boundary at $x=0$ . Uniform regions are depicted by points in the hodograph planes and are labelled $U_{i}$ $(i=1-4)$ , while simple wave regions correspond to dashed straight line segments aligned with the axes and are labelled $S_{i}$ $(i=1-3)$ . Finally, complex wave regions are labelled $C_{i}\,(i=1-3)$ . This figure shows the wave pattern that emerges following flow initiation and the first three reflections from the end wall and free-draining boundary. This pattern of successive reflections will continue indefinitely, approaching the end state in which the fluid is quiescent and of depth equal to the barrier height; this is where the two boundary curves intersect at $\unicode[STIX]{x1D6FC}=-\unicode[STIX]{x1D6FD}=2\sqrt{D}$ .

3.1 Initial motion: $t\leqslant 1$

The drainage condition (2.3) yields Froude numbers in the range $0\leqslant Fr<1$ (provided $D>0$ ) and so at that boundary $(x=1)$ there is one outflowing characteristic direction and one inflowing. This means that it is straightforward to construct the initial solution (valid for $t\leqslant 1$ ) since the characteristic variable $\unicode[STIX]{x1D6FC}$ is determined by the initial conditions. Thus throughout the portion of the domain associated with the initial conditions, we find that $\unicode[STIX]{x1D6FC}=2$ and within this portion there is a rarefaction centred on $x=1$ , which we term simple wave region $S_{1}$ . Within $S_{1}$ , $-2\leqslant \unicode[STIX]{x1D6FD}\leqslant \unicode[STIX]{x1D6FD}_{\ast }$ , where $\unicode[STIX]{x1D6FD}_{\ast }$ is a function of $D$ and $w$ and is the maximum value of $\unicode[STIX]{x1D6FD}$ in this region, determined from (2.3) with $\unicode[STIX]{x1D6FC}=2$ . (Since $0\leqslant Fr<1$ , we can show that $-2\leqslant \unicode[STIX]{x1D6FD}_{\ast }<-2/3$ .) The solution for the velocity and height fields within $S_{1}$ are given by

(3.1a,b ) $$\begin{eqnarray}u=\frac{2}{3}\left(1+\frac{(x-1)}{t}\right)\quad \text{and}\quad h=\frac{1}{9}\left(2-\frac{(x-1)}{t}\right)^{2},\end{eqnarray}$$

for $x_{b1}\equiv 1-t<x<x_{b2}\equiv 1+(2+3\unicode[STIX]{x1D6FD}_{\ast })t/4$ (see figure 2). Upstream of this rarefaction, the fluid remains undisturbed from its initial state with $h=1$ and $u=0$ ; we denoted this uniform region as $U_{1}$ . Downstream of the rarefaction, there is a uniform region attached to the right boundary $(x=1)$ , within which $u=1+\unicode[STIX]{x1D6FD}_{\ast }/2$ and $h=(2-\unicode[STIX]{x1D6FD}_{\ast })^{2}/16$ ; we denote this as region $U_{2}$ . Example profiles of height and velocity are plotted in figure 5. We observe that the fluid slumps gravitationally, progressively mobilising fluid within the reservoir and that its motion is ‘choked’ by draining at the outflow. These solutions for the flow are valid until the rearward propagating characteristic reaches the back wall of the reservoir, which occurs at $t=1$ . Thereafter there are reflections from this wall and the motion becomes more complicated.

3.2 Motion for $t\geqslant 1$

When $t\geqslant 1$ , the structure of the solution for the height and velocity fields is altered due to reflection first from the back wall of the reservoir and then from the free-draining boundary at the front. It is useful to treat this dynamics in the characteristic plane (figure 2) and in the hodograph plane (figure 3), both of which reveal the successive wave reflections. As described above, the initial phase of the motion involves a simple wave region, $S_{1}$ , sandwiched between two uniform regions, $U_{1}$ and $U_{2}$ , which are attached to the boundaries of the reservoir. In the hodograph plane, a simple wave region corresponds to a line segment parallel to the coordinate axes and $S_{1}$ is given by $\unicode[STIX]{x1D6FC}=2$ and $-2\leqslant \unicode[STIX]{x1D6FD}\leqslant \unicode[STIX]{x1D6FD}_{\ast }$ . Uniform regions correspond to points in the hodograph plane; $U_{1}$ is given by $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=(2,-2)$ and $U_{2}$ by $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=(2,\unicode[STIX]{x1D6FD}_{\ast })$ . Complex wave regions correspond to connected domains in the hodograph plane (see figure 3).

The reflection from the back wall generates a complex wave region $C_{1}$ , which is bounded by the characteristic curve on which $\unicode[STIX]{x1D6FC}=2$ . We denote that curve parametrically by $x_{a1}(\unicode[STIX]{x1D6FD})$ and $t_{a1}(\unicode[STIX]{x1D6FD})$ and it is straightforward to show that

(3.2a,b ) $$\begin{eqnarray}t_{a1}=\frac{8}{(2-\unicode[STIX]{x1D6FD})^{3/2}}\quad \text{and}\quad x_{a1}=1+\frac{2(2+3\unicode[STIX]{x1D6FD})}{(2-\unicode[STIX]{x1D6FD})^{3/2}}.\end{eqnarray}$$

The region, $C_{1}$ , is also bounded by the characteristic on which $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{\ast }$ , denoted by $x_{b2}(t)$ , which occurs at the leading edge of the rarefaction fan in the simple wave region, $S_{1}$ (see figure 2). It is given from (3.1) by $x_{b2}=1+(2+3\unicode[STIX]{x1D6FD}_{\ast })t/4$ . The curve $x_{a1}$ intersects $x_{b2}$ at the point $P_{2}$ with coordinates, $(x,t)=(1+2(2+3\unicode[STIX]{x1D6FD}_{\ast })/(2-\unicode[STIX]{x1D6FD}_{\ast })^{3/2},8/(2-\unicode[STIX]{x1D6FD}_{\ast })^{3/2})$ . Thereafter $x_{a1}$ follows a straight line that bounds the uniform region, $U_{2}$ until reaching the outflow at $x=1$ . In this portion the characteristics are given by

(3.3) $$\begin{eqnarray}\frac{x-x_{a1}(\unicode[STIX]{x1D6FD}_{\ast })}{t-t_{a1}(\unicode[STIX]{x1D6FD}_{\ast })}=\frac{6+\unicode[STIX]{x1D6FD}_{\ast }}{4},\end{eqnarray}$$

and the point, $P_{3}$ , at the outflow is $(1,16/[(6+\unicode[STIX]{x1D6FD}_{\ast })(2-\unicode[STIX]{x1D6FD}_{\ast })^{1/2}])$ . The complex wave region, $C_{1}$ , may be calculated using the Riemann function, integrated around a trajectory in the hodograph plane, as demonstrated by Hogg (Reference Hogg2006); it is compactly given by

(3.4) $$\begin{eqnarray}t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=B(2,-2;\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}),\end{eqnarray}$$

and the solution for $x$ is given by integrating along $\unicode[STIX]{x1D6FC}$ -characteristics from the back wall, using (2.6a ) to show

(3.5) $$\begin{eqnarray}x(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=\frac{(3\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD})}{4}t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})-\frac{\unicode[STIX]{x1D6FC}}{2}t(\unicode[STIX]{x1D6FC},-\unicode[STIX]{x1D6FC})-\frac{1}{4}\int _{-\unicode[STIX]{x1D6FC}}^{\unicode[STIX]{x1D6FD}}t(\unicode[STIX]{x1D6FC},b)\,\text{d}b.\end{eqnarray}$$

In the characteristic plane, the continuation of the $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{\ast }$ characteristic beyond $P_{2}$ bounds the region $C_{1}$ (see figure 2). It meets the back wall at $P_{4}=(0,t_{b2}(-\unicode[STIX]{x1D6FD}_{\ast }))$ , where

(3.6) $$\begin{eqnarray}t_{b2}(\unicode[STIX]{x1D6FC})=B(2,-2;\unicode[STIX]{x1D6FC},-\unicode[STIX]{x1D6FD}_{\ast }).\end{eqnarray}$$

The $\unicode[STIX]{x1D6FD}$ -characteristics starting from the uniform region, $U_{2}$ generate the simple wave region, $S_{2}$ , within which $-\unicode[STIX]{x1D6FD}_{\ast }\leqslant \unicode[STIX]{x1D6FC}\leqslant 2$ (see figures 2 and 3). The trajectory of the characteristic from point $P_{3}$ on which $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{\ast }$ is denoted by $(x_{b3},t_{b3})$ and satisfies

(3.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}x_{b3}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}=\frac{(\unicode[STIX]{x1D6FC}+3\unicode[STIX]{x1D6FD}_{\ast })}{4}\frac{\unicode[STIX]{x2202}t_{b3}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}.\end{eqnarray}$$

However on $\unicode[STIX]{x1D6FC}$ -characteristics within $S_{2}$ , we also find that

(3.8) $$\begin{eqnarray}\frac{x-x_{b2}}{t-t_{b2}}=\frac{3\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD}_{\ast }}{4}.\end{eqnarray}$$

Thus we may establish

(3.9) $$\begin{eqnarray}t_{b3}-t_{b2}=\frac{-8(2+3\unicode[STIX]{x1D6FD}_{\ast })}{(6+\unicode[STIX]{x1D6FD}_{\ast })(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FD}_{\ast })^{3/2}},\end{eqnarray}$$

and then $x_{b3}-x_{b2}$ may be deduced from (3.8). This $\unicode[STIX]{x1D6FD}$ -characteristic, $(x_{b3},t_{b3})$ , forms one of the boundaries to the complex wave region at the front, $C_{2}$ , and to the uniform region, $U_{3}$ , adjacent to the back wall. In particular the point $P_{5}$ is determined by the first location along this characteristic on which the velocity vanishes, which is given by $\unicode[STIX]{x1D6FC}=-\unicode[STIX]{x1D6FD}=-\unicode[STIX]{x1D6FD}_{\ast }$ . Beyond $P_{5}$ , the $\unicode[STIX]{x1D6FD}$ -characteristic is a straight line given by

(3.10) $$\begin{eqnarray}\frac{x-x_{b3}(-\unicode[STIX]{x1D6FD}_{\ast })}{t-t_{b3}(-\unicode[STIX]{x1D6FD}_{\ast })}=-\frac{\unicode[STIX]{x1D6FD}_{\ast }}{2}.\end{eqnarray}$$

This straight segment intersects the back wall at point $P_{7}= (\!0,t_{b3}(-\unicode[STIX]{x1D6FD}_{\ast })+2x_{b3}$ $(-\unicode[STIX]{x1D6FD}_{\ast })/\unicode[STIX]{x1D6FD}_{\ast }\!)$ .

Constructing the solution within the complex wave region, $C_{2}$ entails the use of the Riemann function (2.7) and the conversion of the boundary conditions into hodograph variables. In particular the free-draining condition, given by (2.3), may be written as

(3.11) $$\begin{eqnarray}\frac{\text{d}x}{\text{d}\unicode[STIX]{x1D6FC}}=0\quad \text{on }\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D6FC}),\end{eqnarray}$$

We note that given a barrier height, $D$ , that $\unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D6FC})$ determines a curve in the hodograph plane such that $\unicode[STIX]{x1D6FD}_{s}(2)=\unicode[STIX]{x1D6FD}_{\ast }$ and that $\unicode[STIX]{x1D6FD}(2\sqrt{D})=-2\sqrt{D}$ , the latter condition corresponding to no flow when the fluid depth matches the height of the barrier. (This function, $\unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D6FC})$ , is plotted as a solid curve in figure 3 for $D=0.05$ .) From (2.3), we also find that

(3.12) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6FD}_{s}}{\text{d}\unicode[STIX]{x1D6FC}}=\frac{(Fr^{4/3}-1/w^{2/3})(Fr-2)-2DFr^{1/3}/h}{(Fr^{4/3}-1/w^{2/3})(Fr+2)-2DFr^{1/3}/h},\end{eqnarray}$$

and this implies that $\text{d}\unicode[STIX]{x1D6FD}_{s}/\text{d}\unicode[STIX]{x1D6FC}=-1$ when $Fr=0$ . Thus the curve $\unicode[STIX]{x1D6FD}_{s}$ approaches its end point at $\unicode[STIX]{x1D6FD}_{s}=-2\sqrt{D}$ tangentially to the line $\unicode[STIX]{x1D6FD}=-\unicode[STIX]{x1D6FC}$ (see figure 3). (The case without a barrier is different: when $D=0$ , the outflow Froude number, $Fr$ , does not vary with the flow depth, $h$ and so the boundary curve, $\unicode[STIX]{x1D6FD}_{s}$ , is a straight line of gradient $\text{d}\unicode[STIX]{x1D6FD}_{s}/\text{d}\unicode[STIX]{x1D6FC}=(Fr-2)/(Fr+2)$ .) Using the characteristic equations (2.6a ) and (2.6b ), we can show that the boundary condition (3.11) is given by

(3.13) $$\begin{eqnarray}2(\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD}_{s})\left(\frac{\unicode[STIX]{x2202}t}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}+\frac{\text{d}\unicode[STIX]{x1D6FD}_{s}}{\text{d}\unicode[STIX]{x1D6FC}}\frac{\unicode[STIX]{x2202}t}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}\right)+(\unicode[STIX]{x1D6FD}_{s}-\unicode[STIX]{x1D6FC})\left(\frac{\unicode[STIX]{x2202}t}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}}-\frac{\text{d}\unicode[STIX]{x1D6FD}_{s}}{\text{d}\unicode[STIX]{x1D6FC}}\frac{\unicode[STIX]{x2202}t}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}\right)=0\quad \text{on }\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D6FC}).\end{eqnarray}$$

The other perimeter data in the hodograph plane are supplied by $t=t_{b3}(\unicode[STIX]{x1D6FC})$ along $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{\ast }$ . We then integrate around the closed curve $OQRT$ , as depicted in figure 4. The coordinates of these points in the hodograph plane are $O$ , $(2,\unicode[STIX]{x1D6FD}_{\ast })$ ; $Q$ , $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{\ast })$ ; $R$ , $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s})$ ; and $T$ , $(\unicode[STIX]{x1D6FC}_{s},\unicode[STIX]{x1D6FD}_{s})$ . The curves $OQ$ , $QR$ and $RT$ correspond to straight line segments in the hodograph plane aligned with one of the coordinate axes, and the curve $RT$ corresponds to the free-draining boundary condition ( $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{s}$ and $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{s}$ ). Integration of (2.8) around $OQRT$ yields

(3.14) $$\begin{eqnarray}\displaystyle & & \displaystyle t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s})=B(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{\ast };\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s})t_{b3}(\unicode[STIX]{x1D6FC})+B(2,\unicode[STIX]{x1D6FD}_{\ast };\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s})t_{b3}(2)\left(\frac{2+\unicode[STIX]{x1D6FD}_{\ast }}{2-\unicode[STIX]{x1D6FD}_{\ast }}-\frac{1}{2}\right)\nonumber\\ \displaystyle & & \displaystyle \quad +\,B(\unicode[STIX]{x1D6FC}_{s},\unicode[STIX]{x1D6FD}_{s};\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s})t(\unicode[STIX]{x1D6FC}_{s},\unicode[STIX]{x1D6FD}_{s})\left(\frac{1}{2}-\frac{(\unicode[STIX]{x1D6FC}_{s}+\unicode[STIX]{x1D6FD}_{s})}{\unicode[STIX]{x1D6FC}_{s}-\unicode[STIX]{x1D6FD}_{s}}\right)-\int _{2}^{\unicode[STIX]{x1D6FC}}\left(\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}a}-\frac{3B}{2(a-\unicode[STIX]{x1D6FD}_{\ast })}\right)t_{b3}(a)\,\text{d}a\nonumber\\ \displaystyle & & \displaystyle \quad -\,\int _{\unicode[STIX]{x1D6FC}_{s}}^{2}\frac{t(a_{s},b_{s})}{2(a_{s}-b_{s})}\left((3a_{s}+b_{s})\left(\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}a}-\frac{B}{a_{s}-b_{s}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.(a_{s}+3b_{s})\frac{\text{d}b_{s}}{\text{d}a}\left(\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}b}+\frac{B}{a_{s}-b_{s}}\right)\right)\,\text{d}a.\end{eqnarray}$$

In the first integral of (3.14), the Riemann function and its derivative are evaluated at $b=\unicode[STIX]{x1D6FD}_{\ast }$ , whereas in the final integral of (3.14), the Riemann function and its derivatives are evaluated along the curve $TO$ given parametrically by $a=a_{s}$ and $b=b_{s}$ . Before we can use this expression to evaluate $t$ at a general position within the complex wave region, $C_{2}$ , we must first evaluate the time field along the boundary curve $TO$ ( $t_{s}(\unicode[STIX]{x1D6FC})\equiv t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D6FC}))$ ). Setting $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{s}$ in (3.14), we derive an integral equation for $t_{s}$ of the form

(3.15) $$\begin{eqnarray}t_{s}(\unicode[STIX]{x1D6FC})=f_{1}(\unicode[STIX]{x1D6FC})+\int _{2}^{\unicode[STIX]{x1D6FC}}f_{2}(a,\unicode[STIX]{x1D6FC})t_{s}(a)\,\text{d}a,\end{eqnarray}$$

where $f_{1}$ and $f_{2}$ are functions determined from (3.14). This is a Volterra integral equation of the second kind and is readily solved by iteration; we find numerical convergence to a relative tolerance of $10^{-12}$ within 10 iterates. Example solutions are plotted in figures 5 and 6.

Figure 4. Sections of the hodograph plane showing: (a) the closed curve $OQRT$ , which is used to construct the solution in the complex wave region, $C_{2}$ ; and (b) the rectangle $VWXY$ , which is used to construct the solution in $C_{3}$ .

Figure 5. The height, $h(x,t)$ and velocity, $u(x,t)$ , as functions of position at various instances of time, $t=1,2,3,4$ and $5$ (labelled (i)–(v), respectively) for $D=0.05$ and $w=1$ . The solid line depicts the solution calculated using the quasi-analytical hodograph methods, while the dot-dashed line depicts the solution from the direct numerical integration of the shallow-water equations (2.1a )–(2.1b ). We note that there is negligible difference between the curves, which overlie each other; the largest differences occur close to positions in the profiles where the solutions have discontinuous gradient. Such features are captured exactly by the analytical approach, whereas the numerical results are slightly smoothed. The dashed line indicates the height of the barrier $(D=0.05)$ .

Figure 6. The height at the free-draining boundary, $h(1,t)$ , as a function of time for (a) various barrier heights, $D$ , and $w=1$ ; and (b) various constrictions, $w$ , and $D=0.1$ . The solid lines correspond to results from integration in the hodograph plane and the dot-dashed lines correspond to numerical integration of the shallow water equations (2.1a )–(2.1b ). The quasi-analytical hodograph calculations were performed until the end of region $U_{4}$ , but the numerical computations were continued over much longer times. We note the very close agreement between the two approaches, with the only very slight differences occurring at the junctions of the intervals during which the height is uniform or decreasing.

Given $t_{s}(\unicode[STIX]{x1D6FC})$ , the dependent variables can be evaluated immediately by numerical quadrature, using (3.14) for $t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ and integrating along an $\unicode[STIX]{x1D6FC}$ -characteristic using (2.6a ) to determine $x(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ , with condition $x=x_{b3}$ . This enables the boundary to $C_{2}$ to be computed as the $\unicode[STIX]{x1D6FC}$ -characteristic on which $\unicode[STIX]{x1D6FC}=-\unicode[STIX]{x1D6FD}_{\ast }$ , starting from $P_{5}$ . The $\unicode[STIX]{x1D6FC}$ -characteristic reaches the free-draining boundary at $P_{6}=(1,t_{s}(-\unicode[STIX]{x1D6FD}_{\ast }))$ .

Finally we note that we also construct the characteristic on which $\unicode[STIX]{x1D6FC}=-\unicode[STIX]{x1D6FD}_{\ast }$ emanating from the point $P_{7}$ at the back wall. This curve denoted parametrically by $(x_{a3},t_{a3})$ is determined from the characteristic equation (2.6a ), which here yields

(3.16) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}x_{a3}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}=\frac{(-3\unicode[STIX]{x1D6FD}_{\ast }+\unicode[STIX]{x1D6FD})}{4}\frac{\unicode[STIX]{x2202}t_{a3}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FD}}.\end{eqnarray}$$

Also in the simple wave region, $S_{3}$ , in which $\unicode[STIX]{x1D6FC}=-\unicode[STIX]{x1D6FD}_{\ast }$ , we have that

(3.17) $$\begin{eqnarray}\frac{x-x_{a2}}{t-t_{a2}}=\frac{-\unicode[STIX]{x1D6FD}_{\ast }+3\unicode[STIX]{x1D6FD}}{4}.\end{eqnarray}$$

Hence we deduce that

(3.18) $$\begin{eqnarray}t_{a3}-t_{a2}=-\frac{2}{\unicode[STIX]{x1D6FD}_{\ast }}x_{b3}(-\unicode[STIX]{x1D6FD}_{\ast })\left(\frac{-2\unicode[STIX]{x1D6FD}_{\ast }}{-\unicode[STIX]{x1D6FD}_{\ast }-\unicode[STIX]{x1D6FD}}\right)^{3/2},\end{eqnarray}$$

and then $x_{a3}$ follows from (3.17). The uniform region adjacent to the free-draining boundary corresponds to $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=(-\unicode[STIX]{x1D6FD}_{\ast },\unicode[STIX]{x1D6FD}_{m})$ , where $\unicode[STIX]{x1D6FD}_{m}=\unicode[STIX]{x1D6FD}_{s}(-\unicode[STIX]{x1D6FD}_{\ast })$ . Then we may determine $P_{8}$ as $(x_{a3}(\unicode[STIX]{x1D6FD}_{m}),t_{a3}(\unicode[STIX]{x1D6FD}_{m}))$ and $P_{9}$ as $(1,t_{a3}(\unicode[STIX]{x1D6FD}_{m})+4(1-x_{a3}(\unicode[STIX]{x1D6FD}_{m}))/(\unicode[STIX]{x1D6FD}_{m}-3\unicode[STIX]{x1D6FD}_{\ast }))$ .

The construction of the solution may be extended to the complex wave region $C_{3}$ using the boundary data along the $\unicode[STIX]{x1D6FC}$ - characteristic on which $\unicode[STIX]{x1D6FC}=-\unicode[STIX]{x1D6FD}_{\ast }$ ; these data are $x_{a3}(\unicode[STIX]{x1D6FD})$ and $t=t_{a3}(\unicode[STIX]{x1D6FD})$ , and are derived above. In addition we enforce the impermeability of the back wall (2.2) by imposing symmetry on the $t$ -variable, so that $t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=t(-\unicode[STIX]{x1D6FD},-\unicode[STIX]{x1D6FC})$ . We then integrate around a rectangle in the hodograph plane, $VWXY$ , with edges aligned with the coordinate axes (see figure 4). The vertices of the rectangle are at $V=(-\unicode[STIX]{x1D6FD}_{\ast },\unicode[STIX]{x1D6FD}_{\ast })$ , $W=(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{\ast })$ , $X=(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ and $Y=(-\unicode[STIX]{x1D6FD}_{\ast },\unicode[STIX]{x1D6FD})$ . The characteristic data $t=t_{a3}(\unicode[STIX]{x1D6FD})$ provide the solution on the boundary segment $YV$ , while the imposed symmetry determines the solution on the boundary segment $VW$ , $t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{\ast })=t(-\unicode[STIX]{x1D6FD}_{\ast },-\unicode[STIX]{x1D6FC})=t_{a3}(-\unicode[STIX]{x1D6FC})$ . We may then integrate (2.8) around $VWXY$ to yield

(3.19) $$\begin{eqnarray}\displaystyle t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}) & = & \displaystyle B(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{\ast };\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})t_{a3}(-\unicode[STIX]{x1D6FC})+B(-\unicode[STIX]{x1D6FD}_{3},\unicode[STIX]{x1D6FD};\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})t_{a3}(\unicode[STIX]{x1D6FD})\nonumber\\ \displaystyle & & \displaystyle -\,B(-\unicode[STIX]{x1D6FD}_{\ast },\unicode[STIX]{x1D6FD}_{\ast };\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})t_{a3}(\unicode[STIX]{x1D6FD}_{\ast })-I_{1}-I_{2},\end{eqnarray}$$

where the integrals $I_{1}$ and $I_{2}$ are given by

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle I_{1}=\int _{-\unicode[STIX]{x1D6FD}_{\ast }}^{\unicode[STIX]{x1D6FC}}\left(-\frac{3B}{2(a-\unicode[STIX]{x1D6FD}_{\ast })}+\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}a}\right)t_{a3}(-a)\,\text{d}a, & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle I_{2}=\int _{\unicode[STIX]{x1D6FD}}^{\unicode[STIX]{x1D6FD}_{\ast }}\left(-\frac{3B}{2(-\unicode[STIX]{x1D6FD}_{\ast }-b)}-\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}b}\right)t_{a3}(b)\,\text{d}b. & \displaystyle\end{eqnarray}$$

In the integrand of $I_{1}$ , $b=\unicode[STIX]{x1D6FD}_{\ast }$ , whereas in the integrand of $I_{2}$ , $a=-\unicode[STIX]{x1D6FD}_{\ast }$ . This construction is essentially the same as was used to derive (3.4) (Hogg Reference Hogg2006). The evaluation of $x(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ follows immediately by integrating along $\unicode[STIX]{x1D6FC}$ -characteristics (2.6a ) and using the boundary condition that $x(\unicode[STIX]{x1D6FC},-\unicode[STIX]{x1D6FC})=0$ .

The description of the solution presented in the previous paragraphs completes the computation of the flow fields for the first three reflections; subsequent reflections could in principle be computed using an identical formulation to find the complex wave regions next to the free-draining boundary and then to calculate the solutions close to the back wall. In the hodograph plane we can see that this sequence of regions forms an infinite ‘staircase’ towards the final state in which the fluid is quiescent of a flow depth equal to the height of the barrier. Notably, each reflection will generate uniform regions adjacent to the boundaries.

We illustrate the solution by plotting the solutions for the height and velocity fields as function of position at various fixed times for $D=0.05$ (see figure 5). The profiles, plotted at times $t=1,2,3,4$ and $5$ , draw from the uniform, simple and complex wave regions illustrated in figure 2. We note that they feature the uniform regions adjacent to the boundaries and that they illustrate the wave reflections along the reservoir. For example, we see that the flow oscillates so that at times it is thickest close to the end wall, while at other times it is thickest close to the free-draining boundary. We also plot the height at the draining boundary, $h(1,t)$ , as a function of time for various barrier heights and a constant width, $w=1$ , (figure 6 a) and for various constrictions and a constant barrier height, $D=0.1$ (figure 6 b). We note that this height at the end of the channel has intervals within which it is constant, followed by periods within which it declines; this is in accord with the description that we have developed above. The successive wave reflections lead to uniform regions adjacent to the boundaries within which the height and velocity are constant; in figure 6, we observe the height variation within regions, $U_{1}$ , $C_{2}$ and $U_{4}$ (see figure 2) and we have calculated these using the techniques described above. In principle we could extend this further, but in this study we access later times numerically and through using a different analytical approach (§ 3.3).

3.3 Late times: numerical results and asymptotic analysis

From the structure of the hodograph solution we anticipate that at late times we will continue to see a ‘staircase’ of uniform and complex wave regions. This is consistent with our numerical simulations, computed up to time $t=100$ , and presented in figure 7. In this figure we plot the height excess over the barrier and the velocity incident on the barrier as functions of time for a selection of barrier heights. (We have introduced a time offset, $t_{0}$ and have also plotted the asymptotic trends, which are introduced and deduced below.) It is evident from these plots that the solution tends rapidly towards a power law, around which there are strong, but decaying, oscillations, whose period and shape are functions of the barrier height, $D$ .

In what follows we investigate the asymptotic form of the solutions for the height and velocity fields trends at late times after the flow initiation. For this analysis we utilise a multi-scale asymptotic expansion in time to capture the progressive unsteady draining of the reservoir and the more rapid oscillations, the amplitude of which decay temporally and the phase speed of which approach a constant value dependent upon the barrier height. This analysis is based upon embedding a set of time scales into the dependent variables and then examining their evolution at relatively late times. We first construct these appropriate time scales, then present the resulting expansion with the algebraic details confined to appendix B. Finally we compare the asymptotic forms to the numerical simulations.

Figure 7. (a) The elevation of the fluid layer above the barrier relative to the height of the barrier, $h(1,t)/D-1$ ; and (b) and the velocity at the outflow $u(1,t)/\sqrt{D}$ as functions of time for $w=1$ and $D=0.9,0.7,0.5,0.3,0.1$ (denoted (i)–(v), respectively). The solid lines are from numerical integration of the governing equations and reveal strong oscillatory behaviour superimposed upon progressive decay of each variable. The dashed-dot lines are the single-scale asymptotic solution from (3.26a,b ).

We investigate the solution at late times, where $t=O(1/\unicode[STIX]{x1D716})$ and $\unicode[STIX]{x1D716}$ is a small, arbitrary ordering parameter $(0<\unicode[STIX]{x1D716}\ll 1)$ . In the end, the solutions cannot depend upon $\unicode[STIX]{x1D716}$ , but its introduction is useful in organising the magnitude of various terms. We also introduce the following set of time scales, $\hat{T}_{n}(\unicode[STIX]{x1D716},t)$ $(n=0,1,2,\ldots )$ , such that at the long time scale of interest $\hat{T}_{n}=O(\unicode[STIX]{x1D716}^{n-1})$ . Provided that all variables other than $t$ and $\hat{T}_{n}$ are independent of $\unicode[STIX]{x1D716}$ we can always express the time scales as

(3.22) $$\begin{eqnarray}\hat{T}_{n}(\unicode[STIX]{x1D716},t)=\unicode[STIX]{x1D716}^{n-1}\unicode[STIX]{x1D714}_{n}(\unicode[STIX]{x1D716}t),\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{n}$ are functions of order unity when $t=O(1/\unicode[STIX]{x1D716})$ . For the purposes of constructing an asymptotic expansion the time scale $\hat{T}_{1}=\unicode[STIX]{x1D714}_{1}(\unicode[STIX]{x1D716}t)$ is equivalent to $\tilde{T}_{1}=\unicode[STIX]{x1D716}t$ . Therefore, without loss of generality, we take

(3.23a-c ) $$\begin{eqnarray}\hat{T}_{0}=\unicode[STIX]{x1D716}^{-1}\unicode[STIX]{x1D714}_{0}(T_{1}),\quad \hat{T}_{1}=\unicode[STIX]{x1D716}t,\quad \hat{T}_{n}=\unicode[STIX]{x1D716}^{n-1}\unicode[STIX]{x1D714}_{n}(T_{1}),\quad n\in \{2,3,\ldots \},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{n}^{\prime }\neq 0$ , $n\in \{0,2,3,\ldots \}$ .

For most problems multi-scale asymptotics fails to produce a closed system of equations for high-order corrections involving high-order time scales, as mentioned in Bender & Orszag (Reference Bender and Orszag1999, § 11.2). To understand why this is, we consider a problem in $m$ time scales which, at each order, generates $r$ new arbitrary functions dependent on $m-1$ time scales (the dependence on one of the time scales being set at that order), along with $s$ solvability conditions providing restrictions on the functions introduced at lower order. Considering an expansion including $N$ terms, and ignoring contributions independent of $N$ , this means there is a total of $r(m-1)N$ functional dependencies and $sN$ solvability conditions, resulting in $[r(m-1)-s]N$ undetermined dependencies. To stress how pathologically bad this lack of closure is, we observe that if a full set of $O(N)$ time scales is used then we will have $O(N^{2})$ undetermined dependencies. If we are to have a closed expansion, introducing the same number of dependencies and conditions as we increase from order $N$ to $N+1$ , we must have $m=1+s/r$ .

For many problems, such as the one considered here, $s=a$ and thus we require two time scales. We anticipate that, as with many of the problems discussed in Kevorkian & Cole (Reference Kevorkian and Cole1996), the time scale $\hat{T}_{0}$ is the time scale in which oscillations and waves occur, $\hat{T}_{1}$ governs the amplitude of these waves, and higher-order time scales act as a phase shift to $\hat{T}_{0}$ in the form of a strained coordinate (Hinch Reference Hinch1992, § 7.1). This results in the two-scale ansatz: that all time dependencies can be expressed by the time scales

(3.24a ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{0}=k_{0}\hat{T}_{0}+k_{2}\hat{T}_{2}+k_{3}\hat{T}_{3}+\cdots =\unicode[STIX]{x1D716}^{-1}\unicode[STIX]{x1D714}_{0}(T_{1})+\unicode[STIX]{x1D716}\unicode[STIX]{x1D714}_{2}(T_{1})+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D714}_{3}(T_{1})+\cdots \,,\qquad & \displaystyle\end{eqnarray}$$

(3.24b ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{1}=\hat{T}_{1}=\unicode[STIX]{x1D716}t, & \displaystyle\end{eqnarray}$$

$k_{n}$

$\unicode[STIX]{x1D714}_{n}$

Our asymptotic analysis then proceeds by assuming that the height and velocity fields are dependent on both time scales $T_{0}$ and $T_{1}$ as well as the spatial variable  $x$ . We substitute series expansion of the form (B 9b ) into the governing equations (2.1a )–(2.1b ) and deduce expressions that satisfy the boundary conditions and remain bounded as time evolves. These latter ‘secular’ conditions are key to determining the functions, $\unicode[STIX]{x1D714}_{n}$ , and the decaying amplitudes of the wave-like behaviour.

It is convenient to write the solutions for the dependent fields as the sum of two contributions,

(3.25a,b ) $$\begin{eqnarray}h=h_{a}+h_{w}\quad \text{and}\quad u=u_{a}+u_{w},\end{eqnarray}$$

where $h_{a}$ and $u_{a}$ are asymptotic components of the solutions that depend on $T_{1}$ (and will be henceforth termed the ‘single time-scale asymptotics’) and $h_{w}$ and $u_{w}$ are the wave-like components. From appendix B, we find that correct to $O(1/(\sqrt{D}(t-t_{0}))^{5})$

(3.26a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{h_{a}(x,t)}{D}=1+\frac{27}{2w^{2}(\sqrt{D}(t-t_{0}))^{2}}+\frac{27(3x^{2}+5)}{2w^{2}(\sqrt{D}(t-t_{0}))^{4}}, & \displaystyle\end{eqnarray}$$

(3.26b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{a}(x,t)}{\sqrt{D}}=\frac{27x}{w^{2}(\sqrt{D}(t-t_{0}))^{3}}, & \displaystyle\end{eqnarray}$$

(3.26c ) $$\begin{eqnarray}\displaystyle \frac{h_{w}(x,t)}{D} & = & \displaystyle \frac{1}{(\sqrt{D}(t-t_{0}))^{3}}[F(\unicode[STIX]{x1D702})+F(\unicode[STIX]{x1D709})]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{(\sqrt{D}(t-t_{0}))^{4}}\left[\frac{w^{2}}{9}(F(\unicode[STIX]{x1D702})^{2}+F(\unicode[STIX]{x1D709})^{2})+3x(F(\unicode[STIX]{x1D702})-F(\unicode[STIX]{x1D709}))\right],\qquad\end{eqnarray}$$

(3.26d ) $$\begin{eqnarray}\displaystyle \frac{u_{w}(x,t)}{\sqrt{D}} & = & \displaystyle \frac{1}{(\sqrt{D}(t-t_{0}))^{3}}[F(\unicode[STIX]{x1D702})-F(\unicode[STIX]{x1D709})]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{(\sqrt{D}(t-t_{0}))^{4}}\left[\frac{w^{2}}{9}(F(\unicode[STIX]{x1D702})^{2}-F(\unicode[STIX]{x1D709})^{2})+3x(F(\unicode[STIX]{x1D702})+F(\unicode[STIX]{x1D709}))\right],\qquad\end{eqnarray}$$

(3.27a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D709}(x,t)=\sqrt{D}t-\frac{27}{4w^{2}\sqrt{D}(t-t_{0})}+x,\quad \text{and}\quad \unicode[STIX]{x1D702}(x,t)=\sqrt{D}t-\frac{27}{4w^{2}\sqrt{D}(t-t_{0})}-x.\quad\end{eqnarray}$$

The variable $t_{0}$ denotes a temporal offset that we are free to choose. The arbitrary function $F$ is related to the initial conditions that gave rise to the unsteady draining flow. In terms of the variables in the appendix B, it is given by $F(y)=D\hat{F}_{3}(y-\sqrt{D}\unicode[STIX]{x1D719})$ , and is thus of period $2$ . The function is affected by the choice of $t_{0}$ via the transformation

(3.28) $$\begin{eqnarray}t_{0}\mapsto t_{0}+k,\quad \text{leads to }F\mapsto F+\frac{27k}{2\sqrt{D}w^{2}}.\end{eqnarray}$$

Thus we may choose $F$ to have zero integral over each period by choosing $t_{0}$ , which means that the wave component oscillates around the single-scale component. From (3.26) we see that $D$ only affects the scales of $h$ , $u$ and $t$ , which is equivalent to choosing to make the variables dimensionless with respect to $HD$ rather than $H$ . The constant small parameter $\unicode[STIX]{x1D716}$ used to compute the solution at the time $t$ must have magnitude approximately $1/t$ for $\unicode[STIX]{x1D716}t$ to be order unity. We define the time dependent small parameter $\unicode[STIX]{x1D716}_{t}$ to mean the value of $\unicode[STIX]{x1D716}$ which rescales $t$ to exactly unity, taking into account appropriate scales and offset, that is

(3.29) $$\begin{eqnarray}\unicode[STIX]{x1D716}_{t}=\frac{1}{\sqrt{D}(t-t_{0})}.\end{eqnarray}$$

We compare these asymptotic results with the numerical simulations as follows: we suppose that we know $t_{0}$ as well as the dependent variables at some comparison time $t=t_{c}$ , given by $h(x,t_{c})$ , $u(x,t_{c})$ , and we use these fields to determine the function $F$ over a single period. The differences between the numerically computed fields and the single-scale asymptotic components, $h_{wn}=h(x,t_{c})-h_{a}(x,t_{c})$ and $u_{wn}=u(x,t_{c})-u_{a}(x,t_{c})$ are added and subtracted to give

(3.30a ) $$\begin{eqnarray}\displaystyle \frac{h_{wn}(x,t_{c})}{2D}+\frac{u_{wn}(x,t_{c})}{2\sqrt{D}} & = & \displaystyle \frac{F(\unicode[STIX]{x1D702}_{c})}{(\sqrt{D}(t_{c}-t_{0}))^{3}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{(\sqrt{D}(t_{c}-t_{0}))^{4}}\left[\frac{w^{2}F(\unicode[STIX]{x1D702}_{c})^{2}}{9}+3xF(\unicode[STIX]{x1D702}_{c})\right],\end{eqnarray}$$

(3.30b ) $$\begin{eqnarray}\displaystyle \frac{h_{wn}(x,t_{c})}{2D}-\frac{u_{wn}(x,t_{c})}{2\sqrt{D}} & = & \displaystyle \frac{F(\unicode[STIX]{x1D709}_{c})}{(\sqrt{D}(t_{c}-t_{0}))^{3}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{(\sqrt{D}(t_{c}-t_{0}))^{4}}\left[\frac{w^{2}F(\unicode[STIX]{x1D709}_{c})^{2}}{9}-3xF(\unicode[STIX]{x1D709}_{c})\right],\end{eqnarray}$$

$\unicode[STIX]{x1D709}_{c}(x)=\unicode[STIX]{x1D709}(x,t_{c})$

$\unicode[STIX]{x1D702}_{c}(x)=\unicode[STIX]{x1D702}(x,t_{c})$

$0\leqslant x\leqslant 1$

$S(y)$

(3.31) $$\begin{eqnarray}S(y)=(\sqrt{D}(t_{c}-t_{0}))^{3}\left[\frac{h_{wn}(|\unicode[STIX]{x1D712}|,t_{c})}{2D}-\operatorname{sgn}(\unicode[STIX]{x1D712})\frac{u_{wn}(|\unicode[STIX]{x1D712}|,t_{c})}{2\sqrt{D}}\right],\end{eqnarray}$$

where $y=\sqrt{D}t_{c}-27/(4w^{2}\sqrt{D}(t_{c}-t_{0}))+\unicode[STIX]{x1D712}$ and $-1<\unicode[STIX]{x1D712}<1$ . Therefore the function $F$ satisfies the following quadratic equation,

(3.32) $$\begin{eqnarray}F^{2}(y)-2BF(y)-C=0,\end{eqnarray}$$

where

(3.33a,b ) $$\begin{eqnarray}B=\frac{27\unicode[STIX]{x1D712}-9\sqrt{D}(t_{c}-t_{0})}{2w^{2}},\quad \text{and}\quad C=\frac{9\sqrt{D}(t_{c}-t_{0})}{w^{2}}S(y).\end{eqnarray}$$

We expect that $F\sim S$ for $t_{c}\gg t_{0}$ , which selects the solution to (3.32) given by

(3.34) $$\begin{eqnarray}F=B+\sqrt{B^{2}+C}.\end{eqnarray}$$

Finally, we evaluate the temporal offset, $t_{0}$ , by imposing

(3.35) $$\begin{eqnarray}\int _{0}^{2}F(y)\,\text{d}y=0.\end{eqnarray}$$

In table 1 we present the evaluation of the temporal offset, $t_{0}$ , for various barrier heights, $D$ , and uniform width, $w=1$ , using either the initial conditions $(t_{c}=0)$ or the numerically determined profiles at $t_{c}=100$ . We find that the two values of $t_{0}$ for each barrier height are broadly similar; the relative difference is less than $10\,\%$ and systematically decreases with $D$ . The computed value of $\unicode[STIX]{x1D716}_{t}$ given in table 1 according to (3.29) provides an indication of the magnitude of the difference between the numerically computed profiles and their asymptotic forms. Indeed we find that even using the initial conditions $(t_{c}=0)$ to determine $t_{0}$ appears to provide an accurate means of determining the asymptotic form for the larger barrier heights ( $D\gtrsim 0.5$ ).

Table 1. The temporal offset, $t_{0}$ , for various values of the barrier height $D$ with $w=1$ , computed using (3.35) enforced at time $t_{c}=0$ and $t_{c}=100$ . (For the latter, the dependent variables, $h$ and $u$ , were evaluated from the numerical simulations.) The corresponding values of $\unicode[STIX]{x1D716}_{t}$ are also given, along with the relative difference in the two computed values of $t_{0}$ .

In figure 7 we have plotted the single time-scale asymptotics, $h_{a}$ and $u_{a}$ (see (3.26a )–(3.26b )) as functions of time at the end of the domain $(x=1)$ . For these plots we evaluated the temporal offset using the initial data $(t_{c}=0)$ . We observe that these asymptotic forms accurately capture the numerical computations over relatively long computational times, and thus at late times for all barrier heights the outflow velocity decays with the inverse cube of time, whereas the excess height decays as the inverse square of time. We further note that there are oscillations of diminishing amplitude about the progressive decay of the height and velocity fields, which in these figures correspond to periods during which the outflow velocity and height remain constant; these oscillations are captured in the wave-like behaviour, $h_{w}$ and $u_{w}$ and are discussed below.

At late times when the depth of the fluid at the outflow is close to the barrier height $(h-D\ll 1)$ , we find that the velocity field satisfies $u\sim w(h-D)^{3/2}/D$ . Thus the flow speed is weakening and eventually viscous stresses, which are not included in the dynamical model, may become non-negligible. We show in appendix A how to assess their magnitude within the shallow flow over the crest of the weir and for the typical parameter values at environmental scales, viscous effects appear negligible until $h-D\sim 10^{-3}$ .

The rescaled volume flux at the outflow boundary, $u(1,t)h(1,t)/D^{3/2}$ is plotted in figure 8 as a function of rescaled time $\sqrt{D}(t-t_{0})$ for various values of the barrier height, $D$ . In this plot we have calculated the offset time, $t_{0}$ from the initial profile (i.e.  $t_{c}=0$ ), according to the procedure described above; values of $t_{0}$ are given as the second column of table 1. We note that for $\sqrt{D}(t-t_{0})\gtrsim 4$ , the evolution follows a universal behaviour around which each of the separate cases oscillate due to the successive wave reflections between the back wall and the outflow. Furthermore we have also plotted the analytically calculated asymptotic form, $u_{a}(1,t)h_{a}(1,t)/D^{3/2}$ , where $h_{a}$ and $u_{a}$ are given by (3.26a ) and (3.26b ), respectively. We observe that this expression accurately matches the computed flux.

Figure 8. The rescaled volume flux of fluid leaving the domain at the outflow, $u(1,t)h(1,t)/D^{3/2}$ , as a function of rescaled time $\sqrt{D}(t-t_{0})$ with the temporal offset, $t_{0}$ , calculated from the initial profile, for $D=0.1,0.3,0.5,0.7$ and $0.9$ (solid lines) and $w=1$ . The asymptotic form, $u_{a}(1,t)h_{a}(1,t)/D^{3/2}$ , is plotted with a dot-dashed line. In this figure, the curves from the different numerical computations overlay at late times and it is hard to distinguish them, indicating a universal behaviour that is accurately represented by the asymptotic form.

The shape of the waveform of the oscillations can also be determined asymptotically and may in some cases be deduced from the initial conditions. When $t_{c}=0$ , the fluid is stationary $(u(x,t_{c})=0)$ and the reservoir is of unit depth $(h(x,t_{c})=1)$ . In this case using the transform (3.28) and enforcing (3.35) we arrive at

(3.36) $$\begin{eqnarray}F(y)=\frac{27\unicode[STIX]{x1D712}}{2w^{2}}\equiv \frac{27}{2w^{2}}\left(y-\sqrt{D}t_{c}+\frac{27}{4w^{2}\sqrt{D}(t_{c}-t_{0})}\right)\end{eqnarray}$$

for $-1<\unicode[STIX]{x1D712}<1$ . This function, constructed to be periodic, contains a discontinuity at $\unicode[STIX]{x1D712}=\pm 1$ . This is exactly the same as what would be produced if we had truncated our asymptotic expansion at a lower order and included terms up to $(t-t_{0})^{-3}$ . To understand this observation, we note that the rarefaction fan which smooths out the mismatch between the initial conditions and boundary condition at early times only exists because of the nonlinearity of the shallow water equations. A linear system of hyperbolic partial differential equations would simply advect this mismatch as a discontinuity. Thus, we should only expect the asymptotic solution to capture this smoothing behaviour if we expand to sufficiently high order to capture the nonlinear interaction of the waves, and to calculate (3.26) we only expanded to capture the nonlinear effects of the boundary conditions and single-scale component, not the wave component. We therefore find that the asymptotic solution simply advects the initial mismatch as a discontinuity (which is what we see in figure 9). This failure to capture the nonlinear self-interaction of the waves means we fail to capture the correct wave shape over the first few periods. We explore this effect by making comparisons between the simulation and asymptotic solution at a number of early times to find the time at which the nonlinear effects become small and consequentially the simulation results are accurately represented by the asymptotic expressions.

Figure 9. The waveform of the oscillation at long times, $F$ , for barrier of heights (a)  $D=0.1$ ; (b)  $0.5$ ; (c)  $0.7$ ; and (d)  $0.9$ and for $w=1$ . In each figure, the waveform is extracted from the numerical simulations at seven comparison times $t_{c}$ . For each case the waveforms converge to a single profile when $t_{c}$ becomes sufficiently large, in accord with the asymptotic analysis. In (d) the waveform at late times is quite close to the waveform deduced from initial conditions ( $t_{c}=0$ , dot-dashed line), but for (a–c) there are differences due to appreciable nonlinear interactions during the initial phases of the motion.

Figure 9 shows the waveform $F(y)$ as computed from the numerical simulations for a selection of barrier heights, $D$ , and comparison times, $t_{c}$ with $w=1$ . (The curve at $t_{c}=0$ is given by (3.36) in each case.) For all cases, as time passes, the wave seen in the simulation converges on a time-independent shape, which is in accord with the asymptotic analysis. Indeed this confirms that the amplitude of the wave-like part of the dependent fields, $h_{w}$ and $u_{w}$ , decays as $(t-t_{0})^{-3}$ and that the phase speed approaches $\sqrt{D}$ (see (3.26c ) and (3.26d )). In fact we observe for $\sqrt{D}(t_{c}-t_{0})\gtrsim 10$ there is very little change in the computed shape of $F$ , indicating that by that stage the numerical results are accurately approximated by the asymptotic form. We now discuss the figures in order of increasing $\unicode[STIX]{x1D716}_{t}$ (see table 1), i.e. most accurate to least accurate.

Figure 9(d) shows that for relatively large values of $D$ , the asymptotic expansion based upon $t_{c}=0$ accurately predicts the numerically computed phase speed and the wave shape. Figures 9(c) and 9(d), in contrast, reveal some discrepancies. Examining the portions with positive gradient (which correspond to uniform and simple wave regions in the hodograph solution) we see excellent correspondence between different times suggesting that we have both the correct phase speed and wave shape for these portions. However, the portions with negative gradient (corresponding to simple and complex wave regions) are not very accurate, and these are precisely the regions in which the nonlinear wave self-interaction is important. Figure 9(a) reveals that for relatively small values of $D$ , we do not have the correct wave phase speed for early times, because the simulated waveform drifts off of the predicted wave even in the positive gradient regions. Furthermore it is evident that the calculation of the waveform $F$ from the initial conditions $(t_{c}=0)$ does not provide a useful indication of the behaviour at later times. Instead one must deduce the asymptotic form from the profiles at later times, $t_{c}$ , or one must employ a higher-order expansion, to predict the phase speed accurately.

4.1 Characteristic and hodograph methods

The draining boundary at $x=1$ (4.1) imposes a constant Froude number on the flow at this location and leads to a considerable simplification of the solution derived for flow over a barrier. In the hodograph plane this boundary condition (4.1) is represented by

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{s}=-\frac{(2-Fr)}{(2+Fr)}\unicode[STIX]{x1D6FC}.\end{eqnarray}$$

It is a straight line, rather than a curved segment as in § 3, and terminates at the origin.

The solution features a similar pattern of uniform, simple and complex wave regions (see figures 2 and 10). The key value of the characteristic variable is $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{\ast }=\unicode[STIX]{x1D6FD}_{s}(2)=-2(2-Fr)/(2+Fr)$ and the $\unicode[STIX]{x1D6FD}$ -characteristics associated with this value bound the first complex wave region attached to the back wall, $C_{1}$ and the first complex wave region attached to the drainage boundary. The construction of the solution for $t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ and $x(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ in complex wave region $C_{1}$ is identical to that presented in § 3 (see (3.4)). The integral equation for $t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ within region $C_{2}$ is, however, somewhat simpler than (3.14), because the boundary condition is enforced along a straight line in the hodograph plane. For a flow through a constriction with Froude number $Fr$ , it is given by

(4.3) $$\begin{eqnarray}\displaystyle t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s}) & = & \displaystyle B(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{\ast };\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s})t_{b3}(\unicode[STIX]{x1D6FC})+B(2,\unicode[STIX]{x1D6FD}_{\ast };\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s})t_{b3}(2){\textstyle \frac{1}{2}}(Fr-1)\nonumber\\ \displaystyle & & \displaystyle +\,B(\unicode[STIX]{x1D6FC}_{s},\unicode[STIX]{x1D6FD}_{s};\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s})t(\unicode[STIX]{x1D6FC}_{s},\unicode[STIX]{x1D6FD}_{s}){\textstyle \frac{1}{2}}(1-Fr)-\int _{2}^{\unicode[STIX]{x1D6FC}}\left(\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}a}-\frac{3B}{2(a-\unicode[STIX]{x1D6FD}_{\ast })}\right)t_{b3}(a)\,\text{d}a\nonumber\\ \displaystyle & & \displaystyle -\,\int _{\unicode[STIX]{x1D6FC}_{s}}^{2}\frac{t(a_{s},b_{s})}{2}\left((1+Fr)\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}a}+\frac{(Fr-1)(Fr-2)}{(Fr+2)}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}b}-\frac{3Fr}{2a}B\right)\,\text{d}a.\end{eqnarray}$$

As in § 3, we evaluate this expression numerically by first computing $t_{s}(\unicode[STIX]{x1D6FC})=t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D6FC}))$ , which is the $t$ -field along the boundary $(x=1)$ , before employing numerical quadrature to determine both $t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ and $x(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ . Numerical iteration is used to find the solution for $t(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D6FC}))$ – we find that a converged solution to within a tolerance of $10^{-12}$ is typically found within $6$ iterative steps. In the hodograph plane, the complex wave region $C_{2}$ forms a triangle given by $\{(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}):\unicode[STIX]{x1D6FD}_{\ast }\leqslant \unicode[STIX]{x1D6FD}\leqslant \unicode[STIX]{x1D6FD}_{s}(\unicode[STIX]{x1D6FC}),-\unicode[STIX]{x1D6FD}_{\ast }\leqslant \unicode[STIX]{x1D6FC}\leqslant 2\}$ and we note that $\unicode[STIX]{x1D6FD}_{m}=\unicode[STIX]{x1D6FD}_{s}(-\unicode[STIX]{x1D6FD}_{\ast })=2[(2-Fr)/(2+Fr)]^{2}$ .

Having numerically evaluated the solution in region $C_{2}$ we can compute the $\unicode[STIX]{x1D6FC}$ -characteristics, $x=x_{a2}$ and $x=x_{a3}$ on which $\unicode[STIX]{x1D6FC}=-\unicode[STIX]{x1D6FD}_{\ast }$ (see figures 2 and 10). The latter bounds the complex wave region $C_{3}$ and using (3.19) we may evaluate the solution within it. These steps then complete the computation of the solution following the first three reflections from the boundaries.

We illustrate the solution by plotting the temporal variation of the height at the draining boundary, $h(1,t)$ , for various values of the Froude number, $Fr$ , and consequentially, the constricted width, $w$ (see figure 11). We observe that the height at that location features periods during which it is constant, corresponding to the uniform regions in figure 10, and periods when it declines. We have computed the first three phases using quasi-analytical techniques described above; in principle these could be continued for later times, but in this study we compute the evolution using the direct numerical integration of the governing shallow-water equations (2.1a )–(2.1b ). We note the close correspondence between the two approaches, with the only apparent differences visible in figure 11 occurring when the solution changes from one wave region to another (and the temporal gradient of the height is discontinuous).

Figure 11. The height at the free-draining boundary, $h(1,t)$ , as a function of time for various values of the outflow Froude number, $Fr$ (and consequentially varying constriction widths, $w$ ). The solid lines correspond to results from integration in the hodograph plane, where the dot-dashed lines correspond to numerical integration of the shallow water equations (2.1a )–(2.1b ). The quasi-analytical hodograph calculations were performed until the end of region $U_{4}$ , but the numerical computations were continued over much longer times. We note the very close agreement between the two approaches, with the only very slight differences occurring at the junctions of the intervals during which the height is uniform or decreasing.

4.2 Long time evolution: similarity solution

After relatively long times, the flow evolves from its initial conditions to a state in which the motion is well captured by similarity solutions that are constructed below. These extend the solution developed by Momen et al. (Reference Momen, Zheng, Bou-Zeid and Stone2017) for an unconstricted outflow, in which the Froude number is unity, to a more general case where the Froude number is given by (4.1). As shown by Momen et al. (Reference Momen, Zheng, Bou-Zeid and Stone2017), for this class of similarity solution, the spatial coordinate, $x$ , is not geared to time and the height and velocity fields have the following form,

(4.4a,b ) $$\begin{eqnarray}h(x,t)=\frac{{\hat{H}}(x)}{t^{2}}\quad \text{and}\quad u(x,t)=\frac{\hat{U} (x)}{t},\end{eqnarray}$$

where ${\hat{H}}(x)$ and $\hat{U} (x)$ are to be determined. This dependence on $t$ is required in order to satisfy the drainage condition, $u(1,t)=Fr\sqrt{h(1,t)}$ , and to balance the terms in the momentum equation (2.1b ). On substitution into the governing equations (2.1a ) and (2.1b ), we find that

(4.5a,b ) $$\begin{eqnarray}\frac{\text{d}{\hat{H}}}{\text{d}x}=\frac{\hat{U} {\hat{H}}}{\hat{U} ^{2}-{\hat{H}}}\quad \text{and}\quad \frac{\text{d}\hat{U} }{\text{d}x}=\frac{\hat{U} ^{2}-2{\hat{H}}}{\hat{U} ^{2}-{\hat{H}}},\end{eqnarray}$$

together with the boundary conditions at the back wall and at the outflow, respectively given by

(4.6a,b ) $$\begin{eqnarray}\hat{U} (0)=0\quad \text{and}\quad \hat{U} (1)=Fr\sqrt{H_{1}},\end{eqnarray}$$

where $H_{1}={\hat{H}}(1)$ . We progress by introducing $\hat{V}=\hat{U} ^{2}/{\hat{H}}$ and then on the assumption that $\text{d}{\hat{H}}/\text{d}x$ does not vanish, we find that

(4.7) $$\begin{eqnarray}\frac{\text{d}\hat{V}}{\text{d}{\hat{H}}}=\frac{\hat{V}-4}{{\hat{H}}}.\end{eqnarray}$$

Integrating (4.7), and imposing ${\hat{H}}=H_{0}$ when $\hat{V}=\hat{U} =0$ at $x=0$ , we deduce that

(4.8) $$\begin{eqnarray}\hat{V}=4\left(1-\frac{{\hat{H}}}{H_{0}}\right).\end{eqnarray}$$

The latter boundary condition (4.6) demands that $\hat{V}=Fr^{2}$ when ${\hat{H}}=H_{1}$ at $x=1$ , and thus from (4.8), we find that

(4.9) $$\begin{eqnarray}\frac{H_{1}}{H_{0}}=1-\left(\frac{Fr}{2}\right)^{2}.\end{eqnarray}$$

To obtain ${\hat{H}}(x)$ we rewrite (4.8) in terms of $\hat{U}$

(4.10) $$\begin{eqnarray}\hat{U} =2\sqrt{{\hat{H}}}\sqrt{1-\frac{{\hat{H}}}{H_{0}}},\end{eqnarray}$$

and substitute into (4.5a ) to derive the implicit solution

(4.11) $$\begin{eqnarray}\frac{x}{\sqrt{H_{0}}}=-\frac{\unicode[STIX]{x03C0}}{2}+2\sqrt{\frac{{\hat{H}}}{H_{0}}}\sqrt{1-\frac{{\hat{H}}}{H_{0}}}+\arcsin \sqrt{\frac{{\hat{H}}}{H_{0}}}.\end{eqnarray}$$

The value of $H_{0}$ is found by setting $x=1$ and substituting (4.9) into the above equation, which yields

(4.12) $$\begin{eqnarray}H_{0}=\left(-\frac{\unicode[STIX]{x03C0}}{2}+Fr\sqrt{1-\left(\frac{Fr}{2}\right)^{2}}+\arccos \frac{Fr}{2}\right)^{-2}.\end{eqnarray}$$

The velocity field $\hat{U} (x)$ is found from (4.10). The solutions for various values of the Froude number are plotted in figure 12. We note that the fluid decreases in depth monotonically from the rear of the flow to the front. This generates a streamwise pressure gradient that accelerates the flow and drives the fluid through the constriction. Lower values of the Froude number, corresponding to smaller widths at the constriction, lead to reduced velocities and weaker pressure gradients.

Figure 12. The similarity functions (a)  $H(x)-H_{1}$ and (b)  $U(x)$ as a function of distance from the back wall, $x$ , for various Froude numbers at the constriction $(x=1)$ : (i)  $Fr=0.1$ ; (ii)  $Fr=0.3$ ; (iii)  $Fr=0.5$ ; (iv)  $Fr=0.7$ ; and (v)  $Fr=0.9$ .

These similarity solutions are compared with the results from the numerical integration of the governing equations; in figure 13, we compare the numerically determined height at the outflow with the prediction from the similarity solution. In this study we have integrated numerically from lock-release initial conditions, which are evidently not in self-similar form of (4.4); thus we find that the similarity solutions is approached progressively. We also note that the self-similarity functions (4.4) are only defined up to an arbitrary origin in time; in other words, we could replace $t$ with $t-\hat{t}_{0}$ in (4.4) and the same solutions are maintained. Evidently at sufficiently long times, $t\gg |\hat{t}_{0}|$ , this makes negligible difference to the predictions of the similarity solutions, but at short times, we may select $\hat{t}_{0}$ to closely match the evolution from lock-release conditions at relatively short times. In figure 13 we have chosen $\hat{t}_{0}=-\sqrt{H_{1}}$ , so that the similarity solution predicts a height of unity at $x=1$ and $t=0$ and we note that this leads to a close correspondence with the numerically computed value (see figure 13). For relatively small values of the Froude number, $Fr$ , the magnitude of the temporal offset, $\hat{t}_{0}$ is quite large; its effect is that the long time, self-similar form in which $h(1,t)$ is proportional to $1/t^{2}$ is not realised until very late times (and in the case $Fr=0.1$ , beyond the computational range plotted in figure 13).

Figure 13. The height at the outflow, $h(1,t)$ , as a function of time for various values of the Froude number: (i)  $F=0.1$ ; (ii)  $F=0.3$ ; (iii)  $F=0.5$ ; (iv)  $F=0.7$ ; and (v)  $F=0.9$ . The solid line plots the results from the numerical integration of the governing equations initiated with lock-release conditions and the dotted lines plot the similarity solution, $H_{1}/(t+\sqrt{H_{1}})^{2}$ .

The similarity solution captures the leading-order decay of the dependent variables. However it does not include the oscillations that occur as fluid sloshes between the back wall and the drainage boundary. These oscillations are manifest as intervals during which the height at the drainage boundary remains constant (figure 13). Interestingly, and in contrast to the draining flow over a barrier § 3, they appear not to have a fixed period, but rather the duration between them progressively increases. This observation will be studied in the next subsection.

4.3 Linear stability of similarity solution

We may analyse the progressive approach of the draining dynamics through a constriction to the similarity solution (4.4), by computing the linear stability of the latter. Following Grundy & Rottman (Reference Grundy and Rottman1985) and Mathunjwa & Hogg (Reference Mathunjwa and Hogg2006), we introduce a small perturbation to the solutions for the height and velocity fields, denoted by $\tilde{H}(x)$ and $\tilde{U} (x)$ respectively, so that

(4.13a,b ) $$\begin{eqnarray}h=\frac{{\hat{H}}(x)}{t^{2}}+\unicode[STIX]{x1D6FF}\frac{\tilde{H}(x)}{t^{\unicode[STIX]{x1D6FE}+2}}\quad \text{and}\quad u=\frac{\hat{U} (x)}{t}+\unicode[STIX]{x1D6FF}\frac{\tilde{U} (x)}{t^{\unicode[STIX]{x1D6FE}+1}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}$ is a small ordering parameter $(\unicode[STIX]{x1D6FF}\ll 1)$ and $\unicode[STIX]{x1D6FE}$ is to be determined. On substitution into the governing equations (2.1a ) and (2.1b ) and after linearisation, we find that

(4.14) $$\begin{eqnarray}\displaystyle & \displaystyle -(2+\unicode[STIX]{x1D6FE})\tilde{H}+\frac{\text{d}}{\text{d}x}\left(\tilde{H}\hat{U} +\tilde{U} {\hat{H}}\right)=0, & \displaystyle\end{eqnarray}$$

(4.15) $$\begin{eqnarray}\displaystyle & \displaystyle -(1+\unicode[STIX]{x1D6FE})\tilde{U} +\frac{\text{d}}{\text{d}x}\left(\tilde{U} \hat{U} +\tilde{H}\right)=0. & \displaystyle\end{eqnarray}$$

The boundary conditions represent no flow at the back wall, $\tilde{U} (0)=0$ , and the linearised drainage condition is given by

(4.16) $$\begin{eqnarray}2\hat{U} (1)\tilde{U} (1)=Fr^{2}\tilde{H}(1).\end{eqnarray}$$

Straightforwardly, we note that $\unicode[STIX]{x1D6FE}=1$ is a solution, with $\tilde{H}=2{\hat{H}}$ and $\tilde{U} =\hat{U}$ ; this corresponds to the introduction of a shift of the temporal origin, which, as discussed above, does not change the similarity solution. In addition, $\unicode[STIX]{x1D6FE}=0$ is a solution with $\tilde{U} =\text{d}\hat{U} /\text{d}x$ and $\tilde{H}=\text{d}{\hat{H}}/\text{d}x$ ; this corresponds to a translation of the $x$ -axis. Both of the these solutions are inevitably found in the study of the stability of similarity solutions and do not reveal the linear stability properties (Mathunjwa & Hogg Reference Mathunjwa and Hogg2006). Instead we seek non-trivial solutions to (4.14)–(4.15), which determine the exponent, $\unicode[STIX]{x1D6FE}$ , as an eigenvalue.

We determine the eigenvalue, $\unicode[STIX]{x1D6FE}$ , through numerical integration of (4.14)–(4.15), using the solution for ${\hat{H}}(x)$ and $\hat{U} (x)$ given by (4.10) and (4.11); in this computation, we allow the eigenvalue and functions to be complex valued, noting that since the governing equations have real-valued coefficients, they are also satisfied by the complex conjugates of the solutions. We plot the value of $\unicode[STIX]{x1D6FE}$ with the smallest real part in figure 14 as a function of the drainage Froude number, $Fr$ . We note that $\Re (\unicode[STIX]{x1D6FE})>0$ and thus the similarity solutions are linearly stable; in other words, all disturbances decay more rapidly than the solutions itself. In fact, we note that $\Re (\unicode[STIX]{x1D6FE})$ is only relatively weakly dependent upon $Fr$ . However the eigenvalue also has a non-vanishing imaginary part and its relatively strong variation is also plotted in figure 14.

Figure 14. The exponent, $\unicode[STIX]{x1D6FE}$ , as a function of the Froude number, $Fr$ .

We may use this solution to analyse the oscillations noted above in the numerically determined solution to the drainage from lock-release initial conditions (figure 13). Denoting the ratio between successive times at which the height at the end of the domain exhibits a spatially uniform region by ${\mathcal{R}}$ , we find that

(4.17) $$\begin{eqnarray}{\mathcal{R}}=\exp \left(\frac{2\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FE}_{I}}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{I}\equiv \Im (\unicode[STIX]{x1D6FE})$ denotes the imaginary part of the exponent $\unicode[STIX]{x1D6FE}$ . In figure 15 we plot ${\mathcal{R}}$ as a function of the Froude number, $Fr$ ; we observe that the ratio of successive periods increases monotonically with $Fr$ . We also examine this ratio in our numerical results by computing ${\mathcal{R}}_{n}=(t_{n}-t_{0})/(t_{n-1}-t_{0})$ , the ratio of the $n$ th to $(n-1)$ th times at which the period of a constant fluid depth at $x=1$ starts as a function of $n$ for $Fr=0.1$ and $Fr=0.5$ (see inset figures in figure 15). In both cases we find that the numerically determined ratio approaches the asymptotic prediction; initially there are differences because the motion is strongly nonlinear, but these diminish as the flow evolves into a state that is accurately represented by this linearised analysis. We note that for $Fr=0.1$ there are many oscillations observed during the computational period, because the ratio, ${\mathcal{R}}$ is quite close to unity. In contrast, when $Fr=0.5$ , the ratio of the successive oscillations is larger and so fewer can be determined in the numerical results before their amplitude has very substantially diminished.

Figure 15. The ratio of successive times at which the height at the outflow exhibits periods during which it is constant, ${\mathcal{R}}$ , as a function of the Froude number, $Fr$ . The inset figures show the numerically determined ratio, ${\mathcal{R}}_{n}$ , as a function of the oscillation number $n$ (shown with symbols) and the theoretical prediction (dashed line) for $Fr=0.1$ and $Fr=0.5$ .