Trapped free surface waves for a Lamb–Oseen vortex flow

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 modes 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 corotate with the flow.

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 phenomenon, together with a schematic representation, is shown in figure 1; interested readers are encouraged to follow the reference for a video.

Figure 1. The physical scenario motivating this work. Surface waves are trapped within a rotating vortex flow, but propagate in the direction opposite to the direction of rotation. (a) Photograph of two vortices in a swimming pool (a still taken from the video of Skipp (2020)). (b) 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.

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, Coutant & Richartz (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. The term ‘quasibound’ references that temporally weakly damped modes are not strictly bounded, and in fact may have weak exponential growth at infinity in order to satisfy a causal radiation condition (Briggs 1964; Bers 1983; Berti, Cardoso & Starinets 2009, pp. 17–18; Le Dizès & Billant 2009; Riedinger, Le Dizès & Meunier 2010); despite this weak growth, such modes are highly localized to a finite region (as discussed in § 4.1), and here we simply refer to these as trapped modes. 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 do not consider the related wave scattering problem.

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

(1.1)\begin{equation} \boldsymbol{U_0} = \boldsymbol{\hat{\theta}} \frac{\varGamma_0}{2{\rm \pi} r} (1-\exp\{{-r^2/a^2}\}), \end{equation}

where $r$ is the horizontal distance from the vortex centre, $\boldsymbol {\hat {\theta }}$ is a unit vector in the direction of rotation and $\varGamma _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 time scales compared with the viscous diffusion time scale, $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, Fabre & Lacaze 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, Sipp & Jacquin (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, Fabre & Lacaze (2015); Mougel et al. (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 non-dispersive 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) and Hunt & Dutykh (2021), although in the present work we do not consider any electric or magnetic 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 (Engquist & Majda 1977; Bayliss & Turkel 1980; Keller & Givoli 1989; Grote & Keller 1995; Grote 2000). 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, Neta & Giraldo 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 (Chew & Liu 1996; Hastings, Schneider & Broschat 1996) and fluid dynamics (Söderstrom, Karlsson & Museth 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 organized as follows. In § 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 § 3. The results of this numerical solution are described in § 4. Finally, resulting conclusions and possible future directions for continued research are highlighted in § 5.

2. Mathematical model

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

(2.1)\begin{equation} \frac{\partial\boldsymbol{U}}{\partial t} + \boldsymbol{U \cdot \nabla}\boldsymbol{U} + \frac{1}{\rho}\boldsymbol{\nabla} P + g\boldsymbol{\hat{z}} = 0,\quad \boldsymbol{\nabla \cdot U} = 0, \end{equation}

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

(2.2)\begin{equation} \frac{\partial H}{\partial t} + \boldsymbol{U \cdot \nabla} H = \boldsymbol{U\cdot \hat{z}} \quad \text{and}\quad P = \bar{P}\quad \text{on}\ z = H. \end{equation}

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

(2.3a–c)\begin{equation} \boldsymbol{U} = U_0(r)\boldsymbol{\hat{\theta}} + \varepsilon(u_r\boldsymbol{\hat{r}} + u_\theta\boldsymbol{\hat{\theta}} + u_z\boldsymbol{\hat{z}}), \quad P = P_0(r,z) + \varepsilon p,\quad H = h_0(r) + \varepsilon h. \end{equation}

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

(2.4a)$$\begin{gather} P_0(r, z) = \bar{P} + \rho g ( h_0(r) - z), \end{gather}$$

(2.4b)$$\begin{gather}h_{0}(r) = h_\infty - \frac{1}{g}\int_r^\infty \frac{U^{2}_{0}(r')}{r'}\,\mathrm{d} r', \end{gather}$$

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

(2.5)\begin{equation} U_{0}(r) = \frac{\varGamma_{0}}{2{\rm \pi} r}(1-\exp({-}r^2/a^2)), \end{equation}

where $a$ sets the radial size of the core and $\varGamma _{0}$ sets the circulation of the vortex. Note that, for the Lamb–Oseen vortex, for small $r$, we have $U_0(r) \approx r\varGamma _0/2{\rm \pi} a^2$, so that (2.5) is a solid body rotation near the centre of the vortex, while for large $r$ we have $U_0(r)\approx \varGamma _0/2{\rm \pi} r$, so that (2.5) 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.1) about the base solution $(U_{0}, P_{0}, h_{0})$ given above, the governing equations for the perturbation are

(2.6a)$$\begin{gather} D_{t}u_{r} - 2\varOmega_{0}(r)u_\theta + \frac{1}{\rho}\frac{\partial p}{\partial r} = 0, \end{gather}$$

(2.6b)$$\begin{gather}D_{t}u_{\theta} + \frac{1}{r}(rU_0(r))'u_{r} + \frac{1}{\rho r}\frac{\partial p}{\partial \theta} = 0, \end{gather}$$

(2.6c)$$\begin{gather}D_{t}u_{z} + \frac{1}{\rho}\frac{\partial p}{\partial z} = 0, \end{gather}$$

(2.6d)$$\begin{gather}\frac{1}{r}\frac{\partial}{\partial r}(ru_{r}) + \frac{1}{r}\frac{\partial u_{\theta}}{\partial \theta} + \frac{\partial u_{z}}{\partial z} = 0, \end{gather}$$

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

(2.7a)$$\begin{gather} u_{z} = \frac{1}{\rho g}D_{t}p + h'_{0}u_{r} \quad\text{and}\quad h = \frac{p}{\rho g} \quad\text{on}\ z = h_{0}(r), \end{gather}$$

(2.7b)$$\begin{gather}u_{z} = 0 \quad \text{on} \ z = 0. \end{gather}$$

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

(2.8a)$$\begin{gather} \boldsymbol{u} = [u(r, z)\boldsymbol{\hat{r}} + v(r, z)\boldsymbol{\hat{\theta}} + w(r, z)\boldsymbol{\hat{z}}]\exp(-\mathrm{i}\omega t + \mathrm{i} m\theta), \end{gather}$$

(2.8b)$$\begin{gather}p = \phi(r, z)\exp(-\mathrm{i}\omega t + \mathrm{i} m\theta), \end{gather}$$

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

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

(2.9)\begin{equation} \left.\begin{gathered} (r,\theta,z) = (a\tilde{r}, \tilde{\theta}, a\tilde{z}), \quad t = \sqrt{\frac{a}{g}}\tilde{t}, \quad U_{0}(r) = \frac{\varGamma_{0}}{2{\rm \pi} a}\widetilde{U}_{0}(\tilde{r}), \quad \varOmega_0 = \frac{\varGamma_{0}}{2{\rm \pi} a^2}\tilde{\varOmega_0},\\ h = a\tilde{h}, \qquad p = \rho ag\tilde{p}, \quad \omega = \sqrt\frac{g}{a}\tilde{\omega}, \quad \boldsymbol{u} = \sqrt{ag}\boldsymbol{\tilde{u}}. \end{gathered}\right\} \end{equation}

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

(2.10a,b)\begin{equation} F =\frac{\varGamma_0}{2{\rm \pi}\sqrt{ga^3}}, \quad \tilde{h}_\infty = \frac{h_\infty}{a}. \end{equation}

Here $F$ is the Froude number and sets the non-dimensionalized 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 (2.13c) 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 (2.13) 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 non-dimensional eigenvalue problem is

(2.11a)$$\begin{gather} ({-\mathrm{i}}\omega + \mathrm{i} m F \varOmega_{0}(r))u - 2F\varOmega_{0}(r)v + \frac{\partial \phi}{\partial r} = 0, \end{gather}$$

(2.11b)$$\begin{gather}({-\mathrm{i}}\omega + \mathrm{i} m F \varOmega_{0}(r))v + \frac{F}{r}(rU_0(r))'u + \frac{\mathrm{i} m}{r}\phi = 0, \end{gather}$$

(2.11c)$$\begin{gather}({-\mathrm{i}}\omega + \mathrm{i} m F \varOmega_{0}(r))w + \frac{\partial \phi}{\partial z} = 0, \end{gather}$$

(2.11d)$$\begin{gather}\frac{1}{r}\frac{\partial}{\partial r}(ru) + \frac{\mathrm{i} m}{r}v + \frac{\partial w}{\partial z} = 0, \end{gather}$$

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

(2.12a)$$\begin{gather} w = 0, \quad \text{on} \ z = 0, \end{gather}$$

(2.12b)$$\begin{gather}w = ({-\mathrm{i}}\omega + \mathrm{i} m F \varOmega_{0}(r))\phi + F^{2} r\varOmega^{2}_{0}(r)u, \quad \text{on} \ z = h_{0}(r). \end{gather}$$

For the Lamb–Oseen vortex, in dimensionless terms

(2.13a,b)$$\begin{gather} U_0(r) = \frac{1 - \exp({-}r^2)}{r}, \quad \varOmega_0(r) = \frac{U_0(r)}{r}, \end{gather}$$

(2.13c)$$\begin{gather}h_{0}(r) = h_{\infty} - F^2 \int_r^\infty \frac{U^{2}_{0}(r')}{r'}\,\mathrm{d} r'. \end{gather}$$

One consequence of the harmonic assumption (2.8) 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\varOmega _0(r_c) = 0$. It is not immediately obvious from (2.11) 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 have no difficulty in this case. However, in fact the critical layer is related to behaviour that is not of the harmonic form assumed in (2.8); investigating the effect of the critical layer requires a different mathematical and numerical technique (such as Frobenius 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 three-dimensional 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 NRBC at this boundary. However, exact NRBCs are generally derived using the method of characteristics, but this is not available to us since (2.11d) 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 re-enter 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 (2.11d), which then becomes

(3.1)\begin{equation} \xi(r)\phi + \frac{1}{r}\frac{\partial}{\partial r}(ru) + \frac{\mathrm{i} m}{r}v + \frac{\partial w}{\partial z} = 0, \end{equation}

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

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

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$. We show in § 4.2 that this leads to the correct causal behaviour. A representation of the damping function as well as the subdivision of the numerical domain is shown in figure 3.

Figure 3. (a) Profile of the quadratic damping function defined in (3.2) with parameters: $\bar {\xi } = 5$, $R_c = 5$ and $R = 10$. (b) Schematic of the computational domain in the $r$–$z$ plane.

The governing equations (2.11a–c) together with the modified continuity equation (3.1) and boundary conditions (2.12) form an eigenvalue problem to solve for the allowable frequencies $\omega$ permitting a non-zero modal solution.

3.2. Numerical discretization

In order to solve (2.11a–c, 2.12, 3.1), 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

(3.3)\begin{equation} \boldsymbol{\mathsf{A}}\boldsymbol{w} = \omega \boldsymbol{\mathsf{B}}\boldsymbol{w}, \end{equation}

with $\boldsymbol {w} = (u_{ij}, v_{ij}, w_{ij}, \phi _{ij})$ representing our array containing the spectral coefficients of each unknown, $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{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. We sort eigenvalues by imaginary part and focus on the leading eigenvalues, i.e. those with the most positive imaginary part. These correspond to unstable modes, should there be any, or the least damped stable modes, if not.

Not all solutions to the discretized problem (3.3) correspond to solutions to the continuous problem being approximated, however. To remove under-resolved eigenmodes and spurious eigenvalues, the numerical solutions to (3.3) are filtered, as described in the Appendix A.2. The numerical scheme has been validated against three published cases: two with walls at finite radius (Mougel et al. 2014, 2015), and one unbounded but in shallow water (Ford 1994). We obtain unstable modes where they are known to exist. Details are given 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 $-\textrm {Im}(\omega )$ is the decay rate of the mode, variations in damping typically have a strong effect on $\textrm {Im}(\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 (A11) 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 NRBC. For varying magnitudes of damping $\bar {\xi }$, the radiating eigenfunction is plotted in figure 4, and the corresponding eigenvalues are given in table 1. 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 $\textrm {Im}(\omega )$ to variations in damping strength, as expected.

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. (a) A small computational domain $R_{c} = 5$, $R = 10$. (b) A large computational domain $R_c=15$, $R=30$.

Table 1. Eigenvalues as function of the amount of damping for $m=7$ and $F = 0.3$, $R_{c} = 5$ and $R = 10$.

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, and the corresponding eigenvalues in table 2. 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 $\textrm {Im}(\omega )$ being particularly sensitive.

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$ (a) and a larger domain $R=30$ (b).

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$.

Finally, in figure 6 we compare the radiating mode for $m = 7$, $F = 0.3$ computed using two different sizes 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.

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.

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 with $\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 detail 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 wavenumbers 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 wavenumber.

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, p. 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 functions of 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 7(a), and a long tail in the radial direction, as shown in figure 8(a). 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. (An animation of the radiating mode in figures 7(a) and 8(a) is given as Movie 1 in the online Supplementary material available at https://doi.org/10.1017/jfm.2024.645.) In contrast, trapped modes are confined to a bounded region in the radial direction, and without the spiralling 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 7(e) and 8(e). (An animation of the the trapped mode in figures 7(e) and 8(e) is given as Movie 2 in the online Supplementary material.) 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.

Figure 7. Plots of the free surface height $h(r, \theta, t=0) = \textrm {Re}[\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. Animations of (a) and (e) are given in Movies 1 and 2 in the online Supplementary material.

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$. Animations of (a) and (e) are given in Movies 1 and 2 in the online Supplementary material.

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 where the modulus of the free surface height is plotted as a function of $r$ for different Froude numbers.

Figure 9. Modulus of the free surface height along the radius for $m = 7$ in the transition regime $F \in [0.3, 0.36]$ from radiating to trapped modes. The trapped mode at $F = 0.5$ is also shown for comparison.

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. In this way each eigensolution will be indexed by both $m$ and $n$, with corresponding eigenvalues $\omega = \omega _{mn}$. In figure 11(a), the two types of modes can be seen: radiating modes with a significant negative $\textrm {Im}(\omega )$, and trapped modes having an almost null $\textrm {Im}(\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 § 4.3). 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 hypothesize 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 (Fabre et al. 2006), even in the presence of a free surface.

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.

Figure 11. Imaginary (a) and real (b) parts of the eigenvalue $\omega _{mn}$ as function of the Froude number $F$ for $m = 7$ and varying $n$.

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 11(b), which shows how the real part of the eigenvalues varies with the Froude number, again for $m = 7$. For $\textrm {Re}(\omega ) < 0$, waves rotate opposite to the base flow (counter-rotating waves), while for $\textrm {Re}(\omega ) > 0$, waves rotate in the same direction as the base flow (corotating waves). It is clear that for each $n$, there is a value of Froude number separating counter-rotating from corotating surface waves.

4.2. Far-field behaviour and causality

In this section, we investigate the behaviour of our numerical solutions in the far field, and show that they have the correct behaviour expected of causal outgoing waves. Away from the vortex core, in the far field, the behaviour of surface waves is well understood. We first consider the behaviour of deep-water waves in the far-field (see, e.g. Batchelor 2000, p. 371)

(4.1)\begin{equation} h \propto r^{{-}1/2}\exp\{{\pm}\mathrm{i} kr\} \quad \text{with the dispersion relation}\ \omega^2=k. \end{equation}

For a complex frequency $\omega = \omega _r + \mathrm {i}\omega _i$, the dispersion relation gives $\textrm {Im}(k)=2\omega _r\omega _i$. Causality (Briggs 1964; Bers 1983; Berti et al. 2009, pp. 17–18; Le Dizès & Billant 2009; Riedinger et al. 2010) implies that the outgoing wave must decay when $\omega _i>0$, so we must take $\pm = \mathrm {sgn}(\omega _r)$. Thus, for temporally damped modes when $\omega _i<0$ we expect spatially growing behaviour in the far field. This growth in the far field is found in our numerical solutions, as shown in figure 12(a). This behaviour is also seen as $\omega _i \to 0$, but only becomes apparent at extremely large radii (e.g. $r= 10^{10}$ gives $|h(r)| = 0.01$ in figure 12b), hence our assertion that eigenmodes with sufficiently small $\omega _i$ can be deemed trapped.

Figure 12. Comparison of asymptotic and numerical solutions with parameters $m = 7$ and $h_{\infty } = 5$, with $R_c = 40$ and $\bar {\xi } = 10$. (a) A damped radiating mode with $F = 0.3$, giving $\omega =-1.38 - 0.013\mathrm {i}$. (b) A trapped mode with $F = 0.5$, giving $\omega = -0.54 - 6\times 10^{-10}\mathrm {i}$. Solid black lines denote the numerical eigenmode, and blue dashed lines denote the far-field asymptotic solution obtained from the deep-water dispersion relation (4.1). The vertical red dashed lines indicate the edge of the numerical damping region.

The growth of eigenmodes in the far field explains why our numerical method breaks down for radiating modes with large temporal decay rate observed at small $F$. When the temporal decay rate $-\omega _i$ is sufficiently large, the spatial growth of the eigenmode overcomes the artificial numerical damping in the damping region, and our numerical NRBC breaks down. This results in $\omega$ becoming dependent on the amount of damping in the damping region, and the numerical eigenvalue is automatically removed as being unresolved as previously described in § 3.3. A similar result has been obtained in Oliveira, Cardoso & Crispino (2014) studying a model wave equation describing the vortex–waves interaction in the shallow-water limit. 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.

4.3. Extension to other azimuthal wavenumbers.

The results described above for $m = 7$ are found to be typical at larger azimuthal wavenumbers. Figure 13 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 14 for the example cases of $m = 15$ and $m = 20$.

Figure 13. 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 14. Shape of the free surface height $h(r, \theta, t=0) = \textrm {Re}[\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$.

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 15 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 $\textrm {Im}(\omega )=10^{-2}$ to $\textrm {Im}(\omega ) = 10^{-5}$, their close spacing shows that the exact threshold value of $\textrm {Im}(\omega )$ separating trapped and radiating behaviour is not that important, with all contours giving a similar boundary. Also plotted in figure 15 is a dashed line showing $\textrm {Re}(\omega )=0$, which is the boundary between counter-rotating and corotating 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 13, 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$.

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

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

4.4. Effect of the free-surface height on the eigenmodes

In this subsection we investigated the effects associated with a variation in the fluid depth, exploring the regime between shallow and deep water. We again focus on modes with azimuthal mode number $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 16 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 17. Finally, in figure 18, the transition between radiating and trapped is shown as the fluid depth is varied between $h_{\infty } = 0.25$, 0.5, 1. Figure 19 shows the variation in eigenvalue as the Froude number is varied, for different fluid depths. The trend in $\textrm {Im}(\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 $\textrm {Re}(\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.

Figure 16. 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.

Figure 17. 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$.

Figure 18. 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$.

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

4.5. Inertial modes

While the focus of this study is surface gravity waves, our numerics finds modal solutions to the governing equations (2.11) 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 characterized 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 20 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 20.

Figure 20. Inertial modes for $m = 2$, $h_{\infty } = 1$ and $F = 0.2$ (a); $F = 0.8$ (b). Eigenvalues curve as a function of the Froude number for the same parameters (c).

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. Throughout this work, we have assumed a harmonic dependence $\exp \{-\mathrm {i}\omega t + \mathrm {i} m\theta \}$ (2.8) and looked for modal solutions of this form. 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 bathtub vortex), here we have shown that they arise in a fully three-dimensional 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 spatiotemporal 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 15). 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 $\textrm {Im}(\omega )$ for large $F$ was not obvious from the numerical analysis carried out here. Interestingly, our numerics finds no unstable wave modes with $\textrm {Im}(\omega )>0$, unlike in either bounded vortex flows (Tophøj et al. 2013; Mougel et al. 2014, 2017) or unbounded shallow water vortex flows (Ford 1994), both of which possess unstable modes with $\textrm {Im}(\omega )>0$. (As shown in Appendix C, our numerics reproduces such instabilities, and thus would be expected to find instabilities in the present case if they existed.) Finally, we have confirmed what has been seen in the initial motivating experiment of the pool, namely that stable trapped surface wave modes exist and can propagate with or against the base swirling flow (as also summarized in figure 15).

Our assumed harmonic dependence $\exp \{-\mathrm {i}\omega t + \mathrm {i} m\theta \}$ (2.8) means modes with $\textrm {Im}(\omega )>0$ are unambiguously unstable and grow exponentially in time, while modes with $\textrm {Im}(\omega )<0$ are unambiguously stable and decay exponentially in time. In the edge case with $\textrm {Im}(\omega )=0$, a critical layer may be present, defined as a radial location $r=r_c$ for which $D_t = -\mathrm {i}\omega + \mathrm {i} mF\varOmega _0(r_c) = 0$. This possibility has been neglected here. The critical layer is related to the continuous spectrum of the differential operator, and analysing the critical layer requires a different analysis and numerical technique: the temporal Fourier transform implicit in the harmonic dependence must be inverted, and the singular point $r=r_c$ must be avoided by complex deformation of the solution into the complex $r$-plane, with the correct choice of deformation to maintain causality (such as has been done for aeroacoustic waves; e.g. King et al. (2022)). Despite the critical layer being present only for $\textrm {Im}(\omega )=0$, it may result in a solution with algebraic growth in time. Consequently, the lack of solutions with $\textrm {Im}(\omega )>0$ does not prove the stability of perturbations to a Lamb–Oseen vortex with a free surface, but only the lack of unstable solutions with exponential growth in time. However, we note that Fabre et al. (2006) found that the Lamb–Oseen vortex without a free surface is stable to axial-invariant perturbations as disturbances related to the critical layer are necessarily damped in that case.

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 hypothesize, 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, whereas in the shallow water case the waves are non-dispersive 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 NRBC 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 NRBC we implemented while remaining computationally viable. Indeed, introduction of additional unknowns in the mathematical formulation was not needed, thus overcoming the main drawback of ‘perfectly matched layer’ 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 NRBC 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. The assumption of modal solutions with a harmonic dependence $\exp \{-\mathrm {i}\omega t + \mathrm {i} m\theta \}$ (2.8) could also be relaxed, for example through time-domain simulation, which would allow investigation into the critical layer (e.g. as in Riedinger et al. (2010) and 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.