# Trapped Free Surface Waves for a Lamb–Oseen Vortex Flow

E. Zuccoli1  E. J. Brambley12<EMAIL_ADDRESS>D. Barkley1
1Mathematics Institute, University of Warwick, UK 2WMG, University of Warwick,
UK

###### Abstract

Trapped surface waves have been observed in a swimming pool trapped by, and
rotating around, the cores of vortices. To investigate this effect, we have
numerically studied the free-surface response of a Lamb–Oseen vortex to small
perturbations. The fluid has finite depth but is laterally unbounded. The
numerical method used is spectrally accurate, and uses a novel non-reflecting
buffer region to simulate a laterally unbounded fluid. While a variety of
linear waves can arise in this flow, we focus here on surface gravity waves.
We investigate the linear stability of the vortex as a function of the
perturbation azimuthal mode number and the vortex rotation rate. We find that
at low rotation rates, linear modes decay by radiating energy to the far
field, while at higher rotation rates modes become nearly neutrally stable and
trapped in the vicinity of the vortex. While trapped modes have previously
been seen in shallow water surface waves due to small perturbations of a
bathtub vortex, the situation considered here is qualitatively different owing
to the lack of an inward flow and the dispersive nature of non-shallow-water
waves. We also find that for slow vortex rotation rates, trapped waves
propagate in the opposite direction to the vortex rotation, whereas, above a
threshold rotation rate, waves co-rotate with the flow.

###### keywords:

Surface gravity waves; Waves in rotating fluids.

##  1 Introduction

Rotating flows are rather ubiquitous in nature. Examples range from domestic
flows such as the bathtub vortex, to geophysical flows and Jupiter’s Great Red
Spot, to astrophysical flows such as accretion disks. Rotating flows support
the propagation of many different type of waves. These are often viewed as
small perturbations to an equilibrium state, also known as a base flow, and
the behaviour of such waves depends both on the particular velocity field of
the base flow and on the geometry of the problem under consideration. One
important class of such waves are interfacial surface waves. These are waves
forming at the interface between two fluids with different physical
properties, and they often are strongly localized close to the interface.

Our interest here is in free-surface waves on a vortex flow. In fact, our
initial motivation was out of curiosity concerning vortices and waves
interacting in a swimming pool. Changes in the surface height can be easily
visualized by light and dark patterns on the swimming pool floor (for which
the authors recommend a warm sunny climate and an outdoor swimming pool). An
“experiment” (Skipp, 2020) created pairs of vortices by drawing a dinner plate
through the water. The initial waves generated by this disturbance disperse
rapidly, leaving remarkably long-lived vortices with surface waves that appear
trapped in the vortex, but which propagate around the vortex in the opposite
direction to the vortex flow. A photograph of the phenomena, together with a
schematic representation, is shown in figure 1;

Figure 1: The physical scenario motivating this work. Surface waves are
trapped within a rotating vortex flow, but propagate in the direction opposite
the direction of rotation. Left: Photograph of two vortices in a swimming pool
(a still taken from the video of Skipp, 2020). Right: schematic of the flow
and surface waves on the left-hand vortex. Black arrows show the base velocity
field. Blue lines show the wave crests. Red arrows show the direction of
propagation of the waves.

interested readers are encouraged to follow the reference for a video.

The study of surface waves scattered by a vortex flow in a horizontally
unbounded domain in a two-dimensional shallow-water setting with a bathtub
vortex as a background flow was considered by Patrick et al. (2018) and
Patrick (2019), in order to provide an analogy between fluid phenomena and
black holes. Their study revealed the presence of modes that they called
“quasibound states”: normal modes having a very slow decaying rate in time and
a spatial structure which remains trapped within the vortex core region,
similar to standing waves in a potential well. Such states are perhaps
understandable for the bathtub vortex which includes an inflow towards the
vortex. In our work we investigate whether similar states exist also for
purely rotating flows without a shallow-water approximation so that the axial
variation of the perturbations cannot be neglected. We do not consider the
related wave scattering problem.

We will used the term trapped mode to refer to a wave mode localized in space
and with a very slow temporal decay rate. Any slow decay of trapped waves may
be due to dissipation (such as through viscosity), or by the wave not being
perfectly trapped and “leaking” outgoing waves to infinity. Trapped waves
could be excited either by an incoming propagating wave from infinity, or by
any initial disturbance (such as a dinner plate in the video depicted in
figure 1); here, we investigate the existence of trapped waves about a vortex,
being agnostic as to the mechanism creating them in the first place.

We will study free surface waves on rotating fluid whose velocity
$\boldsymbol{U_{0}}$ is that of a Lamb–Oseen vortex

$\boldsymbol{U_{0}}=\boldsymbol{\hat{\theta}}\frac{\Gamma_{0}}{2\upi
r}\Big{(}1-\exp\big{\\{}{-r^{2}\\!/a^{2}}\big{\\}}\Big{)},$ (1)

where $r$ is the horizontal distance from the vortex centre,
$\boldsymbol{\hat{\theta}}$ is a unit vector in the direction of rotation, and
$\Gamma_{0}$ and $a$ are constants (see, e.g. Drazin & Riley, 2006, for
details of the derivation of this flow). The Lamb–Oseen vortex is in fact an
exact solution to the viscous Navier–Stokes equations if we take $a^{2}=4\nu
t$, where $\nu$ is the kinematic viscosity, although for short timescales
compared with the viscous diffusion timescale, $a$ may be approximated as
constant. Unlike the bathtub vortex, this base flow is not the result of
vortex stretching by a downwards outflow. In that regard, it is more similar
to a geophysical vortex, although here we do not consider a stratified fluid,
nor do we assume the depth to be shallow. Unlike other vortices such as the
Rankine vortex (made up of an inner solid body rotation and an outer potential
vortex), which have been found to be unstable (e.g. Ford, 1994; Mougel et al.,
2014), the Lamb–Oseen vortex is usually considered to be a robust and stable
coherent structure. Riedinger et al. (2010) showed that it can become unstable
in a stratified environment, but with no stratification, Fabre et al. (2006)
have shown that the Lamb–Oseen vortex is stable to axial-invariant
perturbations. Their work highlighted the presence of critical layer waves,
arising within a critical layer. However, because the critical layer provides
a damping effect, they find the associated critical layer waves to be damped;
we will not consider the critical layer further here.

Free surface rotating waves have previously been studied computationally by
Mougel et al. (2015, 2017) in a confined geometry for background flows with
particular forms, such as solid-body rotation or a potential vortex flow.
However, since a potential vortex has an infinite velocity at its centre, and
a solid-body rotation has an infinite velocity at infinity, neither solid-body
rotation nor a potential vortex are appropriate for the case we consider here.
Their studies reveal the presence of four main types of waves: surface gravity
waves, inertial waves, Rossby waves, and centrifugal waves. The first three
waves are observed for solid-body rotation, while centrifugal waves and
gravity waves arise in a potential vortex flow. Depending on the problem, each
of these waves can be easily recognized and distinguished by its spatial
structure and the associated frequency of oscillation. Further examples of
problems where surface gravity waves, inertial waves and Rossby waves appear
can be found in Johnson (1997), Greenspan (1969) and McWilliams (2006),
respectively. An interesting feature of the study on the potential vortex flow
is the possibility of an instability of the base solution — known as a Polygon
Instability — due to an interaction between the centrifugal and gravity waves.
This interaction mechanism was initially studied by Tophøj et al. (2013) and
subsequently by Mougel et al. (2017). The stability of vortex flows in bounded
domains is however markedly different to the stability of vortex flows in
unbounded domains, such as considered here, as the instability in bounded
domains found by Tophøj et al. (2013) and Mougel et al. (2017) involves the
interaction of two types of wave modes, one of which is dependent on the outer
boundary. Instabilities of vortices on horizontally unbounded domains have
been found to occur within the shallow-water limit, as studied by Ford (1994).
However, it is worth noting that the study of shallow-water rotating flows is
rather different to the study of finite-depth rotating flows studied here, on
two counts: firstly, shallow-water rotating flows can only rotate very slowly
before the deformation in the free surface at the centre of the vortex touches
the bottom boundary and a dry inner region is formed; and secondly, shallow-
water waves are nondispersive and have a fixed wave speed, while finite-depth
water waves are dispersive such that any wave speed is available to the
system, and for example a wave can exist whose speed matches the flow speed at
a given location. Surface waves have been also studied by Hunt et al. (2015),
particularly considering the effects of an electric field on both linear and
nonlinear inviscid, irrotational waves. Moreover, two-dimensional surface
waves in electrohydrodynamics and magneohydrodynamics have been studied by
Hunt & Vanden-Broeck (2015); Hunt (2019); Hunt & Dutykh (2021), although in
the present work we do not consider any electric or magneic field, but instead
we consider the full three-dimensional linear response to a prescribed
swirling flow.

The present work considers a laterally unbounded domain (i.e. a domain that is
in principle of infinite radial extent), with a finite depth, and thus assumes
neither a shallow water nor a deep water limit. In order to numerically study
wave propagation problems in an unbounded domain, a non-reflecting boundary
condition (NRBC) needs to be imposed at infinity. This is relatively
straightforward in a two-dimensional shallow-water setting, since the
dispersion relation for surface waves is known, and hence a characteristics
analysis (John, 1978) allows the imposition of an exact NRBC. Characteristic
boundary conditions are, however, specific to classical hyperbolic wave
equations (Grote & Keller, 1995; Grote, 2000; Keller & Givoli, 1989; Bayliss &
Turkel, 1980; Engquist & Majda, 1977). Such exact boundary conditions cannot
be directly imposed on our three-dimensional problem, as our problem is not
necessarily hyperbolic, and it is not possible to find real characteristics.
Moreover, surface waves on non-shallow water are dispersive, with the
propagation speed of each mode depending on its frequency, which significantly
complicates a NRBC even for hyperbolic problems (e.g. Lindquist et al., 2012),
and can limit the range of wavenumbers being correctly treated (e.g. Wellens &
Borsboom, 2020). One possible solution to this problem is provided by
absorbing layer methods. The idea behind this class of techniques is to
surround the physical region with computational layers through which waves are
progressively damped such that they reach the edge of the computational domain
with sufficiently small amplitude that any reflections back into the physical
region will be negligible. The absorbing layers need to vary gradually,
however, as otherwise waves can spuriously reflect from the edges of the
absorbing layers. Absorbing layer methods have been applied initially to waves
problems in electromagnetics (e.g. Berenger, 1994) and then extended to other
physical fields like elasticity Hastings et al., 1996; Chew & Liu, 1996 and
fluid dynamics (Söderstrom et al., 2010). A shortcoming of absorbing layer
formulations lies in the number of new unknowns introduced in the problem;
this may be reduced using the formulation of Sim (2010) in the case of second
order hyperbolic equations. Here, a related but different technique is
developed for the absorbing layer formulation in the far field, without
introducing further unknowns in the mathematical problem.

The paper is organised as follows: in section 2, a general description is
given of the mathematical model used to study perturbations to a Lamb–Oseen
swirling flow, including a mathematical description of NRBCs. The resulting
eigenvalue problem is then solved using a numerical procedure described in
section 3. The results of this numerical solution are described in section 4.
Finally, resulting conclusions and possible future directions for continued
research are highlighted in section 5.

##  2 Mathematical Model

Assuming that viscosity is negligible over the timescales of interest here,
the governing equations are the incompressible Euler equations,

$\displaystyle\frac{\partial\boldsymbol{U}}{\partial
t}+\boldsymbol{U\cdot\nabla}\boldsymbol{U}+\frac{1}{\rho}\bnabla
P+g\boldsymbol{\hat{z}}$ $\displaystyle=0,$
$\displaystyle\boldsymbol{\nabla\cdot U}$ $\displaystyle=0,$ (2)

where $\boldsymbol{U}$ is the fluid velocity, $P$ is the fluid pressure,
$\boldsymbol{\hat{z}}$ is a unit vector in the vertical direction, and the
constants $\rho$ and $g$ are the fluid density and the acceleration due to
gravity respectively. The fluid is contained between a bottom boundary at
$z=0$ and an upper free surface at $z=H$. The fluid must satisfy no
penetration through the bottom boundary, giving $\boldsymbol{U\cdot\hat{z}}=0$
at $z=0$. Two boundary conditions, a kinematic and a dynamic boundary
condition, must be satisfied along the free surface itself. Here, we assume
the fluid above the free surface to be dynamically passive, and in particular,
to have a constant pressure $\bar{P}$. Together, these give the boundary
conditions

$\displaystyle\frac{\partial H}{\partial t}+\boldsymbol{U\cdot\nabla}H$
$\displaystyle=\boldsymbol{U\cdot\hat{z}}$ and $\displaystyle P$
$\displaystyle=\bar{P}$ on $\displaystyle z$ $\displaystyle=H.$ (3)

We split the overall velocity and pressure into a steady purely swirling base
flow and a small (magnitude $\varepsilon$) time-dependent perturbation, where

$\displaystyle\boldsymbol{U}$
$\displaystyle=U_{0}(r)\boldsymbol{\hat{\theta}}+\varepsilon(u_{r}\boldsymbol{\hat{r}}+u_{\theta}\boldsymbol{\hat{\theta}}+u_{z}\boldsymbol{\hat{z}}),$
$\displaystyle P$ $\displaystyle=P_{0}(r,z)+\varepsilon p,$ $\displaystyle H$
$\displaystyle=h_{0}(r)+\varepsilon h.$ (4)

###  2.1 The steady base flow solution

A purely swirling steady base flow has a velocity given by
$\boldsymbol{U_{0}}=U_{0}(r)\boldsymbol{\hat{\theta}}$. The governing
equations and boundary conditions for this steady flow are satisfied provided
we take

$\displaystyle P_{0}(r,z)$ $\displaystyle=\bar{P}+\rho
g\left(h_{0}(r)-z\right),$ (5a) $\displaystyle h_{0}(r)$
$\displaystyle=h_{\infty}-\frac{1}{g}\int_{r}^{\infty}\frac{U^{2}_{0}(r^{\prime})}{r^{\prime}}\,\mathrm{d}r^{\prime},$
(5b)

where $h_{\infty}$ is the depth of the fluid at $r=\infty$. This holds for any
velocity profile $U_{0}(r)$. For the specific case of the Lamb–Oseen vortex
considered here, we have

$U_{0}(r)=\frac{\Gamma_{0}}{2\upi
r}\Big{(}1-\exp\left(-r^{2}\\!/a^{2}\right)\Big{)},$ (6)

where $a$ sets the radial size of the core and $\Gamma_{0}$ sets the
circulation of the vortex. Note that, for the Lamb–Oseen vortex, for small
$r$, we have $U_{0}(r)\approx r\Gamma_{0}/2\upi a^{2}$, so that (6) is a solid
body rotation near the centre of the vortex, while for large $r$ we have
$U_{0}(r)\approx\Gamma_{0}/2\upi r$, so that (6) is a potential swirl far from
the centre of the vortex. A typical steady base flow free surface for the
Lamb–Oseen vortex is shown in figure 2.

Figure 2: A typical section through a Lamb–Oseen vortex, showing the steady
base flow free surface.

###  2.2 Perturbation Dynamics

Waves arise when a small perturbation is introduced to the steady base
solution. By linearizing the governing equations (2) about the base solution
$(U_{0},P_{0},h_{0})$ given above, the governing equations for the
perturbation are

$\displaystyle
D_{t}u_{r}-2\Omega_{0}(r)u_{\theta}+\frac{1}{\rho}\frac{\partial p}{\partial
r}$ $\displaystyle=0,$ (7a) $\displaystyle
D_{t}u_{\theta}+\frac{1}{r}\big{(}rU_{0}(r)\big{)}^{\prime}u_{r}+\frac{1}{\rho
r}\frac{\partial p}{\partial\theta}$ $\displaystyle=0,$ (7b) $\displaystyle
D_{t}u_{z}+\frac{1}{\rho}\frac{\partial p}{\partial z}$ $\displaystyle=0,$
(7c) $\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}\left(ru_{r}\right)+\frac{1}{r}\frac{\partial
u_{\theta}}{\partial\theta}+\frac{\partial u_{z}}{\partial z}$
$\displaystyle=0,$ (7d)

where $D_{t}=\partial_{t}+\Omega_{0}(r)\partial_{\theta}$ is the convective
derivative and $\Omega_{0}(r)=U_{0}(r)/r$ is the steady base flow angular
velocity. (Here and throughout, primes denotes derivatives with respect to the
argument for a function of only one variable). Linearizing the boundary
conditions (3) about the steady base flow leads to

$\displaystyle\begin{aligned} u_{z}&=\frac{1}{\rho
g}D_{t}p+h^{\prime}_{0}u_{r}&&\quad\text{and}\quad&h&=\frac{p}{\rho
g}&&\quad\text{on}&z&=h_{0}(r),\end{aligned}$ (8a) $\displaystyle
u_{z}=0\quad\text{on}\quad z=0.$ (8b)

Since our focus is on trapped waves and their formation, in addition to these
boundary conditions there is another implicit condition, which is that there
are no waves entering the domain from $r=\infty$, and thus only outgoing waves
are allowed at $r=\infty$. This rather subtle condition will become more
concrete when we truncate the domain to finite $r$ in order to numerically
solve the above equations.

Before we progress further, we now make two changes to the governing equations
for the perturbation. Firstly, since the governing equations for the
perturbation are linear, we may assume a modal solution of the form

$\displaystyle\boldsymbol{u}$
$\displaystyle=\left[u(r,z)\boldsymbol{\hat{r}}+v(r,z)\boldsymbol{\hat{\theta}}+w(r,z)\boldsymbol{\hat{z}}\right]\exp(-\mathrm{i}\omega
t+\mathrm{i}m\theta),$ (9a) $\displaystyle p$
$\displaystyle=\phi(r,z)\exp(-\mathrm{i}\omega t+\mathrm{i}m\theta),$ (9b)

where $m$ is restricted to integer values, since the solution must be $2\upi$
periodic, while $\omega$ in general will be complex, with
$\Real(\omega)/2\upi$ being the oscillation frequency and $-\Imag(\omega)$
being the decay rate. The eigenvalue problem that will eventually result will
have $\omega$ as the eigenvalue to be found.

Secondly, we re-write the governing equations and boundary conditions in a
non-dimensional form. To do so, we choose the reference lengthscale to be $a$,
the scale of the vortex core (equation (1) and figure 2). The reference time
scale is that given by this lengthscale and gravity: $\sqrt{a/g}$. The
velocity scale is thus $\sqrt{ag}$. All dimensional variable may then be
expressed in terms of a nondimensional variable (denoted by a tilde), as

$\begin{gathered}(r,\theta,z)=(a\tilde{r},\tilde{\theta},a\tilde{z}),\qquad
t=\sqrt{\frac{a}{g}}\tilde{t},\qquad U_{0}(r)=\frac{\Gamma_{0}}{2\upi
a}\widetilde{U}_{0}(\tilde{r}),\qquad\Omega_{0}=\frac{\Gamma_{0}}{2\upi
a^{2}}\tilde{\Omega_{0}},\\\ h=a\tilde{h},\qquad p=\rho
ag\tilde{p},\qquad\omega=\sqrt{\frac{g}{a}}\tilde{\omega},\qquad\boldsymbol{u}=\sqrt{ag}\boldsymbol{\tilde{u}}.\end{gathered}$
(10)

There are two physical parameters that are not scaled to unity by this
nondimensionalization: these may be thought of as the strength of the vortex
and the depth of the fluid at infinity, given respectively as

$\displaystyle F$ $\displaystyle=\frac{\Gamma_{0}}{2\upi\sqrt{ga^{3}}},$
$\displaystyle\tilde{h}_{\infty}$ $\displaystyle=\frac{h_{\infty}}{a}.$ (11)

$F$ is the Froude number and sets the nondimensionalized velocity of the
vortex. Hence, $F\to 0$ corresponds to a slow vortex with negligible steady
surface height variation, while $F\to\infty$ corresponds to a fast vortex with
significant steady surface height variation, as can be seen from ($None$c)
below. The dimensionless depth $\tilde{h}_{\infty}$ is exactly that, so that
the limit $\tilde{h}_{\infty}\to 0$ corresponds to the shallow-water limit and
$\tilde{h}_{\infty}\to\infty$ corresponds to the deep-water limit. Care is
needed, however, in considering the shallow-water limit: in what follows, we
will assume that the steady fluid height never reaches zero, so that the
bottom stays wetted, and consequently the shallow water limit
$\tilde{h}_{\infty}\to 0$ must be taken together with the slow wide vortex
limit $F\to 0$ such that $\tilde{h}_{\infty}/F^{2}$ is bounded away from zero,
as can also be seen from ($None$) below. The azimuthal wavenumber $m$ is also
a dimensionless parameter, and represents the rotational symmetry of the
solution being investigated.

Dropping the tildes, the complete nondimensional eigenvalue problem is

$\displaystyle\big{(}{-\mathrm{i}}\omega+\mathrm{i}mF\Omega_{0}(r)\big{)}u-2F\Omega_{0}(r)v+\frac{\partial\phi}{\partial
r}$ $\displaystyle=0,$ (12a)
$\displaystyle\big{(}{-\mathrm{i}}\omega+\mathrm{i}mF\Omega_{0}(r)\big{)}v+\frac{F}{r}\big{(}rU_{0}(r)\big{)}^{\prime}u+\frac{\mathrm{i}m}{r}\phi$
$\displaystyle=0,$ (12b)
$\displaystyle\big{(}{-\mathrm{i}}\omega+\mathrm{i}mF\Omega_{0}(r)\big{)}w+\frac{\partial\phi}{\partial
z}$ $\displaystyle=0,$ (12c) $\displaystyle\frac{1}{r}\frac{\partial}{\partial
r}\left(ru\right)+\frac{\mathrm{i}m}{r}v+\frac{\partial w}{\partial z}$
$\displaystyle=0,$ (12d)

together with the boundary conditions of no incoming modes at $r=\infty$, and

$\displaystyle w$ $\displaystyle=0,$ $\displaystyle\text{on}\quad z$
$\displaystyle=0,$ (13a) $\displaystyle w$
$\displaystyle=\big{(}{-\mathrm{i}}\omega+\mathrm{i}mF\Omega_{0}(r)\big{)}\phi+F^{2}r\Omega^{2}_{0}(r)u,$
$\displaystyle\text{on}\quad z$ $\displaystyle=h_{0}(r).$ (13b)

For the Lamb–Oseen vortex, in dimensionless terms

$U_{0}(r)=\frac{1-\exp(-r^{2})}{r},\quad\Omega_{0}(r)=\frac{U_{0}(r)}{r},\quad
h_{0}(r)=h_{\infty}-F^{2}\int_{r}^{\infty}\frac{U^{2}_{0}(r^{\prime})}{r^{\prime}}\,\mathrm{d}r^{\prime}.$
$None$

One consequence of the harmonic assumption (9) is the creation of a so-called
“critical layer” at a radial location $r=r_{c}$, given by
$D_{t}=-\mathrm{i}\omega+\mathrm{i}mF\Omega_{0}(r_{c})=0$. It is not
immediately obvious from (12) that anything particularly special occurs at the
critical layer, as in no equation does it cause the highest derivative to
vanish, and therefore our numerics described in §3 has no difficulty in this
case. However, in fact the critical layer is related to behaviour that is not
of the harmonic form assumed in (9); investigating the effect of the critical
layer requires a different mathematical and numerical technique (such as
Frobenious series; see, by way of example, King et al., 2022), which is beyond
the scope of this paper. However, we comment in passing that Fabre et al.
(2006) found that disturbances related to the critical layer are necessarily
damped, suggesting that the undamped trapped wave modes we are interested in
here are not related to the critical layer.

##  3 Numerical methods

###  3.1 Absorbing layer for 3D incompressible Euler equations

In practice, we solve the eigenvalue problem in a finite computational domain
and hence introduce an artificial boundary at finite radius $R$. The far-field
boundary condition of no incoming waves becomes a non-reflecting boundary
condition (NRBC) at this boundary. However, exact NRBCs are generally derived
using the method of characteristics, but this is not available to us since
(12d) is not hyperbolic. A widely used alternative is based on damping layers
or even perfectly matched layers. The idea is to introduce a “damping layer” —
also called an “absorbing layer” — surrounding the physical outer boundary. In
the damping layer, waves are progressively damped until their amplitude is
sufficiently small that reflections from the outer boundary do not reenter the
physical domain. Outgoing waves are effectively absorbed in the layer. The
boundary condition at the outer boundary of the computational domain is then
unimportant, and is typically taken to be either a Dirichlet-type or Neumann-
type boundary condition.

To add damping to our governing equations, we add a “damped compressibility”
term into the continuity equation (12d), which then becomes

$\xi(r)\phi+\frac{1}{r}\frac{\partial}{\partial
r}\left(ru\right)+\frac{\mathrm{i}m}{r}v+\frac{\partial w}{\partial z}=0.$
(15)

where $\xi(r)$ is a sufficiently smooth function that is identically zero
outside the damping region. This is motivated by the introduction of damping
in the acoustic wave equation following Gao et al. (2017, pp. 81-82), and is
explained in further detail in appendix B. In what follows, we find good
results using the simple form for $\xi(r)$ given by

$\xi(r)=\begin{cases}0,&r<R_{c},\\\
\bar{\xi}\left(\dfrac{r-R_{c}}{R-R_{c}}\right)^{2},&R_{c}\leq r\leq
R.\end{cases}$ (16)

This gives a physical region, $r<R_{c}$, with no damping, and computational
region $r\in[R_{c},R]$ with a damping strength governed by the constant
$\bar{\xi}$. We then impose a Dirichlet boundary condition on the pressure at
an artificial boundary $r=R\gg 1$. A representation of the damping function as
well as the subdivision of the numerical domain is shown in figure 3.

Figure 3: Left: profile of the quadratic damping function defined in (16) with
parameters: $\bar{\xi}=5$, $R_{c}=5$ and $R=10$. Right: schematic of the
computational domain in the $r-z$ plane.

The governing equations (12a–c) together with the modified continuity equation
(15) and boundary conditions (13) form an eigenvalue problem to solve for the
allowable frequencies $\omega$ permitting a nonzero modal solution.

###  3.2 Numerical discretization

In order to solve (12a–c,13,15), we used a Galerkin spectral method and expand
with a combination of Legendre polynomials as basis functions in both the
radial and axial coordinate in order to satisfy the Dirichlet boundary
conditions. We also remap the domain from the domain
$D=[0,R]\times[0,h_{0}(r)]$ to the computational square domain
$S=[-1,1]\times[-1,1]$ to account for the shape of the computational domain
with a variable surface height $h_{0}(r)$. We then obtain the weak formulation
of the problem. Full details are given in appendix A. The discretized problem
ends up being of the form

$\displaystyle\mathsfbi{A}\boldsymbol{w}=\omega\mathsfbi{B}\boldsymbol{w},$
(17)

with $\boldsymbol{w}=(u_{ij},v_{ij},w_{ij},\phi_{ij})$ representing our array
containing the spectral coefficients of each unknown, $\mathsfbi{A}$ and
$\mathsfbi{B}$ being matrices of order $4N_{x}N_{y}\times 4N_{x}N_{y}$, where
$N_{x}$ and $N_{y}$ are the number of Legendre polynomials used in the radial
and vertical directions respectively. This discretized problem may then be
solved using any numerical eigenvalue solver; here, we use the eig solver in
Matlab.

Not all solutions to the discretized problem (17) correspond to solutions to
the continuous problem being approximated, however. To remove under-resolved
eigenmodes and spurious eigenvalues, the numerical solutions to (17) are
filtered, as described in the appendix A.2. The numerical scheme is validated
in the case of the simpler Newton’s Bucket problem by comparison against the
results of Mougel et al. (2015) in appendix C.)

###  3.3 Spurious reflected modes

The finite numerical domain and damping region introduces another source of
spurious modes besides those coming from the numerical discretization, namely
those modes which are well-resolved but which include a significant reflection
from either the damping boundary at $r=R_{c}$ or the truncation boundary at
$r=R$. Such spurious eigenmodes are affected by the values of $R_{c}$ and $R$,
and by the amount of damping $\bar{\xi}$, whereas good approximations to the
modes on the infinite domain should be insensitive to these values. Since
-$\Imag(\omega)$ is the decay rate of the mode, variations in damping
typically have a strong effect on $\Imag(\omega)$ for spurious reflected
modes. We may therefore remove these spurious reflected modes by running our
numerical code twice: the first time with a suitably chosen amount of damping
$\bar{\xi}$, and the second time with twice that amount of damping
$2\bar{\xi}$. Only those modes whose eigenvalues do not change significantly
with the change in the damping coefficient are retained (as measured using the
same metric (28) used for the numerical resolution).

###  3.4 Numerical convergence study

The first convergence study involves varying the amount of damping in the
damping layer and checking that the eigenvalues do not change, nor does the
shape of the corresponding eigenfunctions in the physical domain $r<R_{c}$.
Here we present the convergence results for $m=7$ and $F=0.3$, since for these
parameters the system supports radiating (outwardly propagating) modes, and
such modes provide the most stringent test of a non-reflecting boundary
condition. For varying magnitudes of damping $\bar{\xi}$, the radiating
eigenfunction is plotted in figure 4,

Figure 4: Eigenfunctions for $m=7$, $F=0.3$, and different values of
$\bar{\xi}$. A vertical red-dashed line indicates the radial point $R_{c}$
where the damping effect begins. Left: a small computational domain $R_{c}=5$,
$R=10$. Right: a large computational domain $R_{c}=15$, $R=30$.

and the corresponding eigenvalues are given in table 1.

$\bar{\xi}$ | 0.1 | 1 | 5 | 10 | 50
---|---|---|---|---|---
$\Real(\omega)$ | -1.3911 | -1.3894 | -1.3766 | -1.3769 | -1.3688
$\Imag(\omega)$ | $-8\times 10^{-4}$ | -0.0085 | -0.0133 | -0.0138 | \- 0.0123
Table 1: Eigenvalues as function of the amount of damping for $m=7$ and
$F=0.3$, $R_{c}=5$ and $R=10$.

Two domain sizes are shown in figure 4: $R_{c}=5$ and $R=10$ are the values
used for the results in section §4; and $R_{c}=15$ and $R=30$ give an extended
domain so that the effects of damping and resonance can be seen more clearly.
The damping clearly influences the eigenfunction shape, and for
$\bar{\xi}\in\\{0.1,1\\}$ a clear standing wave shape is seen in figure 4. For
$\bar{\xi}\in\\{5,10\\}$, the eigenfunctions are practically identical in the
physical domain $0<r<R_{c}$, and only differ in the damping layer $R_{c}\leq
r\leq R$. For $\bar{\xi}=50$, however, the damping is too strong, and little
oscillations can be seen for $r<R_{c}$, suggesting wave reflection by the edge
of the damping layer. This is also supported by the eigenvalues in table 1,
which show the sensitivity of $\Imag(\omega)$ to variations in damping
strength, as expected.

The second convergence study involves varying the width of the damping layer
whilst maintaining the same size of the computational domain and a fixed
damping strength $\bar{\xi}=5$. The eigenfunctions are displayed in figure 5,

Figure 5: Eigenfunctions for $m=7$, $F=0.3$, $\bar{\xi}=5$, and different
values of $R_{c}$, computed using a small domain $R=10$ (left) and a larger
domain $R=30$ (right).

and the corresponding eigenvalues in table 2.

$R_{c}$ | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---
$\Real(\omega)$ | -1.3766 | -1.3767 | -1.3802 | -1.3788 | -1.3918
$\Imag(\omega)$ | -0.0133 | -0.0131 | -0.0133 | \- 0.0143 | -0.0045
Table 2: Eigenvalues as function of the initial position of the absorbing
layer for $m=7$ and $F=0.3$, $\bar{\xi}=5$ and $R=10$.

For most results in figure 5 the eigenfunctions can be seen to be very similar
in the physical domain $0<r<R_{c}$ and to decay smoothly in $r$ in the damping
layer, while for a damping layer of width $R-R_{c}=1$ a standing wave pattern
can be seen for both domain sizes, implying significant wave reflection from
the domain boundary at $r=R$. Again, this is also seen for the variations in
the eigenvalue in table 2, with again $\Imag(\omega)$ being particularly
sensitive.

Finally, in figure 6

Figure 6: Eigenfunction for $m=7$, $F=0.3$, $\bar{\xi}=10$ computed using two
different domain sizes: $R_{c}=5$ and $R=10$ for the blue solid curve; and
$R_{c}=15$, $R=30$ for the red dashed curve.

we compare the radiating mode for $m=7$, $F=0.3$ computed using two different
size of the domain; $R_{c}=5$, $R=10$ on one hand and $R_{c}=15$, $R=30$ on
the other. The two modes match very well over the physical interval
$r\in[0,5]$. This good match demonstrates that the damping layer successfully
allows our numerical simulation on a small interval to reproduce results that
would have been obtained using a larger interval, thus allowing our numerics
to emulate an infinite unbounded domain using a finite computational domain.

###  3.5 Choice of numerical parameters

Based on numerical convergence studies, for the results that follow we take
$R=10$ and $R_{c}=5$, with $\bar{\xi}=5$ (to which results are compared to
$\bar{\xi}=10$). This choice is motivated by the need for sufficiently high
resolution to resolve all modes of interest in the range of Froude numbers and
azimuthal wavenumbers considered. The eigenvalue tolerance and eigenfunction
resolvedness tolerance are taken to be $\mathrm{tol}=10^{-2}$ and $10^{-1}$
respectively, whereas other numerical parameters are taken to be $N_{x}=50$,
$N_{y}=20$, $b_{x}=12$ and $b_{y}=4$, as explained in more details in appendix
A. For the majority of results presented here, we use a fluid depth of
$h_{\infty}=5$. This is the depth shown in figure 2. This value is large
enough to allow for a wide range of Froude numbers ($0\leq
F<\sqrt{h_{\infty}/\log(2)}\simeq 2.68$) without forming a dry region near
$r=0$, while small enough to exhibit finite-depth effects.

##  4 Results

We compute leading (i.e. least damped) eigenmodes for a range of Froude
numbers and for a range of azimuthal wave numbers as large as $m=20$, and find
that above $m=6$ surface gravity waves dominate. We shall first present the
surface waves for a representative case, wavenumber $m=7$, and discuss in
detail how the modes and eigenvalues depend on Froude number for this case.
Following this, we consider the dependence of modes and eigenvalues on the
azimuthal wave number.

###  4.1 Representative case $m=7$

We begin with a detailed description of the case azimuthal wavenumber $m=7$.
The reason for taking such value of $m$ is dictated both by the explanatory
picture in Patrick et al. (2018, pp.6), who show a picture of the spiral
structure of a normal mode solution for sufficiently high azimuthal
wavenumber, as well as by the reasonable amount of computational time needed
to get the leading surface waves eigenmodes. Indeed, we find numerically that
the spatial structure of surface gravity waves becomes thinner and more
localized close to the free surface as $m$ increases. Following particularly
the latest argument, the case of $m=7$ has been taken as a reference study
case and in the following we are going to show most of the interesting
features applied to this case.

One of the key results of our study is the change in character of modes as
function the rotation rate of the vortex. Specifically, we characterize two
extremes of modes: radiating modes and trapped modes. The spatial structure of
the radiating modes has a characteristic spiral shape, as shown in figure 7a,

Figure 7: Plots of the free surface height
$h(r,\theta,t=0)=\Real[\phi(r,h_{0}(r))\exp\\{\mathrm{i}m\theta\\}]$ for
$m=7$. (a) $F=0.3$. (b) $F=0.32$. (c) $F=0.33$. (d) $F=0.34$. (e) $F=0.36$.
All modes displayed rotate clockwise, i.e. against the vortex flow.

and a long tail in the radial direction, as shown in figure 8a.

Figure 8: Plots of the pressure distribution $|\phi_{m}(r,z)|$ for $m=7$. (a)
$F=0.3$. (b) $F=0.32$. (c) $F=0.33$. (d) $F=0.34$. (e) $F=0.36$.

Such modes behave as travelling waves along the radial direction, and
radiation to infinity is responsible for the dissipation of the initial energy
put into the system to excite the mode. In contrast, trapped modes are
confined to a bounded region in the radial direction, and without the
spiraling structure seen for the radiating modes. Those modes behave as
standing waves along the radial direction and dissipate extremely slowly in
comparison. Examples can be seen in figures 7e and 8e. The distinction between
radiating and trapped modes can further be seen in the eigenvalues. Indeed,
radiating modes have eigenvalues with a significant negative imaginary part,
indicating exponential decay in time, whereas trapped modes have a negligible
imaginary part and so do not appreciably decay in time.

Figures 7 and 8 show representative examples of the continuous, but rapid,
transition from radiating modes to trapped modes as the Froude number (the
dimensionless rotation rate of the vortex), increases. This transition between
radiating and trapped modes is also seen in figure 9

Figure 9: Modulus of the free surface height along the radius for $m=7$ in the
transition regime from radiating to trapped modes.

where the modulus of the free surface height is plotted as a function of $r$
for different Froude numbers.

To look more closely at the trend in the eigensolutions, let us first denote
by $n$ an integer representing the number of peaks in radial direction of the
modulus of the pressure eigenfunctions, as displayed in figure 10.

Figure 10: Modulus of the pressure $|\phi_{m}(r,z)|$ for $m=7$ and $F=0.5$. We
show the four least stable eigenmodes for such parameters, having one, two,
three and four peaks respectively.

In this way each eigensolution will be indexed by both $m$ and $n$, with
corresponding eigenvalues $\omega=\omega_{mn}$. In figure 11a,

(a)

(b)

Figure 11: Imaginary (top) and real (bottom) parts of the eigenvalue
$\omega_{mn}$ as function of the Froude number $F$ for $m=7$ and varying $n$.

the two types of modes can be seen: radiating modes with a significant
negative $\Imag(\omega)$, and trapped modes having an almost null
$\Imag(\omega)$. There is no sharp transition between radiating and trapped
modes. Accordingly, here we set an arbitrary threshold to separate the two
sets of modes by considering a mode to be trapped when its eigenvalue has an
imaginary part smaller than $10^{-5}$ in modulus (although our results are
relatively insensitive to this threshold; see figure 14 discussed below).
Hence, while the trapped modes considered here are almost neutrally stable,
they have eigenvalues with small negative imaginary part and so radiate very
slightly. Interestingly, none of the modes computed here are observed to
become linearly unstable, and we hypothesise that the base solution is at most
marginally stable, but not unstable, to linear perturbations of this type at
large Froude number. This stability can be viewed as a manifestation of the
robustness of the Lamb–Oseen vortex (Fabre et al., 2006), even in the presence
of a free surface.

In our calculations, the eigenmodes for sufficiently low Froude numbers become
dependant on the amount of damping introduced in the system, and are therefore
discarded. A similar result has been obtained in Oliveira et al. (2014)
studying a model wave equation describing the vortex-waves interaction in the
shallow-water limit. Mathematically, this can be justified by considering the
extreme case of $F=0$; in that case, a regularity condition of the modes at
$r=0$ together with a radiating boundary condition at infinity cannot be
satisfied simultaneously, thus leading to the non-existence of solutions in
terms normal modes. We leave to future studies the possibility of tracking the
trend of the spectrum curve as the Froude number goes to zero. It should be
noted, that modes having a large negative decay rate in time do not play a
relevant role in the dynamics of the system.

Another interesting phenomenon concerns the propagation direction of the
surface waves with respect to the rotation of the base vortex flow as its
rotation rate is varied. This is shown in figure 11b, which shows how the real
part of the eigenvalues varies with the Froude number, again for $m=7$. For
$\Real(\omega)<0$, waves rotate opposite to the base flow (counter-rotating
waves), while for $\Real(\omega)>0$, waves rotate in the same direction as the
base flow (co-rotating waves). It is clear that for each $n$, there is a value
of Froude number separating counter-rotating from co-rotating surface waves.

###  4.2 Extension to other azimuthal wavenumbers.

The results described above for $m=7$ are found to be typical at larger
azimuthal wavenumbers.

Figure 12: Imaginary part of the eigenvalues as functions of $F$ for
$m=2,3,4,5,6,7,8,9,10,15,20$, with the arrow indicating the direction of
increasing $m$.

Figure 12 shows the trend of the imaginary part of the eigenvalues as a
function of the Froude number. The results are qualitatively similar to the
$m=7$ case. Eigenvalue branches shift to lower $F$ with increasing $m$, thus
shifting to lower $F$ the value of $F$ separating radiating and trapped modes.
This suggests that for very high azimuthal wavenumber perturbations we expect
to get radiating eigenmodes at lower and lower Froude numbers, as shown in
figure 13

Figure 13: Shape of the free surface height
$h(r,\theta,t=0)=\Real[\phi(r,h_{0}(r))\exp\\{\mathrm{i}m\theta\\}]$ in two
high azimuthal wavenumber cases. (a): $m=15,F=0.19$. (b): $m=20,F=0.16$.

for the example cases of $m=15$ and $m=20$.

We have computed the range of Froude numbers in which modes radiate and
transition towards a neutrally stable state for different values of $m$. This
is summarized in figure 14

Figure 14: Regimes of radiating (large $|\Imag(\omega)|$) and trapped (small
$|\Imag(\omega)|$) modes for different azimuthal wavenumbers $m$ and Froude
numbers $F$. Solid lines are contours of $|\Imag(\omega)|$, with the arrow
indicating the direction of more perfectly trapped behaviour. The dash-dot
line separates counter-rotating modes ($T^{-}$) with $\Real(\omega)<0$ from
co-rotating modes ($T^{+}$) with $\Real(\omega)>0$.

showing the classes of solutions previously described over a wide range of
azimuthal orders $m$ and Froude numbers $F$. The contour lines plotted
separate the regions of parameter space where waves are radiating and where
waves are trapped, with the arrow indicating the direction of better trapping.
While the individual contours range from $\Imag(\omega)=10^{-2}$ to
$\Imag(\omega)=10^{-5}$, their close spacing shows that the exact threshold
value of $\Imag(\omega)$ separating trapped and radiating behaviour is not
that important, with all contours giving a similar boundary. Also plotted in
figure 14 is a dashed line showing $\Real(\omega)=0$, which is the boundary
between counter-rotating and co-rotating trapped waves. Two notable features
of the plot as $m$ is decreased are the shift of contours to larger $F$ and
the widening of the separation between contours. Both of these features are
consistent with figure 12, where one sees not only a shift in the eigenvalue
curves with $m$, but also a steepening of the transition between radiating and
trapped modes with increasing $m$.

While this section considers a wide range of values of $m\geq 3$, we find
$m=0,1$ and $2$ to be dominated by intertial waves rather than surface waves,
which are discussed further in section 4.4.

###  4.3 Effect of the free-surface height on the eigenmodes

In this subsection we investigated the effects associated to a variation in
the fluid depth, exploring the regime between shallow and deep water. We again
focus on modes with azimuthal modenumber $m=7$ and show how the structure of
eigenmodes and corresponding eigenvalues vary with fluid depth, and hence how
the transition between radiating and trapped modes changes with depth.

The variation in the structure of the surface gravity eigenmodes with fluid
depth $h_{\infty}$ is displayed in figure 15

Figure 15: Effect of the fluid height variation on a single peak ($n=1$)
trapped pressure mode computed for $m=7$, and $F=0.5$. (a): $\omega=-0.3058$,
$h_{\infty}=0.25$. (b): $\omega=-0.5316$, $h_{\infty}=0.5$ (c):
$\omega=-0.5366$, $h_{\infty}=1$ (d): $\omega=-0.5370$, $h_{\infty}=5$. The
imaginary part of the eigenvalues reported here is at least of order
$10^{-9}$, hence can be considered as purely real.

for a typical trapped mode computed at Froude number $F=0.5$. The same type of
analysis is carried out for a radiating mode at $F=0.3$, whose spatial
structures are shown in figure 16.

Figure 16: Effect of the fluid height variation on a radiating pressure mode
computed for $m=7$, and $F=0.3$. (a): $\omega=-1.3371-0.0530\mathrm{i}$,
$h_{\infty}=0.4$. (b): $\omega=-1.3633-0.0320\mathrm{i}$, $h_{\infty}=0.5$
(c): $\omega=-1.3772-0.0140\mathrm{i}$, $h_{\infty}=1$ (d):
$\omega=-1.3770-0.0138\mathrm{i}$, $h_{\infty}=5$.

Finally, in figure 17,

Figure 17: Free surface height for $m=7$, $F=0.4$ and different fluid depths.
(a): $h_{\infty}=0.25$. (b): $h_{\infty}=0.5$. (c): $h_{\infty}=1$.

the transition between radiating and trapped mode is shown as the fluid depth
is varied between $h_{\infty}=0.25,0.5,1$. Figure 18

Figure 18: Effect of the fluid height variation on the spectrum curves for the
least damped set of eigenmodes obtained for $m=7$.

shows the variation in eigenvalue as the Froude number is varied, for
different fluid depths. The trend in $\Imag(\omega)$ shown suggests that more
shallow systems require a faster rotating vortex in order to get the same wave
trapping as a deeper water configuration. The variation with depth of
$\Real(\omega)$ is less significant. For $h_{\infty}=1$, the eigenvalues are
essentially identical to those computed for even higher depths, for example,
in this case $h_{\infty}=5$. Hence, in this regard, for a trapped mode a fluid
with $h_{\infty}\in[0.5,1]$ can already be considered deep water.

###  4.4 Inertial modes

While the focus of this study is surface gravity waves, our numerics finds
modal solutions to the governing equations (12) without any assumption about
the modes being surface gravity waves. By way of contrast, therefore, in this
section we briefly discuss another type of mode, namely inertial modes. These
modes persist even at low rotation rates, and are characterised by being
neutrally stable as well as being concentrated within the vortex core and not
being localized close to the free surface. Our results are therefore similar
to the numerical results presented by Mougel et al. (2015) for a solid-body
rotation, since the Lamb–Oseen vortex is very close to solid-body rotation
close to the centre.

Figure 19: Inertial modes for $m=2$, $h_{\infty}=1$ and $F=0.2$ (left);
$F=0.8$ (right). Eigenvalues curve as function of the Froude number for the
same parameters (bottom figure).

Figure 19 shows two inertial modes computed for $m=2$, $h_{\infty}=1$, and
$F=0.2,0.8$. The corresponding eigenvalues are purely real, and are also
plotted as function of the Froude number in figure 19.

##  5 Conclusions and future developments

We have considered the linear response to small perturbations of a free
surface Lamb–Oseen vortex flow in a laterally unbounded domain. We restricted
our focus to surface waves, as it was these waves that were seen in the
motivating experiments; we note that the numerical procedure described here
also calculates other types of modes which we have not investigated in depth
here, and which may in some circumstances be the least damped modes,
especially at low Froude numbers.

Our study classified surface waves into two distinct classes: radiating modes
and trapped modes. While these have previously been seen for shallow water
vortices with inflow (such as a plughole vortex), here we have shown that they
arise in a fully 3D non-shallow-water problem without the need for an inflow
to help trap the modes. Radiating modes are dissipative, and have a spatial
structure which extends in the horizontal plane, therefore behaving as
radially-travelling waves. By contrast, trapped modes, by virtue of their near
zero growth rate, persist for long times with little dissipation and remain
localized within the core region of the vortex and hence could be expected to
be observed on the surface of a rotating vortex after all damped modes have
dissipated. The spatio-temporal behaviour of these modes resembles that of a
radially-standing wave instead of a radially-travelling wave. Our numerical
predictions show the appearance of trapped modes provided the Froude number
$F$ (the dimensionless rotation rate) is above a threshold value, with the
threshold depending on the azimuthal mode number $m$ (as summarized in figure
14). Our computations also suggest that trapped modes will asymptotically
approach a neutrally stable state in the limit of large Froude number without
becoming unstable, although the exact asymptotic scaling of $\Imag(\omega)$
for large $F$ was not obvious from the numerical analysis carried out here.
Interestingly, our numerics finds all wave modes to be stable, unlike in
either bounded vortex flows (Tophøj et al., 2013; Mougel et al., 2014, 2017)
or unbounded shallow water vortex flows (Ford, 1994). Finally, we have
confirmed what has been seen in the initial motivating experiment of the pool,
namely that trapped surface wave modes exist, are stable, and can propagate
with or against the base swirling flow (as also summarized in figure 14).

Based on the results presented here, we may speculate on the underlying
physical mechanism causing this trapping behaviour. One likely candidate might
have been the non-uniform shape of the free surface, with trapped modes
appearing to be localized at a radius close to the largest gradient of the
undeformed free surface. However, in appendix D results are shown for
numerical simulations with an (artificial) flat free surface, and almost
identical results are obtained, suggesting that, at least for the parameters
investigated here, the free surface shape is less important than the base
flow. We hypothesise, therefore, that it is the refraction of the surface
gravity waves by the base flow, causing the waves to curve around the vortex
with exactly the curvature needed to remain trapped, that is the underlying
physical mechanism. This would explain why the non-shallow-water assumption is
important here, as the dispersive nature of the waves in this case means a
frequency can be chosen to obtain the particular wave group velocities needed
for wave trapping, where as in the shallow water case the waves are
nondispersive and have a prescribed group velocity independent of their
frequency.

To numerically simulate a (horizontally-)unbounded fluid on a bounded
numerical domain, a far-field nonreflecting boundary condition or buffer
region is needed. Here, a novel additional term is introduced into the
governing equations, to provide damping of the surface waves in the buffer
region only. This method has proved more accurate than any non-reflecting
boundary condition we implemented while remaining computationally viable.
Indeed, introduction of additional unknowns in the mathematical formulation
was not needed, thus overcoming the main drawback of PML methods. Furthermore,
given that the background vortex flow vanishes in the far-field, the same
absorbing layer formulation can be employed with other similar vortex
distributions.

The numerical expense of our eigenvalue problem (which is a two-dimensional
spatial problem involving both $r$ and $z$ coordinates) might be reduced by
investigating a one-dimensional approximation along the radial coordinate
only; for example in either the deep- or shallow-water limits. Moreover, new
methods and formulations for the imposition of a non-reflecting boundary
condition in the far-field would certainly lead to a saving in the
computational time as it will avoid, for example, the need to numerically
resolve the unphysical buffer region.

Due to the computational expense, we have left to future studies the extension
of our parametric study of waves to extreme Froude numbers. In particular it
would be interesting to investigate the trend of modes for very small and very
large Froude numbers. Also, consideration of lower or higher aspect ratios
$h_{\infty}$ could be investigated. Throughout this work, we have assumed a
harmonic dependence $\exp\\{-\mathrm{i}\omega t+\mathrm{i}m\theta\\}$
(equation 9) and looked for modal solutions, and alternative approaches, such
as time-domain simulation, may give additional insight. This is particularly
in relation to the critical layer, which is where there is a radial location
$r=r_{c}$ for which $D_{t}=-\mathrm{i}\omega+\mathrm{i}mF\Omega_{0}(r_{c})=0$;
this possibility has been neglected here, and investigating it would require a
different numerical and mathematical method (such as has been done for
aeroacoustic waves; e.g. King et al., 2022). Even with the assumption of
harmonic dependence $\exp\\{-\mathrm{i}\omega t+\mathrm{i}m\theta\\}$, a
different but related problem would be to investigate the scattering of an
incoming wave encountering the vortex, which would require a different far-
field boundary condition to introduce the incoming wave as well as to allow
outgoing waves to propagate through the far-field boundary without reflection.
Finally, our model also neglects both nonlinear and surface tension effects;
while this is justified for the swimming pool application we model here, these
assumptions break down for either large amplitudes or short wavelengths.
Therefore, it would be interesting to investigate whether their inclusion
could lead to an instability of the base vortex flow.

###### Acknowledgements.

Acknowledgements. E. Zuccoli was funded through the Warwick Mathematics
Institute Centre for Doctoral Training, and gratefully acknowledges the
support of the University of Warwick and the UK Engineering and Physical
Sciences Research Council (EPSRC grant EP/W523793/1). E.J. Brambley gratefully
acknowledges the support of the UK Engineering and Physical Sciences Research
Council (EPSRC grant EP/V002929/1). D. Barkley gratefully acknowledges
supported from the Simons Foundation (Grant No. 662985, NG). The authors are
grateful to J.M. Skipp for the initial motivation and videos motivating this
project. Declaration of Interests. The authors report no conflict of interest.

##  Appendix A Detailed derivation and analysis of the Numerical Solution

###  A.1 Numerical solution of the eigenvalue problem

In order to solve (12,13,15), we used a Galerkin spectral method and expand
with a combination of Legendre polynomials as basis functions in both the
radial and axial coordinate in order to satisfy the Dirichlet boundary
conditions. First of all, to account for the shape of the computational domain
with a variable surface height $h_{0}(r)$, we remap the domain from the
physical domain $D=[0,R]\times[0,h_{0}(r)]$ to the computational square domain
$S=[-1,1]\times[-1,1]$ using the non-orthogonal transformation

$\displaystyle x$ $\displaystyle=\frac{2}{R}r-1,$ $\displaystyle y$
$\displaystyle=2\frac{z}{h_{0}(r)}-1.$ (18)

In this new domain $S$, the problem becomes

$\displaystyle\mathrm{i}mF\Omega_{0}(x)u-2F\Omega_{0}(x)v+\frac{2}{R}\\!\left[\frac{\partial\phi}{\partial
x}-\frac{h^{\prime}_{0}(x)}{h_{0}(x)}(y+1)\frac{\partial\phi}{\partial
y}\right]\\!=\mathrm{i}\omega u,$ (19a)
$\displaystyle\mathrm{i}mF\Omega_{0}(x)v+\frac{2F}{R(x+1)}\frac{d}{dx}\left[(x+1)U_{0}(x)\right]u+\frac{2\mathrm{i}m}{R(x+1)}\phi=\mathrm{i}\omega
v,$ (19b)
$\displaystyle\mathrm{i}mF\Omega_{0}(x)w+\frac{2}{h_{0}(x)}\frac{\partial\phi}{\partial
y}=\mathrm{i}\omega w,$ (19c) $\displaystyle\begin{aligned}
\frac{2}{R(x+1)}\frac{\partial}{\partial
x}\big{[}(x+1)u\big{]}&+\frac{2}{R}\frac{h^{\prime}_{0}(x)}{h_{0}(x)}u+\frac{2\mathrm{i}m}{R(x+1)}v\\\
&+\frac{2}{h_{0}(x)}\frac{\partial}{\partial
y}\\!\left[w-\frac{h^{\prime}_{0}(x)}{R}(y+1)u\right]\\!+\xi(x)\phi=0,\end{aligned}$
(19d)
$\displaystyle\mathrm{i}mF\Omega_{0}(x)\phi(x,1)+\frac{2h_{0}^{\prime}(x)}{R}u(x,1)-w(x,1)=\mathrm{i}\omega\phi(x,1),$
(19e) $\displaystyle
w(x,-1)=0,\qquad\quad\phi(1,y)=0,\qquad\quad|\boldsymbol{u}(-1,y)|<\infty,\qquad\quad|\phi(-1,y)|<\infty.$
(19f)

For $R_{c}\leq r\leq R$, we include a damping layer with a nonzero $\xi(r)$.
Here, we choose a simple form for $\xi(r)$, which we demonstrate in section
3.4 works well for the parameters we consider here:

$\xi(r)=\begin{cases}0,&r<R_{c},\\\
\bar{\xi}\left(\dfrac{r-R_{c}}{R-R_{c}}\right)^{2},&R_{c}\leq r\leq
R.\end{cases}$ (20)

Since $\xi(r)$ is continuously differentiable at $r=R_{c}$, no special
transmission conditions are needed there (Sim, 2010, p. 46).

We next obtain the weak discretized formulation of this problem. The basis
functions are chosen such that they identically satisfy the homogeneous
Dirichlet boundary data. Also, the specific expansion of each variable in
terms of basis functions will depend on the value of the azimuthal wavenumber,
as this is associated with a singularity at the origin for the pressure and
the aximuthal velocity component. In the following we make the distinction
between the case of axisymmetric perturbations ($m=0$), and generic non-
axisymmetric perturbations ($m\neq 0$).

####  A.1.1 Discrete problem for non-axisymmetric perturbations

In the general case of a nonzero value of $m$, the weak formulation of the
eigenvalue problem can be obtained in the following way. First of all, we
notice the presence of singular terms at $r=0$ for $u,v,\phi$ so we will
require these function to be null at the origin. Consequently, even $w$ will
be so. Thus, we look at the unknowns in the following spaces:

$\displaystyle[u,v,\phi]\in V_{H}(S)=\\{(u,v,\phi)\in
H^{1}(S):(u,v,\phi)=0,\quad at\quad r=0\\},$ (21a) $\displaystyle w\in
V_{v}(S)=\\{w\in H^{1}(S):w=0,\quad at\quad r=0,\quad z=0\\},$ (21b)

where $H^{1}(S)$ is the usual Sobolev space (Quarteroni, 2009, ch. 2). We
multiply the azimuthal component of the momentum equation and the continuity
equation by $(x+1)$. Then, by multiplying each of the equations by suitable
test functions $(v_{x},v_{t},v_{y},q)$ in the same space as the corresponding
unknowns, after integrating over the square $S$ and exploiting the boundary
conditions, the weak formulation of the differential problem reads: find
$([u,v,\phi],w)\in V_{H}(S)\times V_{v}(S)$ such that
$\forall([v_{x},v_{t},q],v_{v})\in V_{H}(S)\times V_{v}(S)$ the following
holds:

$\displaystyle\int_{S}v_{x}\mathrm{i}mF\Omega_{0}(x)u-\int_{S}v_{x}2F\Omega_{0}(x)v-\int_{S}\frac{2}{R}\frac{\partial
v_{x}}{\partial
x}\phi-\int_{S}v_{x}\frac{2}{R}\frac{h^{\prime}_{0}(x)}{h_{0}(x)}(y+1)\frac{\partial\phi}{\partial
y}=\mathrm{i}\omega\int_{S}v_{x}u,$ (22a)
$\displaystyle\int_{S}v_{t}\left[\mathrm{i}mF\Omega_{0}(x)(x+1)v+\frac{2F}{R}\\!\left[(x+1)U_{0}(x)\right]\\!u+\frac{2\mathrm{i}m}{R}\phi\right]=\mathrm{i}\omega\int_{S}v_{t}(x+1)v,$
(22b)
$\displaystyle\int_{S}v_{y}\left[\mathrm{i}mF\Omega_{0}(x)w+\frac{2}{h_{0}(x)}\frac{\partial\phi}{\partial
y}\right]=\mathrm{i}\omega\int_{S}v_{y}w,$ (22c) $\displaystyle\begin{aligned}
\int_{S}&\frac{2}{R}q\frac{\partial}{\partial
x}\left[(x+1)u\right]+\int_{S}\frac{2}{R}\frac{h^{\prime}_{0}(x)}{h_{0}(x)}(x+1)\left[qu+\frac{\partial
q}{\partial y}u\right]+\int_{S}\frac{2\mathrm{i}m}{R}qv\\\
&-\int_{S}\frac{2}{h_{0}(x)}(x+1)\frac{\partial q}{\partial
y}w+\int_{x=-1}^{x=1}\frac{2\mathrm{i}mF(x+1)\Omega_{0}(x)}{h_{0}(x)}q(x,1)\phi(x,1)\\\
&=\mathrm{i}\omega\int_{x=-1}^{x=1}\frac{2(x+1)}{h_{0}(x)}q(x,1)\phi(x,1).\end{aligned}$
(22d)

Let us define the bilinear forms $\mathcal{A}:V_{H}\times
V_{v}\rightarrow\mathbb{R}$ and $\mathcal{B}:V_{H}\times
V_{v}\rightarrow\mathbb{R}$ such that the generalised eigenvalue problem above
can be compactly written as: find $([u,v,\phi],w)\in V_{H}\times V_{v}$ such
that

$\displaystyle\mathcal{A}([u,v,w,\phi],[v_{x},v_{t},v_{y},q])=\omega\mathcal{B}([u,v,w,\phi],[v_{x},v_{t},v_{y},q]),\quad\forall([v_{x},v_{t},q],v_{y})\in
V_{H}\times V_{v}$ (23)

At this point, in order to discretize the problem, we need to expand the
unknowns in terms of proper basis functions. Such basis functions are taken in
such a way the homogeneous Dirichlet boundary conditions are automatically
satisfied in both the axial and radial coordinates. In particular, we define
the following set of polynomials, $P_{n}^{\star}(x)$ as follows:

$P_{n}^{\star}(x)=P_{n}(x)+P_{n-1}(x),\quad n\geq 1,$ (24)

where $P_{n}(x)$ are standard Legendre polynomials.

Thanks to the above definition, a Dirichlet boundary condition at $x=-1$ is
automatically satisfied. Therefore, we expand the velocity components and the
pressure as

$\displaystyle[u,v,\phi](x,y)=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}[u_{ij},v_{ij},\phi_{ij}]P_{i}^{\star}(x)P_{j}(y),$
(25) $\displaystyle
w(x,y)=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}w_{ij}P_{i}^{\star}(x)P_{j}^{\star}(y).$

After substituting the expansion into the weak formulation (23), we obtain the
final discrete eigenvalue problem, which reads

$\displaystyle\mathsfbi{A}\boldsymbol{w}=\omega\mathsfbi{B}\boldsymbol{w},$
(26)

where $\boldsymbol{w}=(u_{ij},v_{ij},w_{ij},\phi_{ij})$ is our array
containing the spectral coefficients previously defined, $\mathsfbi{A}$ and
$\mathsfbi{B}$ are the matrices of order $4N_{x}N_{y}\times 4N_{x}N_{y}$coming
from the evaluation of the integrals appearing in the weak formulation for
each equation.

####  A.1.2 Discrete problem for axisymmetric perturbations

In the axisymmetric case ($m=0$), the differential problem (19) over the
square $S$ simplifies, and it can be noted that only the radial component of
the velocity must go to zero as $r\rightarrow 0$, with the other unknowns
allowed to take any finite value. Thus, studying the axisymmetric perturbation
problem means to look for the eigensolutions in the following spaces: $u\in
V_{H}(S)$, $(v,\phi)\in H^{1}(S)$ and $w\in V_{v0}=\\{w\in H^{1}(S):w=0,\quad
at\quad z=0\\}$. The weak formulation is then obtained as shown previously.
Regarding the disctetization process, we express the four unknowns (and
corresponding test functions) as follows

$\displaystyle
u(x,y)=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}u_{ij}P_{i}^{\star}(x)P_{j}(y),$
(27)
$\displaystyle[v,\phi](x,y)=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}[v_{ij},\phi_{ij}]P_{i}(x)P_{j}(y),$
$\displaystyle
w(x,y)=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}w_{ij}P_{i}(x)P_{j}^{\star}(y).$

After substitution of the expansion above into the weak formulation, the
resulting algebraic generalised eigenproblem (26) is recalled in the case
$m=0$.

###  A.2 Spurious numerical modes and resolvedness conditions

Solving the discretized eigenvalue problem (26) is known to give rise to
spurious eigenvalues and eigenfunctions. These are numerical artifacts caused
by numerical under-resolution of highly-oscillatory modes; that is,
eigenfunctions do satisfy the discretized problem but do not approximate
solutions of the continuous eigenvalue problem. Following Brambley & Peake
(2008) (see also Brambley, 2007, p. 78–81), we implement two tests to remove
such spurious numerical modes.

The first condition requires that the eigenvalues do not move significantly
when the numerical discretization is changed. Suppose eigenvalues
$\\{\omega_{p}\\}$ have been computed using a discretization $(N_{x},N_{y})$,
and suppose $\\{\hat{\omega}_{q}\\}$ have been computed with a slightly
increased resolution, such as $(N_{x}+1,N_{y}+1)$. The first condition is then
given by requiring that, for each eigenvalue $\omega_{p}$,

$\min_{q}|\omega_{p}-\hat{\omega}_{q}|<\text{tol},$ (28)

where tol is a prescribed tolerance. For the results computed in this paper
$\text{tol}=10^{-2}$ was used.

The second resolvedness condition involves the spectral coefficients of the
eigenfunctions. In particular, we require the modulus of the spectral
coefficients to decay smoothly as the number of polynomials used gets larger.
Both the aforementioned conditions have been applied to the pressure
eigenfunctions as this is the unknown of most interest because it contains
information on the shape of the free surface height. So, let $\phi_{ij}$ be a
generic spectral coefficient, the following condition has been implemented

$\displaystyle\frac{\sup_{(i,j)\in\mathcal{B}}|\phi_{ij}|}{\sup_{\mathrm{all}\,i,j}|\phi_{ij}|}<\text{tol},$
(29)

where $\mathcal{B}$ is a part of the domain in the $N_{x}-N_{y}$ plane,
defined as

$\mathcal{B}=\big{\\{}(i,j):N_{x}+1-b_{x}\leq i\leq N_{x}+1\text{ or
}N_{y}+1-b_{y}\leq j\leq N_{y}+1\big{\\}},$ (30)

tol is a prescribed tolerance on the magnitude of the coefficients and
$b_{x},b_{y}$ two prescribed borders widths in the $N_{x}-N_{y}$ plane.

Figure 20: Example of a well-resolved mode (Top) and a poorly-resolved mode
(Bottom).

Figure 20 compares the magnitude of the coefficients for a well resolved mode
and a poorly-resolved mode for $N_{x}=50$, $N_{y}=20$, $b_{x}=12$, and
$b_{y}=4$, $\text{tol}=10^{-1}$. The well-resolved mode can be seen to have
spectral coefficient $|\phi_{ij}|$ decreasing exponentially quickly. By
contrast, for the poorly-resolved mode, the spectral coefficients
$|\phi_{ij}|$ are oscillatory increasing as function of the number of
polynomials used.

##  Appendix B Damping layer formulation

In section 3.1, the continuity equation was modified by the introduction of
the “damped compressibility” term $\xi(r)\phi$ in equation (15) in order to
emulate an infinite domain using our finite computational domain, by damping
out disturbances as they reach the unphysical computational boundary at $r=R$.
Here, we justify the inclusion of that term by analogy with the equations of
acoustics.

Small (linear) homentropic perturbations $(\boldsymbol{u},p,\rho)$ to a static
fluid $(\boldsymbol{0},p_{0},\rho_{0})$ are governed by the linearized Euler
Equations:

$\displaystyle\rho_{0}\frac{\partial\boldsymbol{u}}{\partial t}+\nabla p$
$\displaystyle=0,$ $\displaystyle\frac{1}{c^{2}}\frac{\partial p}{\partial
t}+\rho_{0}\nabla\cdot\boldsymbol{u}$ $\displaystyle=0,$ $\displaystyle
p=c^{2}\rho,$ (31)

where $c$ is the speed of sound. These can be combined into a single wave
equation for the pressure,

$\frac{1}{c^{2}}\frac{\partial^{2}p}{\partial t^{2}}-\nabla^{2}p=0.$ (32)

Following Gao et al. (2017, pp. 81-82), in order to introduce a sponge layer
to damp outgoing waves, we add a damping term of the form
$\xi(r)\partial_{t}p$ into the wave equation (32),

$\frac{1}{c^{2}}\frac{\partial^{2}p}{\partial t^{2}}+\xi(r)\frac{\partial
p}{\partial t}-\nabla^{2}p=0.$ (33)

This wave equation can be split back into the original physical mass- and
momentum-equations in the original physical variables
$(\boldsymbol{u},p,\rho)$ as

$\displaystyle\rho_{0}\frac{\partial\boldsymbol{u}}{\partial t}+\nabla p$
$\displaystyle=0,$ $\displaystyle\frac{1}{c^{2}}\frac{\partial p}{\partial
t}+\xi(r)p+\rho_{0}\nabla\cdot\boldsymbol{u}$ $\displaystyle=0,$
$\displaystyle p=c^{2}\rho,$ (34)

If we now take the incompressible limit $1/c^{2}\to 0$ in system (34), we
obtain the modified continuity equation (15) introduced in section 3.1.

Figure 21

Figure 21: Comparison of the structure of the eigenmode $|h(r,0,0)|$ for
$m=7$, $F=0.3$ and different value of the compressibility $\alpha=1/c^{2}$.

and table 3

$\alpha$ | $0$ | $10^{-4}$ | $10^{-2}$ | $10^{-1}$ | $1$
---|---|---|---|---|---
$\Real(\omega)$ | -1.3769 | -1.3769 | -1.3759 | -1.3662 | -1.2711
$\Imag(\omega)$ | -0.0138 | -0.0138 | -0.0139 | -0.0147 | -0.0241
Table 3: Eigenvalues as function of the compressibility
$\alpha=\frac{1}{c^{2}}$ for a radiating mode, for $m=7$, and $F=0.3$.

show the effect of the compressibility $\alpha=1/c^{2}$ on the structure of
the radiating eigenmode and the eigenvalue respectively, for parameters $m=7$
and $F=0.3$. From both, it is clear that the eigenmode computed does not vary
significantly as $\alpha\to 0$, and so we conclude that the damping properties
of the $\xi(r)p$ term are carried over in the limit $\alpha\to 0$.

##  Appendix C Code validation and comparison with results for a Solid-Body
Rotation background flow

Here, we provide evidence of the correctness of our numerical method by means
of results obtained for a simpler Newton’s bucket problem, consisting of a
solid body rotation background flow of the form
$\boldsymbol{U_{0}}=\bar{\Omega}r\boldsymbol{\hat{\theta}}$, with rigid walls
at $r=R$ and no damping region, as analysed by Mougel et al. (2015). In this
appendix, we adopt the same non-dimensionalization as Mougel et al. (2015), so
that we can directly compare numerical solutions. Results from our numerics
are plotted in figure 22

Figure 22: Eigenmodes comparison with Newton’s Bucket Problem Mougel et al.
(2015) using their same parameters. (a) $\lambda=4.289$. (b): $\lambda=7.435$.
(c): $\lambda=9.043$. (d): $\lambda=-5.201$. (e): $\lambda=-7.458$. (f):
$\lambda=-9.034$. (g): $\lambda=1.407$. (h): $\lambda=1.097$. (i):
$\lambda=1.783$. (j): $\lambda=-1.698$. (k): $\lambda=-1.329$. (l):
$\lambda=-1.914$. (m): $-0.03723$. (n): $\lambda=-0.01422$. (o):
$\lambda=-0.007505$. Here the number of points in the radial and axial
direction is $N_{x}=10$ and $N_{y}=20$ respectively.

and can be compared with figure 3 of Mougel et al. (2015, p. 223). Excellent
agreement is seen, giving further confidence in the correctness of our
numerics here.

##  Appendix D Results without the base free surface deformation

In this appendix, we demonstrate that the dominant contribution to the
trapping of modes by the vortex shown in section §4 is the swirl of the base
flow, and that the deformation of the base free surface has little effect
itself. In the following we artificially assume the base free surface to be
flat, meaning that $h_{0}(r)=h_{\infty}$ and $h^{\prime}_{0}(r)=0$. All the
remaining terms in the governing equations and boundary conditions do not
change. Solving the eigenvalue problem in this case, we compare the
eigenvalues to the complete case where the base free surface deformation has
been taken into account in figure 23.

Figure 23: Imaginary and real part of the eigenvalues as function of the
Froude number for modes with $m=7$ and $n=1$. Blue-circle curve: full
numerics; red-asterisk curve: without the free surface contribution.

Little difference is seen between the two. Indeed, figure 24

Figure 24: Absolute difference between eigenvalues computed with and without
the free surface contribution, $|\omega-\omega_{\mathrm{flat}}|$, for the
first radial mode $n=1$ and $m=7$.

plots the difference in eigenvalues as function of the Froude number; i.e.
$|\omega-\omega_{\mathrm{flat}}|$, where $\omega_{\mathrm{flat}}$ are the
eigenvalues without the free surface variation. It is clear that for the range
of Froude numbers considered, the two agree remarkably closely.

In terms of the eigenfunctions, the free surface perturbation is displayed in
figure 25,

Figure 25: Plots of the perturbation free surface height
$h(r,\theta,t=0)=\Real[\phi(r,h_{0}(r))\exp\\{\mathrm{i}m\theta\\}]$ for $m=7$
without the base free surface deformation, to be compared against figure 7.
(a) $F=0.3$. (b) $F=0.32$. (c) $F=0.33$. (d) $F=0.34$. (e) $F=0.36$. All modes
displayed rotate clockwise, i.e. against the vortex flow.

which should be compared against the free surface perturbation with a base
free surface height variation plotted in figure 7; the two figures can be seen
to be almost identical. The corresponding structure in the $r-z$ plane is
shown in figure 26,

Figure 26: Plots of the pressure distribution $|\phi_{m}(r,z)|$ for $m=7$
without the base free surface deformation, to be compared against figure 8.
(a) $F=0.3$. (b) $F=0.32$. (c) $F=0.33$. (d) $F=0.34$. (e) $F=0.36$.

which should be compared against figure 8 in the results section; again, the
two figures are almost identical apart from the base flow height.

We therefore conclude that the driving mechanism behind the trapping of modes
by the vortex is the rotation of the base flow, and not the free surface
deformation.

## References

  * Bayliss & Turkel (1980) Bayliss, A. & Turkel, E. 1980 Radiation boundary conditions for wave-like equations. Communications on Pure and Applied Mathematics 33, 707–725.  doi: 10.1002/cpa.3160330603.
  * Berenger (1994) Berenger, J. P. 1994 A perfectly matched layer for the absorprion of electromagnetic waves. J. Comput. Phys. 114, 185–200.  doi: 10.1006/jcph.1994.1159.
  * Brambley (2007) Brambley, E. J. 2007 The acoustics of curved and lined cylindrical ducts with mean flow. PhD thesis, University of Cambridge.  doi: 10.17863/CAM.16086.
  * Brambley & Peake (2008) Brambley, E. J. & Peake, N. 2008 Sound transmission in strongly curved slowly varying cylindrical ducts with flow. J. Fluid Mech. 596, 387–412.  doi: 10.1017/S0022112007009603.
  * Chew & Liu (1996) Chew, W. C. & Liu, Q. H. 1996 Perfectly matched layers for elastodynamics: A new absorbing boundary condition. Journal of Computational Acoustics 4, 341–359.  doi: 10.1142/S0218396X96000118.
  * Drazin & Riley (2006) Drazin, P. G. & Riley, N. 2006 The Navier-Stokes equations: a classification of flows and exact solutions. Cambridge University Press.  ISBN 9780521681629.
  * Engquist & Majda (1977) Engquist, B. & Majda, A. 1977 Absorbing boundary conditions for the numerical simulation of waves. Mathematics of Computation 31, 629–651.  doi: 10.1090/S0025-5718-1977-0436612-4.
  * Fabre et al. (2006) Fabre, D., Sipp, D. & Jacquin, L. 2006 Kelvin waves and the singular modes of the lamb-oseen vortex. J. Fluid Mech. 551, 235–274.  doi: 10.1017/S0022112005008463.
  * Ford (1994) Ford, R. 1994 The instability of an axisymmetric vortex with monotonic potential vorticity in rotating shallow water. J. Fluid Mech 280, 303–334.  doi: 10.1017/S0022112094002946.
  * Gao et al. (2017) Gao, Y., Song, H., Zhang, J. & Yao, Z. 2017 Comparison of artificial absorbing boundaries for acoustic wave equation modelling. Exploration Geophysics 48, 76–93.  doi: 10.1071/EG15068.
  * Greenspan (1969) Greenspan, H. P. 1969 The theory of rotating fluids. Cambridge University Press.  ISBN 9780521051477.
  * Grote (2000) Grote, M. J. 2000 Non-reflecting boundary conditions for time dependent wave propagation. SAM Research Report 2000-04. ETH Zurich.  doi: 10.3929/ethz-a-004288373.
  * Grote & Keller (1995) Grote, M. J. & Keller, J. B. 1995 Exact nonreflecting boundary conditions for the time dependent wave equation. SIAM J. Appl. Math. 55(2), 280–297.  doi: 10.1137/S0036139993269266.
  * Hastings et al. (1996) Hastings, F. D., Schneider, J. B. & Broschat, S. L. 1996 Application of the perfectly matched layer (pml) absorbing boundary condition to elastic wave propagation. Journal of Acoustical Society of America 100, 3061–3069.  doi: 10.1121/1.417118.
  * Hunt (2019) Hunt, M. 2019 Two-dimensional surface waves in magnetohydrodynamics. J. Plasma Phys 85.  doi: 10.1017/S0022377819000497.
  * Hunt et al. (2015) Hunt, M., Parau, E., Vanden-Broeck, J.-M. & Papageorgiou, D. 2015 Linear and nonlinear surface waves in electrohydrodynamics, arXiv: 1501.02783.  doi: 10.48550/arXiv.1501.02783.
  * Hunt & Dutykh (2021) Hunt, M. J. & Dutykh, D. 2021 Free surface flows in electrohydrodynamics with a constant vorticity distribution. Water Waves 3, 297–317.  doi: 10.1007/s42286-020-00043-9.
  * Hunt & Vanden-Broeck (2015) Hunt, M. J. & Vanden-Broeck, J.-M. 2015 A study of the effects of electric field on two-dimensional inviscid nonlinear free surface flows generated by moving disturbances. Journal of Engineering Mathematics 92, 1–13.  doi: 10.1007/s10665-014-9766-6.
  * John (1978) John, F. 1978 Partial differential equations. Springer.  ISBN 9781468400618.
  * Johnson (1997) Johnson, R. S. 1997 A Modern Introduction to the Mathematical Theory of Water Waves. Cambridge University Press.  ISBN 9780521598323.
  * Keller & Givoli (1989) Keller, J. B. & Givoli, D. 1989 Exact non-reflecting boundary conditions. Computational Physics 82, 172–192.  doi: 10.1016/0021-9991(89)90041-7.
  * King et al. (2022) King, M. J., Brambley, E. J., Liupekevicius, R., Radia, M., Lafourcade, P. & Shah, T. H. 2022 The critical layer in quadratic flow boundary layers over acoustic linings. J. Fluid Mech. 950, A8.  doi: 10.1017/jfm.2022.753.
  * Lindquist et al. (2012) Lindquist, J. M., Neta, B. & Giraldo, F. X. 2012 High-order non-reflecting boundary conditions for dispersive waves in polar coordinates using spectral elements. Applied Mathematics and Computation 218, 6666–6676.  doi: 10.1016/j.amc.2011.12.023.
  * McWilliams (2006) McWilliams, J. C. 2006 Fundamentals of Geophysical Fluid Dynamics. Cambridge University Press.  ISBN 9780521856379.
  * Mougel et al. (2014) Mougel, J., Fabre, D. & Lacaze, L. 2014 Waves and instabilities in rotating free surface flows. Mechanics and Industry 15, 107–112.  doi: 10.1051/meca/2014007.
  * Mougel et al. (2015) Mougel, J., Fabre, D. & Lacaze, L. 2015 Waves in newton’s bucket. J. Fluid Mech. 783, 211–250.  doi: 10.1017/jfm.2015.527.
  * Mougel et al. (2017) Mougel, J., Fabre, D., Lacaze, L. & Bohr, T. 2017 On the instabilities of a potential vortex with a free surface. J. Fluid Mech. 824, 230–264.  doi: 10.1017/jfm.2017.341.
  * Oliveira et al. (2014) Oliveira, L. A., Cardoso, V. & Crispino, L. C. B. 2014 Ergoregion instability: The hydrodynamic vortex. Phys. Rev. D 89, 124008.  doi: 10.1103/PhysRevD.89.124008.
  * Patrick (2019) Patrick, S. 2019 On the analogy between black holes and bathtub vortices. PhD thesis, University of Nottingham.  doi: 10.48550/arXiv.2009.02133.
  * Patrick et al. (2018) Patrick, S., Coutant, A. & Richartz, M. 2018 Black-hole quasibound states from a draining bathtub vortex flow. Phys. Rev. Lett. 121, 061101.  doi: 10.1103/PhysRevLett.121.061101.
  * Quarteroni (2009) Quarteroni, A. 2009 Numerical models for differential problems. Springer.  ISBN 9788847010710.
  * Riedinger et al. (2010) Riedinger, X., Dizes, S. Le & Meunier, P. 2010 Viscous stability properties of a lamb-oseen vortex in a stratified fluid. J. Fluid Mech. 645, 255–278.  doi: 10.1017/S002211200999262X.
  * Sim (2010) Sim, I. 2010 Non-reflecting boundary conditions for time-dependent wave propagation. PhD thesis, University of Basel.  doi: 10.5451/unibas-005269594.
  * Skipp (2020) Skipp, J. M. 2020 Single vortex ring. YouTube https://youtu.be/ApwCq7UkCfY, accessed 09/11/2021.
  * Söderstrom et al. (2010) Söderstrom, A., Karlsson, M. & Museth, K. 2010 A pml-based nonreflective boundary for free surface fluid animation. ACM Transactions on Graphics 29.  doi: 10.1145/1857907.1857912.
  * Tophøj et al. (2013) Tophøj, L., Mougel, J., Bohr, T. & Fabre, D. 2013 Rotating polygon instability of a swirling free surface flow. Physical Review Letters 110, 194502.  doi: 10.1103/PhysRevLett.110.194502.
  * Wellens & Borsboom (2020) Wellens, P. & Borsboom, M. 2020 A generating and absorbing boundary condition for dispersive waves in detailed simulations of free-surface flow interaction with marine structures. Computers and Fluids 200.  doi: 10.1016/j.compfluid.2019.104387.