Modon solutions in an N-layer quasi-geostrophic model

Abstract

Modons, or dipolar vortices, are common and long-lived features of the upper ocean, consisting of a pair of monopolar vortices moving through self-advection. Such structures remain stable over long times and may be important for fluid transport over large distances. Here we present a semi-analytical method for finding fully nonlinear modon solutions in a multi-layer quasi-geostrophic model with arbitrarily many layers. Our approach works by reducing the problem to a multi-parameter linear eigenvalue problem which can be solved using numerical techniques from linear algebra. The method is shown to replicate previous results for one and two-layer models and is applied to a three-layer model to find a solution describing a mid-depth propagating, topographic vortex.

1 Introduction

Modons, or dipolar vortices, are common coherent structures in the upper ocean (Ni et al., 2020) and are observed to exist for long times and transport fluid over large distances (Nycander & Isichenko, 1990). Modons are similarly relevant to atmospheric dynamics, and have been used to model various atmospheric processes such as atmospheric blocking (McWilliams, 1980) and Madden-Julian Oscillation events (Rostami & Zeitlin, 2021).

Early studies of modon solutions to the two-dimensional Euler equation were carried out independently by Lamb and Chaplygin (Lamb, 1932; Meleshko & van Heijst, 1994) and similar analytical solutions to the equivalent barotropic, quasi-geostrophic problem were identified by Larichev & Reznik (1976). The long-term behaviour of these structures remains an active area of study, with much recent work devoted to understanding modon breakdown processes such as instabilities (Makarov & Kizner, 2011; Flierl & Haines, 1994; Brion et al., 2014; Johnson & Crowe, 2021; Davies et al., 2023; Protas, 2024), energy loss through wave generation (Flierl & Haines, 1994; Johnson & Crowe, 2021), and viscous decay (Flór et al., 1995; Nielsen & Rasmussen, 1997).

Layered quasi-geostrophic models are commonly used in modelling the ocean and atmosphere. These models consist of two-dimensional layers coupled through a vertical pressure gradient and are generally valid on scales where the Earth’s rotation is dominant. Determining modon solutions in such models allows for the study of coherent structures in a more realistic setting than simple 2D dynamics. Kizner et al. (2003) thoroughly explore modon solutions in a two-layer model however calculating modon solutions in models with arbitrarily many layers would be mathematically tedious using the same analytical methods.

Here, we extend a method we developed for finding surface QG dipoles (Johnson & Crowe, 2023; Crowe & Johnson, 2023) to the layered QG system by analytically reducing the problem to a linear eigenvalue problem which can be solved numerically. We begin by presenting the model and solution method in Sections 2 and 3, before applying it to some examples in Section 4. Our method is found to be consistent with previous one-layer and two-layer studies and effective at finding new solutions in a three-layer model. Finally, we discuss our method and results in Section 5.

2 Problem Setup

Consider a dipolar vortex (or modon) of radius $a$ moving in the zonal ($x$) direction with speed $U$. In the frame of the moving vortex, the $N$-layer quasi-geostrophic equations are

$$(\partial_{t}-U\partial_{x})q_{i}+J(\psi_{i},q_{i})+\beta_{i}\partial_{x}\psi_{i}=0,\quad i\in\{1,2,\dots,N\},$$ (1)

for streamfunction $\psi_{i}$ and potential vorticity anomaly $q_{i}$ in each layer. Here $\beta_{i}$ denotes the background vorticity gradient in the meridional ($y$) direction and may vary by layer to incorporate the effects of bottom topography. The potential vorticities in each layer are given by

$$q_{1}= \,\nabla^{2}\psi_{1}+R_{1}^{-2}(\psi_{2}-\psi_{1}),$$ (2)

$$q_{i}= \,\nabla^{2}\psi_{i}+R_{i}^{-2}(\psi_{i-1}-2\,\psi_{i}+\psi_{i+1}),\quad i\in\{2,\dots N-1\},$$ (3)

$$q_{N}= \,\nabla^{2}\psi_{N}+R_{N}^{-2}(\psi_{N-1}-\psi_{N}),$$ (4)

where

$$R_{i}=\frac{\sqrt{g^{\prime}H_{i}}}{f},$$ (5)

is the Rossby radius of deformation in each layer. Here $g^{\prime}$ denotes the buoyancy difference between consecutive layers and is taken to be the same for all layers for simplicity. Layer depth is denoted by $H_{i}$ and $f$ denotes the Coriolis parameter. The effects of a background flow can be easily incorporated into this model but will not be considered here.

Eq. 1 can be written as

$$\partial_{t}q_{i}+J(\psi_{i}+Uy,q_{i}+\beta_{i}y)=0,\quad i\in\{1,2,\dots,N\}$$ (6)

and we now seek steady, nonlinear modon solutions by letting $\partial_{t}=0$. Eq. 6 may be written as

$$q_{i}+\beta_{i}y=F_{i}(\psi_{i}+Uy),\quad i\in\{1,2,\dots,N\}$$ (7)

where $F_{i}$ is an arbitrary, piece-wise differentiable function of a single variable. To proceed, we take $F_{i}$ to be a piece-wise linear function

$$F_{i}(z)=\begin{cases}(\mathcal{K}_{i}/a^{2})\,z,&x^{2}+y^{2}<a^{2},\\ (\beta_{i}/U)\,z,&x^{2}+y^{2}>a^{2},\end{cases}$$ (8)

where $a$ denotes the radius of the vortex. The form of the $F_{i}$ outside the vortex is set by the requirement that vorticity and streamfunction perturbations decay towards infinity. Inside the vortex, we may either choose $\mathcal{K}_{i}$ to match the exterior solution, i.e. $\mathcal{K}_{i}=\beta_{i}a^{2}/U$, or take $\mathcal{K}_{i}=-K_{i}^{2}$ where $K_{i}$ is an eigenvalue which will be determined by continuity of $\psi_{i}$ across the vortex boundary. These cases will be respectively referred to as ‘passive’ and ‘active’ vortex regions. An active region consists of an isolated region of vorticity whereas the vorticity in a passive region is determined only by the advection of the background vorticity, $\beta_{i}y$, by the vortex flow-field.

3 Solution Method

Using the form of $F_{i}$ from Eq. 8, Eq. 7 can be written as

$$(\nabla^{2}-R_{1}^{-2})\psi_{1}+R_{1}^{-2}\psi_{2}+\beta_{1}y =(\mathcal{K}_{1}/a^{2})(\psi_{1}+Uy),$$ (9)

$$(\nabla^{2}-2R_{i}^{-2})\psi_{i}+R_{i}^{-2}(\psi_{i-1}+\psi_{i+1})+\beta_{i}y =(\mathcal{K}_{i}/a^{2})(\psi_{i}+Uy),$$ (10)

$$(\nabla^{2}-R_{N}^{-2})\psi_{N}+R_{N}^{-2}\psi_{N-1}+\beta_{N}y =(\mathcal{K}_{N}/a^{2})(\psi_{N}+Uy),$$ (11)

inside the vortex ($x^{2}+y^{2}<a^{2}$) and

$$(\nabla^{2}-R_{1}^{-2})\psi_{1}+R_{1}^{-2}\psi_{2} =(\beta_{1}/U)\psi_{1},$$ (12)

$$(\nabla^{2}-2R_{i}^{-2})\psi_{i}+R_{i}^{-2}(\psi_{i-1}+\psi_{i+1}) =(\beta_{i}/U)\psi_{i},$$ (13)

$$(\nabla^{2}-R_{N}^{-2})\psi_{N}+R_{N}^{-2}\psi_{N-1} =(\beta_{N}/U)\psi_{N},$$ (14)

outside the vortex ($x^{2}+y^{2}>a^{2}$) for $i\in\{2,\dots N-1\}$.

3.1 Solution using Hankel transforms

To proceed, we move to plane polar coordinates, $(x,y)=(r\cos\theta,r\sin\theta)$, and represent $\psi_{i}$ using the Hankel tranform

$$\psi_{i}=Ua\,\sin\theta\int_{0}^{\infty}\hat{\psi}_{i}(\xi)\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi,$$ (15)

where $s=r/a$ and $\operatorname{J}_{1}$ denotes the Bessel function of the first kind of order $1$. The angular dependence in Eq. 15 is chosen to match that of a typical dipolar vortex. Substituting Eq. 15 into Eqs. (9) to (14) gives

$$\int_{0}^{\infty}\left[(\xi^{2}+\lambda_{1}^{2}+\mathcal{K}_{1})\hat{\psi}_{1}-\lambda_{1}^{2}\hat{\psi}_{2}\right]\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi =(\mu_{1}-\mathcal{K}_{1})s,$$ (16)

$$\int_{0}^{\infty}\left[(\xi^{2}+2\lambda_{i}^{2}+\mathcal{K}_{i})\hat{\psi}_{i}-\lambda_{i}^{2}(\hat{\psi}_{i+1}+\hat{\psi}_{i-1})\right]\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi =(\mu_{i}-\mathcal{K}_{i})s,$$ (17)

$$\int_{0}^{\infty}\left[(\xi^{2}+\lambda_{N}^{2}+\mathcal{K}_{N})\hat{\psi}_{N}-\lambda_{N}^{2}\hat{\psi}_{N-1}\right]\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi =(\mu_{N}-\mathcal{K}_{N})s,$$ (18)

inside the vortex ($s<1$) and

$$\int_{0}^{\infty}\left[(\xi^{2}+\lambda_{1}^{2}+\mu_{1})\hat{\psi}_{1}-\lambda_{1}^{2}\hat{\psi}_{2}\right]\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi =0,$$ (19)

$$\int_{0}^{\infty}\left[(\xi^{2}+2\lambda_{i}^{2}+\mu_{i})\hat{\psi}_{i}-\lambda_{i}^{2}(\hat{\psi}_{i+1}+\hat{\psi}_{i-1})\right]\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi =0,$$ (20)

$$\int_{0}^{\infty}\left[(\xi^{2}+\lambda_{N}^{2}+\mu_{N})\hat{\psi}_{N}-\lambda_{N}^{2}\hat{\psi}_{N-1}\right]\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi =0,$$ (21)

outside the vortex ($s>1$) for $i\in\{2,\dots N-1\}$. Here $\mu_{i}=\beta_{i}a^{2}/U$ and $\lambda_{i}=a/R_{i}$ respectively denote the nondimensional background vorticity gradient and the ratio of the vortex radius to the Rossby radius in each layer.

We now define the matrix function

$$\mathsfbi{K}(\xi)=\begin{pmatrix}\xi^{2}+\lambda_{1}^{2}&-\lambda_{1}^{2}&0&\dots&0&0\\ -\lambda_{2}^{2}&\xi^{2}+2\lambda_{2}^{2}&-\lambda_{2}^{2}&\dots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\dots&-\lambda_{N}^{2}&\xi^{2}+\lambda_{N}^{2}\end{pmatrix},$$ (22)

and diagonal matrix

$$\mathsfbi{D}(\textbf{c})=\begin{pmatrix}c_{1}&0&\dots&0\\ 0&c_{2}&\dots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&c_{N}\end{pmatrix},$$ (23)

for some vector $\textbf{c}=(c_{1},c_{2},\dots c_{N})^{T}$. We also introduce the vectors $\boldsymbol{\mu}=(\mu_{1},\mu_{2},\dots,\mu_{N})^{T}$, $\boldsymbol{\mathcal{K}}=(\mathcal{K}_{1},\mathcal{K}_{2},\dots,\mathcal{K}_{N})^{T}$ and $\hat{\boldsymbol{\psi}}=(\hat{\psi}_{1},\hat{\psi}_{2},\dots,\hat{\psi}_{N})^{T}$ so Eqs. (16) to (21) can be written concisely as

$$\int_{0}^{\infty}\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mathcal{K}})\right]\hat{\boldsymbol{\psi}}(\xi)\,\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi =(\boldsymbol{\mu}-\boldsymbol{\mathcal{K}})s, s<1,$$ (24)

$$\int_{0}^{\infty}\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mu})\right]\hat{\boldsymbol{\psi}}(\xi)\,\operatorname{J}_{1}(s\xi)\,\xi\,\textrm{d}\xi =0, s>1.$$ (25)

To proceed, we follow the approach of Johnson & Crowe (2023) and let

$$\textbf{A}(\xi)=\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mu})\right]\hat{\boldsymbol{\psi}}(\xi)\,\xi,$$ (26)

where A is expanded in terms of Bessel functions as

$$\textbf{A}(\xi)=\sum_{j=0}^{\infty}\textbf{a}_{j}\operatorname{J}_{2j+2}(\xi),$$ (27)

where the $\textbf{a}_{j}$ are vectors of expansion coefficients. This expansion allows us to exploit the integral relation

$$\int_{0}^{\infty}\operatorname{J}_{2j+2}(\xi)\operatorname{J}_{1}(s\xi)\,\textrm{d}\xi=\begin{cases}\operatorname{R}_{j}(s)&\textrm{for}\quad s<1,\\ 0&\textrm{for}\quad s>1,\end{cases}$$ (28)

and its inverse result

$$\operatorname{J}_{2k+2}(\xi)/\xi=\int_{0}^{1}s\operatorname{R}_{k}(s)\operatorname{J}_{1}(s\xi)\,\textrm{d}s,$$ (29)

where

$$\operatorname{R}_{n}(s)=(-1)^{n}s\operatorname{P}^{(0,1)}_{n}(2s^{2}-1),$$ (30)

denotes the Zernike radial function, a set of polynomials of degree $2n+1$ (Born & Wolf, 2019) which are orthogonal over $s\in[0,1]$ with weight $s$, and $\operatorname{P}_{n}^{(\alpha,\beta)}$ denotes the Jacobi polynomial. From Eq. 28 for $s>1$, we observe that Eq. 25 is automatically satisfied by our expansion of $\textbf{A}(\xi)$ in Eq. 27 for all choices of $\textbf{a}_{j}$.

Eq. 24 may now be written as

$$\sum_{j=0}^{\infty}\left[\int_{0}^{\infty}\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mathcal{K}})\right]\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mu})\right]^{-1}\operatorname{J}_{2j+2}(\xi)\operatorname{J}_{1}(s\xi)\,\textrm{d}\xi\right]\textbf{a}_{j}=(\boldsymbol{\mu}-\boldsymbol{\mathcal{K}})s,$$ (31)

for $s<1$. We now multiply Eq. 31 by $s\operatorname{R}_{k}(s)$ and integrate over $s\in[0,1]$ to get

$$\sum_{j=0}^{\infty}\left[\int_{0}^{\infty}\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mathcal{K}})\right]\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mu})\right]^{-1}\xi^{-1}\operatorname{J}_{2j+2}(\xi)\operatorname{J}_{2k+2}(\xi)\,\textrm{d}\xi\right]\textbf{a}_{j}=\frac{\delta_{k0}}{4}(\boldsymbol{\mu}-\boldsymbol{\mathcal{K}}),$$ (32)

where $k\in\{0,1,2,\dots\}$. Defining

$$\mathsfbi{A}_{kj} =\int_{0}^{\infty}\mathsfbi{K}(\xi)\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mu})\right]^{-1}\xi^{-1}\operatorname{J}_{2j+2}(\xi)\operatorname{J}_{2k+2}(\xi)\,\textrm{d}\xi,$$ (33)

$$\mathsfbi{B}_{kj} =\int_{0}^{\infty}\left[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mu})\right]^{-1}\xi^{-1}\operatorname{J}_{2j+2}(\xi)\operatorname{J}_{2k+2}(\xi)\,\textrm{d}\xi,$$ (34)

$$c_{k} =\frac{\delta_{k0}}{4},$$ (35)

allows us to write Eq. 32 as

$$\sum_{j=0}^{\infty}\left[\mathsfbi{A}_{kj}+\mathsfbi{D}(\boldsymbol{\mathcal{K}})\mathsfbi{B}_{kj}\right]\textbf{a}_{j}=(\boldsymbol{\mu}-\boldsymbol{\mathcal{K}})c_{k}.$$ (36)

We may now truncate this sum after a finite number of terms, say $M$, so $j,k\in\{0,1,2,\dots,M-1\}$. The resulting problem is an $NM\times NM$ inhomogeneous, eigenvalue problem with up to $N$ eigenvalues. Up to $N$ additional equations are needed to solve this system and may be obtained by the requirement that the vortex boundary is a streamline in the vortex frame, $\psi_{i}+Uy=0$ on $s=1$. This condition gives $N$ equations in the form of a single equation for the coefficient vectors $\textbf{a}_{j}$:

$$\int_{0}^{\infty}\textbf{A}(\xi)\operatorname{J}_{1}(\xi)\,\xi\,\textrm{d}\xi=\sum_{j=0}^{\infty}\textbf{a}_{j}\!\int_{0}^{\infty}\operatorname{J}_{2j+2}(\xi)\operatorname{J}_{1}(\xi)\,\textrm{d}\xi=0\!\quad\!\implies\!\quad\sum_{j=0}^{\infty}(-1)^{j}\textbf{a}_{j}=0,$$ (37)

which may similarly be truncated after $M$ terms.

We note that modon solutions do not exist if there’s a singularity of $[\mathsfbi{K}(\xi)+\mathsfbi{D}(\boldsymbol{\mu})]^{-1}$ for $\xi\in(0,\infty)$ as the integrals $\mathsfbi{A}_{kj}$ and $\mathsfbi{B}_{jk}$ do not converge. These singularities corresponds to the excitation of linear waves in the layered model. In the 2D case, the condition for the existence of singularities is equivalent to the conditions discussed in Kizner et al. (2003).

3.2 The truncated eigenvalue problem

Our truncated eigenvalue problem from Eqs. 36 and 37 may be written as

$$\left[\mathsfbi{A}+\sum_{i=1}^{N}\mathcal{K}_{i}\mathsfbi{B}_{i}\right]\textbf{a}=\sum_{i=1}^{N}(\mu_{i}-\mathcal{K}_{i})\textbf{c}_{i},\quad\textrm{s.t.}\quad\textbf{d}_{i}^{T}\textbf{a}=0\quad\textrm{for}\quad i\in\{1,\dots,N\},$$ (38)

for

$$\mathsfbi{A}=\begin{pmatrix}\mathsfbi{A}_{0,0}&\dots&\mathsfbi{A}_{0,M-1}\\ \vdots&\ddots&\vdots\\ \mathsfbi{A}_{M-1,0}&\dots&\dots\mathsfbi{A}_{M-1,M-1}\end{pmatrix},\quad\mathsfbi{B}=\begin{pmatrix}\mathsfbi{B}_{0,0}&\dots&\mathsfbi{B}_{0,M-1}\\ \vdots&\ddots&\vdots\\ \mathsfbi{B}_{M-1,0}&\dots&\dots\mathsfbi{B}_{M-1,M-1}\end{pmatrix},$$ (39)

$\textbf{a}=(\textbf{a}_{0}^{T},\textbf{a}_{1}^{T},\dots,\textbf{a}_{M-1}^{T})^{T}$, $\textbf{c}=(\boldsymbol{1},\boldsymbol{0},\dots,\boldsymbol{0})^{T}$ and $\textbf{d}=(\boldsymbol{1},-\boldsymbol{1},\dots,(-1)^{M-1}\boldsymbol{1})^{T}$ where 1 and 0 denote length $N$ row vectors of ones and zeros respectively. The indexed quantities $\mathsfbi{B}_{i}$ and $\textbf{c}_{i}$ are given by $\mathsfbi{B}$ and c respectively where all rows with a row index, $n\neq i\pmod{N}$, are set to zero.

Eq. 38 is a multi-parameter eigenvalue problem (Atkinson, 1972) and may be solved using a range of methods (see Appendix A). For active layers, we must solve for the eigenvalue $\mathcal{K}_{i}=-K_{i}^{2}$ while for passive layers we let $\mathcal{K}_{i}=\mu_{i}$ so the term $\mathcal{K}_{i}\mathsfbi{B}_{i}$ may be merged into $\mathsfbi{A}$ and the right-hand side term $(\mu_{i}-\mathcal{K}_{i})\textbf{c}_{i}$ vanishes. Additional conditions of the form $\textbf{d}_{i}^{T}\textbf{a}=0$ are required for active layers but not for passive layers where continuity is automatically enforced by the fact that the components of a corresponding to a passive layer are zero. Mathematically, this is equivalent to the requirement to include an additional equation each time a new (unknown) eigenvalue is introduced. There will be a discrete spectrum of Eigenvalues for each layer, corresponding to different modes in the radial direction. However, modon studies typically focus on the first order radial mode which has the smallest value of $K_{i}$ and is generally the most stable. For multi-layer models, modon solutions may have different radial mode numbers in each layer.

3.3 Determining streamfunctions and potential vorticies

Once the $K_{i}$ and $\textbf{a}_{j}$ are calculated, we may determine the streamfunctions, $\psi_{i}$, by inverting Eq. 26 for $\hat{\boldsymbol{\psi}}_{i}(\xi)\,\xi$ and evaluating the Hankel transform in Eq. 15 directly. Alternatively, using Eqs. 26, 27 and 28 we can show that

$$\boldsymbol{\Delta}_{N}(\boldsymbol{\beta})\,\boldsymbol{\psi}=-\frac{U}{a}\sin\theta\sum_{j=0}^{M-1}\textbf{a}_{j}\operatorname{R}_{j}(r/a),$$ (40)

where $\boldsymbol{\psi}=(\psi_{1},\psi_{2},\dots,\psi_{N})^{T}$ and

$$\boldsymbol{\Delta}_{N}(\boldsymbol{\beta})=\begin{pmatrix}\nabla^{2}-\frac{1}{R_{1}^{2}}-\frac{\beta_{1}}{U}&\frac{1}{R_{1}^{2}}&0&\dots&0&0\\ \frac{1}{R_{2}^{2}}&\nabla^{2}-\frac{2}{R_{2}^{2}}-\frac{\beta_{2}}{U}&\frac{1}{R_{2}^{2}}&\dots&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&\dots&\frac{1}{R_{N}^{2}}&\nabla^{2}-\frac{1}{R_{N}^{2}}-\frac{\beta_{N}}{U}\end{pmatrix},$$ (41)

where $\boldsymbol{\beta}=(\beta_{1},\beta_{2},\dots,\beta_{N})^{T}$. Eq. 40 may be easily inverted using discrete Fourier transforms so this method is often easier numerically than inverting the Hankel transforms. From Eqs. 2, 3 and 4, potential vorticity anomalies may be similarly calculated using $\textbf{q}=\Delta_{N}(\boldsymbol{0})\,\boldsymbol{\psi}$ where $\textbf{q}=(q_{1},q_{2},\dots,q_{N})^{T}$.

4 Examples

We now present three examples to demonstrate that this method gives consistent results with previous studies and can be used to find new modon solutions in layered QG models. A Matlab script to solve for the general N-layer problem is included as supplementary material.

4.1 The Larichev and Reznik dipole

(a)

(b)

Figure 1: A 1-layer modon solution with $(U,a,R_{1},\beta_{1})=(1,1,1,1)$. We show (a) $\psi_{1}$ and (b) $q_{1}$. This solution corresponds to an eigenvalue of $K_{1}\approx 4.108$.

Consider a 1-layer model and choose parameters $(U,a,R_{1},\beta_{1})=1$ and $\boldsymbol{\mathcal{K}}_{1}=-K_{1}^{2}$. Modon solutions may be found using the method of Larichev & Reznik (1976) which combines a Bessel function solution with root finding to determine $K_{1}$. Using our approach, we construct a linear problem of the form

$$\left[\mathsfbi{A}-K_{1}^{2}\mathsfbi{B}\right]\,\textbf{a}=(\mu_{1}+K_{1}^{2})\,\textbf{c},\quad\textrm{s.t.}\quad\textbf{d}^{T}\textbf{a}=0,$$ (42)

using Eqs. 36 and 37. This system may be solved using the method outlined in Appendix A. We determine a value of $K=4.10787...$ which matches the value calculated via root finding to seven significant figures using only $M=12$ terms. Plots of the streamfunction $\psi_{1}$ and potential vorticity anomaly $q_{1}$ are given in Fig. 1. Higher order radial modes and solutions for different parameters, including the case of $R_{1}=\infty$, may be easily calculated using the same method.

4.2 Beta-plane baroclinic topographic modons

Kizner et al. (2003) provide an extensive discussion of baroclinic beta-plane topographic modons using a two layer quasi-geostrophic model. Here we present a two layer solution with $(U,a,R_{1},R_{2},\beta_{1},\beta_{2})=(1,1,1,1,0,1)$ and $(\mathcal{K}_{1},\mathcal{K}_{2})=-(K_{1}^{2},K_{2}^{2})$ corresponding to two active layers. Using our approach, this problem is straightforwardly reduced to solving the two-parameter eigenvalue problem

$$\left[\mathsfbi{A}-K_{1}^{2}\mathsfbi{B}_{1}-K_{2}^{2}\mathsfbi{B}_{2}\right]\textbf{a}=(\mu_{1}+K_{1}^{2})\textbf{c}_{1}+(\mu_{2}+K_{2}^{2})\textbf{c}_{2},\quad\textrm{s.t.}\quad\textbf{d}_{1}^{T}\textbf{a}=\textbf{d}_{2}^{T}\textbf{a}=0,$$ (43)

which may be approached using the methods outlined in Appendix A. Our 2-layer modon solutions are shown in Fig. 2. Solutions with one passive layer—referred to as ‘modons with one interior domain’ by Kizner et al. (2003)—may be obtained by solving the same problem using $K_{2}^{2}=-\mu_{2}$ and neglecting the condition $\textbf{d}_{2}^{T}\textbf{a}=0$. This demonstrates a major advantage of our approach; different vortex regimes can be considered using the same problem setup. By contrast, when using the Bessel function approach of Kizner et al. (2003), care is needed to ensure that the correct Bessel function is used depending on the sign of $\mathcal{K}_{i}$. As such, switching between active and passive layers is less straightforward as it requires both a solution with a different functional form, and new matching conditions across the vortex boundary.

(a)

(b)

(c)

(d)

Figure 2: A 2-layer modon solution with $(U,a,R_{1},R_{2},\beta_{1},\beta_{2})=(1,1,1,1,0,1)$. We show (a) $\psi_{1}$, (b) $q_{1}$, (c) $\psi_{2}$, and (d) $q_{2}$. This solution corresponds to eigenvalues $(K_{1},K_{2})\approx(3.800,3.950)$.

4.3 A mid-depth topographic modon

To demonstrate the flexibility of our approach, we consider a mid-depth vortex in a 3-layer model. We assume that the top and bottom layers are passive and consider an active middle layer. A topographic slope with gradient in the $y$ direction is modelled by taking $\beta_{1}=0,\beta_{2}=0$ in the top two layers and $\beta_{3}=1$ in the bottom layer. For simplicity we take $R_{i}=1$ for $i\in\{1,2,3\}$, $a=1$ and $U=1$. This problem reduces to the eigenvalue problem

$$\left[\mathsfbi{A}+\mu_{1}\mathsfbi{B}_{1}-K_{2}^{2}\mathsfbi{B}_{2}+\mu_{3}\mathsfbi{B}_{3}\right]\textbf{a}=(\mu_{2}+K_{2}^{2})\textbf{c}_{2},\quad\textrm{s.t.}\quad\textbf{d}_{2}^{T}\textbf{a}=0,$$ (44)

which can be straightforwardly solved using standard methods (see Appendix A). The solution consists of an isolated region of strong vorticity in the middle layer, with velocity fields in all three layers. Weak vorticity is generated in the bottom layer due to the advection of fluid in the up-slope ($y$) direction while no vorticity anomaly exists in the surface layer. Fig. 3 shows $\psi_{i}$ for $i\in\{1,2,3\}$ and $q_{2}$ for this solution, the fields $q_{1}=0$ and $q_{3}=\psi_{3}$ are not shown.

(a)

(b)

(c)

(d)

Figure 3: A 3-layer mid-depth vortex solution with $(U,a,R_{1},R_{2},R_{3},\beta_{1},\beta_{2},\beta_{3})=(1,1,1,1,1,0,0,1)$. We show (a) $\psi_{1}$, (b) $\psi_{2}$, (c) $q_{2}$, and (d) $\psi_{3}$. In this case, $q_{1}=0$ and $q_{3}=\psi_{3}$. Layers 1 and 3 are passive and layer 2 has eigenvalue $K_{2}\approx 4.1835$.

5 Discussion and Conclusions

We have presented a semi-analytical method for finding modon, or dipolar vortex, solutions in a layered quasi-geostrophic model. While the system is, in principle, solvable using combinations of Bessel functions inside and outside the vortex region, the matching conditions required to determine all the coefficients can become extremely complicated for multiple layers. Further, these matching conditions must be solved numerically with solutions depending on the vortex parameters in a complicated, nonlinear way. As such, any additional insight obtained from having a closed form solution is limited. By contrast, our procedural method can be straightforwardly used to convert a model with arbitrarily many layers into a multi-parameter eigenvalue problem (Atkinson, 1972) which can be solved numerically.

Additional advantages of our method include the ability to solve for any combination of active and passive layers without having to recalculate any terms in the linear system or change the form of the solution. Though we have not examined this in detail, we expect conditions for the existence of steady modon solutions—corresponding to a lack of linear wave excitation—to follow from examining the singularities of the quantity $[\mathsfbi{K}+\mathsfbi{D}(\boldsymbol{\mu})]^{-1}$. Finally, cases which would have to be treated separately using the standard analytical approach, such as an infinite Rossby radius, require no special treatment here.

We suggest that methods combining Hankel transforms and Zernike radial functions may be relevant to a range of physical problems where important dynamics occur in an isolated circular region. These methods are not restricted to variants of the Laplace and Helmholtz equations, as shown for the Dirichlet-to-Neumann operators of Crowe & Johnson (2023); Johnson & Crowe (2023). Additional effects—such as variations in density difference between layers—could be easily included into our model. It may also be possible to include the effects of layer-dependent background flow or extend the method to a three-dimensional quasi-geostrophic model.

Code Availability. Matlab scripts implementing this method for the N-layer case and demonstrating the three examples considered are included as supplementary material.

Declaration of Interests. The authors report no conflict of interest.

Appendix A Notes on Inhomogeneous Eigenvalue Problems

Consider the problem

$$(\mathsfbi{A}-\lambda\mathsfbi{B})\,\textbf{a}=\textbf{c}_{0}+\lambda\textbf{c}_{1},\quad\textrm{such that}\quad\textbf{d}^{T}\textbf{a}=0.$$ (45)

Assuming that $\mathsfbi{A}$ is non-singular, we can left-multiply Eq. 45 by $\textbf{d}^{T}\mathsfbi{A}^{-1}$ and multiply the resulting scalar by $\textbf{c}_{0}+\lambda\textbf{c}_{1}$ to get

$$-\frac{\lambda(\textbf{c}_{0}\textbf{d}^{T}+\lambda\textbf{c}_{1}\textbf{d}^{T})\mathsfbi{A}^{-1}\mathsfbi{B}}{(\textbf{d}^{T}\mathsfbi{A}^{-1}\textbf{c}_{0})+\lambda(\textbf{d}^{T}\mathsfbi{A}^{-1}\textbf{c}_{1})}\,\textbf{a}=\textbf{c}_{0}+\lambda\textbf{c}_{1},$$ (46)

which can be used to eliminate the inhomogeneity from Eq. 45. The resulting homogeneous problem is the quadratic eigenvalue problem

$$\left[(\textbf{d}^{T}\mathsfbi{A}^{-1}\textbf{c}_{0})\mathsfbi{A}+\lambda\left((\textbf{d}^{T}\mathsfbi{A}^{-1}\textbf{c}_{1})\mathsfbi{A}-(\textbf{d}^{T}\mathsfbi{A}^{-1}\textbf{c}_{0})\mathsfbi{B}-(\textbf{c}_{0}\textbf{d}^{T})\mathsfbi{A}^{-1}\mathsfbi{B}\right)\right.\\ \left.-\lambda^{2}\left((\textbf{d}^{T}\mathsfbi{A}^{-1}\textbf{c}_{1})\mathsfbi{B}+(\textbf{c}_{1}\textbf{d}^{T})\mathsfbi{A}^{-1}\mathsfbi{B}\right)\right]\textbf{a}=0.$$ (47)

Eq. 47 may be put into canonical form and solved as a linear eigenvalue problem for $\lambda$. The corresponding eigenvector, a, should be directly calculated from Eq. 45 as the eigenvectors from Eq. 47 will not necessarily have the correct magnitude since the quadratic problem is linear in a.

Consider the general problem

$$\left[\mathsfbi{A}-\sum_{i=1}^{N}\lambda_{i}\mathsfbi{B}_{i}\right]\textbf{a}=\textbf{c}_{0}+\sum_{i=1}^{N}\lambda_{i}\textbf{c}_{i},\quad\textrm{s.t.}\quad\textbf{d}_{j}^{T}\textbf{a}=0\quad\textrm{for}\quad j\in\{1,\dots,N\},$$ (48)

where $\mathsfbi{A}$ and the $\mathsfbi{B}_{i}$ are $M\times M$ matrices and $M>N$. We can begin by extending the $\textbf{d}_{j}$ to a basis of $\mathrm{I\!R}^{M}$ by introducing $M-N$ orthonormal vectors $\textbf{e}_{k}$ for $k\in\{1,\dots,M-N\}$ such that each $\textbf{e}_{k}$ is orthogonal to the $\textbf{d}_{j}$. Therefore the vector a lies in the $M-N$ dimensional space spanned by the $\textbf{e}_{k}$ and can be written as

$$\textbf{a}=\sum_{k=1}^{M-N}a^{\prime}_{k}\textbf{e}_{k}.$$ (49)

We now define the vector valued function

$$\textbf{f}(\textbf{x})=\left[\mathsfbi{A}-\sum_{i=1}^{N}\lambda_{i}\mathsfbi{B}_{i}\right]\textbf{a}-\textbf{c}_{0}-\sum_{i=1}^{N}\lambda_{i}\textbf{c}_{i},$$ (50)

for $\textbf{x}=[\lambda_{1},\dots,\lambda_{N},a_{1}^{\prime},\dots,a_{M-N}^{\prime}]^{T}$ so solving for our eigenvalues and eigenvectors is equivalent to finding roots of $\textbf{f}(\textbf{x})=\textbf{0}$. Gradients of f with respect to the $\lambda_{i}$ and $a^{\prime}_{k}$ can easily be determined analytically and standard gradient based techniques can be used to find these roots. Alternatively, numerical techniques developed specifically for multi-parameter problems (Atkinson, 1972; Muhic̆ & Plestenjak, 2010) may prove effective, however, these rely on converting the problem to a large single parameter problem which scales poorly with $N$ and $M$.