1. Introduction
Mesoscale and submesoscale ocean vortices play an important role in the distribution of heat, salt and other tracers in the global ocean circulation. These vortices proliferate in both the oceans and atmosphere, and they frequently interact with currents and other vortices (Yasuda, Okuda & Hirai Reference Yasuda, Okuda and Hirai1992; Schultz Tokos, Hinrichsen & Zenk Reference Schultz Tokos, Hinrichsen and Zenk1994; Paireau, Tabeling & Legras Reference Paireau, Tabeling and Legras1997; Carton et al. Reference Carton, Daniault, Alves, Cherubin and Ambar2010; Rodríguez-Marroyo et al. Reference Rodríguez-Marroyo, Viúdez and Ruiz2011; Chouksey Reference Chouksey2023; Zoeller & Viúdez Reference Zoeller and Viúdez2024). In the Mediterranean Sea, for example, Carton et al. (Reference Carton, Daniault, Alves, Cherubin and Ambar2010) compared experimental observations with the results of a three-layer quasi-geostrophic numerical model and concluded that mesoscale vortices interact via their mutual advection by partner vortices. Arguably, understanding the dynamics of mesoscale and submesoscale vortices, and in particular their interactions, is essential to understanding the dynamics of the whole ocean.
Early experimental studies of Brown & Roshko (Reference Brown and Roshko1974) and Winant & Browand (Reference Winant and Browand1974) showed that, in a turbulent mixing layer, like-sign vortices can merge, thereby forming a larger vortex. Vortices may merge if the distance separating them is less than a threshold known as the critical merger distance, see Dritschel (Reference Dritschel1995), Reinaud & Dritschel (Reference Reinaud and Dritschel2005) and references therein. Vortex merger, also known as coalescence or fusion, has been exhaustively described in many experimental and numerical studies (Christiansen & Zabusky Reference Christiansen and Zabusky1973; Dritschel Reference Dritschel1985; Melander, Zabusky & Mcwilliams Reference Melander, Zabusky and Mcwilliams1988; Dritschel & Waugh Reference Dritschel and Waugh1992; Waugh Reference Waugh1992; Ritchie & Holland Reference Ritchie and Holland1993; Perrot & Carton Reference Perrot and Carton2010; Ciani, Carton & Verron Reference Ciani, Carton and Verron2016). The merger of shielded vortices, i.e. vortices consisting of a core surrounded by an opposite-signed vorticity annulus may, however, be inhibited when the vorticity distributions are not overlapping (Ritchie & Holland Reference Ritchie and Holland1993; Valcke & Verron Reference Valcke and Verron1997; Carton Reference Carton2001). In this respect Carton (Reference Carton1992) concluded that shielding generally reduces the critical separation distance for merging, and moreover, in some cases, barotropic instability within each vortex may inhibit their merging. Indeed, in some shielded vortices, the vorticity of the annulus re-organises to generate dipolar moments which make the vortices move apart before they can merge.
The interaction of two shielded vortices depends on their initial separation, their vorticity distribution and other properties such as the baroclinicity and stratification (Carton Reference Carton1992; Chan & Law Reference Chan and Law1995; Valcke & Verron Reference Valcke and Verron1997). In particular, it has been observed that ocean vortices remain approximately isolated and present a shield of relative vorticity, that is, vertical vorticity relative to the Earth’s rotating reference frame (Valcke & Verron Reference Valcke and Verron1997; Schmidt et al. Reference Schmidt, Beckers, Nielsen, Rasmussen and van Heijst1998; Beckers et al. Reference Beckers, Verzicco, Clercx and Van Heijst2001; Viúdez Reference Viúdez2021; Chouksey Reference Chouksey2023). The investigation of how such ocean vortices interact is therefore merited.
It is well known that vortices may either completely merge, partially merge or separate (Ritchie & Holland Reference Ritchie and Holland1993; Holland & Dietachmayer Reference Holland and Dietachmayer1993; Wang & Holland Reference Wang and Holland1995; Valcke & Verron Reference Valcke and Verron1997; Ciani et al. Reference Ciani, Carton and Verron2016), but the detailed physical processes responsible for the formation of robust near-equilibrium multipolar structures are still largely unknown. This study focuses on the interaction of two like-signed, continuous, shielded neutral vortices, both in two-dimensional (2-D) and in three-dimensional (3-D) quasi-geostrophic (QG) flows relevant to the atmosphere/ocean dynamics. We find that vortices can reach an oscillating near-equilibrium state where the vortices attract and repel each other through the exchange of peripheral vorticity. We verify that this new form of interaction is a robust process using simulations conducted by two very different numerical methods, namely a pseudo-spectral code and a contour-advection code, both in 2-D and in 3-D QG flows.
This study is structured as follows. Section 2 describes the mathematical and numerical formulation, and explains the initial conditions used, both in 2-D as well as in 3-D QG. Next, the results of our numerical simulations are presented in § 3. Finally, conclusions are offered in § 4.
2. Mathematical and numerical formulation
2.1. Basic equations
In the 3-D QG model, the dynamical equation describing the flow evolution under inviscid adiabatic conditions is the material conservation of potential vorticity anomaly (hereafter PVA), denoted by
$\varpi^q({\boldsymbol{x}},t)$
, where the superscript
$q$
indicates quasi-geostrophic dynamics. The PVA is advected by the horizontal geostrophic flow
\begin{equation} \frac {\partial {\varpi ^q}}{\partial {t}} + \boldsymbol{u}_{h}^{g} \boldsymbol{\cdot }\boldsymbol{\nabla }_h{\varpi ^q} = 0 . \end{equation}
Here,
$\boldsymbol{u}^g_h(\boldsymbol{x},t)$
is the horizontal geostrophic velocity, where
$h$
indicates horizontal components and
$g$
indicates geostrophic balance, scaled by
$f^{-1}_0$
, where
$f_0$
is the constant background planetary vorticity, or Coriolis frequency. The geostrophic velocity
$\boldsymbol{u}^g_h$
is defined from the geopotential anomaly field
$\phi (\boldsymbol{x}, t)$
by
\begin{equation} \boldsymbol{u}^g_h\equiv -\boldsymbol{\nabla } \times (\phi \, {\hat {\boldsymbol{z}}}) . \end{equation}
The QG PVA,
\begin{equation} \varpi ^q \equiv \zeta ^g + \mathcal{S} = \hat {{\nabla} }^2 \phi , \end{equation}
is defined as the sum of the vertical component of geostrophic vorticity
$\zeta ^g(\boldsymbol{x}, t)$
and the dimensionless vertical stratification anomaly
$\mathcal{S}(\boldsymbol{x},t)$
. The operator
$\hat {\boldsymbol{\nabla }}$
is the 3-D gradient in the rescaled coordinates
$(x,y,\hat {z})$
. Under the leading-order geostrophic approximation, the vertical vorticity
\begin{equation} \zeta ^g \equiv \hat {\boldsymbol{z}} \boldsymbol{\cdot }\boldsymbol{\nabla } \times \boldsymbol{u}_{h}^{g} , \end{equation}
becomes the horizontal Laplacian of the geopotential
\begin{equation} \zeta ^g = {\nabla} ^2_h{\phi } . \end{equation}
The dimensionless vertical stratification anomaly takes the form
\begin{equation} \mathcal{S} \equiv - \frac {\partial {D }}{\partial {z}} = \frac {\partial ^2{\phi }}{\partial {\hat {z}}^2} , \end{equation}
where
$D(\boldsymbol{x},t)$
is the vertical isopycnal displacement, while
$ \hat {z} \equiv (N_0 /f_0 )z$
is the stretched QG vertical coordinate, with
$N_0$
being a constant background Brunt–Väisälä frequency. The vertical isopycnal displacement
$D(\boldsymbol{x},t)\equiv z - \eta (\boldsymbol{x},t)$
, where
$\eta (\boldsymbol{x},t)$
is the vertical location that an isopycnal located at
$\boldsymbol{x}$
at time
$t$
has in the reference configuration defined by
$\varrho _Zz+\rho _0$
. The density anomaly is
$\rho '(\boldsymbol{x},t)\equiv \rho (\boldsymbol{x},t)-\varrho _Z z-\rho _0$
, where
$\rho (\boldsymbol{x},t)$
is the mass density, and
$\rho _0$
and
$\varrho _Z$
are a constant background density and a constant background density stratification. Thus, the density field is expressed in terms of displacements.
In the 2-D model of non–divergent Euler flows, the governing equation reduces to the material conservation of vertical vorticity
$\zeta (\boldsymbol{x}_h,t)$
\begin{equation} \frac {{\rm d}\zeta }{{\rm d}t}\equiv \frac {\partial {\zeta }}{\partial {t}} + \boldsymbol{u}_{h} \boldsymbol{\cdot }\boldsymbol{\boldsymbol{\nabla }}_h{\zeta } = 0 , \end{equation}
where
$\boldsymbol{u}_{h}(\boldsymbol{x}_h,t)$
is the divergence-free horizontal velocity and
$\boldsymbol{\boldsymbol{\nabla }}_h$
is the 2-D gradient operator.
2.2. Initial conditions
This subsection describes the initial conditions used in the 2-D and 3-D QG numerical simulations.
2.2.1. Two-dimensional flows
Following Viúdez (Reference Viúdez2021), the initial vorticity field consists of a finite linear combination of shifted, normalised and truncated cylindrical Bessel modes of order 0,
${J}_0(\kappa \rho )$
, where
$\kappa$
is the cylindrical radial wavenumber (
$\kappa =1$
in the results below), and
$\rho$
is the cylindrical radius of the polar coordinate system
$(\rho ,\varphi )$
. The Bessel modes are truncated at the
$p$
th zero of the cylindrical Bessel function of order 1, namely
$j_{1,p}$
, where
$p \in \mathbb{N}$
, arranged so that
${J}_1(j_{1,p})=0$
. Each vorticity layer mode,
$\zeta _{p}(\rho )$
, is defined piecewise as
\begin{equation} \zeta _p(\rho ) \equiv 2c_p \left \{ \!\!\begin{array}{cc} {J}_0(\kappa \rho ) - {J}_0 (j_{1,p}) & \kappa \rho \leq j_{1,p} \\ 0 & \kappa \rho \gt j_{1,p} , \end{array} \right . \end{equation}
where
$c_{p}$
is a normalisation constant defined by
\begin{equation} c_p \equiv - \frac {\kappa ^2}{{J}_0(j_{1,p}) \, j^2_{1,p}} . \end{equation}
The value of
$c_p$
ensures that each layer-mode
$\zeta _p$
has the same area integral of vorticity (circulation), without loss of generality taken to be
$2\pi$
\begin{equation} \int _0^{2\pi }\int _0^{j_{1,p}/\kappa } \zeta _{p} (\rho ) \, \rho \, {\rm d}\rho \, {\rm d}\varphi = 2\pi . \end{equation}
A vortex is thus constructed as a linear superposition of
$N$
layer modes
\begin{equation} \zeta _{a_1, a_2, \ldots a_N}(\rho ) \equiv \sum _{p=1}^{N} a_p \zeta _p(\rho ) , \end{equation}
where
$a_p$
is the amplitude of mode
$p$
. Since each layer-mode
$\zeta _p$
has the same area integral of vorticity, neutral vortices satisfy the relation
\begin{equation} \sum _{p=1}^{N} a_p = 0 . \end{equation}
2.2.2. Three-dimensional QG flows
In the 3-D QG model, the initial PVA field consists of a finite superposition of radial PVA layer modes
$\varpi _{p}(r)$
(see Zoeller & Viúdez Reference Zoeller and Viúdez2023b
). The layer modes are spherically symmetric, with a PVA distribution given by a truncated, shifted and normalised spherical Bessel mode of order 0,
${j}_0(kr)$
, where
$k$
is the spherical radial wavenumber (
$k=1$
in the results below). Its boundary radius coincides with a zero
$j_{3/2,p}$
of the first-order spherical Bessel function
${j}_1(kr)$
, so that
${j}_1(j_{3/2,p})=0$
. Each PVA layer-mode
$\varpi _{p}(r)$
is defined piecewise as
\begin{equation} \varpi _{p} (r) \equiv 2C_p k^3 \left \{\!\! \begin{array}{cc} {j}_{0}\left ( j_{3/2,p} \right ) - {j}_{0}\left ( kr \right ) & kr \leq j_{3/2,p} \\\\ 0 & kr \gt j_{3/2,p} , \end{array} \right . \end{equation}
where
$C_p$
is a normalisation constant given by
\begin{equation} C_p \equiv \frac {3}{4 {j}_0 \left ( {j}_{3/2,p}\right ) \, {j}_{3/2,p}^3} . \end{equation}
This constant is chosen so that each layer mode for
$p=1,2,\ldots$
has the same volume integral of PVA, namely
\begin{equation} \int _0^{2\pi }\int _0^{\pi }\int _0^{j_{3/2,p}/k} \varpi _{p} (r) \, r^{2} \sin \theta \, {\rm d}r \, {\rm d}\theta \, {\rm d}\varphi = 2\pi , \end{equation}
in spherical coordinates
$(r,\theta ,\varphi )$
.
As in the 2-D case, we define the initial PVA distribution
$\varpi (r)$
as a superposition of
$N$
PVA layer modes
\begin{equation} \varpi (r) \equiv \sum _{p=1}^N a_p \varpi _p(r) . \end{equation}
In the case of neutral vortices, defined as vortices having a vanishing volume integral of PVA, the sum of the layer-mode coefficients
$a_p$
vanishes
\begin{equation} \sum _{p}^{N} a_p = 0 . \end{equation}
Although in general
$a_p\in \mathbb{R}$
, we restrict these values to
$a_p\in \mathbb{Z}$
. Thus, every vortex may be defined by its initial PVA configuration, given by an
$N$
-tuple of layer-mode coefficients
$\boldsymbol{a}\equiv \{ {a_1, \ldots, a_N} \}$
.
Vorticity (
$\zeta , \varpi$
) and azimuthal velocity (
$v$
) profiles for both (a) 2-D and (b) 3-D QG vortices.
In this study we use two of the simplest neutral vortices one may construct in either the 2-D or 3-D QG cases. The first is a two-layer vortex, which is known to evolve spontaneously into a stable tripolar structure (Viúdez Reference Viúdez2021; Zoeller & Viúdez Reference Zoeller and Viúdez2023b ), and has vorticity and PVA distributions
\begin{equation} \zeta _{1,-1} (\rho ) \equiv \zeta _1 (\rho ) - \zeta _2 (\rho ) \quad \textrm {and} \quad \varpi _{1,-1} (r) \equiv \varpi _1 (r) - \varpi _2 (r) . \end{equation}
The second vortex is a stable, axisymmetric three-layer vortex with vorticity and PVA distributions given by
\begin{align} \zeta _{2,-1,-1} (\rho ) \equiv 2\zeta _1 (\rho ) - \zeta _2 (\rho ) - \zeta _3 (\rho ) \quad \textrm {and} \quad \varpi _{2,-1,-1} (r) \equiv 2\varpi _1 (r) - \varpi _2 (r) - \varpi _3 (r).\\[-16pt]\nonumber \end{align}
For simplicity, we refer to these initial vortex profiles as
$\zeta _u$
or
$\varpi _u$
for the unstable two-layer vortex, and
$\zeta _s$
or
$\varpi _s$
for the stable three-layer vortex. The initial velocity and vorticity profiles for the 2-D and 3-D QG cases are illustrated in figure 1.
2.3. Numerical algorithms
We use two different numerical methods to run our simulations and to ensure our main conclusions are robust and do not depend on details of the numerical method. The first one, a purely Eulerian pseudo-spectral method, is detailed in Appendix A. The second one, a hybrid Lagrangian/Eulerian contour-advection method, is detailed in Appendix B. The latter method allows for a particularly fine resolution of the vorticity/potential vorticity field, thereby achieving high accuracy at low computational cost.
Vorticity field
$\zeta \times 10$
of two initially touching
$\zeta _u$
vortices at times (a)
$t=0$
, (b)
$t=2190$
, (c)
$t=2580$
, (d)
$t=3000$
, (e)
$t=4610$
and ( f)
$t=5650$
, using the 2-D pseudo-spectral code.
3. Numerical results
3.1. Two-dimensional flows
This section considers the interaction between two neutral vortices, either with the same initial profile,
$\zeta _u$
, or with different initial profiles,
$\zeta _s$
and
$\zeta _u$
. We investigate the interaction between the vortices for various vortex separations, from touching vortices to vortices separated by
$\delta _0 = j_{1,2}$
.
3.1.1. Interaction between
$\zeta _u$
vortices
It is known that an isolated (neutral) vortex may destabilise and evolve into a tripole, mainly due to the growth of azimuthal mode
$m=2$
(Carton & Legras Reference Carton and Legras1994; Morel & Carton Reference Morel and Carton1994; Kizner & Khvoles Reference Kizner and Khvoles2004; Trieling et al. Reference Trieling, van Heijst and Kizner2010; Viúdez Reference Viúdez2021). The time evolution of two
$\zeta _u$
vortices whose boundaries are touching initially is shown in figure 2. At the initial time, when the two neutral vortices are placed close together (but not overlapping), the vortices do not interact, as the flow field due to each vortex vanishes beyond their boundary. Afterwards, as expected, an azimuthal mode 2 develops on the vortices from numerical noise, eventually causing them to reorganise into robust (stable) tripoles. These tripolar structures now produce an external potential flow (absent for the original vortices). This causes small anomalies in the negative vorticity poles which are no longer perfectly aligned with the positive centre of the vortex core, causing the vortices to move, specifically to attract and eventually exchange vorticity at two different locations symmetrically (figure 3). When in contact, the vortices touch along a curve. Vorticity is then exchanged across the contact curve, redistributing negative vorticity between the two vortices. This exchange generates a dipolar moment within each vortex which (by advection) transfers negative vorticity from one vortex to the other. In this sense, the left vortex transfers negative vorticity at a location
$y\lt 0$
to the right vortex, but amplifies negative vorticity at a location
$y\gt 0$
(figure 3). This slightly changed vorticity anomaly distribution, or dipolar moment, causes each vortex to move in the opposite direction, so that the vortices repel slightly. These vorticity changes are, however, rapidly redistributed azimuthally along the vortex periphery so that their dipolar moment weakens, and the vortices start approaching again figure 2. In this manner, the tripolar structures reach an oscillating near-equilibrium state, where vorticity exchange and vorticity redistribution take place nearly periodically.
Zoom of the interaction shown in figure 2 (c), (d), (e) of the vorticity
$\zeta \times 10$
at times (a)
$t=2580$
, (b)
$t=3000$
, (c)
$t=4610$
.
The pseudo-spectral and the contour-advection simulations produce similar qualitative results, as shown in figures 2 and 4. In both simulations the
$\zeta _u$
vortex evolves into a stable tripole and both capture the vortex interaction arising from vorticity exchange, and an alternating attraction and repulsion between the vortices.
Evolution of the vorticity
$\zeta \times 10$
of two initially touching
$\zeta _u$
vortices, here using the 2-D contour-advection code, combined Lagrangian advection method (CLAM).
This mechanism of vorticity exchange leading to an oscillating near-equilibrium state is new. It has already been shown that two shielded vortices that are in close proximity may evolve through different near-equilibrium states, attracting and/or repelling each other, depending on the changes in the vorticity distribution of their outer shield (Carton Reference Carton1992; Holland & Dietachmayer Reference Holland and Dietachmayer1993; Valcke & Verron Reference Valcke and Verron1997; Carton et al. Reference Carton, Ciani, Verron, Reinaud and Sokolovskiy2016). However, in most cases, excepting very weak interactions, the interaction results in either the formation of a larger merged vortex, a partial merger or straining out where smaller vortices are formed or the repulsion of the two vortices (Carton Reference Carton1992; Ritchie & Holland Reference Ritchie and Holland1993; Holland & Dietachmayer Reference Holland and Dietachmayer1993; Chan & Law Reference Chan and Law1995; Dritschel Reference Dritschel1995; Wang & Holland Reference Wang and Holland1995; Valcke & Verron Reference Valcke and Verron1997; Reinaud & Dritschel Reference Reinaud and Dritschel2005).
The observed differences between the pseudo-spectral and contour-advection simulations mainly concern the precise location of the vortices and the time it takes for the vortices to interact. Such differences are expected because of the very long integration times required for these interaction processes to develop, and the accumulated impact of hyperviscosity in the pseudo-spectral model. As a result, the vorticity leakage in the contour-advection simulation is larger, thus the tripolar structures rotate faster and their displacement is larger than in the pseudo-spectral simulation, cf. Figure 4.
The motion of each vortex has two components, one due to the rotating vorticity poles (anomalies) within the tripole, and the other due to the potential flow generated by its tripole partner. An analysis of the relative motion of the two touching
$\zeta _u$
vortices for the contour-advection simulation follows. The centre of each vortex
$\boldsymbol{X}_j$
(for
$j = 1,\, 2$
) over the time interval
$80 \leq t \lt 592$
is computed every data-save interval
$\Delta t = 1$
. The vortex centres are defined as the geometric centre of the two contiguous regions where the vorticity exceeds the root-mean-square value of the vorticity field in the domain. The vorticity field is split into two parts: the inner vorticity field,
$q_i^j(\boldsymbol{x},t)$
, and the outer vorticity field,
$q_o^j(\boldsymbol{x},t)$
, for
$j = 1, 2$
. These are defined by
\begin{equation} q_{i}^j(\boldsymbol{x},t) \equiv \begin{cases} q(\boldsymbol{x},t), & |\boldsymbol{x} - \boldsymbol{X}_j | \leq r_{c}, \\ 0, & |\boldsymbol{x} - \boldsymbol{X}_j | \gt r_{c} \end{cases}, \qquad q_{o}^j(\boldsymbol{x},t) \equiv \begin{cases} 0, & | \boldsymbol{x} - \boldsymbol{X}_j | \leq r_{c}, \\ q(\boldsymbol{x},t) & | \boldsymbol{x} - \boldsymbol{X}_j | \gt r_{c}, \end{cases} \end{equation}
where
$r_{c} \equiv \alpha R_0$
, with
$R_0$
the initial vortex radius. In practice, we set
$\alpha =1.05$
to take into account the deformation of the vortices. We then evaluate the velocities
$\boldsymbol{u}_{i,o}^j$
induced by
$q_{i,o}^j$
and define
$w_{i,o}^j \equiv \boldsymbol{u}_{i,o}^j \boldsymbol{\cdot }\boldsymbol{t}^j$
as the projection of the induced velocity parallel to the line joining the two vortex centres. Here,
$\boldsymbol{t}^j$
is the unit vector in the direction
$\boldsymbol{X}_{3-j}-\boldsymbol{X}_j$
. Therefore a velocity
$w_{i,o}^j \gt 0$
makes vortex
$j$
move towards vortex
$3-j$
.
Results for
$d=|\boldsymbol{X}_2-\boldsymbol{X}_1|$
and
$w^1_{i,o}$
are presented in figure 5 (results for
$w^2_{ i,o}$
are virtually identical and are not shown). They show that the distance between the two vortex centres is characterised by fast, small-amplitude oscillations superimposed on a slow evolution. These fast oscillations dominate the evolution of the projected velocity
$w_o^j$
but are also noticeable in the evolution of
$w_i^j$
. They are caused by the rotation of the tripoles. To extract the slow motion, the fast oscillations were filtered out from
$w_{i,o}^j$
by applying a low-pass spectral filter over the time interval
$80\leq t \lt 592$
. To do this, we perform a Fourier transform of
$w_{i,o}^j$
and set to zero all Fourier coefficients having wavenumbers exceeding 30; a subsequent inverse transform defines the filtered velocities
$\bar {w}_{i,o}^j$
. The results show that the full projected velocity
$w_i^j+w_o^j$
, and hence the overall relative motion, is dominated by the self-induced velocity field
$w_i^j$
(figure 5). This velocity is due to the dipolar moment that arises from the misalignment of the vorticity poles relative to the tripole centre (i.e. they are not co-linear), a finding which is consistent with previous studies (Holland & Dietachmayer Reference Holland and Dietachmayer1993; Chan & Law Reference Chan and Law1995).
Interaction of two initially touching
$\zeta _u$
vortices simulated using the 2-D contour-advection model, CLAM. The panels show the evolution of various diagnostics: (a) the distance
$d$
between the two vortex centres; (b) the projected velocity
$w_i^1$
(red) and the low-pass filtered projected velocity
$\bar {w}_i^1$
(black); (c) the projected velocity
$w_o^1$
(red) and the low-pass filtered projected velocity
$\bar {w}_o^1$
(black); (d) the full projected velocity
$w_i^1 + w_o^1$
;.
3.1.2. Interaction between
$\zeta _s$
and
$\zeta _u$
vortices
The interaction of different neutral vortices, a
$\zeta _u$
and a
$\zeta _s$
vortex, is discussed next. We illustrate one case where the vortices are initially separated by a distance
$\delta _0 = 0.5j_{1,2}$
between their boundaries (figure 6). The evolution is qualitatively similar to the evolution described in § 3.1.1. The
$\zeta _u$
vortex evolves, as expected, into a robust compact tripole, while the
$\zeta _s$
vortex hardly deforms at all (figure 6
b). The linear stability analysis for both a
$\zeta _u$
and a
$\zeta _s$
vortex when isolated may be found in Appendix C. Only the tripole then generates an exterior potential flow and the vortices attract until they exchange negative vorticity at their peripheries (figure 6
c). This vorticity exchange again generates a dipolar moment in each vortex, which makes them move apart slightly, until the vorticity anomalies are redistributed along the vortex peripheries and the vortices start approaching again. The two vortices maintain an oscillating, near-equilibrium state due to peripheral vorticity exchange and vorticity redistribution. During this interaction the
$\zeta _s$
vortex is only very weakly disturbed.
Evolution of the vorticity
$\zeta \times 10$
of the
$\zeta _u$
and
$\zeta _s$
vortices at times (a)
$t=0$
, (b)
$t=5000$
, (c)
$t=38\,720$
using the 2-D pseudo-spectral code.
The vortex evolution depends weakly on the initial distance
$\delta _0$
between vortices. When the vortex outer layers are initially touching, the near-equilibrium state exhibits oscillations of small amplitude. When
$\delta _0$
increases the vortices also evolve towards a near-equilibrium state, exchanging vorticity, alternatively attracting and repelling each other, but the oscillations become less periodic and the separation distance between them oscillates with larger amplitude. In these cases, the rotation rate of the whole vortex structure is larger, which is associated with the greater erosion of negative vorticity into the exterior flow in form of thin vortex filaments.
3.2. Three-dimensional quasi-geostrophic flows
As in 2-D flows, isolated (neutral) vortices may destabilise and evolve into tripoles in 3-D QG flows (Zoeller & Viúdez Reference Zoeller and Viúdez2023b
; Barabinot et al. Reference Barabinot, Reinaud, Carton, de Marez and Meunier2024). The
$\varpi _u$
vortices evolve into spheroidal stable compact tripoles, and therefore generate an exterior potential flow. The re-organisation of the vortices into tripoles is due to the growth of an azimuthal mode with
$m=2$
as found previously for 2-D vortices. However, once the vortices closely approach, higher modes develop similar to those seen in previous surface QG numerical simulations (Carton et al. Reference Carton, Ciani, Verron, Reinaud and Sokolovskiy2016).
Notably, 3-D QG vortices having comparable maximum vorticity approach much more slowly than their 2-D counterparts. This is due to the facts that (i) the turnover period of a 2-D circular patch of uniform vorticity
$\zeta _0$
is
$T_a=4\pi /\zeta _0$
while it is
$T_a=6\pi /\varpi _0$
for a 3-D QG spherical vortex of uniform PVA
$\varpi _0$
, and (ii) the exterior velocity field decays as
$1/r$
in two dimensions, while decays as
$1/r^2$
in 3-D QG, where
$r$
is the distance from the vortex centre. In fact, it decays even more strongly for neutral vortices having a tripolar shape, by a further factor of
$r$
.
The interaction of two
$\varpi _u$
vortices is shown in figure 7, both from a vertical and a horizontal perspective (top and bottom rows, respectively). Initially, an
$m=2$
perturbation with amplitude 100 times smaller than the peak PVA is added to speed up the early stages of evolution. The spheroidal shape of the QG vortices makes the interaction between the two vortices weaker, since the exchange of PVA takes place mainly near the mid-plane
$z=0$
where the vortices touch. By contrast, close non-neutral 3-D QG vortices strongly interact over a height range comparable to that of the original vortices (Reinaud & Dritschel Reference Reinaud and Dritschel2005). Therefore shielding greatly diminishes the vortex interaction, in particular the exchange of material. The exchange generates a small dipolar moment (smaller than in two dimensions), predominantly in layers close to the vortex mid-plane. This causes the vortices to move slightly apart and the azimuthal PVA variations are subsequently redistributed along the vortex peripheries.
Images of the PVA
$\varpi \times 10$
of two initially touching
$\varpi _u$
vortices at times (a)
$t=0$
and (b)
$t=3960$
using the 3-D QG pseudo-spectral code. The top row shows a side view, while the bottom row shows a top view.
Notably, the interaction and PVA exchange between the vortices, in both the contour advection and pseudo-spectral simulations, does not occur when the vortex boundaries are at least one diameter apart from each other. Then, the vortices evolve into compact tripoles, moving apart from each other at late times. When the initial separation distance is reduced to half of a vortex diameter, PVA exchange and subsequent vortex oscillations again occur. Nolan, Montgomery & Grasso (Reference Nolan, Montgomery and Grasso2001) observed that in 3-D baroclinic hurricane-like vortices an asymmetric instability of azimuthal wavenumber
$m=1$
grows which leads to inner-core vorticity redistribution and the displacement of the vortex centre, so that the vortex acquires a small-amplitude trochoidal wobble.
4. Concluding remarks
The present idealised study has focused on the interaction of neutral, shielded vortices in both 2-D and 3-D QG flows. We have found that some pairs of neutral vortices reach an oscillating near-equilibrium state due to a (potential) vorticity exchange mechanism. This oscillating near-equilibrium mechanism is new. It involves a periodic exchange of vorticity and the generation of dipolar moments within the vortices. These dipolar moments separate the vortices. However, the formation of an exterior potential flow arising from the breaking of circular symmetry, and the subsequent vorticity advection and redistribution of peripheral vorticity causes the vortices to attract. Again a dipolar moment grows, and again the vortices move apart, resulting in an oscillating near-equilibrium state between the vortices.
Our results have been verified using independent and markedly different numerical models, namely a pseudo-spectral model and a high-resolution contour-advection model, both in two and three dimensions. The models exhibit a qualitatively similar behaviour, with differences arising predominantly from numerical diffusion. This shows that the solutions given by Viúdez (Reference Viúdez2021) and Zoeller & Viúdez (Reference Zoeller and Viúdez2023b ) are robust. The idealised initial conditions used in this work suggest that the particular vorticity and potential vorticity distributions assumed could provide a plausible explanation for the persistence of some vortices in the ocean, even when they are interacting with other vortices. Furthermore, we have also seen that shielding, occurring when the area or volume integrated vorticity vanishes within each vortex, greatly weakens and retards the interaction. This is particularly evident in 3-D QG flows, since there potential vorticity exchange only occurs near the vortex mid-planes.
The vortex configurations are not all equally robust. In some cases (and in many additional ones not reported in this study), there is a small vorticity leakage from the vortex peripheries to the exterior flow. This vorticity leakage is seen for example in the contour-advection simulation (figure 4). This vorticity transport to the exterior flow is small but can accumulate with time, thereby partially breaking down the shielding effect, which in turn causes the whole configuration to start rotating. The sign of rotation is the same as that of the exterior flow near the vortices. Eventually, this exterior vorticity can become large enough to prevent vorticity exchange and may arrest the oscillations.
We acknowledge that such vortices, which are no longer entirely shielded, may be widespread in geophysical flows. While we do not address their interactions in the present study, our findings for shielded vortices nonetheless help us understand how shielding can increase the robustness and longevity of vortices. For example, the stable equilibrium state found in this study could be an example of what Chouksey, Gula & Carton (Reference Chouksey, Gula and Carton2023) identified as a ‘doublet’ of two same-signed, deep-water coherent vortices, in their study of the coexistence of two anticyclonic vortices over their lifetimes in the North Atlantic. A similar oscillating behaviour to that seen between neutral vortices in the present study (including small vorticity exchange) has been observed for two finite-area, uniform-temperature vortices when separated at a certain distance by Carton et al. (Reference Carton, Ciani, Verron, Reinaud and Sokolovskiy2016) in a surface QG model. However, the separation distance between the vortices where the near equilibrium occurs has a much wider range in the present study, from vortices touching at their boundaries to vortices being separated between their boundaries by
$\delta _0 = j_{3/2,2}/2 = 3.86$
. Additionally, none of the interactions investigated here led to merging, because of shielding, even when the vortices at the beginning were touching.
In future research, it would be worth carrying out a weakly nonlinear stability analysis on these 3-D QG vortices to better understand the formation of vortex tripoles. Furthermore, one could investigate similar vorticity exchange interactions but extended to a set of multiple (more than 2) vortices similar to the ones described in Viúdez (Reference Viúdez2021). Groups of two or more coherent vortices travelling together have been observed in the oceans in deep water (Chouksey Reference Chouksey2023). Additionally, other parameter variations should be examined, such as shielded vortices but of opposite sign, or two vortices of different sizes, or offset vertically in the 3-D QG model.
Acknowledgements
We are very grateful to three anonymous reviewers for their very helpful comments on this manuscript. We also acknowledge the ‘Severo Ochoa Centre of Excellence’ accreditation (CEX2019-000928-S).
Funding
This work has been funded by the Spanish Government through the project SAGA (Ministerio de Ciencia, Innovacion y Universidades, Ref. No. RTI2018-100844-B-C33) and the FPI PhD grant given to V.Z. (Ministerio de Ciencia e Innovacion PRE2019-090309). Furthermore, the exchange between the ICM and St Andrews University was facilitated by the Erasmus+ Programme of the European Union through a short-term doctoral mobility grant.
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Pseudo-spectral numerical model
The numerical simulations have been carried out using a 2-D and a 3-D pseudo-spectral code in doubly and triply periodic domains, respectively, using a leap-frog time-stepping method together with a weak Robert–Asselin time filter to avoid the decoupling of even and odd time levels. This follows the approach of Dritschel & Viúdez (Reference Dritschel and Viúdez2003), where the equation of material conservation of vorticity (2.7) or QG PVA (2.1) is numerically integrated in time, as detailed in Zoeller & Viúdez (Reference Zoeller and Viúdez2023a
). The numerical hyperviscosity e-folding rate
$e_f$
was chosen as small as possible (
$e_f = 0.2$
in 3-D QG and
$e_f=1$
in two dimensions) to prevent the growth of numerical noise.
In the 2-D case, the spatial domain uses a numerical grid of
$n_x \times n_y = 1024\times 1024$
grid points, and the horizontal extent is
$L_x=L_y= 40$
. In the 3-D QG case, the number of grid points is
$n_x \times n_y \times n_z = 256^3$
, which is relatively coarse in order to enable the analysis of a large number of long numerical simulations. The time step is fixed at
$\delta t=0.01$
as there is no growth of vorticity or PVA – the flow does not speed up. The horizontal and (stretched) vertical extent is
$L_x = L_y = c_0 L_z = 40$
(where
$c_0=N_0/f_0$
), which is large compared with the radii of the vortices,
$r_p=j_{1,p}$
with
$\{ j_{1,2}, j_{1,3}\}\simeq \{ 7.01 , 10.17 \}$
in the 2-D case and
$\{ j_{3/2,2}, j_{3/2,3}\}\simeq \{ 7.73 , 10.9 \}$
in the 3-D case. The 3-D numerical results are independent of the ratio
$c_0 \equiv N_0 /f_0$
. Typical geophysical (atmospheric and oceanographic) applications of the 3-D QG model assume values
$c_0\in [10,100]$
.
Appendix B. Contour-advection numerical model
In contour advection, the vorticity field in two dimensions, or the PVA field on each isopycnal surface in three dimensions, is discretised into a finite set of contours, each of which consists of a variable number of nodes connected together by local cubic splines. Across each contour, the vorticity or PVA jumps by a fixed, prescribed amount, equal to
$(\zeta _{max }-\zeta _{min })/n_l$
or
$(\varpi _{max }-\varpi _{min })/n_l$
, with
$n_l=100$
(a large value which well represents the actual continuous distributions).
The 2-D model uses the CLAM developed by Dritschel & Fontane (Reference Dritschel and Fontane2010), which is an extension of the contour-advective semi-Lagrangian (CASL) algorithm developed by Dritschel & Ambaum (Reference Dritschel and Ambaum1997). The CLAM combines a pseudo-spectral method, which is especially accurate for large scales and energy conservation, with Lagrangian contours, which accurately and efficiently represent intermediate to small scales, well below those available to the pseudo-spectral method. The CLAM allows for the representation of general forcing and dissipation (non-conservation), although this feature is not used here. The basic (inversion) grid resolution is
$1024^2$
, but studies show that the effective grid resolution is approximately 16 times higher in each dimension (Dritschel & Tobias Reference Dritschel and Tobias2012). The range of the spatial numerical domain is
$[-\pi , \pi ] \times [-\pi , \pi ]$
.
The 3-D QG model extends the 2-D CASL model (Dritschel & Ambaum Reference Dritschel and Ambaum1997), to permit one to study the dynamics of complex, fine-scale conservative potential vorticity distributions in three dimensions (Dritschel, de la Torre Juárez & Ambaum Reference Dritschel, de la Torre Juárez and Ambaum1999). The basic (inversion) grid resolution is
$256^3$
, although the effective resolution is much higher. In the vertical, 4 times as many isopycnal surfaces are used than grid points to better represent the vertical variation of the potential vorticity field. We fix the extent of the (triply periodic) physical domain to be
$(2\pi )^3$
and rescale the radii of the vortices and the distances between them to ensure that all vortex configurations have identical total circulation and are not strongly affected by the flow periodicity. The method is essentially inviscid, with ‘contour surgery’ being the only source of pseudo-diffusion (Dritschel Reference Dritschel1988). The initial conditions are unperturbed.
Appendix C. Linear stability of a two-dimensional shielded vortex
We analyse the linear stability of a single 2-D shielded vortex. The vorticity field
$\zeta$
is discretised using
$N$
equal jumps
$\Delta q = (\zeta _{max } -\zeta _{min })/N$
spanning the range of vorticity
$(\zeta _{min },\zeta _{max })$
within the vortex. We then discretise the vortex by
$N_c \gt N$
circular contours
${\mathcal C}_j$
of radii
$\bar {\rho }_j$
, for
$1\leq j \leq N_c$
. The radii are determined from
$\zeta (\bar {\rho }_j) = \zeta _{min } + (i-0.5)\Delta q$
, for
$1\leq i \leq N$
(
$N_c$
exceeds
$N$
because
$\zeta$
is non-monotonic in a shielded vortex). We associate the signed vorticity jump
$\Delta q_j = \pm \Delta q$
with the contour
${\mathcal C}_j$
depending on whether the vorticity increases or decreases across the contour in the direction of decreasing radius.
To examine linear stability, we displace each contour radially
\begin{equation} \rho _i(\varphi ,t) = \bar {\rho }_i + \epsilon \eta _i {e}^{{i} (m\varphi -\sigma t)}, \end{equation}
where
$\epsilon \ll 1$
,
$\epsilon \eta _i$
is the complex amplitude (the real part of
$\rho _i$
is intended), and
$\sigma (m)=\sigma _r(m) +{i} \sigma _i(m)$
is the eigenvalue to be determined (
$\sigma _r$
gives the growth rate while
$\sigma _i$
gives the frequency of wave-like disturbances). The circular symmetry of the basic-state flow permits us to study each azimuthal wavenumber
$m=1,\,2\,\ldots$
separately.
The kinematic condition on contour
${\mathcal C}_i$
is
\begin{equation} \dfrac {\textrm{D} \rho _i}{\textrm{D} t} = {i}(m\bar \varOmega _i-\sigma )\eta _i = v_i, \end{equation}
where
$\bar \varOmega _i=\bar {u}_i/\rho _i$
is the basic-state angular velocity, while
$v_i$
is the radial velocity. The latter is obtained by contour integration over the
$N_c$
contours
\begin{equation} v_i(\varphi ,t) = \sum _{j=1}^{N_c} \Delta q_j \oint _{{\mathcal C}_j} G(\rho _i(\varphi ,t),\rho _j(\varphi ',t)) {\rm d}\varphi ', \end{equation}
where
$G$
is the Green’s function giving the contribution of contour
${\mathcal C}_j$
to contour
${\mathcal C}_i$
. Linearising (C3) for
$\epsilon \ll 1$
leads to
\begin{equation} (m\bar \varOmega _i-\sigma )\eta _i= \sum _{i=1}^{N_c} \Delta q_j {\mathcal K}(\bar {\rho }_i,\bar {\rho }_j) \eta _j , \end{equation}
where
\begin{equation} {\mathcal K}(\bar {\rho }_i,\bar {\rho }_j) = \dfrac {1}{2} \begin{cases} \left ( \dfrac {\bar {\rho }_i}{\bar {\rho }_j}\right )^{1+m}, & \textrm {for} \quad \bar {\rho }_i\lt \bar {\rho }_j, \\[12pt] \left ( \dfrac {\bar {\rho }_i}{\bar {\rho }_j}\right )^{1-m}, & \textrm {for} \quad \bar {\rho }_i\geq \bar {\rho }_j,\end{cases} \end{equation}
see Flierl (Reference Flierl1988). For each
$m$
, this results in an
$N_c\times N_c$
eigenvalue problem for
$\sigma$
and the associated eigenvector
$\eta _j$
(
$1\leq j \leq N_c$
). In general, there are
$N_c$
eigenvalues
$\sigma$
.
For the
$\zeta _u$
vortex, only azimuthal mode
$m=2$
is linearly unstable, with a growth rate
$\sigma _i=0.042$
and
$\sigma _r=0.25$
(to two digits of accuracy). Convergence of the results has been tested by varying the number of vorticity jumps
$N$
from 100 to 5000. Figure 8 shows the contour deformation of the eigenmode for
$N=100$
. The greatest deformation occurs near the centre of the vorticity profile, in the region of weak negative vorticity. By contrast, the
$\zeta _s$
vortex has been found to be linearly stable for all
$m$
.
Contours (vorticity jumps)
$C_j$
,
$1\leq j\leq N_c$
, deformed by the unstable eigenmode with azimuthal number
$m=2$
, using
$N=100$
, for the
$\zeta _u$
vortex. The red contours correspond to the region where
$\zeta \gt 0$
, while the blue contours correspond to the region where
$\zeta \lt 0$
.
Open access
