2. Set-up

We consider a semi-infinite layer of rotating, linearly stratified fluid of constant buoyancy frequency over a slope with constant gradient. We assume that the slope is sufficiently small for the motion to be governed by the three-dimensional QG equations, and introduce coordinates moving at speed $U$ in the direction ($x$) perpendicular to the down-slope direction ($y$). Finally, we make the SQG approximation that the interior potential vorticity (anomaly) in the layer vanishes. Our system then satisfies

(2.1)\begin{equation} \nabla^2 \psi = 0\quad\textrm{for}\ z>0, \end{equation}

where $\nabla ^2 = \partial _x^2+\partial _y^2+\partial _z^2$, with $(x,y)$ scaled on a horizontal length scale $L$, and $z$ scaled on $fL/N$ for Coriolis parameter $f$ and buoyancy frequency $N$. The evolution equation for the buoyancy on the bottom boundary gives

(2.2)\begin{equation} \left(\frac{\partial }{\partial t}-U\,\frac{\partial }{\partial x}+J\left[\psi,*\right]\right)\frac{\partial \psi}{\partial z} + w = 0\quad\textrm{on}\ z = 0, \end{equation}

which may be combined with the no-normal flow condition

(2.3)\begin{equation} w = \lambda\,\frac{\partial \psi}{\partial x} \quad\textrm{on}\ z = 0, \end{equation}

to obtain a governing equation for $\psi$ on the boundary:

(2.4)\begin{equation} \left(\frac{\partial }{\partial t}-U\,\frac{\partial }{\partial x}+J\left[\psi,*\right]\right)\frac{\partial \psi}{\partial z}+\lambda\,\frac{\partial \psi}{\partial x} = 0\quad\textrm{on} \ z = 0. \end{equation}

Here, $J[g_1,g_2] = \partial _x g_1\, \partial _y g_2 - \partial _x g_2\, \partial _y g_1$ is the Jacobian derivative, and $\lambda = (f^2L/UN)\, s$, for geometric slope $s$, is the slope in the QG limit. Details of this derivation are given in Rodney & Johnson (Reference Rodney and Johnson2014). Motion is assumed to be confined to a region near the bottom boundary, so we impose the far-field condition that $\psi \to 0$ as $z \to \infty$. This system is equivalent mathematically to the case of an SQG layer in a background buoyancy gradient, with $\lambda$ representing the constant horizontal buoyancy gradient (Held et al. Reference Held, Pierrehumbert, Garner and Swanson1995). The results below thus apply to modons moving along a large-scale density gradient in the atmosphere or ocean.

Throughout this work, we will consider two cases, referred to as ‘retrograde’ and ‘prograde’. Here, retrograde describes modons that are moving in the direction opposite to all topographic Rossby waves and are therefore not expected to generate a wave wake. Conversely, prograde describes modons moving with the Rossby waves, which may generate a wake.

3. Steady retrograde modons, $\mu =\lambda a/U\ge 0$

We now seek steady solutions in the form of dipolar vortices (or modons) travelling at speed $U$ in the cross-slope ($x$) direction. We consider vortices that are circular in the plane $z=0$ with (cylindrical) radius $r = a$. Taking $\partial _t = 0$, (2.4) may be written as

(3.1)\begin{equation} J\left[\psi+Uy,\frac{\partial \psi}{\partial z}+\lambda y\right] = 0\quad \textrm{on}\ z = 0, \end{equation}

which gives that

(3.2)\begin{equation} \frac{\partial \psi}{\partial z}+\lambda y = F\left(\psi+Uy\right)\quad \textrm{on}\ z = 0, \end{equation}

for some arbitrary function $F$. We will consider the class of solutions where $F$ is a linear function to obtain the Long's model (Long Reference Long1955) system

(3.3)\begin{equation} \frac{\partial \psi}{\partial z}+\lambda y ={-}\frac{K}{a}\left(\psi+Uy\right)\quad \textrm{on}\ z = 0, \end{equation}

where the value of $K$ may be different inside and outside the vortex. Far from the vortex, we require that both $\psi$ and $\partial _z \psi$ vanish, hence $K/a = -\lambda /U$ outside the vortex ($r > a$). Conversely, inside the vortex ($r < a$), $K$ must be positive to ensure that $\psi$ decays away from the boundary but will otherwise remain arbitrary until determined during the solution. From (3.3), continuity of ${\partial \psi }/{\partial z}$ and the change in the value of $K$ across the vortex boundary require that $\psi +Uy = 0$ on $r = a$. This is equivalent to the vortex boundary being a streamline in coordinates moving with the vortex. Therefore, (3.3) becomes

(3.4)\begin{equation} \frac{\partial \psi}{\partial z} =\begin{cases} -\dfrac{K}{a}\,\psi-\left(\dfrac{UK}{a}+\lambda\right) y & \textrm{for}\ z = 0,\, r < a, \\[8pt] \dfrac{\lambda}{U}\,\psi & \textrm{for}\ z = 0,\ r > a. \end{cases} \end{equation}

Steady modon solutions are now obtained by seeking solutions that satisfy (2.1) and (3.4). Following Muraki & Snyder (Reference Muraki and Snyder2007) and Johnson & Crowe (Reference Johnson and Crowe2023), we write

(3.5)\begin{equation} \psi(x,y,z) = Ua \sin\theta \int_0^\infty \hat{\psi}(\xi) J_1(\xi r/a)\exp(-\xi z/a) \,\xi \,\textrm{d} \xi, \end{equation}

where $\theta$ is the (cylindrical) polar angle, $J_1$ denotes the Bessel function of the first kind of order $1$, and $\hat {\psi }$ describes the Hankel transform of $\psi$ for $\theta ={\rm \pi} /2$ and $z = 0$. Substituting (3.5) into (3.4) gives

(3.6a)\begin{gather} \int_0^\infty \hat{\psi}(\xi) J_1(\xi s) \,\xi^2 \,\textrm{d} \xi - K\int_0^\infty \hat{\psi}(\xi) J_1(\xi s) \,\xi \,\textrm{d} \xi= \left(K+\mu\right) s \quad\textrm{for}\ s < 1, \end{gather}

(3.6b)\begin{gather}\int_0^\infty \hat{\psi}(\xi) J_1(\xi s) \,\xi^2 \,\textrm{d} \xi + \mu\int_0^\infty \hat{\psi}(\xi) J_1(\xi s) \,\xi \,\textrm{d} \xi= 0 \quad\textrm{for}\ s > 1, \end{gather}

where $s = r/a$ is the rescaled vortex radius, and $\mu = \lambda a/U$ is the rescaled slope. We now substitute

(3.7)\begin{equation} A(\xi) = \left(\xi^2+\mu\xi\right)\hat{\psi}(\xi), \end{equation}

so (3.6b) becomes

(3.8)\begin{equation} \int_0^\infty A(\xi) J_1(\xi s) \,\textrm{d} \xi = 0\quad\textrm{for}\ s > 1. \end{equation}

At this point, it is convenient to expand $A(\xi )$ in terms of Bessel functions as

(3.9)\begin{equation} A(\xi) = \sum_{n = 0}^\infty a_n J_{2n+2}(\xi), \end{equation}

for some undetermined coefficients $a_n$ (Tranter Reference Tranter1971; Johnson & Crowe Reference Johnson and Crowe2023). This form allows us to exploit the integral relationship

(3.10)\begin{equation} \int_0^\infty J_{2n+2}(\xi)J_1(\xi s) \,\textrm{d}\xi = \begin{cases} \textit{R}_n(s) & \textrm{for}\ s<1, \\ 0 & \textrm{for}\ s>1, \end{cases} \end{equation}

such that (3.6b) is satisfied automatically for all choices of $a_n$. Here, $\textit {R}_n(s)$ denotes the Zernike radial function (Born & Wolf Reference Born and Wolf2019), a set of degree $2n+1$ polynomials orthogonal over $s\in [0,1]$ with weight $s$. These functions are defined in Appendix B, where their properties and integral relationships discussed. Substituting for $A$ in (3.6a) gives

(3.11)\begin{equation} \int_0^\infty \frac{\xi-K}{\xi+\mu}\,A(\xi) J_1(\xi s) \,\textrm{d} \xi = (K+\mu) s \quad\textrm{for}\ s < 1, \end{equation}

hence

(3.12)\begin{equation} \sum_{n=0}^\infty a_n \int_0^\infty \frac{\xi-K}{\xi+\mu}J_{2n+2}(\xi) J_1(\xi s) \,\textrm{d} \xi = (K+\mu) s \quad\textrm{for}\ s < 1. \end{equation}

We now multiply (3.12) through by $s \textit {R}_m(s)$ and integrate over $s \in [0,1]$. Note that from (3.10), $\textit {R}_n(s)$ is the Hankel transform of $J_{2n+2}(\xi )/\xi$, so inverting this relation gives

(3.13)\begin{equation} J_{2m+2}(\xi)/\xi = \int_0^1 s \textit{R}_m(s) J_1(\xi s) \,\textrm{d} s. \end{equation}

Using (3.13), we obtain

(3.14)\begin{equation} \sum_{n = 0}^\infty a_n \int_0^\infty \frac{\xi-K}{\xi(\xi+\mu)} J_{2n+2}(\xi)J_{2m+2}(\xi) \,\textrm{d} \xi = \frac{K+\mu}{4}\,\delta_{m 0}, \end{equation}

where $m \in \{0,1,2,\dots \}$, and $\delta _{ij}$ is the Kronecker delta. Equation (3.14) takes the form of an infinite, inhomogeneous, generalised eigenvalue problem of the form

(3.15)\begin{equation} \left(\boldsymbol{\mathsf{A}}-K\boldsymbol{\mathsf{B}}\right){\boldsymbol{a}} = {\boldsymbol{c}}. \end{equation}

Here, ${\boldsymbol {a}}$ denotes the $a_n$,

(3.16)\begin{gather} {}[\boldsymbol{\mathsf{A}}]_{mn} = A_{mn} = \int_0^\infty \frac{1}{\xi+\mu}J_{2m+2}(\xi)J_{2n+2}(\xi) \,\textrm{d} \xi, \end{gather}

(3.17)\begin{gather} {}[\boldsymbol{\mathsf{B}}]_{mn} = B_{mn} = \int_0^\infty \frac{1}{\xi(\xi+\mu)} J_{2m+2}(\xi) J_{2n+2}(\xi) \,\textrm{d} \xi \end{gather}

and

(3.18)\begin{equation} [{\boldsymbol{c}}]_m = c_m = \frac{K+\mu}{4}\,\delta_{m0}. \end{equation}

The integrals in (3.16) and (3.17) are oscillatory and hence difficult to calculate directly. A method of rewriting these integrals in a non-oscillatory form that can be calculated easily is discussed in Appendix A.

The inhomogeneity of (3.15) is confined to the row $m = 0$ and may be removed by considering the boundary conditions on $s = 1$. By (3.4), we require

(3.19)\begin{equation} \frac{\partial \psi}{\partial z} = \frac{\lambda}{U}\,\psi \quad\textrm{on}\ s = 1, \end{equation}

therefore

(3.20)\begin{equation} \sum_{n = 0}^\infty a_n \int_0^\infty J_{2n+2}(\xi)J_1(\xi)\,\textrm{d}\xi = 0, \end{equation}

and hence

(3.21)\begin{equation} a_0 = \sum_{n=1}^\infty ({-}1)^{n+1}a_n. \end{equation}

We may therefore remove the first row of (3.15) and replace all instances of $a_0$ using (3.21). The resulting linear system allows us to determine the $a_n$ up to a multiplicative constant, which may be determined subsequently using the first row of (3.15). In the case $\mu = 0$, the integrals $A_{mn}$ and $B_{mn}$ may be found analytically as shown in Johnson & Crowe (Reference Johnson and Crowe2023). Our resulting linear system is

(3.22)\begin{equation} \left[ \boldsymbol{\mathsf{B}} - \frac{1}{K}\, \boldsymbol{\mathsf{A}} \right] {\boldsymbol{a}} = \boldsymbol{0}, \end{equation}

where

(3.23)\begin{equation} \tilde{A}_{mn} = A_{mn}+({-}1)^{n+1} A_{m0} \end{equation}

and

(3.24)\begin{equation} \tilde{B}_{mn} = B_{mn}+({-}1)^{n+1} B_{m0}, \end{equation}

for $m,n \in \{1,2,\dots \}$, and the eigenvectors $a_n$ are scaled using the condition

(3.25)\begin{equation} \left(A_{00}-K B_{00}\right)a_0 + \sum_{n=1}^\infty\left(A_{0n}-K B_{0n}\right)a_n = \frac{K+\mu}{4}, \end{equation}

where $a_0$ is given by (3.21). This linear system may be truncated to a finite number of terms, $m,n \in \{1,2,\dots,N\}$, to obtain numerical results using any eigensolver. We solve for $N$ eigenvalues $K$, and sets of coefficients $a_n$, with the smallest value of $K$ corresponding to the required dipolar vortex. Solutions with higher $K$ are valid and correspond to vortex solutions with a higher-order radial structure. Johnson & Crowe (Reference Johnson and Crowe2023) show that for $\lambda =0$, the higher-order solutions are unstable through pairing of the inner vortices. Higher-order solutions for $\lambda >0$ can be expected to be susceptible to the same instability, so will not be considered further. This method allows modon solutions to be found quickly and accurately, and a value $N = 19$ is sufficient to ensure that the largest neglected coefficient is less than $10^{-6}$. Our solution satisfies

(3.26)\begin{equation} \left[\frac{\partial \psi}{\partial z}-\frac{\lambda}{U}\,\psi\right]_{z = 0} = \begin{cases} -U \sin\theta \sum_{n=0}^\infty a_n \textit{R}_n(r/a) & \textrm{for}\ r< a,\\ 0 & \textrm{for}\ r>a, \end{cases} \end{equation}

and therefore, numerically, once the $a_n$ are determined, it is often easier to solve $\nabla ^2 \psi = 0$ subject to (3.26) than evaluate the Hankel transforms in (3.5). Note also that the Zernike sum in (3.26) can be evaluated from a three-term recurrence relation without requiring the form of the Zernike polynomials (Johnson & Crowe Reference Johnson and Crowe2023). The solutions for $\lambda \ge 0$ are the same as the solutions of Muraki & Snyder (Reference Muraki and Snyder2007); however, the approach here allows solutions to be found using a compact semi-analytical approach rather than the nonlinear root-finding method of Muraki & Snyder (Reference Muraki and Snyder2007).

Figure 1 shows $\psi$ and ${\partial \psi }/{\partial z}$ for a range of non-negative values of $\lambda$. We observe that the maximum values of both $\psi$ and ${\partial \psi }/{\partial z}$ increase with increasing $\lambda$. Further, in the plots of ${\partial \psi }/{\partial z}$, we observe regions just outside $r = a$ of values that are signed oppositely to those within $r < a$. This is due to the relation ${\partial \psi }/{\partial z} = \lambda \psi /U$ for $r > a$, and corresponds to the vortex being surrounded by regions of fluid with buoyancy signed oppositely to the peaks. Evolving these solutions forward in time using the method discussed in § 5 shows that they are both steady and stable over long times.

Figure 1. Plots for (a,c,e) $\psi$ and (b,d,f) $\psi _z$ on $z = 0$ for a range of values of $\lambda$: (a,b) $\lambda = 0$, (c,d) $\lambda = 0.5$, and (e,f) $\lambda = 2$. We use $(U,a) = (1,1)$ throughout. The dashed line denotes the vortex boundary $r = a$.

We note that the maximum values of $\psi /(Ua)$ and $[{\partial \psi }/{\partial z}]/U$ on $z = 0$ depend only on $U$ and $a$ through the parameter $\mu$. This parameter dependence occurs as the integral in (3.5) depends only on $r/a$, $z/a$ and $\mu$, as $\hat {\psi }$ depends only on $\mu$. Figure 2 shows $K$, $\max [\psi ]_{z = 0}$ and $\max [{\partial \psi }/{\partial z}]_{z = 0}$ as functions of non-negative $\mu$. All quantities are observed to increase with $\mu$.

Figure 2. Plots of (a) $K$, (b) $\max [\psi ]_{z = 0}$, and (c) $\max [{\partial \psi }/{\partial z}]_{z = 0}$ as functions of $\mu = \lambda a/U$.

4. Decaying prograde modons, $\mu =\lambda a/U<0$

If a vortex has $\mu < 0$, then the integrals in (3.16) and (3.17) are singular at $\xi = |\mu |$ and do not converge. This corresponds to the case where the moving vortex generates a topographic Rossby wave wake, stating that the phase speed in the negative $x$ direction of a Rossby wave of wavenumber $|\mu |/a = |\lambda /U|$ coincides with the modon speed. Here, solutions to the governing equations are not unique unless a radiation condition is applied (Crowe & Johnson Reference Crowe and Johnson2021; Crowe, Kemp & Johnson Reference Crowe, Kemp and Johnson2021) to ensure that no wave energy is entering the flow from the far field.

To study the case $\mu < 0$, we follow the approach of Flierl & Haines (Reference Flierl and Haines1994) and Johnson & Crowe (Reference Johnson and Crowe2021). We begin by determining the wave field generated by a moving dipolar vortex. The energy flux into this wave field is determined and then equated to the loss of vortex energy to obtain a prediction for the vortex evolution in the case of small $\mu$. We describe this limit as the ‘weak slope’ (or equivalently, strong vortex) limit corresponding to geometric slopes satisfying $s \ll |U|\,N/f^2L$.

4.1. The wave field solution

In the far field, the vortex resembles a point dipole at the origin, so the wave field satisfies

(4.1)\begin{equation} \frac{\partial \psi}{\partial z}-\frac{\lambda}{U}\,\psi ={-}m\, \delta(x)\,\delta'(y) \quad\textrm{on}\ z = 0 \end{equation}

and

(4.2)\begin{equation} \nabla^2 \psi = 0 \quad\textrm{in}\ z > 0, \end{equation}

where $m$ denotes the dipole impulse. We now write $\kappa = -\lambda /U$, or equivalently, $\kappa = -\mu /a$, so $\kappa$ will correspond to the wavenumber of the generated wave field. Equations (4.1) and (4.2) can be solved using a two-dimensional Fourier transform in $x$ and $y$ to obtain

(4.3)\begin{equation} \psi = \frac{m}{4{\rm \pi}^2}\iint \frac{{\rm i}l}{\sqrt{k^2+l^2}-\kappa} \exp\left[ {\rm i}k x + {\rm i}ly - z\sqrt{k^2+l^2}\right] {\rm d}k \,\textrm{d} l. \end{equation}

For $\kappa >0$ (and hence $\mu < 0$), there is a singularity in the integrand of (4.3) at $k^2+l^2 = \kappa ^2$. This singularity corresponds to the appearance of topographic Rossby waves, and the integral must be evaluated such that the solution obeys the radiation condition of no disturbances propagating inwards from the far field. Our dipole solution may be written as the $y$ derivative of the monopole solution of Johnson (Reference Johnson1978a), so

(4.4)\begin{equation} \psi = m\,\frac{\partial G}{\partial y}, \end{equation}

where

(4.5)\begin{align} G(r,\theta,z) &=\frac{1}{2{\rm \pi}\sqrt{r^2+z^2}} \nonumber\\ &\quad - \frac{\kappa}{4}\exp(-\kappa z)\left[ \mathcal{H}_0(\kappa r) + Y_0(\kappa r) + \frac{2}{\rm \pi}\int_0^{\kappa z} \frac{\exp(t)}{\sqrt{t^2+\kappa^2 r^2}} \,\textrm{d} t + \frac{8}{\rm \pi}\,S(\kappa r,\theta)\right], \end{align}

with $\mathcal {H}_0$ the zero-order Struve function, and $Y_0$ a Bessel function of the second kind. The final term, $S$, cancels the wave field ahead of any disturbance and reinforces the wave field behind the disturbance. It also appears in the treatment of barotropic flow past a cylinder on a beta-plane (McCartney Reference McCartney1975), and can be written as

(4.6)\begin{equation} S(\kappa r,\theta) = \sum_{k=0}^\infty \frac{\cos[(2k+1)\theta]\,{\rm J}_{2k+1}(\kappa r)}{2k+1}. \end{equation}

4.2. The wave energy flux

We now reduce to the case of small $\kappa$ and consider a semi-infinite cylinder of radius $R\gg 1$ centred on the vortex. As $\kappa$ is small, we take $R$ such that $\kappa R \ll 1$, and note that the outward energy flux through this cylinder must be independent of $R$. We are therefore able to calculate this energy flux using the streamfunction in the limit of small $\kappa r$. Using polar coordinates,

(4.7)\begin{equation} \psi = m\left[\sin\theta\,\frac{\partial G}{\partial r}+\frac{\cos\theta}{r}\,\frac{\partial G}{\partial \theta}\right], \end{equation}

hence

(4.8)\begin{align} \psi &\sim \frac{m}{\rm \pi}\left( \sin\theta\left[-\frac{r}{2(r^2+z^2)^{3/2}} + O(\kappa)\right] \right.\nonumber\\ &\quad \left.{}+ \sin\theta\cos\theta\left[ \frac{\kappa^4}{3}\,r^2 \exp(-\kappa z) + O(\kappa^6)\right] + O(\kappa^6) \right). \end{align}

The rate of change of vortex energy is given by the sum of the work done by the pressure on the cylinder and the energy advected through the cylinder. This energy loss is calculated in Appendix E and given by (E11). Using our asymptotic result from (4.8), we can calculate the rate of energy loss in the limit of small $\mu$ as

(4.9)\begin{equation} \frac{\mathrm{d} E}{\mathrm{d} t} ={-}\frac{m^2\,|U|\,\kappa^4}{3{\rm \pi}} ={-}\frac{{\rm \pi}\,|a_0|^2\,\lambda^4 a^6}{48\,|U|}, \end{equation}

hence

(4.10)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t}\left[ U^2 a^3\right] ={-}\frac{ |a_0|\,\lambda^4 a^6}{24\,|U|}, \end{equation}

where $a_0$ is defined in (3.21), and the results for the vortex energy $E$ and momentum $m$ are given in the limit of small $\mu$ in Appendix D (D3 and D9) as

(4.11a,b)\begin{equation} m = \frac{{\rm \pi}\,|a_0|\,U a^3}{4}, \quad E = \frac{{\rm \pi}\,|a_0|\,U^2 a^3}{2}. \end{equation}

4.3. Asymptotic predictions for the vortex decay

To predict the vortex decay, we will assume that the vortex form remains self-similar throughout the evolution and hence depends only on its speed, $U$ and radius $a$. Both $U$ and $a$ may vary with time, hence we need two equations involving these parameters in order to determine the evolution. Equation (4.10) corresponds to conservation of energy and is one such equation. Conservation of total momentum or potential vorticity might be expected to hold and give further equations linking $U$ and $a$. However, our previous work (Johnson & Crowe Reference Johnson and Crowe2021) has shown that these quantities are not well conserved throughout the evolution as some fluid escapes the vortex as its radius decreases, carrying significant amounts of momentum and potential vorticity with it. This issue was resolved by following Flierl & Haines (Reference Flierl and Haines1994) and instead considering the value of a tracer within the innermost streamline. Here, the natural choice is buoyancy, and we therefore assume that the buoyancy at the maximum of $\psi$ within the vortex is conserved. This gives

(4.12)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t}\left. \frac{\partial \psi}{\partial z}\right|_{\psi = \psi_{max}} = 0, \end{equation}

and is shown in § 5 to be a reasonable assumption. Noting that in the limit of small $\mu$, ${\partial \psi }/{\partial z} \propto U$ independently of $a$, we have

(4.13)\begin{equation} \frac{\mathrm{d} U}{\mathrm{d} t} = 0 \end{equation}

throughout the vortex evolution. Combining this result with (4.10) gives our predictions for the evolution of the vortex speed and radius as

(4.14a,b)\begin{equation} U(t) = U(t_0), \quad a(t) = a(t_0)\left[ 1+\frac{|a_0|\,a(t_0)^3 \lambda^4}{72\,|U|^3}\,(t-t_0) \right]^{{-}1/3}, \end{equation}

for some initial time $t = t_0$.

5. Numerical solutions for decaying prograde modons

We now test the predictions of § 4 using numerical simulations of (2.1) and (2.4). Solving this system directly requires an inversion of the Laplacian at each time step in order to calculate ${\partial \psi }/{\partial z}$. The value of ${\partial \psi }/{\partial z}$ on the boundary is then stepped forward using (2.4). Since $\psi$ satisfies the Laplace equation, it is standard practice to determine ${\partial \psi }/{\partial z}$ from $\psi$ using a Dirichlet-to-Neumann operator (Held et al. Reference Held, Pierrehumbert, Garner and Swanson1995). This avoids solving in the full three-dimensional domain and reduces our problem to two spatial dimensions. Working in Fourier space in $x$ and $y$, (2.1) becomes

(5.1)\begin{equation} \frac{\partial ^2\hat\psi}{\partial z^2} - (k^2+l^2)\, \hat\psi = 0, \end{equation}

hence

(5.2)\begin{equation} \left.\frac{\partial \hat\psi}{\partial z}\right|_{z = 0} ={-}\sqrt{k^2+l^2}\left.\hat\psi\right|_{z = 0}, \end{equation}

where the choice of sign comes from the requirement that $\psi \to 0$ as $z\to \infty$. We now define $\mathcal {D}$ to be a linear operator that acts in Fourier space as $\hat {\mathcal {D}} = -\sqrt {k^2+l^2}$. Therefore, $\mathcal {D}$ may be thought of as both the fractional Laplacian operator $\mathcal {D} = -\sqrt {-\nabla ^2}$ and the required Dirichlet-to-Neumann map ${\partial \psi }/{\partial z} = \mathcal {D}\,\psi$. We therefore solve the two-dimensional system

(5.3)\begin{equation} \left(\frac{\partial }{\partial t}-C\,\frac{\partial }{\partial x}+J\left[\psi,*\right]\right)\mathcal{D}\psi+\lambda\,\frac{\partial \psi}{\partial x} ={-}\nu\,\nabla^4 \mathcal{D}\psi \quad\textrm{on}\ z = 0, \end{equation}

where hyperdiffusion is included for numerical stability. Here, $C$ denotes the speed of the moving coordinates, hence taking $C = U$ results in the vortex remaining centred on the origin.

Equation (5.3) is solved using spectral methods via the Dedalus package (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). We use a doubly periodic two-dimensional domain of size $102.4 \times 102.4$, with $4096$ grid points in each direction, and expand $\psi$ in a Fourier basis in both directions. Time stepping is performed using a $3-\epsilon$ order implicit–explicit Runge–Kutta scheme. Our simulations are initialised by placing an SQG modon at the centre of the domain using the solutions from § 3 with $\lambda = 0$ and $(U,a) = (-1,1)$. Simulations are run for $\lambda \in \{0,0.2,0.4,0.6,0.8,1.0\}$ over the time interval $t \in [0,50]$. This time interval is chosen as it ensures that the generated waves do not have sufficient time to loop around the domain and interact with the vortex.

Initially, over a short time interval $t \in [0,t_0]$, the vortex solution adjusts to the non-zero slope parameter $\lambda$ and begins to generate a wave field. To compare our simulations with the theoretical predictions of (4.14a,b), we take a value $t_0 = 5$ and calculate $a(t_0)$ and $U(t_0)$ from our numerical data. This value of $t_0$ is found to be sufficient such that the transient motions arising from the initial adjustment have decayed.

Figure 3 shows $\psi$ from our numerical simulations at times $t = t_0 = 5$ and $t = 40$ for slope $\lambda = 0.4$. The formation of a wave wake can be seen clearly. Results are shown in coordinates moving with speed $C = -1$. Since the vortex speed $U$ changes slightly during the initial adjustment phase, the vortex does not remain at the origin throughout the evolution. Figure 4 shows the maximum values of $\psi$ and ${\partial \psi }/{\partial z}$ as functions of time for each simulation run. Our asymptotic decay predictions from (4.14a,b) are shown as dashed lines and found to give accurate predictions for the decay, even for the cases where $\lambda$ is order 1. The value of $\max [{\partial \psi }/{\partial z}]$ is found to be well conserved throughout the evolution – varying by less than $3\,\%$ for most cases – justifying the assumptions made in § 4.3. The curves for $\lambda \ge 0.6$ in figure 4 are shown over the limited interval $t \in [t_0, t_\lambda ]$, where $t_\lambda = 40, 25, 20$, respectively, for $\lambda = 0.6, 0.8, 1.0$. Outside these intervals, the dipolar vortices break down and the assumption of self-similar vortex structure is no longer valid. This breakdown is explored in § 6.

Figure 3. The streamfunction $\psi$ as a function of $(x,y)$ for $\lambda = 0.4$ at times (a) $t = 5$ and (b) $t = 40$.

Figure 4. Plots of (a) $\max [\psi ]$ and (b) $\max [{\partial \psi }/{\partial z}]$ from numerical simulations for $\lambda \in \{0, 0.2, 0.4, 0.6, 0.8, 1\}$. Results are shown as functions of time $t$ from $t_0 = 5$ onwards, and all curves are normalised by their initial value. Dashed lines denote the asymptotic predictions of (4.14a,b).

6. Vortex breakdown

Numerical solutions reveal that over long times, decaying dipolar vortices with $\mu < 0$ will break down into a pair of monopolar vortices. When this breakdown occurs, the structures cease to move purely in the east–west direction and instead follow complex trajectories. This breakdown occurs sooner for steeper slopes, i.e. sooner for larger values of $\lambda$ and $|\mu |$.

Examination of simulation results suggests that the breakdown occurs through the development of waves inside the vortex. These waves lead to an asymmetry that causes the two halves of the dipole to split. The wave amplitude is larger for larger values of $|\mu |$, hence the breakdown occurs sooner. Figure 5 shows the breakdown of a dipolar vortex for $(U,a,\lambda ) = (-1,1,0.8)$ that occurs from approximately $t_\lambda = 25$. The breakdown is seen most clearly in the buoyancy field ${\partial \psi }/{\partial z}$, hence we plot buoyancy rather than streamfunction $\psi$. We observe that after an asymmetry develops in the dipolar structure, streamers of buoyancy (and potential vorticity) are stripped away from each side. These streamers later roll up into lines of small monopolar vortices that interact with each other through a series of merging and splitting events, as predicted by the theory of SQG turbulence (Held et al. Reference Held, Pierrehumbert, Garner and Swanson1995). As the initial asymmetry results in the two halves having different vortex strengths, they cease to move together in a straight line. Instead, the modified vortex pair deflects south before following a complicated path as the vortex strength changes due to subsequent interactions. As expected, waves are still generated post-breakdown if there is a significant disturbance moving westwards. Supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.607 shows the full evolution of the modon for the case $\lambda = 0.8$ depicted in figure 5. A similar breakdown due to asymmetric instability is observed by Makarov & Kizner (Reference Makarov and Kizner2011) and Jalali & Dritschel (Reference Jalali and Dritschel2020) for the case of oppositely signed vortex patches.

Figure 5. Plots of the surface buoyancy ${\partial \psi }/{\partial z}$ for $\lambda = 0.8$ and (a) $t = 25$, (b) $t = 27$, (c) $t = 29$, (d) $t = 30$, (e) $t = 35$, and (f) $t = 40$. The breakdown of a dipolar vortex into two monopolar vortices can be seen clearly.

Snyder et al. (Reference Snyder, Muraki, Plougonven and Zhang2007) considered the evolution of SQG modons in a primitive equation model and observed that waves could form within the modon. These inertia-gravity waves were found to result from both initial adjustment and imbalance, and were stronger for larger Rossby numbers. In our case, we are considering a balanced model and, as such, expect that waves occur through a combination of initial adjustment and topographic effects. A small disturbance advected in the up-slope $y$ direction by the interior flow can increase in vorticity due to the conservation of potential vorticity. This magnification of small disturbances may in turn lead to vortex breakdown.

7. Stabilisation by a coastal boundary

By the image effect, we would expect a monopolar vortex moving along a coastal boundary to behave in the same way as a dipolar modon and be described by the same solutions as discussed in §§ 2–4. In particular, we expect steady solutions for retrograde vortices, and a decay rate for prograde vortices that follows (4.14a,b) in the limit of small $\lambda$.

To test our decay prediction, we rerun our Dedalus simulations on the half-domain $(x,y) \in [-51.2,51.2]\times [-51.2,0]$ using a sine basis in the $y$ direction to enforce the condition that $\psi = 0$ along the wall $y = 0$. All other parameters and analysis techniques remain unchanged. Figure 6 shows the maximum values of $\psi$ and ${\partial \psi }/{\partial z}$ as functions of time for each simulation run. Our asymptotic prediction from (4.14a,b) is included as a dashed line and found to be consistent with the observed decay, even for order $1$ values of $\lambda$. Figure 7 shows the surface streamfunction and buoyancy for the case $\lambda = 0.8$ at two different times. The formation of the wave field can be seen clearly in the plots of $\psi$, while the close-up plots of ${\partial \psi }/{\partial z}$ show the decrease in vortex size at later times. We observe that at late times and larger values of $\lambda$, the vortex shape becomes elliptical, with a shorter semi-radius in the $y$ direction. This will reduce the validity of the assumption that the vortex remains self-similar throughout the evolution, and may explain partially why the decay prediction is less accurate in these cases. Supplementary movie 2 shows the full evolution of the coastal vortex for the case $\lambda = 0.8$ depicted in figure 7.

Figure 6. As figure 4 for a monopolar vortex moving along a coastal boundary.

Figure 7. The time evolution of a monopolar vortex on a wall for $\lambda = 0.8$. We plot the surface streamfunction $\psi$ for (a) $t = 10$ and (b) $t = 40$, and the surface buoyancy ${\partial \psi }/{\partial z}$ for (c) $t = 10$ and (d) $t = 40$.

A major difference observed between the dipolar vortices and the coastal monopoles is that the vortex breakdown discussed in § 6 does not occur during the monopole evolution. This is likely a consequence of the enforced symmetry of the image vortex across the line $y = 0$, and suggests that the breakdown of the dipoles occurs due to the growth of an asymmetry between the regions of positive and negative vorticity.