# On the stability of fully nonlinear hydraulic-fall solutions to the forced
water-wave problem

J. S. Keeler<EMAIL_ADDRESS>School of Mathematics, University of East
Anglia, Norwich, NR4 7TJ, UK  M. G. Blyth<EMAIL_ADDRESS>School of
Mathematics, University of East Anglia, Norwich, NR4 7TJ, UK

###### Abstract

Two-dimensional free-surface flow over localised topography is examined with
the emphasis on the stability of hydraulic-fall solutions. A Gaussian
topography profile is assumed with a positive or negative amplitude modelling
a bump or a dip, respectively. Steady hydraulic-fall solutions to the full
incompressible, irrotational Euler equations are computed, and their linear
and nonlinear stability is analysed by computing eigenspectra of the pertinent
linearised operator and by solving an initial value problem. The computations
are carried out numerically using a specially developed computational
framework based on the finite element method. The Hamiltonian structure of the
problem is demonstrated and stability is determined by computing eigenspectra
of the pertinent linearised operator. It is found that a hydraulic-fall flow
over a bump is spectrally stable. The corresponding flow over a dip is found
to be linearly unstable. In the latter case, time-dependent simulations show
that the flow ultimately settles into a time-periodic motion that corresponds
to an invariant solution in an appropriately defined phase space. Physically,
the solution consists of a localised large amplitude wave that pulsates above
the dip while simultaneously emitting nonlinear cnoidal waves in the upstream
direction and multi-harmonic linear waves in the downstream direction.

††preprint: APS/123-QED

##  I Introduction

The incompressible flow of an inviscid liquid with a free surface over
topography and/or subject to a surface pressure distribution is a classical
problem in fluid dynamics that has received considerable attention over the
last two hundred years [see, for example 42, 2, 16]. The principal aim is to
determine the free-surface profile and to explore how it changes as the
relevant parameters, such as the oncoming flow speed or the amplitude of the
forcing (either the topography or pressure distribution) are varied. For a
localised forcing, it is well known that in one possible steady flow
configuration, herein referred to as a hydraulic fall, the fluid level drops
from a higher uniform level upstream to a lower uniform level downstream
(figure 1). Such flows have been observed experimentally [e.g. 19, 37], and
have been computed as solutions to the fully nonlinear irrotational Euler
(hereinafter FNL) equations [e.g. 18, 15], and as solutions to reduced-order
models such as the forced Korteweg-de Vries (hereinafter fKdV) equation [see,
for example, 19, 16]. A recent review is provided by [5].

Most of the studies in the literature that concern hydraulic-fall solutions
focus on steady flow, and there appears to have been very little effort
devoted to gaining a theoretical understanding of the stability of these
flows. Experimental studies have tended to concentrate on flow over a bump as
depicted in figure 1(b), and, indeed, we have been unable to identify any
experimental work on hydraulic falls over a dip (figure 1(a)). This may simply
be due to the fact that presumably a wave tank with a bump is easier to
engineer than one with a dip. In this article we discuss the stability
properties of hydraulic-fall solutions for both positive and negative
forcings; (i) by carrying out a linear stability analysis, and (ii) by
developing a computational framework based on the finite-element method that
we use to analyse nonlinear stability. We show that over a bump the hydraulic-
fall solution is spectrally stable, but over a dip it is linearly unstable.
Most intriguingly, from our suite of time-dependent calculations we identify a
new time-dependent invariant solution of the FNL system which corresponds to a
stable time-periodic orbit in a suitably defined phase space.

Figure 1: Hydraulic-fall solutions. Sketch of the two basic types in the non-
dimensional domain, the arrow indicates the flow direction. (a) Bump forcing,
$Fr<1$ (b) Dip forcing, $Fr<1$ where $Fr$ is the Froude number, defined in
(1). The downstream Froude number, $Fr_{\mathrm{d}}$, is of opposite
criticality to $Fr$.

The fKdV and FNL steady solution spaces for various topographic forcings are
by now well established [e.g. 19, 16, 40, 41]. For the fKdV system a number of
different solution types can be identified [2, 6, 16]; for a recent review see
[5]. Broadly speaking, these can be classified as cnoidal wave solutions which
are flat upstream and wavy downstream, solitary-wave solutions which are flat
both upstream and downstream, and transcritical solutions which we are herein
referring to as hydraulic falls. These solutions are characterised by the
upstream depth-based Froude number given by

$Fr=\frac{U}{\sqrt{gH}},$ (1)

where $g$ is the gravitational constant, and $U$ and $H$ are the speed of the
flow and the fluid depth far upstream. Cnoidal-wave solutions occur when
$Fr<1$ (subcritical flow), and solitary-wave solutions occur when $Fr>1$
(supercritical flow). A distinguishing feature of hydraulic-fall solutions is
that the flow upstream is subcritical and the flow downstream is
supercritical. For a given topographic forcing, a steady hydraulic-fall
solution exists only for a particular value of $Fr$. Furthermore, for
hydraulic-fall solutions, if the forcing is negative (a dip) then the wave
profile has a local maximum above the forcing, and if the forcing is positive
(a bump) the fluid depth decreases monotonically from upstream to downstream
(see figure 1).

Using numerical simulations of the fKdV system carried out for a number of
different types of perturbations to the base steady state, [17], [12] and [13]
all concluded that a hydraulic fall over a positive forcing is stable. [32]
showed via time-dependent calculations starting from a small wavy Gaussian
perturbation, that the hydraulic fall solution appeared to be stable although
they did not track the long-term behaviour. None of these studies attempted a
formal linear stability analysis.

The article proceeds as follows. In § II we state the mathematical problem and
the weak formulation of the problem that forms the basis of our numerical
framework. In § III.1 we briefly describe the steady bifurcation structure of
the FNL system. We carry out a linear stability analysis in § III.2 and a
nonlinear stability analysis via direct numerical simulations of the FNL model
in § III.3. Finally in § III.4 and § IV we discuss our results and consider
possible avenues for future research.

##  II Problem formulation

We consider inviscid, irrotational, incompressible flow over a bottom
topography under the influence of gravity and neglecting surface tension. The
general flow scenario is depicted in figure 2. The fluid flows from left to
right over a localised topographic structure seen as a positive bump in the
figure. In the absence of the topography so that the bottom is flat, or
sufficiently far upstream of the obstacle, the flow consists of a uniform
stream of strength $U$ with uniform fluid depth $H$. Our concern is with the
generally unsteady disruption that is provoked at the free surface by the
topography.

We solve for the fluid flow and the deformation of the free surface using a
numerical finite-element method. In the following sub-sections we first
describe the mathematical problem to be solved before discussing the truncated
computational problem that is required for the numerical implementation. To
aid the discussion, we use primes to denote regions and boundaries in the
mathematical problem. For example, $\Omega^{\prime}(t)$ represents the fluid
domain in the mathematical problem, and $\Omega(t)$ represents the
corresponding truncated computational domain.

Figure 2: Sketch of the non-dimensional time-dependent problem domain for the
hydraulic-fall problem. (a) The mathematical domain. The fluid domain is
$\Omega^{\prime}(t)$, the bottom boundary is $\Gamma^{\prime}_{0}$ and the
free-surface is $\Gamma^{\prime}_{2}$. As $x\to-\infty$, we impose uniform
flow; $\phi\sim x$. (b) The computational domain. The inflow boundary is
denoted $\Gamma_{3}$ and the outflow boundary is $\Gamma_{0}$. The upstream
and downstream sponges are indicated by dotted lines on the free-surface
$\Gamma_{2}$. In the computational domain, far downstream we impose that
$\phi=Vx$. In both panels, the normal vectors to each boundary $\Gamma_{i}$
are given as $\textbf{n}_{i}$.

###  II.1 Mathematical formulation

We non-dimensionalise the problem by scaling all lengths and velocities by the
height $H$ and strength $U$ of the oncoming uniform stream, respectively. The
corresponding time scale is $H/U$. Figure 2(a) shows the domain for the
mathematical problem expressed in terms of dimensionless variables and with
reference to a Cartesian set of axes $(x,y)$. The free-surface is described by
$\textbf{r}=(x_{f},y_{f})$ and the bottom topography by
$\boldsymbol{\sigma}=(x_{b},y_{b})$. Written in terms of the velocity
potential, $\phi$, the dimensionless governing equation and boundary
conditions are

$\displaystyle\nabla^{2}\phi=0,\quad\textbf{x}\in\Omega^{\prime}(t)\quad$
(Conservation of mass) (2)
$\displaystyle\nabla\phi\cdot\textbf{n}_{0}=0,\quad\textbf{x}\in\Gamma^{\prime}_{0},\quad$
(No penetration on bottom) (3)
$\displaystyle{\color[rgb]{0,0,0}\frac{\partial\textbf{r}}{\partial
t}\cdot\textbf{n}_{2}=\nabla\phi\cdot\textbf{n}_{2},\quad\textbf{x}\in\Gamma^{\prime}_{2}(t),}\quad$
(Kinematic condition) (4)
$\displaystyle\phi_{t}+\frac{1}{2}|\nabla\phi|^{2}+\frac{1}{Fr^{2}}\left(y_{f}-1\right)-\frac{1}{2}=0\quad\textbf{x}\in\Gamma^{\prime}_{2}(t).\>$
(Dynamic condition) (5)

The Froude number $Fr$ was defined in (1). The dynamic boundary condition
follows from applying Bernoulli’s equation at the free surface and utilising
the upstream uniform flow conditions, namely

$\displaystyle y_{f}\to 1,\>\phi\to x\qquad\mbox{as}\qquad x\to-\infty.$ (6)

Throughout we shall adopt the Gaussian topographic forcing function

$y_{b}=a\mbox{e}^{-b^{2}x_{b}^{2}},$ (7)

for some $a,b\in\mathbb{R}$.

###  II.2 Truncated computational formulation

To prepare for the numerical discretisation, as shown in figure 2(b) we
truncate the infinite domain so that ${x}\in[-L,L]$, with $L$ to be chosen to
be sufficiently large. The problem on the truncated domain is given by (2)–(5)
(with the primes on the flow domain and boundary symbols removed). Boundary
conditions must be imposed at the artificial outflow and inflow boundaries,
$\Gamma_{1}$ and $\Gamma_{3}$, respectively, that arise due to the truncation.

At the inflow boundary $\Gamma_{3}$, corresponding to ${x}=-L$, the velocity
potential and free surface level are set to mirror the condition (6) in the
mathematical formulation. Accordingly, we impose

${y}_{f}=1,\quad{\phi}={x},\qquad{\textbf{x}}\in\Gamma_{3},\qquad(\mbox{Inflow
condition).}$ (8)

which imposes a Dirichlet condition on $\phi$. The conditions to be imposed at
the outflow boundary are more subtle. Anticipating the presence of waves
propagating downstream, for our unsteady calculations we will make use of a
sponge layer to force the free surface to become locally flat downstream
allowing us to impose the outflow Neumann boundary condition

$\nabla{\phi}\cdot\textbf{n}_{1}=V\qquad\mbox{on}\qquad\textbf{x}\in\Gamma_{1},\qquad(\mbox{Ouflow
condition)}$ (9)

where $\textbf{n}_{1}=\textbf{i}$ is the unit vector in the $x$ direction, and
$V$ is set to ensure conservation of mass. In this way we expect to capture
the genuine flow features across the main part of the computational domain
(the ‘computational window’ – see figure 2b) except in a short region close to
the outflow boundary where the flow is artificially forced into becoming a
locally uniform stream. The choice of $V$ together with details of the sponge-
layer will be discussed below.

###  II.3 Steady calculations

In presenting the steady solution space we include both types for completeness
and to provide greater context. Since hydraulic rise solutions are subcritical
at the downstream end (contrast the hydraulic fall solutions sketched in
figure 1), we found that in practice it is convenient to replace the
conditions on $\phi$ in (8) and (9) with a Neumann condition on $\phi$ at the
inflow boundary $\Gamma_{3}$ and a Dirichlet condition on $\phi$ at the
outflow boundary. In this way we can prevent the occurrence of cnoidal waves
downstream,

As was mentioned above, the downstream speed $V$ in (9) is chosen to ensure
conservation of mass. If $\gamma_{s}=y_{f}(x=L)$ is the a priori unknown free
surface level local to $\Gamma_{1}$, then we must choose $V=1/\gamma_{s}$. The
steady form of the dynamic boundary condition (5) then requires

$\displaystyle\frac{1}{2}\gamma_{s}^{-2}+\frac{1}{Fr^{2}}\left(\gamma_{s}-1\right)-\frac{1}{2}=0.$
(10)

This is satisfied for any $Fr$ if $\gamma_{s}=1$, but if $\gamma_{s}\neq 1$ it
imposes a constraint on the Froude number, and thus the Froude number must
come as part of the solution.

While our focus in this paper is on hydraulic fall solutions, we can also
calculate hydraulic-rise solutions for which the fluid depth increases
monotonically from its upstream level, where $Fr>1$, to its downstream level,
where $Fr<1$. These solutions are not observed in experiments and so the main
focus of our stability analysis is hydraulic falls, although we do present a
brief description of the nonlinear stability of hydraulic rises for
completeness.

###  II.4 Time-dependent calculations

In the time-dependent problem when waves reach the artificial boundaries in
the truncated computational domain at $x=\pm L$ they will either be reflected
back into the domain or the simulation will fail. A number of different
approaches have been proposed for dealing with this issue (see [35] for a
review). These include imposing appropriate radiation conditions [10, 39] or
introducing perfectly-matched layers (see, for example, [4] for acoustic
waves). We choose a simpler approach by introducing so-called sponge-layers
adjacent to the inflow and outflow boundaries. The sponge-layer was described
by [9] and [3], and it has been implemented, for example, by [21] and [27] in
the case of the forced KdV equation.

In the sponge layers, disturbances near to the inflow and outflow boundaries
are exponentially damped so that

$y_{f}-y_{\pm}\propto\mbox{e}^{-A_{\pm}t},$ (11)

where $y_{-}$ and $y_{+}$ are, respectively, target upstream and downstream
free surface levels that will be discussed below, and $A_{+}$ and $A_{-}$ are
specified decay rate constants. To achieve this we modify the kinematic
condition, (4), to become

$\frac{\partial\textbf{r}}{\partial
t}\cdot\textbf{n}_{2}-\underbrace{S_{-}(x)(y_{f}-y_{-})}_{\mbox{Upstream
sponge}}-\underbrace{S_{+}(x)(y_{f}-y_{+})}_{\mbox{Downstream
sponge}}=\nabla\phi\cdot\textbf{n}_{2}.$ (12)

The sponge functions $S_{\pm}(x)$ take the form

$S_{+}(x)=A_{+}[1+\tanh(B_{+}(x-C_{+}))],\quad
S_{-}(x)=A_{-}[1+\tanh(-B_{-}(x-C_{-}))],$ (13)

so that the exponential damping only occurs over specified regions dictated by
the parameters $B_{\pm}$ and $C_{\pm}$. Care must be exercised when choosing
values for $A_{\pm},B_{\pm},C_{\pm}$ to ensure that free surface waves do not
get reflected back into the main part of the computational domain. After
extensive experimentation, we found that taking $A_{\pm}=5$ is appropriate,
and setting $B_{+}=1$ ensures that small amplitude waves are absorbed
downstream, and taking $B_{-}=0.01$ avoids large amplitude waves being
reflected back into the domain at the upstream end. We set the locations of
the sponge layers by taking $C_{\pm}=\pm L\mp 10$, and typically we choose
$L>100$ – see figure 2(b) where the sponge region is marked approximately by
dashed lines on $\Gamma_{2}$.

The upstream target level we set as $y_{-}=1$, whereas for the downstream
target level, we either set $y_{+}=\gamma_{s}$ or $y_{+}=1$ depending on the
initial condition. This choice will be discussed in § III.3 and § III.4.
Regardless of the choice of $\gamma_{s}$, we fix $V=1/y_{f}(L,t)$ so that the
outflow flux is equal to unity. We stress that this technique of damping the
waves upstream and downstream is independent of the underlying numerical
discretisation of the system. It could in principle be applied to the
analogous 3D problem as an alternative to the radiation conditions used in
[10, 39].

###  II.5 Weak formulation

The governing equations, (2)–(5) are highly nonlinear and numerical methods
are typically required to solve them. Boundary-integral methods are popular
[see, for example 8, 16, 7, 19, 40, 41], although finite-difference schemes
have been used [23] and more recently a spectral method has been implemented
[20]. We choose a different approach and develop a numerical framework based
on a weak formulation of the governing equations over the computational domain
(figure 2b). The advantages of this approach shall be discussed in § II.6 but
first we describe the mathematical weak formulation.

We multiply (2) by a test function, $\psi(\textbf{x})$, which is required to
vanish on boundaries where Dirichlet conditions in $\phi$ are imposed.
Integrating (2) over the computational domain and integrating by parts yields

$\iint_{\Omega(t)}\left(\nabla^{2}\phi\right)\,\psi\,\mbox{d}V\equiv\int_{\partial\Omega(t)}\left(\textbf{n}\cdot\nabla\phi\right)\psi\,\mbox{d}S-\iint_{\Omega(t)}\nabla\phi\cdot\nabla\psi\,\mbox{d}V=0,$
(14)

where $\partial\Omega(t)$ represents the boundary of $\Omega(t)$, n is the
outwards-pointing unit normal vector on each part of the boundary and
$\mbox{d}V$, $\mbox{d}S$ are differential area and line elements,
respectively. The domain boundary can be decomposed into
$\partial\Omega(t)=\Gamma_{0}+\Gamma_{1}+\Gamma_{2}(t)+\Gamma_{3}$ and so,
imposing the Neumann conditions (3), (4) and (10), means that (14) becomes

$\mathcal{R}_{\mathrm{Bulk}}(\textbf{x},\phi(\textbf{x}))\equiv\int_{\Gamma_{2}(t)}\underbrace{\left(\frac{\partial\textbf{r}}{\partial
t}\cdot\textbf{n}_{2}\right)}_{\mbox{Kinematic
Cond.}}\psi\,\mbox{d}S+\int_{\Gamma_{1}}\underbrace{V}_{\mbox{Outflow}}\psi\,\mbox{d}S-\iint_{\Omega(t)}\nabla\phi\cdot\nabla\psi\,\mbox{d}V=0,$
(15)

with the unit normal $\textbf{n}_{2}$ identified in figure 2. To satisfy the
dynamic boundary condition we multiply (5) by a test function and integrate
over the free-surface to obtain

$\mathcal{R}_{\mathrm{Dyn}}(\textbf{x},\phi(\textbf{x}))\equiv\int_{\Gamma_{2}(t)}\left(\phi_{t}+\frac{1}{2}|\nabla\phi|^{2}+\frac{1}{Fr^{2}}\left(y_{f}-1\right)-\frac{1}{2}\right)\psi\,\mbox{d}S=0.$
(16)

This finalises the weak formulation of the problem. Equations (15) and (16)
are to be solved to determine for the free-surface location, $y_{f}(x,t)$, and
the velocity potential in the fluid, $\phi(\textbf{x},t)$.

###  II.6 Numerical discretisation

To solve (15) and (16) numerically, we utilise the open-source C++ package
oomph-lib [24] that discretises the governing equations using a FEM Galerkin
method and allows us to take advantage of state-of-the-art linear solvers,
mesh-update techniques and provides flexibility to switch between steady,
time-dependent and linear stability calculations. We utilise an isoparametric
representation, wherein the same interpolating shape functions are used for
$\psi$ and for the position variables. For this problem we used piecewise
cubic shape functions and typically chose 400 $\times$ 10 elements in the
$x\times y$ directions. We use a structured quadrilateral mesh with either a
spine-node update strategy (for the stability analysis and time-dependent
calculations) or a pseudo-solid elastic mesh update strategy (for the steady
calculations), although in principle other mesh geometries, such as an
unstructured triangular mesh, can be used. Newton iterations are used to
calculate steady states and a backwards-difference Euler order 2 (BDF2)
implicit method is used for unsteady time-stepping; after extensive
experimentation we set the timestep $\Delta t=1.0$, which is small enough to
resolve the temporal features of the solutions to be discussed. An attractive
feature of this formulation is that time-dependent calculations only require a
very minor (and easy) augmentation of the Jacobian matrix used in the Newton
iterations and so it is straightforward to switch between steady-state
calculations and time-dependent calculations. Another notable feature is that
despite calculating quantities in the bulk fluid as well as the free-surface,
the corresponding Jacobian matrix is sparse, in contrast to the boundary-
integral approaches, and the linear inversion is quick based on the open-
source SuperLu package [31] which performs LU factorisation. Finally, in the
weak formulation, the numerical generalised eigenvalue problem that results
from the linear stability analysis is highly rank deficient as time-
derivatives do not explicitly appear in the bulk fluid. To compute the
eigenvalues and eigenmodes accurately and efficiently, we implement the
eigensolver from the Anasazi linear algebra library [25] that is based on
Arnoldi-iteration and has been used successfully in other rank-deficient
generalised eigenproblems [38, 29].

##  III Steady solutions, stability analysis, and time-dependent simulations

In this section we first describe the steady hydraulic-fall solution
structure. Next we establish the linear stability properties of the hydraulic-
fall solutions by solving a generalised eigenvalue problem and probe the
nonlinear stability properties by performing time-dependent simulations.

###  III.1 Steady bifurcation structure

Figure 3: Steady solution space for $b=0.3$. The amplitude of the forcing,
$a$, is the horizontal measure and the vertical measure is the Froude number,
$Fr$. The inset panels show FNL profiles on the branch indicated by arrows,
Here $a^{*}\approx-0.016$, $(a_{F},\,Fr_{F})\approx(-0.0165,1.006)$.

We briefly describe the steady bifurcation structure for the FNL system. In
figure 3 we show the hydraulic-fall solution space in the $(a,Fr)$ plane. The
hydraulic-rise steady states are included for completeness.

We focus on the difference between the free-surface profiles for $a<0$ and
$a>0$. For $a>0$ the hydraulic-fall/rise profiles are monotonic
decreasing/increasing functions of $x$, respectively, and there is a unique
hydraulic fall solution for each $a$. For $a<0$ the profiles are non-monotonic
and in each case the wave height reaches a maximum close to $x=0$ before
monotonically decreasing (vice-versa for the hydraulic-rise profiles).

Unlike the corresponding curve for the delta-function-forced KdV equation
[see, for example 5], the solution curve in figure 3 is neither symmetric
about $a=0$ nor about $Fr=1$ (although the latter appears true, close
inspection reveals this not to be the case – see the zoomed inset diagram in
the figure). Again in contrast to what is observed for the forced KdV
equation, the solution curve has a turning point which is located at
$(a_{F},Fr_{F})$, with $Fr_{F}>1$, so that there is a parameter window,
$a_{F}<a<a^{*}$, in which there is a multiplicity of hydraulic-rise solutions.

There are no hydraulic fall solutions for $a\leq a^{*}$. We encounter
numerical difficulties in reaching the point $a=a^{*}$. As the branch
approaches this point the wave-profile, as shown in the top left inset of
figure 3 (labelled A) appears to approach a solitary-wave type solution in
which the downstream height $\gamma_{s}\to 1$. We believe that the profile at
$a=a^{*}$ corresponds to that at the termination point of the $Fr=1$ branch of
solution (not computed here), as described by [28] for the fKdV model.

###  III.2 Linear stability analysis

To analyse the stability of the steady solutions discussed in the previous
section, we first introduce the velocity potential on the free-surface,
$\varphi(x,t)\equiv\phi(x,y_{f},t)$, and the surface elevation $\eta\equiv
y_{f}-1$. We then write

$\varphi=\varphi_{s}(x)+\hat{\varphi}(x,t),\qquad\eta=\eta_{s}(x)+\hat{\eta}(x,t),$
(17)

where $\varphi_{s}(x)$ and $\eta_{s}(x)$ are the velocity potential and
elevation for a steady hydraulic fall solution, and the hatted variables are
time-dependent perturbations which are assumed to vanish at infinity. We
emphasise that the perturbations are arbitrary and not at this point assumed
to be small.

It is helpful for the subsequent linear stability analysis to write the system
of equations in Hamiltonian form. However, since we perturb about steady
hydraulic-fall solutions for which the surface potential and elevation
functions, $\varphi_{s}$ and $\eta_{s}$, do not vanish as $x\to\infty$, some
modification to the well-known Craig-Sulem-Zakharov (CSZ) Hamiltonian
formulation of the water wave problem [44, 14] is necessary, which we now
describe. We note that the base state in (17) does not necessarily have to
represent a hydraulic-fall solution; the following analysis is also valid if
the underlying steady state is a hydraulic-rise or solitary-wave solution, for
example.

In constructing the Hamiltonian it is important to realise that (17) reflects
a nonlinear perturbation to free-surface quantities, and that some care must
be exercised when defining corresponding bulk velocity potentials. In
particular we define $\xi(x,y,t)$ and $\hat{\phi}(x,y,t)$ to be the solutions
to the following Dirichlet-Neumann problems, denoted A and B, which are both
defined on the time-dependent domain $\Omega^{\prime}(t)$:

$\mbox{Problem
A}:\left\\{\begin{split}\nabla^{2}\xi=0\qquad\mbox{in}\qquad\Omega^{\prime}(t),\\\
\xi=\varphi_{s}\qquad\mbox{on}\qquad\Gamma^{\prime}_{2}(t),\\\
\textbf{n}_{0}\cdot\nabla\xi=0\qquad\mbox{on}\qquad\Gamma^{\prime}_{0},\\\
\xi\sim x\quad\mbox{as}\quad x\to-\infty,\\\
\xi\sim\gamma_{s}^{-1}x\quad\mbox{as}\quad x\to\infty,\\\
\end{split}\>\right.\quad\mbox{Problem
B}:\left\\{\begin{split}\nabla^{2}\hat{\phi}=0\qquad\mbox{in}\qquad\Omega^{\prime}(t),\\\
\hat{\phi}=\hat{\varphi}\qquad\mbox{on}\qquad\Gamma^{\prime}_{2}(t),\\\
\textbf{n}_{0}\cdot\nabla\hat{\phi}=0\qquad\mbox{on}\qquad\Gamma^{\prime}_{0},\\\
\hat{\phi}\to 0\quad\mbox{as}\quad x\to-\infty,\\\ \hat{\phi}\to
0\quad\mbox{as}\quad x\to\infty.\end{split}\>\right.$ (18)

We assume $\hat{\eta}$ and higher-order corrections to $\xi$ decay
sufficiently fast as $x\to\pm\infty$ for the integrals that we will state
below (e.g. (20)) to be convergent. A subtle point to note is that $\xi$ does
not correspond to the velocity potential of the steady state as it is defined
on the time-dependent domain, which in general does not coincide with that for
the steady state. However, the trace of $\xi$ on the free-surface coincides
with that for the steady problem (we note that this construction relies on the
assumption that the free-surface is a graph).

The combined velocity potential $\phi=\xi+\hat{\phi}$ satisfies Laplace’s
equation, (2), in the time-dependent domain. Taking inspiration from the CSZ
construction, Problems A and B have been formulated in such a way that we may
write down a Hamiltonian for our problem in the form

$\begin{split}\hat{\mathcal{H}}(\hat{\varphi},\hat{\eta})=\hat{\mathcal{H}}_{0}(\xi,\eta_{s},\hat{\eta})+\underbrace{\frac{1}{2}\iint_{\Omega^{\prime}(t)}|\nabla\hat{\phi}|^{2}\,\mbox{d}V+\frac{1}{2Fr^{2}}\int_{-\infty}^{\infty}\hat{\eta}^{2}\,\mbox{d}x}_{\mathrm{CSZ}}\\\
+\underbrace{\iint_{\Omega^{\prime}(t)}\nabla\xi\cdot\nabla\hat{\phi}\,\mbox{d}V+\frac{1}{Fr^{2}}\int_{-\infty}^{\infty}\eta_{s}\hat{\eta}\,\mbox{d}x}_{\mathrm{Perturbation}}.\end{split}$
(19)

We have split $\hat{\mathcal{H}}$ into a term that replicates the original CSV
Hamiltonian, a perturbation term that arises from the nonlinear perturbation
in (17), and a term that can be considered as the ‘base’ energy,
$\hat{\mathcal{H}}_{0}$, defined as

$\begin{split}\hat{\mathcal{H}}_{0}(\xi,\eta_{s},\hat{\eta})=\int_{-\infty}^{\infty}\left[\frac{1}{2}\int_{0}^{1+\eta_{s}+\hat{\eta}}|\nabla\xi|^{2}-1\,\mbox{d}y+\frac{1}{2Fr^{2}}\eta_{s}^{2}+Q\right]\,\mbox{d}x,\quad
Q=\frac{1}{2Fr^{2}}\eta_{s}\left(\eta_{s}+2\right).\end{split}$ (20)

It is important to note that the integral in (20) is convergent due to the
addition of the term, $Q$ (Since $Q$ does not depend on $\hat{\phi}$ or
$\hat{\eta}$ it does not contribute on taking variations of
$\hat{\mathcal{H}}$ with respect to these variables). As $x\to-\infty$ the
terms inside the square brackets in (20) vanish due to the inflow conditions
(8). As $x\to\infty$ they vanish due to the downstream condition (10), posed
on the mathematical domain as $x\to\infty$, assuming that we restrict the flow
to approach a uniform stream with free-surface height either 1 or
$\gamma_{s}$. This restriction as $x\to\infty$ is satisfied when discussing
computations for nonlinear stability in § III.3.

Finally, by using variational arguments as described in [44] and [14] it can
be shown that the system in (2)–(5), together with (17) can be written as

$\frac{\partial\boldsymbol{\zeta}}{\partial
t}=\textbf{K}\frac{\delta\hat{\mathcal{H}}}{\delta{\boldsymbol{\zeta}}},\qquad\textbf{K}=\begin{pmatrix}0&I\\\
-I&0\end{pmatrix},\qquad\frac{\delta}{\delta\boldsymbol{\zeta}}=\left(\frac{\delta}{\delta\hat{\eta}},\frac{\delta}{\delta\hat{\varphi}}\right)^{T},$
(21)

where $\boldsymbol{\zeta}=(\hat{\eta},\hat{\varphi})^{T}$, where
$\delta/\delta\boldsymbol{\zeta}$ is the variational derivative.

Although the stability theory of Hamiltonian systems is well known [see, for
example 26], it is helpful to repeat the salient details. First, we linearise
(21) about $\hat{\varphi}\equiv 0\equiv\hat{\eta}$, that is we set
$\boldsymbol{\zeta}=\varepsilon\bar{\boldsymbol{\zeta}}$, with
$0<\varepsilon\ll 1$ and
$\bar{\boldsymbol{\zeta}}=(\bar{\eta},\bar{\varphi})^{T}$. At $O(\varepsilon)$
we find

$\frac{\partial\bar{\boldsymbol{\zeta}}}{\partial
t}=\textbf{K}\textbf{L}\>\bar{\boldsymbol{\zeta}},\qquad\textbf{L}=\frac{1}{2}\frac{\delta}{\delta\bar{\boldsymbol{\zeta}}}(\delta^{2}\hat{\mathcal{H}}),$
(22)

where $\delta^{2}\hat{\mathcal{H}}$ is the symmetric Hessian matrix (see, for
example [26]). Proceeding further, we write
$\bar{\boldsymbol{\zeta}}=\textbf{g}_{s}\mbox{e}^{st}$, where
$\textbf{g}_{s}=(g_{\hat{\eta},s}(x),g_{\hat{\varphi},s}(x))^{T}$ so that (22)
becomes the eigenvalue problem

$s\textbf{g}_{s}=\textbf{K}\textbf{L}\>\textbf{g}_{s}$ (23)

for eigenvalues $s$ and eigenmodes $\textbf{g}_{s}$.

We distinguish between two classes of solutions to (23) which correspond to i)
a continuous essential spectrum, $s_{\mathrm{ess}}$, and ii) a discrete point
spectrum, $s_{\mathrm{p}}$. The eigenmodes of $s_{\mathrm{ess}}$ are bounded
as $|x|\to\infty$ whilst the eigenmodes of $s_{\mathrm{p}}$ decay to zero as
$|x|\to\infty$. With K skew-symmetric and assuming L self-adjoint, it is easy
to show that if $s=\rho\in\mathbb{C}$ satisfies (23), then so does $s=-\rho$
and $s=\pm\rho^{*}$ (stars here indicate complex conjugates) and therefore
there is a four-fold symmetry in the complex plane. This is a standard
distinguishing feature of Hamiltonian systems [see, for example 26].

We calculate $s_{\mathrm{ess}}$ and $s_{\mathrm{p}}$ numerically using our
finite-element framework (this implementation is identical to the procedure
described in [38]). The calculation is delicate and requires care to ensure
convergence. In addition, we remark that to calculate $s_{\mathrm{ess}}$ and
$s_{\mathrm{p}}$ we have to solve three separate numerical problems (one for
$s_{\mathrm{ess}}$ and two for $s_{\mathrm{p}}$); the precise details will be
discussed in parallel with the results below.

####  III.2.1 The essential spectrum, $s_{\mathrm{ess}}$

The essential spectrum, $s_{\mathrm{ess}}$, can be calculated by examining the
form of (23) in the limit as $|x|\to\infty$ [see, for example 36]. An
interesting aspect of this problem is that the operator L takes different
forms as $x\to-\infty$ and as $x\to\infty$, and consequently different
dispersion relations are obtained in these limits. This has implications for
the spatial wave number of the eigenmodes.

Figure 4: The dispersion curves for $Fr=0.9659$ and $\gamma_{s}=0.9550$
(corresponding to the hydraulic-fall solution when $a=-0.01$, $b=0.3$). Panel
(a): the upstream wave frequency $\omega$ as a function of the upstream wave
number $k_{\mathrm{u}}$. Panel (b): the downstream wave frequency $\omega$ as
a function of the upstream wave number $k_{\mathrm{d}}$. The minimum,
$\omega=\tau=-0.006217$, of the upstream dispersion curve is shown with a
dotted line in panel (a).

Sufficiently far downstream the flow corresponds to a uniform stream of height
$\gamma_{s}$ with speed $V=\gamma_{s}^{-1}$. The dispersion relation relating
the frequency, $\omega$, to the wave number, $k_{d}$, of small amplitude waves
is

$\omega=Vk_{\mathrm{d}}\pm\sqrt{\frac{k_{\mathrm{d}}\tanh{k_{\mathrm{d}}\gamma_{s}}}{Fr^{2}}},\qquad\mbox{(Downstream
dispersion relation)}$ (24)

corresponding to waves which travel faster/slower, respectively, than the
downstream speed, $V$. Upstream where the fluid depth and speed are both
unity, with an obvious shift in notation the dispersion relation is

$\omega=k_{\mathrm{u}}\pm\sqrt{\frac{k_{\mathrm{u}}\tanh{k_{\mathrm{u}}}}{Fr^{2}}}.\qquad\mbox{(Upstream
dispersion relation)}$ (25)

In both cases the essential spectrum is such that

$s_{\mathrm{ess}}=\mathrm{i}\omega,\qquad\textbf{g}_{\mathrm{s}}\propto(\mathrm{e}^{\mathrm{i}kx},\mathrm{i}\mathrm{e}^{\mathrm{i}kx})^{T},$
(26)

with $k=k_{u}$ or $k_{d}$, and $\omega\in\mathbb{R}$. The profile of
$g_{\hat{\eta},s}$ will be such that it connects a spatially oscillatory wave
comprising wave numbers $k_{\mathrm{d}}$ that satisfy (24) as $x\to\infty$, to
a different spatially oscillatory wave consisting of wave numbers
$k_{\mathrm{u}}$ that satisfy (25) as $x\to-\infty$.

The dispersion relations in (24) and (25) are shown in figure 4 for
$a=-0.01,Fr=0.9659$. We plot the dispersion curves for positive $k$ and for
the minus sign in (24) and (25), but note that the dispersion curve for the
plus sign can be obtained by simply replacing $k$ and $\omega$ with $-k$ and
$-\omega$. We emphasise an important distinction between the upstream and
downstream dispersion curves. As can be seen in panel (a) the upstream
dispersion curve has a stationary point at a frequency which we denote $\tau$,
whereas the downstream dispersion curve monotonically increases with
$k_{\mathrm{d}}$, as can be seen in panel (b). This observation will be
important when we describe the construction of $s_{\mathrm{ess}}$ below.

To calculate $s_{\mathrm{ess}}$ numerically, we set $\phi=-L$ at the upstream
boundary $\Gamma_{3}$ and allow the upstream height to come as part of the
solution (in so doing we are able to capture waves upstream), and we ‘pin’ the
downstream height at the outflow boundary $\Gamma_{1}$. Since we do not
constrain $\phi$ at the outflow boundary this allows us to capture waves
downstream. We filter out spurious essential modes by removing modes which
have $|\mbox{Re}(s_{\mathrm{ess}})|>10^{-7}$. We note that we do not apply the
sponge layer in this calculation, i.e. we set $A_{\pm}=0$.

The essential spectrum comprises two types of modes, denoted types i and ii.
Type i modes occur when $|\mathrm{Im}(s_{\mathrm{ess}})|>|\tau|$. A set of
calculations associated with a typical member of this class is shown in figure
5. In panel (a) $s_{\mathrm{ess}}$ is shown on the imaginary axis in the
complex $s$-plane for the underlying steady-state shown in panel (b) for
$a=0.01,Fr=0.8823$. In panel (a) we highlight a particular member of
$s_{\mathrm{ess}}$, denoted $\lambda$, which satisfies
$|\mathrm{Im}(\lambda)|>|\tau|$ (values quoted in the caption). The
corresponding mode (real part) is shown in panel (c) and comprises two
distinct wave patterns upstream and downstream, denoted $g_{\mathrm{u}}$ and
$g_{\mathrm{d}}$ and highlighted in panels (d) and (e) respectively. The
dominant wave numbers of the type i modes can be determined, numerically, by
calculating the power spectra of $g_{\mathrm{u}}$ and $g_{\mathrm{d}}$ as a
function of the wave number $k$, as shown in panel (g). We calculate the power
spectra using the fast Fourier transform (fft) of $g_{\mathrm{u}}$ and
$g_{\mathrm{d}}$. For type i modes there is a single dominant wave number for
each of the upstream and downstream signals as shown by the peak in each of
the power spectra in panel (g). The wave numbers associated with these peaks
correspond precisely to the intersections of the upstream and downstream
dispersion curves with the horizontal line $\omega=+|\mathrm{Im}(\lambda)|$,
as shown in panel (f), which has the same horizontal axis and scale as panel
(g) to aid this comparison. This illustrates the excellent agreement between
the theory and numerics and gives us confidence in the fidelity of the
numerical eigensolver.

Type ii modes occur when $|\mathrm{Im}(s_{\mathrm{ess}})|<|\tau|$. A set of
calculations associated with a typical member of this class, denoted by
$s_{\mathrm{ess}}=\mu$, is shown in figure 6. This figure follows an identical
structure to figure 5. In this class of modes, the upstream section of the
mode, shown in panel (d), is clearly multi-harmonic. This is a direct
consequence of the fact that in the upstream dispersion curve, for a given
$|\omega|<|\tau|$, multiple values of $k_{\mathrm{u}}$ satisfy (25). For type
ii modes there are three dominant wave numbers upstream and one dominant wave
number downstream. As for the type i modes, the wave numbers associated with
the peaks in the power spectra correspond precisely to the intersections of
the upstream and downstream dispersion curves, but this time with the
horizontal lines $\omega=\pm|\mathrm{Im}(\mu)|$, as shown in panels (f) and
(g). For type ii modes the wave numbers associated with
$\omega=-|\mathrm{Im}(\mu)|$ will travel (as $t$ progresses) slower than the
uniform upstream speed while wave numbers associated with
$\omega=+|\mathrm{Im}(\mu)|$ will travel faster.

Figure 5: Type i modes of the numerically calculated $s_{\mathrm{ess}}$ for
the hydraulic-fall solution, $Fr=0.8823$, $a=0.01$, $b=0.3$. Panel (a): The
numerically computed $s_{\mathrm{ess}}$ is shown with blue markers on the
imaginary $s$ axis, and a particular element of
$s_{\mathrm{ess}}=\lambda=0.1309\mathrm{i}$ is highlighted with a solid red
marker. The horizontal dotted lines indicate the levels
$s=\pm\mathrm{i}|\tau|=\pm 0.04489\mathrm{i}$; in this calculation
$|\mathrm{Im}(\lambda)|>|\tau|$. Panel (b): The underlying steady state. Panel
(c): The real part of the eigenmode associated with $\lambda$. Panels (d),
(e): The upstream/downstream portion of the real part of the eigenmode,
denoted $g_{\mathrm{u}}/g_{\mathrm{d}}$, respectively. Panel (f): The
downstream and upstream dispersion relation given in (24), (25) respectively
with the minus sign taken in both cases. The dashed horizontal lines indicate
$\omega=\pm|\mathrm{Im}(\lambda)|$. Panel (g): The power spectrum, i.e.
$|\texttt{fft}(g)|^{2}$ the square of the Fourier transform. The horizontal
axes of panels (f) and (g) are identical so a direct comparison between the
peaks of the power spectrum and the intersection of the dispersion curves with
$\omega=\mathrm{Im}(\lambda)$ can be made. We also note that in panels (d)-(g)
all calculations corresponding to the upstream section are shown as solid blue
lines while the downstream section are dotted yellow lines. Figure 6: Type ii
modes of the numerically calculated $s_{\mathrm{ess}}$ for the hydraulic-fall
solution with $Fr=0.8823$, $a=0.01$, $b=0.3$. The figure follows the same
description as figure 5 with the exception that in this calculation we
highlight a different member of $s_{\mathrm{ess}}=\mu=0.03027\mathrm{i}$ such
that $|\mathrm{Im}(\mu)|<|\tau|$. Here $\tau$ has the same value as in figure
5.

####  III.2.2 The point spectrum, $s_{\mathrm{p}}$

Throughout the following discussion we will use $\nu$ to represent the
eigenvalue in the point spectrum, $s_{\mathrm{p}}$, that lies in the first
quadrant of the complex plane. As was mentioned earlier, the point spectrum
has a four-fold symmetry in the complex plane. In contrast to the forced
solitary-wave fKdV calculations of [11] and [27], where the point spectrum
eigenmodes are even-symmetric so that
$g_{\hat{\eta},\nu}(-x)=g_{\hat{\eta},-\nu}(x)$, hydraulic-fall solutions are
not even-symmetric and so we expect that $g_{\hat{\eta},\nu}(-x)\neq
g_{\hat{\eta},-\nu}(x)$.

The numerical computation of the eigenmodes requires a great deal of care to
filter out grid-dependent spurious modes. A very wide computational domain is
needed because evanescent waves on the upstream or downstream side of the
eigenmodes decay extremely slowly (this is quantified below). We identify two
different types of mode in the point spectrum which we term type iii
corresponding to eigenvalues in the left-half plane (here
$s_{\mathrm{p}}=-\nu,\,-\nu^{*}$), and type iv corresponding to eigenvalues in
the right-half plane (here $s_{\mathrm{p}}=\nu,\,\nu^{*}$). To compute type
iii modes, we impose $y_{f}=\gamma_{s}$ and $\phi(L,y)=L$ at the outflow
boundary $\Gamma_{1}$ so that the profile is flat downstream, but there are
evanescent waves upstream. For type iv modes, we set $y_{f}=1$ and
$\phi(-L,y)=-L$, at the inflow boundary $\Gamma_{3}$ to capture eigenmodes
that are flat upstream and have evanescent waves downstream.

In figure 7 we show results for the negative forcing case $a=-0.01$,
$Fr=0.9659$. In panel (a) the point spectrum $s_{\mathrm{p}}$ (green stars) is
shown together with the essential spectrum $s_{\mathrm{ess}}$ (blue circles)
in the complex $s$-plane. The underlying steady state is shown in panel (b).
Panels (c) and (d) show the corresponding iii modes, and panels (e) and (f)
show the corresponding iv modes. Note that the sponge layer is not included
for either the type iii or the type iv calculations. The exponential decay
rate of the type iv downstream waves is scrutinised in panel (g) by plotting
the logarithm of the local wave maxima against $x$ and then estimating its
slope. The decay rate is computed to be $-3.3\times 10^{-4}$.

We only find a non-empty point spectrum when $a<0$ (we calculated the
essential and point spectrum for values of $a,b$ in the range $a^{*}<a<0.1$
and $0.3\leq b\leq 3.0$) and, moreover, that the point spectrum contains only
four eigenvalues. We can investigate the linear instability for $a<0$ further
by computing the point spectrum as $a$ is varied continuously. Figure 8 shows
that as $a\to 0^{-}$ the quartet of eigenvalues appears to approach the
origin. This is consistent with the fact that when $a=0$, i.e. for a flat
bottom, the spectrum contains an isolated zero of multiplicity four [see, for
example 1]. As $a\to a^{*}\approx-0.016$, the unstable eigenmodes approach
those found at the ‘termination point’ solutions discussed in § III.1. The
fact that we have not found a point spectrum when $a>0$ appears to suggest
that hydraulic-fall solutions over positive bumps are linearly stable.

Although we have analysed the stability properties of the linearised system,
as is well known linear stability does not imply nonlinear stability [see, for
example 26]. In the next section we shall perform numerical simulations to
provide evidence for nonlinear stability.

Figure 7: Type iii, iv modes of the numerically calculated $s_{\mathrm{p}}$
for the hydraulic-fall solution with $Fr=0.9659$, $a=-0.01$, $b=0.3$. Panel
(a): $s_{\mathrm{ess}}$ is marked on the imaginary axis and $s_{\mathrm{p}}$
with $\nu=2.375\times 10^{-5}+0.0235\mathrm{i}$ is indicated with green stars.
Panel (b): The underlying steady state. Panels (c)-(d): Type iii eigenmodes
corresponding to $s_{\mathrm{p}}=-\nu,-\nu^{*}$ (real and imaginary parts).
Panels (e)-(f): Type iv eigenmodes corresponding to
$s_{\mathrm{p}}=\nu,\nu^{*}$. Panel (g)
$\log(\max|\mathrm{Re}(g_{\hat{\eta},\nu})|)$ plotted against $x$ indicating
the small spatial decay rate (estimated in the panel) of the downstream wave
in panel (e). Figure 8: Hydraulic-fall solutions for $b=0.3$. The numerically
calculated $s_{\mathrm{p}}$ as $a$ is varied (the corresponding values of $Fr$
can be inferred from figure 3). The solid black lines indicate
$s_{\mathrm{p}}$ and the arrows indicate how $s_{\mathrm{p}}$ evolves as $a$
decreases from zero. $s_{\mathrm{ess}}$ is shown by hollow circular markers.
We note that the eigenmodes of $s_{\mathrm{ess}}$ change as $a$ is varied but
$s_{\mathrm{ess}}$ itself remains on the imaginary axis. The inset diagrams
show eigenmodes corresponding to $s_{\mathrm{p}}$ as indicated by arrows for
the values of $a=-10^{-5}$ (bottom inset panels) and $a=-0.016\approx a^{*}$
(top inset panels).

###  III.3 Nonlinear Stability

In this section we solve an initial value problem in order to confirm the
results of the previous section and to probe the nonlinear stability
properties of the hydraulic-fall solutions. We introduce a localised
perturbation to the steady solution, $y_{s}(x)$, taking

$y_{f}(x,t=0)=y_{s}(x)+\varepsilon x\mbox{e}^{-a_{2}^{2}(x-a_{3})^{2}}.$ (27)

for some choice of the perturbation amplitude $\varepsilon$ and the constants
$a_{2}$ and $a_{3}$. In the results to be presented below we choose
$a_{2}=0.5$ in (7) and set $a_{3}=0$ but the general conclusion on the
solution’s stability remains the same for other values of $a_{3}$. The
computation is carried out as follows. Having deformed the flow domain to
incorporate the perturbation (27), we solve Laplace’s equation over this
domain together with the steady kinematic condition on the free surface
$\Gamma_{2}$, (4), the no-penetration condition on the bottom $\Gamma_{0}$,
(3), and the inflow and outflow conditions on $\Gamma_{3}$ and $\Gamma_{1}$,
(8) and (9), respectively. This provides us with the initial velocity
potential over the perturbed domain, namely $\phi(\textbf{x},t=0)$. We set the
downstream target level for the sponge layer as $y_{+}=\gamma_{s}$.

Our stability analysis of the hydraulic-fall solutions depends on two
parameters: the amplitude of the topographic forcing, $a$, and the size of the
perturbation, $\varepsilon$, in (27). To standardise the size of the
perturbation as $a$ is varied we set

$\varepsilon=l\Delta_{s},\qquad\Delta_{s}\equiv|1-\gamma_{s}|,$ (28)

with $l>0$. This provides a convenient way to compare results for different
values of $a$.

As the system is Hamiltonian the underlying dynamical system cannot contain
attractors or repellors as the Hamiltonian is invariant along any trajectory
in the system’s phase space. Therefore, in this system, the stability
properties of the steady states are interpreted by a local argument; the
system can evolve towards/away from a steady state as energy is emitted far
upstream/downstream.

When discussing the results of our numerical simulations, it will be useful to
refer to the value of the ‘wave resistance coefficient’ that was introduced by
[11] for the fKdV equation, and which is defined to be

$C_{r}=-\int_{\Gamma_{2}}\frac{\partial y_{f}(x,t)}{\partial
x}y_{b}(x)\,\mbox{d}S.$ (29)

In the fKdV model $C_{r}$ is the power supplied by the bottom forcing, $y_{b}$
[11].

####  III.3.1 Positive forcing

Figure 9:  Perturbation of a hydraulic fall over positive forcing with
$a=0.01$ and $Fr=0.88228$ and with initial condition (27) with $l=2.0$,
$b=0.3$, $a_{2}=0.5$, and $a_{3}=0$. (a) Waterfall plot of free-surface
evolution. (b) The initial perturbation (solid red line) and the steady state
(dotted black line). (c) The free-surface (solid red line) at $t=1247$ and the
steady state (dotted black line).

We start by examining hydraulic falls over a positive forcing with $a>0$.
Figure 9 shows the results of a time-dependent calculation using the initial
condition (27) with $l=2.0$ (other parameter values are quoted in the
caption). An animation of the time evolution of the free surface can be found
in the supplementary material (movie_1.mp4). Evidently the system quickly
returns to the original unperturbed stable state. It was found that this
behaviour was typical for forcing amplitudes in the range $a=(0,0.1]$ and
perturbations in the range $l\in(0,4]$, providing evidence that the solution
is stable to localised perturbations as prescribed in (27) for all $a>0$
sampled.

####  III.3.2 Negative Forcing

Figure 10 shows the result of a time-dependent calculation for a topographical
dip with $a=-0.01$, $b=0.3$ and initial condition (27) (the other parameter
values are shown in the captions). We strongly recommend the reader to watch
movie_2.mp4 in the supplementary material to support the following discussion.
A waterfall plot of the time-evolution of the free surface is shown in panel
(g). Panel (c) shows the time-trace of $y_{f}(0,t)$ over the range
$0<t<10^{5}$ (zoomed-in regions of this time-trace are shown in panels (a) and
(e)), and panel (d) shows the $(C_{r},y_{f}(0,t))$ phase plane projection. The
thicker lines in panel (d) indicate regions in phase space that are visited
more frequently. Wave profiles at times $t=5000$ and $t=10^{5}$ are shown in
panels (f),(b), respectively.

Figure 10: Perturbation of the hydraulic fall solution with $a=-0.01$,
$Fr=0.9659$ (and $\Delta_{s}=0.0450$) with the initial condition (27) taking
$l=2.0$, $b=0.3$, $a_{2}=0.5$, and $a_{3}=0$. The underlying type iv eigenmode
of the steady state has $\nu=2.375\times 10^{-5}+0.0235\mathrm{i}$. (a) Time
trace of $y_{f}(0,t)$ over the window $8.6\times 10^{4}<t<9.15\times 10^{5}$.
(b) Free-surface profile at $t=10^{5}$. The estimation of the wave number,
$k$, is indicated. (c) Time trace of $y_{f}(0,t)$ over $0<t<10^{5}$. The
dotted line indicates $y_{s,f}(0)$. We note the frequency of the oscillations
cannot be resolved on this time-scale and zoomed regions are shown in panels
(a) and (e). (d) The phase plane projection $(C_{r},y_{f}(0,t))$ with the
vertical axis on the same scale as panel (c). Thicker lines indicate locations
in phase space more frequently visited. (e) Time trace of $y_{f}(0,t)$ over
$0<t<5000$. (f) Free-surface profile at $t=5000$. The arrows indicate the
direction of wave propagation. (g) Waterfall plot of free-surface evolution.

Considering these results, the time-evolution can be broadly divided into two
distinct stages. The initial stage (over $0<t<5000$, say, and shown in panel
(e)) is characterised by a roughly sinusoidal temporal oscillation of
$y_{f}(0,t)$ about the steady state value, the latter shown as a dotted line
in panel (c). In this stage, the transient growth is well matched by the real
part of the unstable eigenvalue, $s=\nu$, of the underlying state as shown by
the dashed line in panel (c). The frequency of the oscillation is in excellent
agreement with the imaginary part of $\nu$. In fact, the same observation
holds for a broad range of topography amplitudes $a$. This is illustrated in
figure 11(a), where $\mbox{Im}(\nu)$ is plotted against $a$. Also shown is the
frequency of $y_{f}(0,t)$ estimated from the time-trace during the initial
stage of the evolution. Evidently there is excellent agreement between the
two. We also found that this frequency is independent of the perturbation
amplitude $l$.

Figure 11: (a) The frequency of oscillations in the initial stage of the
evolution for the case studied in figure 10 but over a range of $a$ values,
$b=0.3$. The solid lines are the imaginary part of $\nu$ and the line with
circular markers is the estimate of the frequency from the time-dependent
simulations. (b) The linear dispersion relation $\omega=\omega(k)$, with
$k=k_{d}$, according to (24), shown with a dot-dashed curve. The estimated
frequency-wave number relation for the small amplitude waves emitted in panel
(f) of figure 10 is illustrated with red circular markers. The (positive)
imaginary part of the eigenvalue in the point spectrum, $\mbox{Im}(\nu)$, is
shown with a solid blue line.

During this initial stage of the evolution, as the amplitude of the
oscillations increases and small amplitude waves are emitted and travel
downstream, as can be seen in panel (f) of figure 10. We would expect that the
wave number of these emitted waves conforms with the linear dispersion
relation (24). We can test this by estimating the frequency and wave number of
the emitted waves from the time signals and from the instantaneous free-
surface profile, respectively. To approximate the frequency we calculated the
mean separation of the local maxima of the time-signal $y_{f}(0,t)$, and to
estimate the wave number we extracted the mean separation of the local maxima
in the spatial wave profile $y_{f}(x,t)$ at $t=10^{4}$. The results are shown
in figure 11(b), where we see good agreement with the dispersion relation
$\omega=\omega(k)$, with $k=k_{d}$, in (24), although the disparity between
the two increases as the wave number $k$ increases. Also shown in the figure
11(b) is the (positive) imaginary part of the eigenvalue in the corresponding
point spectrum (see section III.2.2).

Eventually the amplitude of the oscillations saturate heralding the second
evolutionary stage, which we refer to as the limiting stage. Here the
oscillations are less sinusoidal in character. The limiting stage is
characterised by a regular stream of larger-amplitude cnoidal waves
propagating in the upstream direction and away from the topographical dip,
and, in contrast to the initial stage, the waves travelling downstream, as
shown in panel (b) of figure 10, now consist of multi-harmonic waves. The
phase-plane projection $(C_{r},y_{f}(0,t))$ shown in panel (d) of the same
figure has a distinctive guitar-pick shape. At large time $t$ the trajectory
appears to have converged to a closed loop (shown in white with arrows), which
represents an invariant solution of the system. We remark that this limiting
behaviour was found to be typical for parameter values $a^{*}<a<0$, $0.3\leq
b\leq 3.0$. We also recommend the reader to watch the animation in movie_3.mp4
focussing on how the limiting stage for figure 10 evolves over a temporal
period.

In conclusion the steady state hydraulic-fall solution with negative forcing
is unstable and perturbations from it settle down to a time-dependent,
nonlinear stable invariant solution irrespective of size of the perturbation
in (27). We discuss this in more detail in § III.4.

####  III.3.3 Hydraulic-rises

For completeness, we briefly describe the nonlinear stability of the steady
hydraulic-rise solutions. Hydraulic rises are unstable for both positive and
negative $a$ as demonstrated in figures 12 and 13, respectively. For positive
forcings, as the initial disturbance propagates downstream a solitary-wave
emanates from the origin and moves upstream whilst a steady cnoidal-wave
pattern emerges downstream. Ignoring the solitary wave, this structure is
known as a generalised steady hydraulic rise (see figure 12). For negative
forcings, similar to what was found for the hydraulic fall, the system settles
into a time-periodic orbit as indicated in panel (c) of figure 13. We observe
that the temporal signal of $y_{f}(0)$ in the limiting stage, as shown in
panel (b), is more sinusoidal than that for the hydraulic fall (see panel (a)
of figure 10), and that a similar cnoidal wave pattern to that seen for the
positively forced hydraulic rise emerges downstream.

Figure 12: Waterfall plot of free-surface evolution of a perturbation of a
hydraulic rise over positive forcing with $a=0.01$ and $Fr=1.1173$ and with
initial condition (27) with $l=2.0$, $b=0.3$, $a_{2}=0.5$, and $a_{3}=0$.
Figure 13: Hydraulic rise. The case $a=-0.01$, $b=0.3$, $Fr=1.0436$ and with
initial condition (27) with $l=2.0$, $a_{2}=0.5$, and $a_{3}=0$. (a) Waterfall
plot for the system starting from a state of near uniform flow. (b) Time-trace
of $y_{f}(0,t)$ in the limiting stage. (c) The projected phase plane,
$(C_{r},y_{f}(0,t))$. The initial condition is shown as a circular marker, and
the trajectory is in orange.

###  III.4 Time-dependent invariant solution

Figure 14: Hydraulic fall. The case $a=-0.01$, $b=0.3$, $Fr=0.9659$. (a)
Waterfall plot for the system starting from a state of near uniform flow. (b)
Time-trace of $y_{f}(0,t)$. (c) Time-trace in the period $9.5\times
10^{4}<t<10^{5}$. (d) The projected phase plane, $(C_{r},y_{f}(0,t))$. The
initial condition is shown as a circular marker, the system trajectory is in
red, and the limiting orbit from figure 10 is shown in white.

To further explore the existence of a time-dependent invariant solution for
negative forcing, $a<0$, we calculate the free-surface response when the
system starts from a condition of near uniform flow with a flat free surface.
Specifically, to provide the initial condition we solve Laplace’s equation
over the flow domain shown in figure 2 but with a flat free surface
$\Gamma_{2}$, enforcing the steady form of the free-surface kinematic
condition, (4), the no-penetration condition on the bottom $\Gamma_{0}$, (3),
and the inflow and outflow conditions on $\Gamma_{3}$ and $\Gamma_{1}$, (8)
and (9), respectively. This yields the initial velocity potential
$\phi(\textbf{x},t=0)$. This starting condition is artificial as the pressure
at the free surface does not match the constant ambient pressure at $t=0$, but
we overlook this slight inconsistency. During the calculation we set the
downstream sponge target level as $y_{+}=1$.

The results are shown in figure 14 for $a=-0.01$ and $Fr=0.9659$, which are
the same values as in figure 10. The waterfall plot shown in panel (a)
indicates that the region left in the wake of the downstream travelling waves
is itself wavy. The signal $y_{f}(0,t)$ shown in panels (b),(c) of figure 14
appears to be periodic and it is similar in character to the oscillations
observed in figure 10(a). Furthermore, in panel (d) of figure 14, the
trajectory in the $(C_{r},y_{f}(0,t))$ plane is shown and it is clear that the
system is evolving towards the guitar-pick shaped limiting orbit observed in
figure 10, which is also shown in the panel (coloured white and labelled as
the ‘limiting orbit’). This calculation provides further evidence of what
appears to be a time-periodic orbit and that it is the same solution as in
figure 10.

In figure 15 we examine the time-dependent invariant solution further by
analysing the calculation in figure 10 for $a=-0.01$. We begin by noting that,
as can be seen in figure 15(a), the mean level of the downstream spatial
signal, denoted $\overline{y_{\mathrm{d}}}$, decreases as time progresses.
Consequently, a modified downstream dispersion relation obtains for the
limiting stage, and this is found by replacing $\gamma_{s}$ and $V$ in (24)
with $\overline{y_{\infty}}$ and $1/\overline{y_{\infty}}$ respectively.

Next we observe that the average temporal frequency, denoted
$\overline{\omega}$, estimated using the local maxima of the time signal
$y_{f}(0,t)$ in figure 15(c), varies in time from a value close to
$\overline{\omega}\approx\mathrm{Im}(\nu)$ in the initial stage of evolution
to a smaller value, denoted $\omega_{1}$, in the limiting stage (see panel (b)
of figure 15). This frequency modulation can be explored further by
calculating the power spectrum of the temporal signal of $y_{f}(0,t)$ in the
limiting stage, i.e. over $10^{5}<t<1.5\times 10^{5}$. The result is shown in
panel (e) of figure 15. This shows that there are actually two main temporal
frequencies in the limiting stage as indicated by the peaks in the temporal
power spectra shown in panel (e) of figure 15. We denote these two dominant
frequencies by $\omega_{1}\approx 0.018$ and $\omega_{2}\approx 0.036$.

In addition, we find there are two dominant wave numbers of the spatial
downstream wave signal, $y_{\mathrm{d}}$ (see panel (h) of figure 15), as
indicated by the spatial power spectrum of $y_{\mathrm{d}}$ at $t=10^{5}$
shown in panel (d) of figure 15. We denote these two dominant wave numbers by
$k_{1}\approx 0.2275$ and $k_{2}\approx 0.3850$.

Bringing this all altogether, in panel (f) of figure 15 we show the limiting
dispersion curve together with two horizontal lines indicating the levels
$\omega=\omega_{1},\omega_{2}$. The values of $k$ where the curve intersects
these lines correspond precisely to the two dominant wave numbers $k_{1}$,
$k_{2}$ seen in panel (d) of figure 15. To aid comparison across panels in the
figure the horizontal axes in panel (d) and panel (f) of figure 15 are
identical, as are the vertical axes of panel (e) and panel (f). These results
indicate that the wave numbers of the downstream disturbance are linked to the
temporal frequencies via the downstream dispersion curve.

In figure 16 we show the power spectrum of the time-signal of $y_{f}(0,t)$ at
the limiting stage for a slightly different case with $a=-0.005$, $b=0.3$. The
energy peak at the dominant frequency $\omega_{1}$ is evident, with other
peaks of diminishing strength located approximately at $\omega=n\omega_{1}$,
where $n=2,3,4,\cdots$, presumably generated via nonlinear effects. This is
strong numerical evidence that the invariant solution is a periodic orbit. The
downstream wave numbers associated with the subharmonic frequencies,
$n\omega_{1}$, must satisfy the linear dispersion relation (24) and, in
general, will be incommensurate. This explains the rather irregular appearance
of the downstream waves in panel (h) of figure 15.

Figure 15: Frequency and wave number analysis of the downstream portion of the
wave, $y_{\mathrm{d}}$, and the time-signal $y_{f}(0,t)$, for the calculation
in figure 10. Panel (a): the downstream mean level,
$\overline{y_{\mathrm{d}}}$ as a function of time. Panel (b) The average
frequency as $t$ varies, estimated from the local maxima of the time signal of
$y_{f}(0,t)$. Panel (c): the time signal of $y_{f}(0,t)$. Panel (d): the power
spectrum of the down stream wave in the limiting stage of the evolution with
the two peaks occurring at $k_{1}\approx 0.2275,k_{2}\approx 0.3850$. Panel
(e): the power spectra of $y_{f}(0,t)$, note the vertical axes is the
frequency, $\omega$, and the horizontal axes is the power. Panel (f) The
limiting stage dispersion curve together with the levels
$\omega=\omega_{1}\approx 0.018,\omega_{2}\approx 0.036$. Note the horizontal
axis and scale are identical in panels (d) and (f) and the vertical axis and
scale are identical in panels (e) and (f) to aid comparison. Panels (g) and
(h): the wave profile at the given time slot and the downstream portion of the
free-surface. Figure 16: The power spectrum of $y_{f}(0,t)$ in the limiting
stage of the evolution for the solution when $a=-0.005,b=0.3$, $l=2$. The
primary harmonic, $\omega_{1}\approx 0.01257$, and the next four subharmonic
frequencies are identified.

In summary, in this subsection we have obtained the most important result of
the present article, namely the apparent existence of a temporally periodic
orbit of the fully nonlinear system. Time-periodic behaviour in free-surface
flow over topography has been observed before, notably for the forced KdV
equation [11], [43] and [22]. Through unsteady simulations these authors found
a similar oscillatory behaviour of the free-surface over the topography. In
this case the period of the oscillation was linked to the frequency at which
solitary waves were emitted in the upstream direction. However, in none of
these studies was the long time behaviour explored.

##  IV Conclusion

We have examined the stability of a hydraulic-fall flow over a localised
bottom topography. The topography is in the form of a gaussian whose amplitude
is either positive (a bump) or negative (a dip). In the initial part of our
analysis we studied the linear stability of these solutions by first
demonstrating that the full equations can be cast in the form of a canonical
Hamiltonian system with a Hamiltonian that is a modification of that derived
by [44] and [14] for the classical water-wave problem. As a consequence the
spectrum of the system linearised about a steady hydraulic-fall solution has a
four-fold symmetry in the complex plane.

In the case of a bump, our numerical computations showed that the point
spectrum is either empty or else all of its eigenvalues are located on the
imaginary axis so that the underlying flow is spectrally stable. For a dip we
identified a single eigenvalue in the point spectrum (to within the four-fold
symmetry) that lies off the imaginary axis with a positive real part, meaning
that the flow is linearly unstable. For both the bump and the dip the
essential spectrum occupies the whole of the imaginary axis. For both the bump
and the dip we carefully analysed the associated eigenfunctions and we
demonstrated that their behaviour in the far-field, away from the localised
topography, could be reconciled with the dispersion relation for small
amplitude waves over a flat bottom.

In the second part of our work we examined the nonlinear stability properties
by solving an initial value problem for the full Euler equations, with the
initial condition taken to represent a small perturbation about a steady
hydraulic fall. For a bump the simulations showed that the disturbance to the
free surface disperses so that the steady solution is ultimately recaptured.
This is in keeping with the prediction of the linear stability analysis.
Consequently, we concentrated most of our attention on the dip. In this case
our numerical simulations showed that the steady hydraulic fall is also
nonlinearly unstable and that, ultimately, the flow approaches a periodic
state in which the free surface exhibits a localised peak that pulsates up-
and-down over the dip, while small amplitude waves are emitted continuously in
both the upstream and the downstream directions. This periodic state may be
interpreted as an invariant solution of the underlying dynamical system. We
have shown numerically that the periodic waves that are emitted downstream
have a dominant temporal frequency that can be accurately linked to the linear
dispersion relation based on the mean downstream depth. The subharmonic
frequencies of these waves were found to be rational multiples of the dominant
frequency, thereby reinforcing our conclusion that a periodic state has been
achieved.

Interesting questions remain. For example, beside periodic solutions as
identified here, there exists the possibility of quasiperiodic solutions. Such
objects exist in phase space as so-called invariant tori [e.g. 30] . However,
directly calculating time-dependent invariant solutions, such as periodic-
orbits or invariant tori, remains a highly non-trivial numerical task. In
other applications in fluid dynamics, periodic orbits have been computed semi-
analytically using a weakly nonlinear approach, e.g. for air bubbles in a
Hele-Shaw channel [29]. [33, 34] have recently applied Koopman analysis and
dynamic mode decomposition to approximate periodic orbits for Burger’s
equation and the Navier-Stokes equations. It would be of interest to apply
such methods to the free-boundary problem considered here in order to further
probe the landscape of invariant solutions in phase space.

Finally, as was noted in the Introduction, there is an apparent lack of
experimental results in the literature concerning hydraulic falls over a dip.
It would be of interest to see if the numerical observations in the present
work, and in particular the time-periodic flow that appears at large time, can
be realised in the laboratory.

## References

  * Akers & Nicholls [2012] Akers, B. & Nicholls, D. P. 2012 Spectral stability of deep two-dimensional gravity water waves: repeated eigenvalues. SIAM J. Appl. Math. 72 (2), 689–711.
  * Akylas [1984] Akylas, T. R. 1984 On the excitation of long nonlinear water waves by a moving pressure distribution. J. Fluid Mech. 141, 455–466.
  * Alias [2014] Alias, A. B. 2014 Mathematical Modelling of Nonlinear Internal Waves in a Rotating Fluid. PhD Thesis.
  * Bermúdez et al. [2007] Bermúdez, A., Hervella-Nieto, L., Prieto, A. & Rodrı, R. et al. 2007 An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems. J. Comp. Phys. 223 (2), 469–488.
  * Binder [2019] Binder, B. J. 2019 Steady two-dimensional free-surface flow past disturbances in an open channel: Solutions of the Korteweg-de Vries equation and analysis of the weakly nonlinear phase space. Fluids 4, 1–24.
  * Binder et al. [2013] Binder, B. J., Blyth, M. G. & McCue, S. W. 2013 Free-surface flow past arbitrary topography and an inverse approach to wave-free solutions. IMA J. Appl. Math. 78, 685–696.
  * Binder et al. [2008] Binder, B. J., Dias, F. & Vanden-Broeck, J-M. 2008 Influence of rapid changes in a channel bottom on free-surface flows. IMA J. Appl. Math. 73, 254–273.
  * Binder et al. [2005] Binder, B. J., Vanden-Broeck, J-M. & Dias, F. 2005 Forced solitary waves and fronts past submerged obstacles. CHAOS 15, 037106 1–13.
  * Boyd [2000] Boyd, J. P. 2000 Chebyshev and Fourier Spectral Methods. Dover Publications.
  * Buttle et al. [2018] Buttle, N. R., Pethiyagoda, R., Moroney, T. J. & McCue, S. W. 2018 Three-dimensional free-surface flow over arbitrary bottom topography. J. Fluid Mech. 846, 166–189.
  * Camassa & Wu [1991] Camassa, R. & Wu, T. Y-T. 1991 Stability of forced steady solitary waves. Phil. Trans. R. Soc. Lond. 10, 429–466.
  * Chardard et al. [2011] Chardard, R., Dias, F., Nyguyen, H. Y. & Vanden-Broeck, J-M. 2011 Stability of some stationary soutions to the forced KdV equation with one or two bumps. J. Eng. Math. 70, 175–189.
  * Choi & Kim [2016] Choi, H. & Kim, H. 2016 The hyperbolic relaxation systems for the forced KdV equations with hydraulic falls. Euro. J. Mech.-B/Fluids 58, 20–28.
  * Craig & Sulem [1993] Craig, W. & Sulem, C. 1993 Numerical simulation of gravity waves. J. Comp. Phys. 108 (1), 73–83.
  * Dias & Vanden-Broeck [1989] Dias, F. & Vanden-Broeck, J.-M. 1989 Open channel flows with submerged obstructions. J. Fluid Mech. 206, 155–170.
  * Dias & Vanden-Broeck [2002] Dias, F. & Vanden-Broeck, J. M. 2002 Generalised critical free-surface flows. J. Eng. Math. 42, 291–302.
  * Donahue & Shen [2010] Donahue, A. S. & Shen, S. S-P. 2010 Stability of hydraulic falls and sub-critical cnoidal waves in water flows over a bump. J. Eng Maths 10, 1007/s.
  * Forbes [1988] Forbes, L. K. 1988 Critical free-surface flow over a semi-circular obstruction. J. Eng. Math. 22 (1), 3–13.
  * Forbes & Schwartz [1982] Forbes, L. K. & Schwartz, L. W. 1982 Free-surface flow over a semicircular obstruction. J. Fluid Mech. 114, 209–314.
  * Forbes et al. [2021] Forbes, L. K., Walters, S. J. & Hocking, G. C. 2021 Ideal planar fluid flow over a submerged obstacle: review and extension. ANZIAM 63 (4), 377–419.
  * Grimshaw & Maleewong [2016] Grimshaw, R. H. J. & Maleewong, M. 2016 Transcritical flow over two obstacles: Forced Korteweg–de Vries framework. J. Fluid Mech. 809, 918–940.
  * Grimshaw & Smyth [1986] Grimshaw, R. H. J. & Smyth, N. 1986 Resonant flow of a stratified fluid over topography. J. Fluid. Mech. 169, 429–464.
  * Grimshaw et al. [2007] Grimshaw, R. H. J., Zhang, D-H. & Chow, K. W. 2007 Generation of solitary waves by transcritical flow over a step. J. Fluid Mech. 537, 235–254.
  * Heil & Hazel [2006] Heil, M. & Hazel, A. L. 2006 oomph-lib–an object-oriented multi-physics finite-element library. In Fluid-structure interaction, pp. 19–49. Springer.
  * Heroux et al. [2003] Heroux, M., Bartlett, R., Howle, V., Hoekstra, R., Hu, J., Kolda, T., Lehoucq, R., Long, K., Pawlowski, R., Phipps, E., Salinger, A., Thornquist, H., Tuminaro, R., Willenbring, J. & Williams, A. 2003 An overview of trilinos. Tech. Rep.. SAND2003-2927. Sandia National Laboratories.
  * Holm et al. [1985] Holm, D. D., Marsden, J. E., Ratiu, T. & Weinstein, A. 1985 Nonlinear stability of fluid and plasma equilibria. Phys. Rep. 123 (1-2), 1–116.
  * Keeler et al. [2017] Keeler, J. S., Binder, B. J. & Blyth, M. G. 2017 On the critical free-surface flow over localised topography. J. Fluid Mech. 832, 73–96.
  * Keeler et al. [2021] Keeler, J. S., Blyth, M. G. & King, J. R. 2021 Termination points and homoclinic glueing for a class of inhomogeneous nonlinear ordinary differential equations. Nonlinearity 34, 532–561.
  * Keeler et al. [2019] Keeler, J. S., Thompson, A. B., Lemoult, G., Juel, A. & Hazel, A. L. 2019 The influence of invariant solutions on the transient behaviour of an air bubble in a hele-shaw channel. Proc. R. Soc. Lond. A 879, 1–27.
  * Kuznetsov [1998] Kuznetsov, Y. A. 1998 Elements of Applied Bifurcation Theory (3rd Ed.). Springer-Verlag.
  * Li [2005] Li, X. S.x 2005 An overview of superlu: Algorithms, implementation, and user interface. ACM Transactions on Mathematical Software (TOMS) 31 (3), 302–325.
  * Page & Părău [2014] Page, C. & Părău, E. I. 2014 Time dependent hydraulic falls and trapped waves over submerged obstacles. Phys. Fluids 26, 126604.
  * Page & Kerswell [2018] Page, J. & Kerswell, R. R. 2018 Koopman analysis of Burgers equation. Phys. Rev. Fluids 3 (7), 071901.
  * Page & Kerswell [2019] Page, J. & Kerswell, R. R. 2019 Koopman mode expansions between simple invariant solutions. J. Fluid Mech. 879, 1–27.
  * Romate [1991] Romate, J. E. 1991 Absorbing boundary conditions for free surface waves. J. Comp. Phys. 99, 135–145.
  * Sandstede & Scheel [2000] Sandstede, B. & Scheel, A. 2000 Absolute and convective instabilities of waves on unbounded and large bounded domains. Physica D 145 (3-4), 233–277.
  * Tam et al. [2015] Tam, A. T., Yu, Z., Kelso, R. M. & Binder, B. J. 2015 Predicting channel bed topography in hydraulic falls. Phys. Fluids 27, 112–106.
  * Thompson et al. [2014] Thompson, A. B., Juel, A. & Hazel, A. L. 2014 Multiple finger propagation modes in hele-shaw channels of variable depth. J. Fluid Mech. 746, 123–164.
  * Ţugulan et al. [2022] Ţugulan, C., Trichtchenko, O. & Părău, E. I. 2022 Three-dimensional waves under ice computed with novel preconditioning methods. J. Comp. Phys. 459, 111129.
  * Wade et al. [2014] Wade, S. L., Binder, B. J., Mattner, T. W. & Denier, J. P. 2014 On the free surface flow of very steep forced solitary waves. J. Fluid Mech. 739, 1–21.
  * Wade et al. [2017] Wade, S. L., Binder, B. J., Mattner, T. W. & Denier, J. P. 2017 Steep waves in free-surface flow past narrow topography. Phys. Fluids 29, 062107.
  * Whitham [1974] Whitham, G. B 1974 Linear and nonlinear waves. John Wiley.
  * Wu [1987] Wu, T. Y-T. 1987 Generation of upstream advancing solitons by moving disturbances. J. Fluid. Mech. 184, 75–99.
  * Zakharov [1968] Zakharov, V. E. 1968 Stability of periodic waves of finite amplitude on the surface of a deep fluid. J. Appl. Mech. Tech. Phys. 9 (2), 190–194.