2. Set-up

Our starting point is the two-dimensional rotating shallow water equations under the rigid lid assumption. Let $Ox'y$ be Cartesian coordinates, fixed in the topography, where $x'$ describes the distance along a straight coastline and $y$ describes the distance in the offshore direction. We consider a near-semicircular vortex moving along the coastline with speed $U(t)$ and introduce coordinates following the vortex centre (Johnson & Crowe Reference Johnson and Crowe2021) by defining

(2.1)\begin{equation} x = x'-\int_0^t U(t') \,\text{d} t', \end{equation}

and working in the coordinate system $Oxy$. All quantities are specified by their value in the frame of the topography unless stated otherwise. We thus simply work with topographic frame variables expressed as functions of the moving coordinates $(x,y)$ rather than working in variables relative to a frame translating with speed $U(t)$. The advantage of this formulation is that fictitious forces, that would arise from treating quantities relative to the accelerating vortex frame, are absent. Throughout, whenever considering vortices which are semicircular, or asymptotically close to semicircular, we will denote the radius by $a(t)$.

The non-dimensional equations in terms of $(x,y)$ governing the horizontal velocity components $(u,v)$ relative to the topography are thus

(2.2a)\begin{gather} \frac{\partial u}{\partial t}-U\frac{\partial u}{\partial x}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y} - \epsilon v ={-}\frac{\partial p}{\partial x}, \end{gather}

(2.2b)\begin{gather}\frac{\partial v}{\partial t}-U\frac{\partial v}{\partial x}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y} + \epsilon u ={-}\frac{\partial p}{\partial y}, \end{gather}

(2.2c)\begin{gather}\frac{\partial }{\partial x}\left(u H\right) + \frac{\partial }{\partial y}\left(v H \right) = 0, \end{gather}

where $\epsilon = 1/{Ro}$ is the inverse Rossby number and $H$ is the layer depth. We take $H$ to represent a shelf with the depth varying only in the offshore, $y$, direction so $H = H(\kern0.06em y)$. Our set-up is shown in figure 1 for a vortex of radius $a(t)$ and a depth profile, $H(\kern0.06em y)$, consisting of a shelf of width $D$ joined to a region of constant depth.

Figure 1.

Our non-dimensional set-up showing a vortex of radius $a(t)$ moving along a coastal boundary with speed $U(t)$. The layer depth, $H(\kern0.06em y)$, is shown as a shelf region of increasing depth of width $D$ matched to a constant depth ocean. The system is rotating with inverse Rossby number of $\epsilon$.

Equation (2.2) can be combined to give a single evolution equation for the potential vorticity

(2.3)\begin{equation} \left(\frac{\partial }{\partial t}-U\frac{\partial }{\partial x}\right)q + \frac{1}{H}J[\psi,q] = 0, \end{equation}

where the velocity can be expressed using a volume flux stream function

(2.4)\begin{equation} (u,v) = \frac{1}{H}\left(-\frac{\partial \psi}{\partial y},\frac{\partial \psi}{\partial x}\right), \end{equation}

the vorticity is given by

(2.5)\begin{equation} \zeta = \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} = \frac{1}{H}\frac{\partial ^2\psi}{\partial x^2}+\frac{\partial }{\partial y}\left[\frac{1}{H}\frac{\partial \psi}{\partial y}\right], \end{equation}

and $q$ denotes the potential vorticity (PV) in the layer

(2.6)\begin{equation} q = \frac{\zeta+\epsilon}{H}. \end{equation}

We impose a wall at $y = 0$ with the boundary condition $v(x,0) = 0$ or equivalently $\psi (x,0) = 0$. Far from the wall, disturbances are assumed to decay so $(u,v) \to 0$ as $y \to \infty$.

Throughout this study we will consider a shelf profile for $H$ consisting of a shelf region with exponentially increasing depth matched to a constant depth ocean. This profile is chosen for simplicity of calculation, however, Huthnance (Reference Huthnance1974) and Gill & Schumann (Reference Gill and Schumann1974) showed that both the form of the dispersion relation remains unchanged and the inner product exists for general topography, so similar results are expected to hold for arbitrary $H(\kern0.06em y)$. Where possible, general results will also be given in terms of the arbitrary profile $H = H(\kern0.06em y)$. Our shelf profile is given by

(2.7)\begin{equation} H(\kern0.06em y) = \begin{cases} \exp[\beta y], & y \leq D,\\ \exp[\beta D], & y \geq D, \end{cases} \end{equation}

where $D$ describes the shelf width and $\beta$ describes the shelf slope.

We take our vortex solution to be dipolar and centred at $(x,y) = (0,0)$ in our vortex-following coordinates. This corresponds to a single vortex in $y>0$ moving by the image effect. The vortex boundary is taken to be close to semicircular with both the vorticity and stream function being continuous across this boundary. In the far field, the vortex appears as an irrotational source doublet (or equivalently, vortex doublet) of strength $\mu$ directed in the positive $x$ direction. Therefore, far from the vortex we have

(2.8)\begin{equation} \psi = \frac{\mu}{2} \delta(x), \quad \text{at} \ y = 0^+. \end{equation}

In the limit of $\beta \to 0$ (with $\epsilon$ such that $\beta \epsilon \to 0$), our vortex solution is given by the classical Lamb–Chaplygin dipole (Meleshko & van Heijst Reference Meleshko and van Heijst1994) and hence

(2.9)\begin{equation} \psi= \begin{cases} -Uy + \dfrac{2 \textrm{J}_1 (Kr)y}{\textrm{J}_0 (Ka) Kr}, & r< a, \\ -Ua^2y/r^2, & r>a, \end{cases} \end{equation}

for $r^2 = x^2+y^2$. Here $a$ is the semicircular vortex radius, $\textrm {J}_0$ and $\textrm {J}_1$ are Bessel functions of the first kind and $K = j_1/a$ where $j_1 \approx 3.8317$ is the first non-zero root of $\textrm {J}_1$. In this case the dipole strength may be calculated as $\mu = 2{\rm \pi} Ua^2$. A numerical approach for finding these vortex solutions for a general depth profile, $H$, and rotation rate, $\epsilon$, is given in § 5.

3. Vortex decay

This shelf system admits shelf wave solutions, so if the vortex is travelling in the same direction as the phase speed of these waves (requiring $\epsilon U > 0$), it may generate a wave field. These waves will remove energy from the vortex resulting in vortex decay (Flierl & Haines Reference Flierl and Haines1994; Crowe et al. Reference Crowe, Kemp and Johnson2021; Johnson & Crowe Reference Johnson and Crowe2021). Here, we determine conditions for the existence of a wave field and hence determine the parameter values for which a vortex will decay. We then use an asymptotic approach to derive an equation for the decay rate of a vortex under the assumption of small rotation rate and shallow slope.

3.1. The topographic wave field

We begin by determining the linear topographic wave solutions admitted by the system in the absence of a vortex. Our wave equation is obtained by linearising (2.3) and setting $U = 0$ to get

(3.1)\begin{equation} \frac{\partial }{\partial t}\left[ \nabla^2 \psi - \frac{H_y}{H}\frac{\partial \psi}{\partial y}\right] - \frac{\epsilon H_y}{H}\frac{\partial \psi}{\partial x} = 0. \end{equation}

We now take $H$ to be our shelf profile from (2.7) and assume wave-like solutions of the form $\psi = \sqrt {H}\tilde {\phi }(\kern0.06em y)\exp (\textrm {i} \omega t - \textrm {i} k x)$ to obtain

(3.2)\begin{equation} \left[\frac{\partial ^2}{\partial y^2} - k^2 \right] \tilde{\phi} = \begin{cases} -\kappa^2\tilde{\phi} & y \leq D,\\ 0 & y \geq D, \end{cases} \end{equation}

where

(3.3)\begin{equation} \kappa^2(k,\omega) = \frac{\epsilon\beta k}{\omega}-\frac{\beta^2}{4}. \end{equation}

We impose boundary conditions of $\tilde {\phi } = 0$ on $y = 0$ and $\tilde {\phi } \to 0$ as $y\to \infty$. Two boundary conditions are also required at $y = D$ so we take the velocity, $(u,v)$, to be continuous here giving that $\tilde {\phi }$ is continuous across $y = D$ and

(3.4)\begin{equation} \left[ \frac{H_y}{2H} \tilde{\phi} + \tilde{\phi}_y \right]_{D^-}^{D^+} = 0. \end{equation}

For $y > D$ our solution is of the form

(3.5)\begin{equation} \tilde{\phi} = C_1(k) \exp\left(-|k|y\right), \end{equation}

for some $C_1$, hence, using (3.4), we can impose the boundary condition

(3.6)\begin{equation} \tilde{\phi}_y + \left(\frac{\beta}{2}+|k|\right)\tilde{\phi} = 0 \quad \text{on} \quad y = D, \end{equation}

and only consider the shelf region $y \in [0,D]$. In the shelf region we have a solution given by

(3.7)\begin{equation} \tilde{\phi} = C_2(k) \sin \left[\sqrt{\kappa^2-k^2} y\right], \end{equation}

where the square root term may be complex. Using the boundary condition at $y = D$, (3.6), we obtain the dispersion relation

(3.8)\begin{equation} \tan\left[\sqrt{\kappa^2-k^2}D\right] ={-}\frac{\sqrt{\kappa^2-k^2}}{|k|+\beta/2}, \end{equation}

with solutions describing a countably infinite set of modes with differing wavenumber and offshore structure. To proceed we define

(3.9)\begin{equation} l = \sqrt{\kappa^2-k^2}, \end{equation}

so $l$ can be thought of as the offshore wavenumber discretised by the shelf boundary at $y = D$. We can now solve (3.8) numerically for $l(k)$ for each mode and plot the frequency, $\omega$, and phase speed, $c_p = \omega /k$, by combining (3.3) and (3.9) to get

(3.10a,b)\begin{equation} \omega = \frac{\epsilon\beta k}{k^2+l^2+\beta^2/4} \quad \text{and} \quad c_p = \frac{\epsilon\beta}{k^2+l^2+\beta^2/4}. \end{equation}

For $\epsilon > 0$, the phase speed of the waves is positive for all wavenumbers and conversely for $\epsilon < 0$, the topographic wave phase speed is always negative. Additionally, the frequency is odd in $k$, so we only need to consider waves with $k \geq 0$. From (3.8), we note that for a given mode, the offshore wavenumber, $l(k)$, is an increasing function of $k$ and lies within the interval

(3.11)\begin{equation} l \in \left( \left[n-\frac{1}{2}\right]\!\frac{\rm \pi}{D}, \frac{n{\rm \pi}}{D} \right), \end{equation}

for mode number $n = \{1,2,3,\ldots \}$ with

(3.12a,b)\begin{equation} l(0) \quad \text{satisfying}\ \tan{Dl} ={-}\frac{2l}{\beta} \quad\text{and} \quad l(k) \to \frac{n{\rm \pi}}{D} \quad \text{as} \ k\to \infty. \end{equation}

Since $\omega = c_p k$, the group velocity is given by

(3.13)\begin{equation} c_g = \frac{\partial \omega}{\partial k} = c_p + k\frac{\partial c_p}{\partial k}, \end{equation}

where the second term may be shown to be negative for $k \neq 0$ hence $c_g < c_p$ for $k \neq 0$.

Figure 2 shows the frequency and phase speed for the first five modes with $\epsilon = 0.2$, $\beta = 0.1$ and $D = 25.6$. These curves are consistent with classical results for topographic Rossby waves. For a vortex to generate a wave field, the speed of the vortex must match the phase speed of one or more waves, hence a vortex cannot generate any waves if it moves faster than the fastest mode or moves in the opposite direction to the topographic waves. Therefore, for our choice of topography, a vortex moving with speed $U$ will only generate waves if

(3.14)\begin{equation} 0 < \epsilon U < \frac{\epsilon^2 \beta}{l_1^2(0) + \beta^2/4}, \end{equation}

where $l_1(0)$ is the smallest solution to $2l+\beta \tan Dl = 0$ and corresponds to the offshore wavenumber of the lowest mode for $k = 0$. Note that ${\rm \pi} /(2D) < l_1(0) < {\rm \pi}/D$ and in the case of $\beta \to 0$ we have $l_1(0) \to {\rm \pi}/(2D)$. The condition in (3.14) has been multiplied through by $\epsilon$ to ensure it holds for both positive and negative $\epsilon$. We note that $c_g \leq c_p$ for all modes, hence for a vortex moving with speed $U = c_p$, energy will be emitted from the rear of the vortex. This radiation condition will be required later.

Figure 2.

Plots of the frequency, $\omega$, (a) and phase speed, $c_p$, (b) for the first five modes with $\epsilon = 0.2$, $\beta = 0.1$ and $D = 25.6$.

If the vortex does not generate waves we expect that steady vortex solutions will exist, these solutions can be found using the method outlined in § 5.

3.2. Waves generated by a moving vortex

We now determine the amplitude of the vortex generated wave field. Working in coordinates following the vortex ($U\neq 0$) and looking for steady wave solutions, we obtain the linearised wave equation

(3.15)\begin{equation} -U\frac{\partial }{\partial x}\left[\frac{\zeta+\epsilon}{H}\right]-\frac{\epsilon H_y}{H^3}\frac{\partial \psi}{\partial x} = 0. \end{equation}

Substituting for $\psi = \sqrt {H} \phi$ gives

(3.16)\begin{equation} \nabla^2 \phi = \begin{cases} \left(\dfrac{\beta^2}{4} - \dfrac{\epsilon\beta}{U}\right)\phi & y \leq D,\\ 0 & y \geq D, \end{cases} \end{equation}

where we note that

(3.17)\begin{equation} \kappa^2(k,Uk) = \frac{\epsilon \beta}{U}-\frac{\beta^2}{4}. \end{equation}

Here $\kappa ^2$ must be positive if any waves are generated by the vortex by (3.14). As described in (2.8), far from the vortex, the vortex appears as a point dipole of strength $\mu$, hence we impose the boundary condition

(3.18)\begin{equation} \phi = \frac{\mu}{2} \delta(x) \quad \text{at} \ y = 0^+. \end{equation}

We also take $\phi \to 0$ as $y \to \infty$ and impose continuity of $\phi$ and $\phi _y+[H_y/(2H)]\phi$ across $y = D$ so that the velocity $(u,v)$ is continuous here. We note that this approach is equivalent to the matching step between an interior vortex and an exterior wave field in a full asymptotic expansion in small $\epsilon ,\,\beta$ (Flierl & Haines Reference Flierl and Haines1994; Crowe et al. Reference Crowe, Kemp and Johnson2021; Johnson & Crowe Reference Johnson and Crowe2021).

We now express $\phi$ using a Fourier transform as

(3.19)\begin{equation} \phi(x,y) = \frac{1}{2{\rm \pi}}\int_{-\infty}^\infty\hat{\phi}(k,y)\exp(\textrm{i} kx)\,\text{d} k, \end{equation}

where $\hat {\phi }$ satisfies the system

(3.20)\begin{equation} \left[\frac{\partial ^2}{\partial y^2} - k^2\right] \hat{\phi} = \begin{cases} -\kappa^2\hat{\phi} & y \leq D,\\ 0 & y \geq D, \end{cases} \end{equation}

subject to

(3.21)\begin{equation} \left. \begin{array}{ll@{}l} \hat{\phi} = \dfrac{\mu}{2} & \text{at} \ y = 0,\\ \hat{\phi} \to 0 & \text{as} \ y \to \infty,\\ \left[\hat{\phi}\right]^+_-{=} 0 & \text{at} \ y = D, \\ \left[\hat{\phi}_y + \dfrac{H_y}{2H} \hat{\phi} \right]^+_-{=} 0 & \text{at} \ y = D, \end{array}\right\} \end{equation}

where we note that $H_y/H = \beta$ for $y< D$ and $H_y/H = 0$ for $y > D$. The solution for $\hat {\phi }$ is given by

(3.22)\begin{equation} \hat{\phi} = \frac{\mu}{2}\begin{cases} \left(\hat{C} \sin\left[ \sqrt{\kappa^2-k^2} \,y\right] + \cos\left[ \sqrt{\kappa^2-k^2} \,y\right]\right) & y \leq D,\\ \left(\hat{C} \sin\left[ \sqrt{\kappa^2-k^2} D\right] + \cos\left[ \sqrt{\kappa^2-k^2} D\right]\right) \textrm{e}^{-|k|\left(\kern0.06em y-D\right)} & y \geq D, \end{cases} \end{equation}

where

(3.23)\begin{equation} \hat{C} = \frac{\sqrt{\kappa^2-k^2}\sin\left[ \sqrt{\kappa^2-k^2} D\right] - \left(|k|+\dfrac{\beta}{2}\right)\cos\left[ \sqrt{\kappa^2-k^2} D\right]}{\sqrt{\kappa^2-k^2}\cos\left[ \sqrt{\kappa^2-k^2} D\right]+\left(|k|+\dfrac{\beta}{2}\right)\sin\left[ \sqrt{\kappa^2-k^2} D\right]}. \end{equation}

We note that the denominator of $\hat {C}$ vanishes if the dispersion relation in (3.8) is satisfied, therefore the wave modes correspond to the residues of these poles.

By the radiation condition that $c_g - U < 0$, we do not expect any waves generated upstream of the vortex and hence we only consider the solution far downstream where $x<0$ and $|x| \gg 1$. For large $x$, the exponential term in (3.19) is strongly oscillatory and the terms without poles decay as $1/x$. We therefore consider only the terms containing $\hat {C}$ and form a closed contour by including the arc $|k| = R$ with $R\to \infty$ in the lower half of the complex $k$ plane. All poles occur along the real line and are taken to lie within the contour. Finally, the contribution from the arc vanishes as $R\to \infty$ giving that

(3.24)\begin{equation} \phi \sim{-}\frac{\textrm{i}\mu}{2}\sum \text{Res}[\hat{\phi},k_n] \exp(\textrm{i} k_n x), \end{equation}

where $\text {Res}[f,x]$ denotes the residue of $f$ and $x$. Here we sum over all poles of $\hat {\phi }$ and have gained an additional factor of $-1$ due to the orientation of the contour.

Since $\hat {\phi }$ is even in $k$ there will be a pole at $-k = k_n$ for each pole at $k = k_n$ with residue of the opposite sign. We therefore have

(3.25)\begin{equation} \phi \sim \frac{\textrm{i}\mu}{2}\sum_{n=1}^N \text{Res}[\hat{\phi},k_n] \left[ \exp(-\textrm{i} k_n x)-\exp(\textrm{i} k_n x)\right] = \mu\sum_{n=1}^N \text{Res}[\hat{\phi},k_n] \sin(k_n x), \end{equation}

where the $k_n$ are the positive poles of $\hat {C}$ and hence are the solutions to the dispersion relation corresponding to a mode of phase speed $U$. Here $N$ describes the number of modes for which $U = c_p$ and will be determined later. By differentiating the denominator of $\hat {C}$ we find that the residues are given by

(3.26)\begin{equation} \text{Res}[\hat{\phi},k_n] = A_n \begin{cases} \sin\left[ \sqrt{\kappa^2-k_n^2} \,y\right] & y\leq D,\\ \sin\left[ \sqrt{\kappa^2-k_n^2} D\right]\textrm{e}^{{-}k_n\left(\kern0.06em y-D\right)} & y \geq D, \end{cases} \end{equation}

where

(3.27)\begin{equation} A_n = \frac{l_n \left[ \epsilon + U k_n \right]}{k_n D \left[\epsilon+ U k_n\right] + \left[ \epsilon + \dfrac{U}{2}\left(k_n-\dfrac{\beta}{2}\right)\right]}, \end{equation}

for offshore wavenumber $l_n = \sqrt {\kappa ^2-k_n^2}$ satisfying

(3.28)\begin{equation} \tan(l_n D) ={-}\frac{l_n}{k_n+\dfrac{\beta}{2}}. \end{equation}

The total wave field for large, negative $x$ is now given by

(3.29)\begin{equation} \phi \sim \mu \begin{cases} \displaystyle\sum_{n=1}^N \left[A_n \sin(l_n y) \sin(k_n x)\right] & y\leq D,\\ \displaystyle\sum_{n=1}^N \left[A_n \sin(l_n D)\sin(k_n x) \,\textrm{e}^{{-}k_n\left(\kern0.06em y-D\right)}\right] & y \geq D. \end{cases} \end{equation}

The values of $(k_n,l_n)$ can be determined numerically by finding the roots of $c_p = U$ using (3.10a,b) (with $l(k)$ given by (3.8)) or by solving (3.28) directly as a function of $k_n$. It may be shown that the modal components of $\phi$ are mutually orthogonal in the $y$ direction.

Finally, by calculating the maximum value of $c_p$ for each mode and comparing this with $U$, we may determine the total number modes, N, as the greatest integer such that

(3.30)\begin{equation} l_N \leq \kappa = \sqrt{\frac{\epsilon\beta}{U}-\frac{\beta^2}{4}}. \end{equation}

Using the bounds on $l$ from (3.11) we have

(3.31)\begin{equation} \left\lfloor \frac{D}{\rm \pi}\sqrt{\frac{\epsilon\beta}{U}-\frac{\beta^2}{4}} \right\rfloor \leq N \leq \left\lfloor \frac{1}{2}+\frac{D}{\rm \pi}\sqrt{\frac{\epsilon\beta}{U}-\frac{\beta^2}{4}} \right\rfloor, \end{equation}

where $\lfloor * \rfloor$ denotes the ‘floor’ function. By examining the form of (3.28), we observe that if the value of $T = \tan (\kappa D) + 2\kappa /\beta$ is non-negative then equality holds in the upper bound of (3.31) whereas if $T < 0$ then equality holds in the lower bound. The case of $N = 0$ corresponds to the vortex moving faster than the fastest wave and is equivalent to the second inequality in (3.14) not being satisfied.

Figure 3 shows the solutions of $c_p = U$ for $\beta = 0.1$, $D = 25.6$, $U = 1$ and $\epsilon \in \{ 0.2, 0.8\}$. The alongshore wavenumber, $k_n$, for which a given mode has $c_p = U$ is shown by an open circle. The value of $l_n$ can be easily determined using $l_n = \sqrt {\kappa ^2-k_n^2}$. We note that if a vortex slows down, it will generate an increased number of modes, as $\kappa$ will increase.

Figure 3.

Plots of the phase speed, $c_p$, for the first five modes with $\beta = 0.1$, $D = 25.6$ and (a) $\epsilon = 0.2$, (b) $\epsilon = 0.8$. The dotted line denotes $U = 1$ and the $k_n$ are determined as the intersections $c_p = U$ and denoted by the open circles. We observe one mode for $\epsilon = 0.2$ and two modes for $\epsilon = 0.8$ for this choice of parameters.

3.3. Wave energy flux and vortex decay

As the vortex generates waves, it loses energy to the wave field and decays. Since the group velocity for all waves is negative in the frame of the vortex, all energy emitted will cross the line $x = -L$ for $L \gg 1$ where the wave field is small in amplitude and hence linear to leading order. Therefore, this energy flux is given to leading order by the quadratic pressure work plus the transport of energy across the line due to the moving coordinates (Crowe et al. Reference Crowe, Kemp and Johnson2021); so

(3.32)\begin{equation} F = \int_0^\infty H\left[{-}pu + \tfrac{1}{2}\left(u^2+v^2\right)U\right] \,\text{d} y, \end{equation}

where the factor of $H$ is obtained by integrating over the layer depth. For linear waves, the pressure may be determined from (2.2) as

(3.33)\begin{equation} p = Uu+\frac{\epsilon}{H}\psi, \end{equation}

and hence

(3.34)\begin{equation} F =\int_0^\infty \left\{\frac{U}{2H}\left[\left(\frac{\partial \psi}{\partial x}\right)^2-\left(\frac{\partial \psi}{\partial y}\right)^2\right]+\frac{\epsilon\psi}{H}\frac{\partial \psi}{\partial y} \right\} \,\text{d} y. \end{equation}

Equation (3.34) may be written in terms of $\phi = \psi /\sqrt {H}$ as

(3.35)\begin{equation} F ={-}\int_0^\infty \left[ \left(\frac{U H_y^2}{8H^2}-\frac{\epsilon H_y}{2 H}\right) \phi^2 + \left(\frac{U H_y}{2H} - \epsilon \right) \phi \phi_y + \frac{U}{2}\phi_y^2 - \frac{U}{2} \phi_x^2\right] \,\text{d} y, \end{equation}

and calculating $F$ using (3.29) gives

(3.36)\begin{equation} F = \frac{U\mu^2}{4}\sum_{n=1}^N \frac{k_n A_n^2}{\epsilon+U k_n} \left( \left[ \epsilon + \frac{U}{2}\left(k_n-\frac{\beta}{2}\right)\right] + k_n D \left[\epsilon+ U k_n\right]\right). \end{equation}

which we note is independent of $x$. Substituting for $A_n$ using (3.27) we have

(3.37)\begin{equation} F = \frac{U\mu^2}{4}\sum_{n=1}^N \frac{k_n l_n^2 \left[\epsilon+ U k_n\right]}{\left[ \epsilon + \dfrac{U}{2}\left(k_n-\dfrac{\beta}{2}\right)\right] + k_n D \left[\epsilon+ U k_n\right]}, \end{equation}

which we note is always positive, corresponding to a loss of vortex energy. Equating this flux with the loss of vortex energy, $E$, gives

(3.38)\begin{equation} \frac{\textrm{d}E}{\textrm{d}t} ={-}F. \end{equation}

If we now assume that the vortex remains self-similar throughout the evolution, the energy, $E$, and dipole strength, $\mu$, may be determined in terms of the vortex speed $U$ and radius $a$. We now have two quantities, $U$ and $a$, with a single evolution equation, (3.38), so a second equation is required to close the system. Following Flierl & Haines (Reference Flierl and Haines1994), Johnson & Crowe (Reference Johnson and Crowe2021) and Crowe et al. (Reference Crowe, Kemp and Johnson2021), we choose conservation of centre vorticity so that the maximum vorticity of the vortex remains constant throughout the evolution. This gives a second equation

(3.39)\begin{equation} \frac{\textrm{d}\eta_c}{\textrm{d}t} = 0, \end{equation}

where $\eta _c = \eta _c(U,a)$ is the maximum vorticity within the vortex and can be determined from the vortex solution, either numerically or analytically.

In the case of asymptotically small $\beta$ and $\epsilon$, the vortex solution reduces to the classical Lamb–Chaplygin dipole and the quantities may be determined analytically from (2.9) as

(3.40a–c)\begin{equation} E(U,a) = {\rm \pi}U^2 a^2,\quad \mu(U,a) = 2{\rm \pi} U a^2 \quad \text{and} \quad \eta_c(U,a) \propto \frac{U}{a}. \end{equation}

Therefore, (3.38) and (3.39) give

(3.41)\begin{equation} \frac{\textrm{d}}{\textrm{d}t}(U^2 a^2) ={-}{\rm \pi} U^3 a^4 \sum_{n = 1}^N \frac{k_n l_n^2 \left[\epsilon+ U k_n\right]}{\left[ \epsilon + \dfrac{U}{2}\left(k_n-\dfrac{\beta}{2}\right)\right] + k_n D \left[\epsilon+ U k_n\right]} \end{equation}

and

(3.42)\begin{equation} \frac{\textrm{d}}{\textrm{d}t}\left(\frac{U}{a}\right) = 0. \end{equation}

Therefor,e $a(t) \propto U(t)$ so (3.41) may be solved as an equation for $U(t)$ subject to some initial condition

(3.43)\begin{equation} (U,a) = (U_0,a_0)\quad \text{at}\ t = t_0. \end{equation}

We note that since the wavevector, $(k_n,l_n)$, and number of modes, $N$, both have a complicated dependence on $U$, this equation would have to be solved numerically.

Figure 4 shows the wave energy flux, $F$, as a function of $U$ for $(\epsilon ,\beta ) = (0.2,0.1)$ and $(1,0.4)$ with $D = 25.6$ and $\mu = 2{\rm \pi} U$. For large values of $U$ the wave energy flux vanishes as there are no modes which match the vortex speed. As we decrease $U$, an increasing number of modes can satisfy $c_p = U$ so new modes appear in our solution. The dashed lines in figure 4 denote the values of $U$ at which these modes appear (disappear) as $U$ is decreased (increased). Whenever a new mode appears, a peak corresponding to this mode is seen in the energy flux similar to the results of Johnson (Reference Johnson1979). New modes appear with $k_n = 0$ (where $c_p$ is maximal), then $k_n$ increases as $U$ decreases, with the associated energy flux moving through a maximum and dropping off. We now consider the limit of small $U$ where the number of modes, $N$ becomes large.

Figure 4.

Plots of the wave energy flux $F$ (blue) and the limit of $F$ as $N \to \infty$, denoted $F_N$ (red). Here $F$ and $F_N$ are shown as a function of $U$ for $(\epsilon ,\beta ) = (0.2,0.1)$ (a) and $(1,0.4)$ (b) with $D = 25.6$ and $\mu = 2{\rm \pi} U$. If the vortex speed, $U$, exceeds the fastest wave there is no wave energy flux, this occurs for $U > 2.14$ for panel (a) and for $U>7.87$ for panel (b).

3.3.1. The large $N$ limit

In the case of a large number of modes, we must have that

(3.44)\begin{equation} N \approx \frac{D}{\rm \pi} \sqrt{\frac{\epsilon\beta}{U}-\frac{\beta^2}{4}} \gg 1. \end{equation}

Therefore, this limit occurs if the velocity is small ($U \ll 1$) or the shelf width is large ($D \gg 1$). Noting that $N \approx D\kappa /{\rm \pi}$ and $k_n = O(\kappa )$ we have that $k_n D$ is large and hence, from (3.37), we have

(3.45)\begin{equation} F \approx \frac{U\mu^2}{4} \sum_{n=1}^N \frac{l_n^2}{D}. \end{equation}

For large $n$, the offshore wavenumber $l_n$ may be approximated using (3.11) as

(3.46)\begin{equation} l_n \approx \frac{n{\rm \pi}}{D}, \end{equation}

so

(3.47)\begin{equation} F \approx \frac{U\mu^2}{4} \sum_{n=1}^N \frac{n^2{\rm \pi}^2}{D^3} \approx \frac{U\mu^2}{12{\rm \pi}} \left[\frac{\epsilon\beta}{U}-\frac{\beta^2}{4}\right]^{{3}/{2}}, \end{equation}

where we have used

(3.48)\begin{equation} \sum_{n=1}^N n^2 \approx \frac{1}{3} N^3, \end{equation}

for large $N$. We now define $F_N$ to be the asymptotic form of $F$ for large $N$ so

(3.49)\begin{equation} F_N = \frac{U\mu^2}{12{\rm \pi}} \left[\frac{\epsilon\beta}{U}-\frac{\beta^2}{4}\right]^{{3}/{2}}, \end{equation}

where $F_N$ is plotted in figure 4 and can be seen to well describe $F$ for small $U$. Additionally, we observe that $F_N$ provides a fairly good approximation to $F$ for order-one values of $U$.

Taking $\epsilon$ and $\beta$ to be small we may use (3.40a–c) to obtain the approximate evolution equations

(3.50)\begin{equation} \frac{\textrm{d}}{\textrm{d}t}\left(U^2 a^2\right) ={-}\frac{U^3a^4}{3} \left[\frac{\epsilon\beta}{U}-\frac{\beta^2}{4}\right]^{{3}/{2}} \end{equation}

and

(3.51)\begin{equation} \frac{\textrm{d}}{\textrm{d}t}\left(\frac{U}{a}\right) = 0, \end{equation}

which may be easily solved in the case of $4\epsilon /U \gg \beta$ for vortex speed and radius

(3.52)\begin{equation} \left( U, a\right) = \left(U_0, a_0 \right) \left[ 1 + \frac{1}{8}\sqrt{\epsilon^3\beta^3 a_0^4/U_0}\, (t-t_0) \right]^{-{2}/{3}}. \end{equation}

This solution describes a polynomial decay of the vortex speed and radius similar to the case of a beta plane modon considered by Flierl & Haines (Reference Flierl and Haines1994) and Johnson & Crowe (Reference Johnson and Crowe2021). Further, we note that (3.50) and (3.51) exactly correspond to the vortex decay in the continuous limit of an unbounded shelf, $D = \infty$, using the method of Johnson & Crowe (Reference Johnson and Crowe2021) and Crowe et al. (Reference Crowe, Kemp and Johnson2021). While the solution in (3.52) does require $4\epsilon /U \gg \beta$, it can be seen that if this condition is initially satisfied then it will remain true as $U$ decreases.

4. Numerical simulations

To test our predictions we perform numerical simulations using Dedalus (Burns et al. (Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), set-up file available as supplementary material). We solve the full nonlinear, rotating shallow water equations under the rigid lid assumption (see (2.2)) in a frame moving with constant speed, $U_f$, in the along-shore ($x$) direction. We use the numerical domain $(x, y) \in [-51.2, 51.2]\times [0,51.2]$ with $1024$ gridpoints in each direction and decompose fields in terms of a Fourier basis in the $x$ direction and a compound Chebyshev basis in the $y$ direction with separate Chebyshev expansions on and off the shelf. Solutions are integrated for $t \in [0,50]$ using a second-order semi-implicit backward differentiation formula (known as BDF) scheme with a timestep of $10^{-3}$. We take boundary conditions of no flow through the walls at $y = 0$ and $y = 51.2$ and include small viscous terms with a viscosity of $\nu = 1.8\times 10^{-5}$ for numerical stability. The inclusion of viscosity requires additional boundary conditions so we impose free-slip conditions on the walls, $\partial _y u = 0$ on $y = 0,\, 51.2$, and note that the leading-order vortex solution, (2.9), also satisfies these conditions. Therefore, there is unlikely to be significant vorticity generation at the boundaries, something which can lead to a modification or breakdown of the vortex and is particularly prevalent using no-slip boundary conditions.

The use of a Fourier basis in the $x$ direction results in a periodic boundary and hence waves may loop around the domain and interfere with the vortex. However, the stop time, $t = 50$, is found to be sufficiently early that these waves do not interact with the vortex. Similarly, the solid wall at $y = 51.2$ differs from the semi-infinite domain used in our theoretical calculations. Since wave-like disturbances will decay exponentially off the shelf (in the region $y > D$), we would expect any effects of this rigid wall to be exponentially small.

For all simulations, the shelf slope, $\beta$, and shelf width, $D$, are chosen as $(\beta ,D) = (0.1, 25.6)$. Simulations are initialised using the velocity fields corresponding to a Lamb–Chaplygin dipolar vortex (see (2.9)) with initial speed $|U(0)| = 1$ and radius $a(0) = 1$. The frame speed, $U_f$, is set to match the speed of this initial vortex. Therefore, we expect the vortex to remain close to $x = 0$ throughout the evolution with deviations occurring as the vortex speed changes.

The effects of non-zero $\beta$ and $\epsilon$ are to modify the initial vortex leading to a transient adjustment phase at the beginning of the simulation where the vortex adjusts to the effects of rotation and shelf slope, and the wave field begins to develop. For small $\epsilon$ and $\beta$, this adjustment is small and the vortex remains approximately a Lamb–Chaplygin dipole with a modified speed and radius. In order to compare with our theoretical predictions, we take the values of $U_0$ and $a_0$ to be the speed and radius after this adjustment phase, with $t = t_0$ describing the time taken for this adjustment to occur. A value of $t_0 = 2$ is found to be sufficient and comparison is made with the theory for $t \geq t_0$. We note that accurately determining the values of $U_0$ and $a_0$ from the numerical data is difficult. This can present issues when comparing with our theoretical predictions due to the sensitive dependence of (3.41) on these quantities. The vortex speed, $U_0$, is determined by tracking the position of the vorticity maximum and the vortex radius, $a_0$, is estimated using the point at which the vorticity becomes $2\,\%$ of its maximum value.

Figure 5 shows the stream function, $\psi$, for the final timestep, $t = 50$, of our numerical simulations for a range of parameters. Figures 5(a) and 5(b) show vortices which are, respectively, moving faster than and in the opposite direction to all shelf wave modes. For these simulations, the value of $\psi _c$ is conserved to within the error expected due to viscous effects and while a very weak wave signature is observed, this is likely the result of transient waves generated during the initial adjustment. Figure 5(c) shows $\psi (x,y,50)$ for $U(0) = 1$ and $\epsilon = 0.2$ which initially matches the speed of a single wave with predicted wavelength of $\lambda _n = 2{\rm \pi} /k_n = 71.2$. This wavenumber approximately matches the observed wave which we note is likely to be restricted by the length of the domain. Finally figure 5(d) shows $\psi (x,y,50)$ for $U(0) = 1$ and $\epsilon = 1$. The initial speed, $U(0) = 1$, is very close to matching the phase speed of the first three modes, however, the vortex undergoes significant adjustment due to the fairly large value of $\epsilon$ and adjusts to a value of $U_0 \approx 1.1$. This value of $U_0$ only matches the speed of the first two modes and, according to our theoretical predictions, corresponds to modes with wavelengths of $\lambda _n = 22.7$ and $\lambda _n = 31.0$ and offshore wavenumbers of $l_n = 0.11$ and $l_n = 0.22$, respectively. This prediction appears consistent with our simulation where both these modes are seen. We note that as the vortex speeds change throughout the evolution, so too will the wavenumbers of the generated mode. The number of generated modes may also change if the vortex speed slows sufficiently to excite a new mode.

Figure 5.

The stream function, $\psi$, as a function of position in a frame moving with the speed of the initial vortex, $U_f = U(0)$. Results are shown for $\beta = 0.1$, $D = 25.6$ and $t = 50$, for various inverse Rossby numbers, $\epsilon$, and initial vortex speeds, $U(0)$. The vortex adjusts slightly due to finite $\epsilon$ and $\beta$ effects so the value of $U = U_0$ taken at $t = t_0 = 1$ can differ from $U(0)$ by up to 10 %–15 %. (a) Here we have $(\epsilon ,U(0)) = (0.05,1)$, a vortex travelling faster than all shelf waves, (b) $(\epsilon ,U(0)) = (0.6,-1)$, a vortex moving in the opposite direction to all shelf waves, (c) $(\epsilon ,U(0)) = (0.2,1)$, the vortex speed matches a single wave and (d) $(\epsilon ,U(0)) = (1,1)$, the vortex speed matches two waves.

For small $\epsilon$ and $\beta$ we can show using the Lamb–Chaplygin solution that the total stream function at the position of the maximum vorticity is proportional to $Ua$. Therefore, assuming that (3.42) holds, we have

(4.1)\begin{equation} \frac{\psi_c(t)}{\psi_0} = \frac{U(t)\, a(t)}{U_0 a_0} = \frac{U^2(t)}{U_0^2}, \end{equation}

where $\psi _c$ is the value of the stream function at the position of maximum vorticity and $\psi _0$ is the value of $\psi _c$ at $t = t_0 = 2$. The normalised value of $\psi _c$ can be easily determined from our simulations and used to test our decay predictions by comparing with solutions of (3.41) and (3.42) as well as the polynomial decay prediction in (3.52).

For very small values of $\epsilon$ and $\beta$ the wave energy flux, $F$, is small such that the vortex decay, and hence the decrease in $\psi _c$, is slow. Since the effect of viscosity is to decrease the domain averaged energy by around 1 %–2 % over the time interval $t \in [0,50]$, it is not possible to accurately determine how the wave energy flux affects the evolution of the vortex energy when the energy lost to the wave field is similar to the viscous dissipation. Conversely, for values of $\epsilon$ greater that $1$, while the wave field remains small due to small $\beta$, the vortex is no longer well described by the Lamb–Chaplygin solution and the asymptotic expressions for $E$, $\mu$ and $\eta _c$ in (3.40a–c) begin to deviate from the true values. Though these deviations are fairly small despite an order $1$ value of $\epsilon$, we find that (3.38) is very sensitive to the values of $E$ and $\mu$ and our prediction gives only the order of magnitude of the decay scale rather than an accurate result. As a compromise between these limits we consider here the cases of $\epsilon \in \{0.4, 0.6, 0.8\}$, $\beta = 0.1$ and $(U(0),a(0)) = (1,1)$ which we observe are well described by the Lamb–Chaplygin solution. As described above, the values of $U_0$ and $a_0$ are predicted from the vortex speed and radius at $t = t_0$.

Figure 6 shows the stream function, $\psi$, from our numerical simulation with $(\epsilon ,\beta ,D,U_0,a_0) = (0.8,0.1,25.6,1.15,1)$. For these parameters we expect two shelf wave modes with alongshore wavelengths of $\lambda _n = 26.7$ and $\lambda _n = 42.7$. We observe evidence of these two modes and note that the mode with the shorter alongshore wavelength (and hence larger alongshore wavenumber $k_n$ and smaller offshore wavenumber $l_n$), propagates slower in the $x$ direction due to a more negative value of the group velocity, $c_g-U_f$, and can be seen looping around the domain ahead of the slower mode.

Figure 6.

Plots of the stream function, $\psi$, showing the formation of the wave field for $\epsilon = 0.8$, $U_0 = 1.15$, $a_0 = 1$, $\beta = 0.1$ and $D = 25.6$ at two values of $t$, $t = 25$ (a) and $t = 50$ (b). Similarly to figure 5, the solution is shown in a frame moving with speed $U_f = 1$ in the $x$ direction.

Figure 7 shows the values of $\psi _c$ and $\eta _c$ from three simulations in which our theoretical predictions are seen to be accurate. We plot two predictions: firstly, the numerical solution to (3.41) and (3.42) using a fourth-order Runge–Kutta scheme is shown with dashed lines; secondly, the polynomial approximation from (3.52) is shown with dot–dashed lines. In figure 7(a), close agreement is observed between our numerical results and both predictions, and due to the sensitive dependence of these predictions on $U_0$ and $a_0$ it is difficult to determine which is closest to the numerical results. The accuracy of (3.52) is particularly interesting given that while the condition of $4\epsilon /U \gg \beta$ holds for the parameters we consider, the number of modes, $N$, is not large. The numerical value of $\psi _c$ appears to slowly oscillate relative to the theoretical prediction; this is likely a consequence of some higher-order wave-like behaviour within the vortex.

Figure 7.

(a) The normalised value of $\psi _c$ as a function of $t$, the solid line shows our numerical results, the dashed line gives our analytical prediction from solving (3.41) and (3.42) and the dot–dashed lines give our asymptotic prediction from (3.52). (b) The normalised value of $\eta _c$ from our numerical simulation as a function of time. Results in both panels are shown for $(\epsilon ,U_0,a_0) = (0.4,0.98,1)$ (blue), $(\epsilon ,U_0,a_0) = (0.6,1.1,0.9)$ (green) and $(\epsilon ,U_0,a_0) = (0.8,1.15,1)$ (red) with $\beta = 0.1$, $D = 25.6$.

In figure 7(b) we plot the value of the maximum vortex vorticity, $\eta _c$, as a function of time in order to test our assumption (see (3.39)) that this quantity is conserved over the decay scale of the vortex. We observe that $\eta _c$ decreases by around $2\,\%$ over the course of the simulations, and since this is much smaller than the decrease in $\psi _c$ and similar in magnitude to the effects of viscous dissipation, we believe that this assumption is valid.

A supplementary movie, available at https://doi.org/10.1017/jfm.2021.790, shows the time evolution of the wave field depicted in figure 6 and 7 with parameters $(\epsilon ,\beta ,D,U_0,a_0) = (0.8,0.1,25.6,1.15,1)$. We plot the mass fluxes, $(Hu,Hv)$, as functions of time, $t$, and position, $(x,y)$, in a frame moving with constant speed, $U_f = 1$. The formation of two modes with different structure and speed can be clearly seen. Additionally, we observe that the waves looping around the domain due to the periodic $x$ direction are unlikely to have a significant effect on the vortex for $t \lesssim 50$. Further, the use of rigid boundary at $y = 51.2$ rather than the semi-infinite domain considered in our theoretical model is justified as there is no noticeable disturbance near $y = 51.2$.

5. The steady nonlinear problem

We have shown that if the speed of the vortex matches the topographic wave speed for some wavenumber, then the vortex will generate waves and decay due to the transfer of energy from the vortex to the wave field. If, however, there is no wave speed which matches the vortex speed, we expect steady vortex solutions to exist. This can occur if the vortex moves in the opposite direction to the waves or moves faster than the fastest topographic wave. The conditions for a decaying vortex are given in (3.14) and we will focus here on how to determine steady vortex solutions to the full nonlinear system when these conditions are not satisfied.

Neglecting the time derivative in (2.3) and combining the constant advection term with the Jacobian gives

(5.1)\begin{equation} J\left[\psi+U\,\int\, H \,\text{d} y,\frac{\zeta+\epsilon}{H}\right] = 0, \end{equation}

hence we have that the potential vorticity can be written as a function of the total stream function as

(5.2)\begin{equation} \frac{\zeta+\epsilon}{H} = F\left(\psi+U\,\int\, H \,\text{d} y\right). \end{equation}

The function $F$ may now be determined outside the vortex using the far-field condition that $\zeta ,\,\psi \to 0$ as $x \to \infty$ and hence

(5.3)\begin{equation} F\left(U\,\int\, H(\kern0.06em y) \,\text{d} y\right) = \frac{\epsilon}{H(\kern0.06em y)}, \end{equation}

for all $y$. Defining

(5.4)\begin{equation} A(\kern0.06em y) = \int_0^y H(\kern0.06em y') \,\text{d} y', \end{equation}

as the cross-sectional area in the offshore region $[0,y]$ gives that

(5.5)\begin{equation} F(z) = \frac{\epsilon}{H\left(A^{{-}1}\left(z/U\right)\right)}. \end{equation}

The full nonlinear problem outside the vortex can now be written as

(5.6)\begin{equation} \frac{1}{H}\frac{\partial ^2\psi}{\partial x^2}+\frac{\partial }{\partial y}\left[\frac{1}{H}\frac{\partial \psi}{\partial y}\right] + \epsilon = \frac{\epsilon H(\kern0.06em y)}{H\left(A^{{-}1}\left(\psi/U+A(\kern0.06em y)\right)\right)} \end{equation}

and solved subject to

(5.7)\begin{equation} \left.\begin{array}{ll@{}l} \psi = 0 & \text{on} \ y = 0,\\ \psi\to 0 & \text{as} \ x^2+y^2 \to \infty,\\ \psi+UA(\kern0.06em y) = 0 & \text{on} \ \mathcal{C},\\ (u,v)\cdot \hat{\boldsymbol{t}} = u_t & \quad \text{on} \quad \mathcal{C}, \end{array}\right\} \end{equation}

where $\mathcal {C}$ is the vortex boundary and $\hat {\boldsymbol {t}}$ is the tangent vector on $\mathcal {C}$. Here $u_t$ is the tangential velocity inside the vortex which is unknown at this stage. The third boundary condition is the no-normal flow condition which states that the vortex boundary is a streamline of the total stream function $\varPsi = \psi +UA(\kern0.06em y)$. Note that this solution is only valid outside the vortex since there are no streamlines which leave the vortex. Therefore, inside the vortex we must instead impose $F$. Inside the vortex, $\psi$ satisfies

(5.8)\begin{equation} \frac{1}{H}\frac{\partial ^2\psi}{\partial x^2}+\frac{\partial }{\partial y}\left[\frac{1}{H}\frac{\partial \psi}{\partial y}\right] + \epsilon = H F(\psi+UA(\kern0.06em y)) \end{equation}

subject to

(5.9)\begin{equation} \left.\begin{array}{ll@{}l} \psi = 0 & \text{on}\ y = 0,\\ \psi+UA(\kern0.06em y) = 0 & \text{on}\ \mathcal{C},\\ (u,v)\cdot \hat{\boldsymbol{t}} = u_t & \text{on} \ \mathcal{C}. \end{array}\right\} \end{equation}

To obtain a full solution we need to determine the vortex boundary, $\mathcal {C}$, and the function, $F$. This can be achieved by fixing $F$ and then determining $\mathcal {C}$ from the requirement that $\psi +UA = 0$ on $\mathcal {C}$ and $u_t$ is the continuous across $\mathcal {C}$. The typical analytical approach would be instead to impose a boundary and then use $u_t$ from the exterior solution to set a parameter in $F$ (Moffatt Reference Moffatt1969; Meleshko & van Heijst Reference Meleshko and van Heijst1994), however, due to the complicated functional dependence on $H$ and nonlinear nature of this problem, there is unlikely to exist a simple expression for the boundary.

To proceed we take

(5.10)\begin{equation} F(z) ={-}K^2 z + \epsilon, \end{equation}

inside the vortex. In the case of $H = \text {const.}$, the system can be solved exactly to obtain the circular Lamb–Chaplygin dipole (Meleshko & van Heijst Reference Meleshko and van Heijst1994) vortex solution centred at a point on the wall. For the Lamb–Chaplygin dipole we have $K = j_1/a$ where $j_1 \approx 3.83$ is the first non-zero root of the Bessel function, $J_1(z)$, and $a$ is the vortex radius. For the case of arbitrary $\epsilon$ and $H$ we expect $K = K(\epsilon ,\langle H\rangle ,a)$ where $\langle H\rangle$ denotes some list of parameters of $H$ and does not depend on $y$. Therefore, imposing a value for $K$ sets the size of the vortex as a function of $(\epsilon ,\langle H\rangle )$. We note that the vortex boundary, $\mathcal {C}$, will not necessarily be a semicircle as in the Lamb–Chaplygin dipole case, hence $a$ here is a parameter describing the vortex size rather than its radius. We do, however, expect that the vortex will be close to a circle if $\epsilon \lesssim 1$ and $H(\kern0.06em y)$ varies slowly inside the vortex.

We can now seek numerical solutions by choosing a value for $K$ and solving (5.6) and (5.8) to obtain the stream function, $\psi$, and hence the vortex size and velocity fields. The numerical method is described in the Appendix (A).

5.1. Results

To illustrate our results we consider $H$ given by (2.7). The cross-sectional area, $A$, is now given by

(5.11)\begin{equation} A(\kern0.06em y) = \begin{cases} \dfrac{1}{\beta}\left[\exp[\beta y] - 1 \right], & y \leq D,\\ \dfrac{1}{\beta}\left[\exp[\beta D] - 1 \right]+(\kern0.06em y-D)\exp[\beta D], & y \geq D, \end{cases} \end{equation}

which can be easily inverted for $A^{-1}(z)$.

We will consider here the case where the vortex radius is less than the shelf width, $a < D$, however, our theory and numerical method is also valid for $a > D$. Here, the linear operator, $\mathcal {L}$, from (A10) can be written as

(5.12)\begin{equation} \mathcal{L} = \nabla^2 + K^2 \exp[2\beta y]\, \theta(a-r) + \left(\frac{\epsilon \beta}{U}\theta(r-a)-\frac{\beta^2}{4} \right)\theta(D-y), \end{equation}

which reduces to the linear operator for the Lamb–Chaplygin dipole problem

(5.13)\begin{equation} \mathcal{L} = \nabla^2 + K^2 \theta(a-r), \end{equation}

in the case of small $\epsilon ,\,\beta$. The terms $C$ and $N(\phi )$ from (A10) can be similarly expressed in terms of $\beta$ and $D$.

Figures 8 and 9 show our nonlinear solutions for $\psi$ for a range of parameters. Figure 8 shows vortices which travel in the opposite direction to the wave field ($U = -1$) for parameter values $(\epsilon ,\beta ) = \{0.25, 1, 4\}\times \{0.1, 1\}$ while figure 9 shows vortices which travel faster than the fastest topographic mode for $\epsilon = 0.1$ and $\beta \in \{0.1, 1\}$. The solutions are calculated on the numerical domain $(x,y) \in [-51.2,\, 51.2]\times [0,\, 51.2]$ using $(Nx,Ny) = (2048,1024)$ grid points and the solution is assumed to have converged if $\delta = 10^{-10}$ in (A13). We plot the stream function in the frame of the vortex, $\psi +UA(\kern0.06em y)$, so the streamline of height $0$ denotes the vortex boundary. While our system depends on five parameters, $\epsilon$, $\beta$, $D$, $U$ and $K$, we can set $|U| = 1$ and fix $K$. Setting $|U| = 1$ is equivalent to setting the velocity scale used for calculating the inverse Rossby number, $\epsilon$, whereas fixing $K$ determines the size of the vortex and hence we can measure the slope, $\beta$, and shelf width, $D$, in units of the Lamb–Chaplygin dipole radius, $a = j_1/K$. We therefore vary only $\epsilon$, $\beta$ and $D$ and note that the dependence of the vortex structure on $D$ is weak; it is not shown here. We note, however, that $D$ can play an important role in setting the speed of the modes and hence is chosen such that the vortices in figure 9 do not generate topographic modes. For the decaying vortex problem where energy is emitted towards the edge of the shelf we expect that $D$ will play a more important role.

Figure 8.

Plots of the stream function in the vortex frame, $\psi + UA(\kern0.06em y)$, for six pairs of parameters, $(\epsilon ,\beta )$, given by (a) $(0.25,0.1)$, (b) $(0.25,1)$, (c) $(1,0.1)$, (d) $(1,1)$, (e) $(4,0.1)$ and (f) $(4,1)$. Solutions are shown for $K = j_1$, $U = -1$ and $D = 12.5$. Dependence on $D$ is weak and hence not shown here. The dotted line shows the vortex boundary for the Lamb–Chaplygin dipole case of $(\epsilon ,\beta ) = (0,0)$ using the same value of $K$.

Figure 9.

Plots of the stream function in the vortex frame, $\psi + UA(\kern0.06em y)$, for $\epsilon = 0.1$, $U = 1$ and $\beta = 0.1$ (a) and $\beta = 1$ (b). Solutions are shown for $K = j_1$ and $D = 12.5$ and the dotted line shows the vortex boundary for the Lamb–Chaplygin dipole case of $(\epsilon ,\beta ) = (0,0)$ using the same value of $K$. These solution correspond to vortices travelling faster than the fastest topographic mode. While the dependence of the vortex structure on $D$ is weak, the value of $D = 12.5$ does ensure that there are no modes matching the vortex speed.

From figures 8 and 9 we observe that the effect of increasing $\beta$ is to reduce the vortex size and to slightly alter its aspect ratio. Conversely, increasing $\epsilon$ has no significant effect on the vortex size and shape though it does increase the peak value of the stream function, corresponding to an increase in peak vorticity. The vortex shape can be discussed in terms of aspect ratio: the ratio of the offshore radius, $a_y$, to the alongshore radius, $a_x$, given by $a_r = a_y/a_x$. Here $a_x$ and $a_y$ are defined as the distance from the origin to the curve $\psi +UA(\kern0.06em y) = 0$ in the $x$ and $y$ directions, respectively. From (5.1) we note that $\psi$ and $\zeta$ may be scaled on $U$ and hence $\epsilon$ enters the system only through the quantity $\epsilon /U$. Therefore, on dimensional grounds we may write the maximum vorticity as

(5.14)\begin{equation} \zeta_{max} = \frac{U}{a} G\left(\frac{\epsilon}{U},\beta,D\right), \end{equation}

where $G$ is some function describing the dependence of the maximum vorticity on the remaining parameters and $a$ is some parameter describing the vortex size, taken here as the offshore radius, $a = a_y$.

In figure 10 we plot $G$ and $a_r$ as functions of $\epsilon /U$ and $\beta$ for $D = 12.5$. Dependence on $D$ is weak and is not discussed here. Wave generation occurs for positive values of $\epsilon /U$ (see (3.14)) and we find that our iterative method fails to converge for vortices close to the decaying regime. Therefore, only a small region of parameter space is shown for small $\epsilon /U > 0$ corresponding to vortices moving faster than all wave modes as in figure 9. The value of $G$ is seen to increase with both $\beta$ and $-\epsilon /U$ corresponding to a higher vorticity within the vortex centre. Finally, we observe that the aspect ratio, $a_r$, increases slowly with $\beta$ and has very weak dependence on $\epsilon /U$. In the case of $\beta = 0$, the results are independent of $\epsilon$ and correspond to the Lamb–Chaplygin case given in (2.9).

Figure 10.

Plots of $G$ (a) and $a_r$ (b) as functions of $\epsilon /U$ and $\beta$ for $D = 12.5$. The range of parameters shown is limited for $\epsilon /U > 0$ due to the appearance of decaying vortices in this region (we plot results for $\epsilon /U \in [-4, 0.2]$). The dashed line denotes $\epsilon /U = 0$.