We systematically study the evolution of Larichev–Reznik dipoles in an equivalent-barotropic quasi-geostrophic beta-plane model in high-resolution numerical simulations. Our results shed light on the self-organization and rich dynamics of dipolar vortices, which are ubiquitous in geophysical flows. By varying both dipole strength and initial angle α0 of dipole tilt to the zonal direction, we discover new breakdown mechanisms of the dipole evolution. The dipoles are quickly destroyed by Rossby wave radiation, if initial tilt is too large or the dipole is too weak; otherwise, via damped oscillations, the dipoles tend to adjust themselves to different states drifting eastward. Two competing physical mechanisms that govern dipole transformations are found: (1) spontaneous dipole instability due to a growing critical linear mode and (2) meridional separation of dipole partners that accumulates over the adjustment period and prevents the above instability. Which mechanism prevails depends on the initial tilt and dipole strength, and the details of this are discussed.

The dynamics of isolated and coherent vortices has attracted the attention of many researchers trying to understand long-lived geophysical vortices, which are ubiquitous in the ocean1–3 and atmospheres.4–6 Within this broad class of phenomena, we focus on vortex dipoles (i.e., pairs of opposite-sign vortices) and their long-term evolution. Dipolar couples transporting fluid inside their cores across large distances are widespread at the ocean surface,7,8 stimulating growing interest in the underlying vortex dynamics.9–11 

A classical example exhibiting steady propagation is the Lamb–Chaplygin dipole (hereafter, LCD),12,13 which is a solution of the two-dimensional (2D) Euler equations. As indicated in Ref. 14, Sec. 4.3, an LCD with a circular core tends to slowly evolve toward a state with smoother vorticity distribution in a slightly elliptical core. While no instability was detected in the inviscid discretized version of the LCD with circular separatrix,15 the authors noted that the lowest vorticity levels (not represented in their isovortical discretization) might be rapidly stripped from the configuration through the neighborhood of the rear stagnation point. Indeed, linearly unstable modes were found when explicit viscosity was included in the system.16 

Building from this, a 2D stationary solution (so-called modon) on the β-plane (taking into account planetary rotation and sphericity) was obtained because Stern17 was found to be unstable.18 Furthermore, a more general solution of the classical quasigeostrophic (QG) equivalent-barotropic model (given by rotating top-layer shallow-water dynamics with deformable lower interface above the motionless deep layer) for propagating dipoles with a circular separatrix was obtained by Larichev and Reznik (hereafter, LRD).19 

When a finite internal Rossby radius of deformation is assumed, the LRD solutions are capable of zonal drift in both eastward and westward directions. The eastward propagation speed is unbounded, whereas the propagation speed to the west must exceed the maximum Rossby wave phase speed20 to prevent the vortex losing energy to generated waves.21 In the past, westward propagating LRDs were suggested as a paradigm for atmospheric blocking,22,23 but were later found unstable,24 hence not suitable for this due to rapid destruction. On the other hand, eastward propagating LRDs were found to be remarkably robust in numerical simulations, even in the presence of weak friction,25 short-wave perturbations,26 and topographic perturbations.27 Proof of their stability, however, remains evasive.24,28,29

Recently, the phenomenon of spontaneous symmetry breaking was discovered in eastward weak-intensity LRDs leading to their ultimate destruction.30 The corresponding unstable growing perturbation has even symmetry about the zonal axis (similar to the unstable modes obtained for the evolving viscous LCD16) despite the odd flow symmetry of the LRD. A normal mode representation of this growing perturbation accurately replicates the time evolution, confirming that the eastward LRD is linearly unstable. This growing mode was initially referred to as the A-mode,30 in accord with some earlier notation for a similar mode; however, in this study, we refer to it as the Davies mode (D-mode), to emphasize its discovery and differentiate it from the A-mode for the LCD when viscosity is imposed in the system.

The effect of dipole tilt relative to the zonal direction was examined in Ref. 31, where two regimes were reported. For weak dipoles with a large tilt, the dipole rapidly disintegrated. However, for most of the cases studied, the dipole experienced damped oscillations along the zonal axis before adjusting to new eastward propagating steady states. Here, we include the effects of dipole tilt using the numerical framework of Ref. 30. We argue that the instabilities were not captured in Ref. 31 due to computational constraints on the spatial resolution of the simulations.

The structure of this paper is the following: in Sec. II, we formulate the nondimensional model, derive the tilted dipole state, and briefly introduce an asymptotic theory for strong dipoles; in Sec. III, we detail our findings for both strong and weak dipoles; and in Sec. IV, we summarize and discuss the results.

To connect to previous studies, we adopt the equivalent-barotropic quasigeostrophic (QG) model on the beta-plane with flat bottom relief,32,33
D Π ̂ D t ̂ = 0 , Π ̂ = q ̂ + β ̂ y ̂ ,
(1)
where β ̂ represents the meridional gradient of the Coriolis parameter; Π ̂ = Π ̂ ( x ̂ , y ̂ , t ̂ ) is the absolute potential vorticity (PV)—a materially conserved quantity on the Lagrangian fluid elements—and the potential vorticity anomaly (PVA) is defined as
q ̂ = ̂ 2 ψ ̂ R ̂ d 2 ψ ̂ ,
(2)
with R ̂ d denoting the internal Rossby radius of deformation and ψ ̂ = ψ ̂ ( x ̂ , y ̂ , t ̂ ) denoting the velocity streamfunction. For compactness, we make use of the differential operators,
D D t ̂ : = t ̂ ψ ̂ y ̂ x ̂ + ψ ̂ x ̂ y ̂ , ̂ 2 : = 2 x ̂ 2 + 2 y ̂ 2 ,
(3)
and hat notation to indicate dimensional physical quantities.
We impose a vortex centered at coordinate ( x ̂ c , 0 ) and perform the change of coordinates x ̂ = X ̂ + x ̂ c ( t ̂ ) and y ̂ = Y ̂ to reposition the origin of our reference frame to the vortex center. Next, Eq. (1) is nondimensionalized by introducing the length and velocity scales, L ̂ and U ̂, respectively,
q t d x c d t q X + J X , Y ( ψ , q ) + β ψ X = 0 ,
(4)
where the Jacobian operator is defined as
J X , Y ( A , B ) : = A X B Y A Y B X ,
(5)
for functions A and B, and the absence of hat notation indicates the following nondimensional quantities:
( X , Y , x c ) = L ̂ 1 ( X ̂ , Y ̂ , x ̂ c ) , t = L ̂ 1 U ̂ t ̂ , β = L ̂ 2 U β ̂ ̂ , ψ = ( L ̂ U ̂ ) 1 ψ ̂ , q = L ̂ U ̂ 1 q ̂ .
(6)
Focusing on steady vortices with purely zonal drift yields the following governing equation:
J X , Y ( ψ + c Y , 2 ψ + β Y ) = 0 ,
(7)
where c = d x c / d t , β = β + c S and S = ( L ̂ / R ̂ d ) 2.
Following Ref. 30, the PVA field can be decomposed into two uniquely defined fields,
q A = q ( X , Y , t ) + q ( X , Y , t ) 2 , q S = q ( X , Y , t ) q ( X , Y , t ) 2 ,
(8)
where qA denotes the A-component, which has even PVA relative to the zonal axis, and qS denotes the S-component, which has odd PVA relative to the zonal axis. The benefit of this presentation is that we can interpret each component as follows: the S-component perturbation leads to symmetric deformation of the vortices around the zonal axis, whereas the A-component perturbation leads to the antisymmetric deformation.
Given that ψ A , S denotes the A- and S-component of the velocity streamfunction, respectively, the quadratic invariants of energy and enstrophy for each component are obtained by integrating over the doubly periodic domain Σ,
E A , S = 1 2 Σ [ ( ψ A , S ) 2 + S ψ A , S 2 ] d Σ , Z A , S = 1 2 Σ q A , S 2 d Σ ,
(9)
where d Σ = d X d Y; and the summations E = E A + E S and Z = Z A + Z S are conserved for the inviscid case, as described by Eq. (7).

Explicit Newtonian viscosity can be incorporated into the equivalent-barotropic model by setting the left-hand side of Eq. (7) equal to R e 1 4 ψ, where R e = ( L ̂ U ̂ ) 1 ν ̂ is the Reynolds number, and ν ̂ denotes eddy viscosity. However, since we are interested in solutions to the inviscid problem (7), we make use of very small implicit numerical viscosity in our simulations, which we discuss later. This decision is also motivated by the results of the sensitivity analysis carried out in Ref. 30, where consistency is found between simulations using explicit Newtonian viscosity and implicit numerical viscosity.

A class of exact solutions to Eq. (7) over an infinite domain are the LRD steady states, which are derived from the relative vorticity relation,
2 ψ = { k 2 ψ c ( k 2 + p 2 ) r sin ( ϑ α 0 ) , r 1 , p 2 ψ , r > 1 ,
(10)
where ( r , ϑ ) are standard polar coordinates. Here, r = 1 defines a circular separatrix ( L ̂ is chosen to be the vortex radius); p 2 = β / c + S > 0; α 0 = 0 ° , 180 ° (eastward and westward propagation, respectively); and k is a positive constant satisfying the nonlinear equation,
k J 1 ( k ) J 2 ( k ) = p K 1 ( p ) K 2 ( p ) ,
(11)
where J μ and K μ are the order-μ Bessel and modified Bessel functions of the first kind, respectively. This equation has an infinite number of solutions for k, but we considered only dipoles with the lowest k. Based on the work of Ref. 18 for shielded vortex structures on a QG beta-plane and Ref. 11 for surface-QG vortex solutions with higher radial mode numbers, we anticipate that higher order modes will be unstable and rapidly break down.
We further assume that the dipole is spatially localized, i.e., ψ 0, as r (more specifically, this is exponential decay), and the solution is continuous and continuously differentiable across the separatrix. The streamfunction field can be expressed in the following form:
ψ ( r , ϑ ) = crF ( r ) sin ( ϑ α 0 ) ,
(12)
where
F ( r ) = { ( p / k ) 2 ( J 1 ( k r ) / r J 1 ( k ) 1 ) 1 , r 1 , K 1 ( p r ) / r K 1 ( p ) , r > 1 ,
(13)
and the corresponding PVA components are as follows:
q A ( r , ϑ ) = X Q ( r ) sin α 0 , q S ( r , ϑ ) = Y Q ( r ) cos α 0 ,
(14)
where
Q ( r ) = { c ( F ( r ) + 1 ) ( S + k 2 ) β r 1 , β F ( r ) , r > 1 .
(15)

In this work, we investigate tilted dipoles for S = 1, i.e., L = Rd, launched at angles 0 < α 0 90 ° to the zonal axis. As a consequence of this, qA is nonzero and the tilted LRDs do not satisfy (7), i.e., the dipoles we consider are no longer steady states. Note that in this case, the A- and S-components are computed about the zonal coordinate axis, rather than the dipole translation axis, which are equivalent only for the eastward LRD. Figure 1 contains visualizations of the A- and S-components and will help the readers to understand the story, as it evolves. Hereafter, the A-component as depicted in Fig. 1(a) is referred to as the tilted mode (T-mode).

FIG. 1.

Physical fields corresponding to the dipole with β = 2 c and tilted at α 0 = 5 °: (a) A-component of PVA, qA (i.e., the T-mode); (b) S-component of PVA, qS; (c) PVA; (d) PV.

FIG. 1.

Physical fields corresponding to the dipole with β = 2 c and tilted at α 0 = 5 °: (a) A-component of PVA, qA (i.e., the T-mode); (b) S-component of PVA, qS; (c) PVA; (d) PV.

Close modal

In the (eastward propagating) case α 0 = 0 ° and for dipoles with β / c 2, the component qA (i.e., D-mode) was initially zero but appeared due to spontaneous symmetry breaking and shown to oscillate and grow exponentially over time.30 This growth caused the dipole to develop asymmetries through elongation and compression of the vortex pair (Fig. 2). The ultimate result was the destruction of the dipole, with the vortex partners decelerating and drifting apart, before propagating in the (opposite) westward direction. At later times, they continued drifting apart and eventually disintegrated into the background flow as in Ref. 34. A combination of dynamical evolution associated with T-mode depicted in Fig. 1(a) and D-mode (Fig. 2) is expected to be seen in the case of α 0 0 °.

FIG. 2.

Upper panels show snapshots of the PVA of the D-mode that develops during the intermediate stage of linear growth for a weak eastward dipole, while the bottom panels show the associated elongations and compressions that the dipole pair experiences, in response to the growing D-mode [e.g., (a) causes deformations so that the dipole tilts, as seen in (e), whereas, (b) causes zonal elongation and zonal compression in the anticyclone and cyclone, respectively, as seen in (f)].

FIG. 2.

Upper panels show snapshots of the PVA of the D-mode that develops during the intermediate stage of linear growth for a weak eastward dipole, while the bottom panels show the associated elongations and compressions that the dipole pair experiences, in response to the growing D-mode [e.g., (a) causes deformations so that the dipole tilts, as seen in (e), whereas, (b) causes zonal elongation and zonal compression in the anticyclone and cyclone, respectively, as seen in (f)].

Close modal
For dipoles, let us consider the coordinates of the positive and negative PVA extrema (X1, Y 1) and (X2, Y 2), respectively. The dipole tilt is, then, characterized by
Δ X = X 1 X 2 = D sin α , Δ Y = Y 2 Y 1 = D cos α ,
(16)
where D is the distance between extrema,
D = ( Δ X ) 2 + ( Δ Y ) 2 ,
(17)
and α is the evolving angle, such that α ( 0 ) = α 0. The dipole center is defined as (Xc, Y c), with
X c = X 1 + X 2 2 , Y c = Y 1 + Y 2 2 ,
(18)
while the zonal and meridional drift components of the dipole center,
u c = d x c d t , v c = d y c d t ,
(19)
respectively, give the propagation velocity. These components allow us to obtain the corresponding dipole speed as
V = u c 2 + v c 2 .
(20)
For very strong dipoles ( c β), a simple kinematic theory described in Ref. 35 predicts the center trajectory, ( x c ( t ) , y c ( t ) ), by neglecting changes in the distance between the dipole extrema (i.e., D = D0) and in the dipole center velocity V = c. This theory tells us that the evolution of the dipole angle, α ( t ), is given by the following equation:
d α d t = β V M d y c ,
(21)
where M d = Q ( r ) r 3 d r 3 V is the magnitude of the LRD moment, which is described by Eqs. (13)–(15) with β = 0 and normalized by the area of fluid trapped inside the separatrix. Combining Eq. (21) with
d y c d t = V sin α
(22)
allowed the full trajectory to be predicted by the equation for the physical pendulum. Thus, for the north-eastward LRDs, we expect to see oscillatory behavior in the displacement of the dipole center around the equilibrium latitude. The maximum excursion Ym is given in Ref. 35 as
Y m = 2 M d β sin ( α 0 2 ) .
(23)
Both the kinematic theory (21) and (22) and numerical simulations in Ref. 31 showed that initially tilted LRDs evolve by damped oscillations, and the corresponding mechanism was identified by Ref. 36, where a PV conservation argument showed that the adjustment is achieved through the dipole losing enstrophy to the surrounding background flow. This mechanism can be incorporated into Eq. (21) as
d α d t = β V ( y c + λ sin α ) M d .
(24)
In this case, the dipole evolution calculated from Eqs. (22)–(24) displays decaying oscillations, which depend on the damping parameter λ. Neither this theory nor the numerical solutions in Ref. 31 take into account potential instability of the eastward LRDs,30 which motivates more in-depth studies. Now that we have reviewed some theoretical work pertaining to tilted dipoles, and we proceed in the following section to extract new information for tilted LRD dynamics using high-resolution numerical simulations.
We solve Eq. (1) on a doubly periodic computational domain of size ( 4 L Y , L Y ) = ( 60 , 15 ) with 4 N × N = 8192 × 2048 nodes. Motivated by mesoscale observations, we consider dipoles on the scale of the internal Rossby radius of deformation, that is, L ̂ = R ̂ d (consequently, S = 1). Furthermore, to draw comparisons with Ref. 31, we nondimensionalize so that c = 0.1 (so that U ̂ = 10 c ̂, where c ̂ is the dimensional equivalent of c) and explore the following parameter space:
( α 0 , β ) = [ 5 ° , 30 ° , 45 ° , 60 ° , 90 ° ] × [ 1 , 2 , 3.5 , 4 ] c ,
(25)
where α 0 45 ° is referred to as small tilt, and α 0 > 45 ° is referred to as large tilt. We also use the terms weak and strong to correspond with c < β 4 c and β c, respectively. Benchmark solution animations for this parameter space are shown in Fig. 3 (Multimedia view).
FIG. 3.

Dipole evolution animations corresponding to the parameter space described in Eq. (25) for 0 t 200. We present animations for this parameter space with values of α0 increasing from top to bottom. The strongest LRD with c / β = 1 is shown in the left column, while other columns show progressively weaker LRDs, with the right-most column corresponding to c / β = 0.25. Multimedia available online.

FIG. 3.

Dipole evolution animations corresponding to the parameter space described in Eq. (25) for 0 t 200. We present animations for this parameter space with values of α0 increasing from top to bottom. The strongest LRD with c / β = 1 is shown in the left column, while other columns show progressively weaker LRDs, with the right-most column corresponding to c / β = 0.25. Multimedia available online.

Close modal

To verify our results, we made use of two distinctly different numerical models. For consistency with Ref. 30, we use finite differences with the CABARET advection scheme; the algorithm is described in detail in Ref. 37, and the numerical convergence properties of the model in relevant flow regimes with rich mesoscale dynamics are discussed in Refs. 38 and 39.

For consistency with Ref. 31 and to validate our results, we also performed simulations using the Dedalus Python package, which combines a pseudo-spectral approach in space with a fourth-order implicit–explicit Runge–Kutta scheme for time integration.40 This approach necessarily includes explicit dissipation represented through a hyperdiffusion term, whereas the CABARET approach has implicit numerical diffusion, which is minimized at each time step through the advection scheme.

We estimated a numerical viscosity value corresponding to dissipation in our CABARET simulations and used it to set equivalent hyperdiffusion parameter in the pseudo-spectral model. We found that both methods yield similar solutions (see the supplementary material); therefore, we primarily present CABARET results, though supplementary animations in Fig. 3 are obtained with the Dedalus package.

The dipole evolution is analyzed by considering both the trajectory of the dipole center (xc, yc) and changes in the internal dipole structure, as characterized by the propagation speed, V, and the distance between extrema, D, defined in Sec. II C. The coordinates of the positive and negative dipole extrema (X1, Y 1) and (X2, Y 2) are evaluated using local 2D-parabolic interpolation.

Comparing with the kinematic theory, we found new behavior in the evolution of dipoles with β / c = 1 (strong dipoles) and with β / c 2 (weak dipoles). These modifications to the dipole dynamics are mostly related to the growing D-mode on top of decaying T-mode, which become clearly visible when the vortex crosses the zonal axis.

Strong dipoles adjust themselves along the zonal axis and appear to continue eastward propagation. Nevertheless, over the course of this adjustment phase, one can see moderate increase in partner separation and associated deceleration in the upper panels in Fig. 4 and Table I. Note that the originally D-mode was described only for weak eastward LRD solution.30 Here, growing D-mode is revealed also for the strong dipoles though they are not destroyed during the time of integration and their structure still resembles the initial LRD.

FIG. 4.

( V / c , D )-plots. Upper rows with β = c and: (a) α 0 = 5 ° (blue), α 0 = 30 ° (red), (b) α 0 = 45 ° (blue), and α 0 = 60 ° (red); and lower rows with β = 2 c and the same values of α0 in corresponding subplots. The data are for 0 t 400.

FIG. 4.

( V / c , D )-plots. Upper rows with β = c and: (a) α 0 = 5 ° (blue), α 0 = 30 ° (red), (b) α 0 = 45 ° (blue), and α 0 = 60 ° (red); and lower rows with β = 2 c and the same values of α0 in corresponding subplots. The data are for 0 t 400.

Close modal
TABLE I.

Percentage change in partner separation and dipole deceleration for β = c, with first entries corresponding to t = 200 and second entries corresponding to t = 400.

α0 D ( t ) / D 0 1 V ( t ) / c 1
5 °  (0.5, 1.1)%  −(1.0, 3.0)% 
30 °  (1.0, 6.0)%  −(5.2, 20.1)% 
45 °  (0.71, 5.1)%  −(7.1, 21.0)% 
60 °  (0.14, 1.5)%  −(6.7, 12.2)% 
90 °  (1.7, 3.2)%  −(16.7, 25.5)% 
α0 D ( t ) / D 0 1 V ( t ) / c 1
5 °  (0.5, 1.1)%  −(1.0, 3.0)% 
30 °  (1.0, 6.0)%  −(5.2, 20.1)% 
45 °  (0.71, 5.1)%  −(7.1, 21.0)% 
60 °  (0.14, 1.5)%  −(6.7, 12.2)% 
90 °  (1.7, 3.2)%  −(16.7, 25.5)% 

Weak dipoles for α 45 ° decelerate with partners separating by more than 30% and a reversal in the direction of zonal propagation (lower panels in Fig. 4), similar to the effects of a growing D-mode in the eastward LRD.30 For strong and weak dipoles with a range of tilts, our findings are summarized in Fig. 5, where there are both similarities with and notable differences from Ref. 31.

FIG. 5.

( α 0 , β )-parameter space summarizing possible outcomes for the dipole evolution for 0 t 400: open circles indicate destruction due to early dipole-Rossby wave interaction when initial tilt is large; diamonds indicate slow D-mode-induced destruction with t s > 100 (where ts is the time at which the dipole changes direction); and filled circles indicate slow decay of dipole intensity and increased meridional partner separation that accumulates over adjustment period and prevents instability.

FIG. 5.

( α 0 , β )-parameter space summarizing possible outcomes for the dipole evolution for 0 t 400: open circles indicate destruction due to early dipole-Rossby wave interaction when initial tilt is large; diamonds indicate slow D-mode-induced destruction with t s > 100 (where ts is the time at which the dipole changes direction); and filled circles indicate slow decay of dipole intensity and increased meridional partner separation that accumulates over adjustment period and prevents instability.

Close modal
Hence, to clarify novelties of the observed phenomena, we make a detailed comparison of the various regimes over the following sections where the modification in the spatial structure of the dipoles is illustrated by snapshots of the PVA A-component, while the relative intensity function,
q m ( t ) = q A ( X 1 , Y 1 , t ) + | q A ( X 2 , Y 2 , t ) | 2 max X , Y q A ( X , Y , 0 ) ,
(26)
describes the time evolution of the A-component amplitude in the dipole normalized by its initial intensity.

Figure 6 shows the evolution of the meridional position of the dipole center, yc, for the case of strong dipoles. Decaying oscillations can clearly be seen. It should be noted that these oscillations persist for all cases studied—see the supplementary material for the case of weak dipoles.

FIG. 6.

Comparative evolution of meridional position of the dipole center, yc, for an initially tilted LRD with β = c: blue curve corresponds to α 0 = 5 °, red curve—to α 0 = 30 °, magenta curve—to α 0 = 60 °, and green curve—to α 0 = 90 °. Solutions were obtained for 0 t 400.

FIG. 6.

Comparative evolution of meridional position of the dipole center, yc, for an initially tilted LRD with β = c: blue curve corresponds to α 0 = 5 °, red curve—to α 0 = 30 °, magenta curve—to α 0 = 60 °, and green curve—to α 0 = 90 °. Solutions were obtained for 0 t 400.

Close modal

In the case of strong dipoles with α 0 = 5 °, the modeled dipole trajectory is compared with the kinematic theory trajectory (Fig. 7), where we used λ = 0.2 and normalized Δ X D sin α and yc with D 0 sin α and Ym [as defined in Eq. (17)], respectively. The similarity between numerical and predicted trajectories suggests that the theory is valid within the strong dipole regime for small initial tilts as the vortex structure remains approximately unchanged throughout the evolution. However, significant deviations in the numerical trajectory in Fig. 7 become visible when t 300, suggesting that a transition in dynamics results in a discrepancy with the kinematic theory to be investigated further.

FIG. 7.

( Δ X , y c )-spiral for ( α 0 , β ) = ( 5 ° , c ) and over the domain 0 t 400: the thick black spiral corresponds to the kinematic theory for strong dipoles, while the dotted red spiral represents our numerical modeling results.

FIG. 7.

( Δ X , y c )-spiral for ( α 0 , β ) = ( 5 ° , c ) and over the domain 0 t 400: the thick black spiral corresponds to the kinematic theory for strong dipoles, while the dotted red spiral represents our numerical modeling results.

Close modal

By t = 400, the vortex has undergone a small deceleration and an increased partner separation of V / c 1 3 % and D / D 0 1 1.1 %, respectively. These alterations are consistent with the dynamics of symmetric eastward drifting LRDs studied in Ref. 30, which continued translating as a steady state, with deceleration and small partner separation resulting from implicit numerical viscosity.

To better understand the spiral deviations in Fig. 7, we consider times t corresponding to intersections of the dipole trajectory with the zonal axis [i.e., y c ( t ) = 0, with = 1 , 2 , 3 ]. Eight intersections are observed in Fig. 6 for α 0 = 5 ° , 30 ° , 60 ° over the time of integration 0 t 400 and seven intersections for α 0 = 90 °. The simple tilted mode (or T-mode) structure of the A-component at initialization [see Fig. 1(a)] is assumed to remain at each intersection according to the kinematic theory, while the amplitude of the A-component, q A sin α, is expected to decay at subsequent crossings together with the angle α ( t ), similar to | Δ X ( t ) | / D sin α 0 in Fig. 7. To evaluate the predictions of the kinematic dipole theory for strong dipoles, we compare the predicted amplitude decay with the values of q m ( t ) in our solutions for different values of α0.

Values of q m ( t ) are plotted in Fig. 8 for different initial angles and correspond well to the theoretical T-mode prediction of α ( t ) for the earlier crossings. However, a significant increase in the value of qm can be seen in Fig. 8 at = 7 , 8 (t > 300) for α 0 = 5 ° and α 0 = 30 °. Such increases indicate the dominance of a growing D-mode (compare with Fig. 2) that becomes clearly visible at the final two crossings (see the supplementary material and Fig. 9). This earlier D-mode appearance in comparison with strong non-tilted eastward LRDs in Ref. 30 indicates that the dipole perturbations due to tilt play a catalytic role and favor the instability, which appears only in weak eastward LRDs over the time of integration.

FIG. 8.

Values of qm at zonal axis crossings, , for β = c and α 0 = 5 ° (blue), α 0 = 30 ° (red), α 0 = 60 ° (magenta), and α 0 = 90 ° (green). Thick black line corresponds to theoretical decay when α 0 = 5 °, as described by Δ X / D 0 sin α 0, where Δ X D sin α. Note that = 0 corresponds to initialization.

FIG. 8.

Values of qm at zonal axis crossings, , for β = c and α 0 = 5 ° (blue), α 0 = 30 ° (red), α 0 = 60 ° (magenta), and α 0 = 90 ° (green). Thick black line corresponds to theoretical decay when α 0 = 5 °, as described by Δ X / D 0 sin α 0, where Δ X D sin α. Note that = 0 corresponds to initialization.

Close modal
FIG. 9.

qA snapshots for ( α 0 , β ) = ( 30 ° , c ) at times where yc = 0. (a) t 5 240, (b) t 6 290, (c) t 7 340, and (d) t 8 390.

FIG. 9.

qA snapshots for ( α 0 , β ) = ( 30 ° , c ) at times where yc = 0. (a) t 5 240, (b) t 6 290, (c) t 7 340, and (d) t 8 390.

Close modal
Moreover, the growth rate between the final zonal axis crossings is
σ = ( q m ( t 8 ) q m ( t 7 ) ) / q m ( t 7 ) t 8 t 7 0.006 ,
(27)
for α 0 = 5 °, which is smaller than D-mode growth for weak eastward LRDs.30 Note that the increase in qm between the last crossings is weaker, σ 0.003, for α 0 = 30 °, despite there being a greater partner separation of 6% and dipole deceleration of 20% at t = 400 [see Fig. 4(a)]. This increased distance between dipole partners contributes to slower D-mode development as seen also for α 0 = 45 ° [Fig. 4(b)].

For an initial launch of α 0 = 60 °, the growing D-mode is not seen; instead, we observed a monotonic and approximately linear decrease in q m ( t ) at consecutive values of (see Fig. 8) and find that the A-component evolves as a decaying T-mode of alternating sign (see Fig. 10). However, at later times, one can see meridional splitting of the T-mode [Fig. 10(d)], reflecting another mechanism to cause the separation between partners, though the apparent partner separation appears less than that caused by the growing D-mode for smaller tilt α 0 = 45 ° [see Fig. 4(b)]. Such adjustment of the dipole to an elliptical (meridionally elongated) shape is accompanied by a decrease in the deceleration magnitude for t > 250. Further investigations in substantially longer domains, over much longer times, and with much higher resolution (to suppress implicit viscosity), go beyond the scope of this paper.

FIG. 10.

qA snapshots for ( α 0 , β ) = ( 60 ° , c ) at times where yc = 0. (a) t 5 245, (b) t 6 285, (c) t 7 335, and (d) t 8 380.

FIG. 10.

qA snapshots for ( α 0 , β ) = ( 60 ° , c ) at times where yc = 0. (a) t 5 245, (b) t 6 285, (c) t 7 335, and (d) t 8 380.

Close modal

Similar T-mode persistence with alternating sign value is observed in the dynamics when α 0 = 90 ° (see the supplementary material). However, despite approximately linear depreciation in qm up until = 5, this decay becomes significantly weaker at = 6 , 7 (see Fig. 8). Even though this is still a monotonically decreasing pattern, the change in decay suggests that there might be some transition in dynamics to be investigated further.

In summary, our analysis for strong dipoles with β = c revealed two distinct scenarios, which are as follows:

  1. The development of a slowly growing D-mode when initial tilt is small, i.e., 0 < α 0 45 °,

  2. The gradual splitting of a decaying T-mode (without the appearance of a D-mode) when the initial tilt is large, i.e., 45 < α 0 90.

Further investigations are needed over a much larger interval of time to clarify the fate of D- and T-modes in strong dipoles.

Our numerical simulations indicate that weak dipoles released at an angle α 0 [ 5 ° , 90 ° ] initially follow oscillatory paths that decay over time (see the supplementary material), comparable to the case of strong dipoles. In Ref. 31, weak dipoles only disintegrate at large α0 (corresponding to α 0 = 90 ° when β = 0.2 and α 0 = 60 ° , 90 ° when β = 0.35), otherwise orienting themselves along the zonal axis, displaying approximate steady translation in the eastward direction. Given the destruction of purely eastward weak dipoles was attributed to the growth of the D-mode in Ref. 30, we anticipate this anomaly to develop on an initially tilted dipole after the adjustment phase.

For weak dipoles with ( α 0 , β ) = ( [ 5 ° , 30 ° , 45 ° ] , 2 c ), we find a phase transition in qA at intersections with the zonal axis (yc = 0). In particular, when α 0 = 5 °, the T-mode at zonal axis crossings rapidly evolves into a growing oscillatory D-mode (see Fig. 11). Consequently, the eddy partners are driven much further apart than for strong dipoles, with D / D 0 1 6.3 % , 36.6 % at t = 200, 400, respectively, [compare Figs. 4(a) and 4(c)]. This is also reflected in their deceleration, where V / c 1 21.6 % , 89.2 % when t = 200, 400, which is consistent with the dynamics of eastward LRDs with the growing D-mode.30 

FIG. 11.

qA snapshots for ( α 0 , β ) = ( 5 ° , 2 c ) at times where yc = 0. (a) t 2 60, (b) t 4 120, (c) t 7 180, and (d) t 15 240.

FIG. 11.

qA snapshots for ( α 0 , β ) = ( 5 ° , 2 c ) at times where yc = 0. (a) t 2 60, (b) t 4 120, (c) t 7 180, and (d) t 15 240.

Close modal

Observed destruction of tilted weak dipoles is consistent with previous numerical simulations using a barotropic model ( L R d) (Ref. 14, Sec. 4.2.2). The similar weak-intensity dipole ( β / c = 2), launched at α 0 = 5 °, relaxed along the zonal axis and proceeded to disintegrate. This disintegration was attributed to the filamentation process and the emission of Rossby waves, suppressed by the use of a cutting filter in the far field, which encouraged drastic deceleration and allowed the separation distance between partners to grow.

In our simulations, similar results are obtained when α 0 = 30 , 45 ° (see Fig. 12), though with longer T-mode persistence before D-mode formation. The dipole partners experience increased separation and greater dipole deceleration [compare Figs. 4(c) and 4(d)], while the zonal dipole drift becomes westward at t S 400. Therefore, even though such opposite-sign couples remain self-advecting eastward until t = 400, the amplifying D-mode will inevitably disintegrate the dipole in longer simulations.

FIG. 12.

qA snapshots for ( α 0 , β ) = ( 45 ° , 2 c ) at times where yc = 0. (a) t 2 60, (b) t 4 130, (c) t 6 200, and (d) t 8 270.

FIG. 12.

qA snapshots for ( α 0 , β ) = ( 45 ° , 2 c ) at times where yc = 0. (a) t 2 60, (b) t 4 130, (c) t 6 200, and (d) t 8 270.

Close modal

For even weaker dipoles with β = 0.35, a growing D-mode emerged for α 0 = 5 ° , 30 °, which results in the subsequent disintegration of the dipole for t < 200 (see the supplementary material and Fig. 3). Such disintegration is consistent with Ref. 30; however, our results throughout Sec. III D differ from those obtained in Ref. 31, where adjustment to a seemingly steady eastward propagating state was observed for 0 t 400. We believe that these discrepancies are due to the 16 times coarser numerical resolution adopted in their 30 years old study, which meant they were unable to capture the instability of weak eastward propagating dipoles.

To summarize our findings in this subsection, we have observed the following:

  1. Weak dipoles developed a D-mode over time when 0 < α 0 < 45 ° (and for α 0 = 45 ° when β = 2 c ).

  2. A growing D-mode decelerates the weak dipoles and increases the core distance between the opposite-sign pair, encouraging the spontaneous symmetry breaking phenomena, as found for symmetric weak eastward LRDs.30 

  3. β = 3.5 c dipoles completely disintegrate within the time interval of 0 t 400, whereas β = 2 c dipoles begin to propagate westward much later and would completely disintegrate if simulated for longer time.

Weak dipoles with ( α 0 , β ) = ( 60 ° , 2 c ) have qA evolve as a decaying T-mode (see Fig. 13) comparable to what we observed with strong dipoles (see Fig. 10). However, a notable difference is that D / D 0 1 62.2 % at t = 200, corresponding to a much greater separation than observed for β = c. This deformation is similar to the meridional separatrix stretching captured in simulations for dipoles with initially circular separatrix in Ref. 14. Despite the large meridional separation at early times, this elliptical deformation does not change significantly in the range of 200 t 400 (see Table II), and these small alterations are likely due to implicit numerical viscosity.41 Furthermore, the much lower deceleration experienced compared to cases with smaller values of α0 suggests that the dipole transforms to a non-circular steady state. This agrees with the steady adjustment discussed in Ref. 31 and the theoretical predictions for strong dipoles in Ref. 36, though these studies did not discuss elliptical deformation.

FIG. 13.

qA snapshots for ( α 0 , β ) = ( 60 ° , 2 c ) at times where yc = 0. (a) t 2 60, (b) t 5 170, (c) t 8 280, and (d) t 11 390.

FIG. 13.

qA snapshots for ( α 0 , β ) = ( 60 ° , 2 c ) at times where yc = 0. (a) t 2 60, (b) t 5 170, (c) t 8 280, and (d) t 11 390.

Close modal
TABLE II.

Percentage change in partner separation and dipole deceleration for β = 2 c, with first entries corresponding to t = 200 and second entries corresponding to t = 400. When α 0 = 90 °, the dipole disintegrates into the background completely at t 280.

α0 D ( t ) / D 0 1 V ( t ) / c 1
5 °  (6.3, 36.6)%  −(21.6, 89.2)% 
30 °  (18.1, 37.8)%  −(60.5, 100.1)% 
45 °  (14.9, 39.8)%  −(52.4, 102.4)% 
60 °  (62.2, 65.3)%  −(58.8, 67.0)% 
90 °  (50.1, 300)%  −(−99.1, 57.4)% 
α0 D ( t ) / D 0 1 V ( t ) / c 1
5 °  (6.3, 36.6)%  −(21.6, 89.2)% 
30 °  (18.1, 37.8)%  −(60.5, 100.1)% 
45 °  (14.9, 39.8)%  −(52.4, 102.4)% 
60 °  (62.2, 65.3)%  −(58.8, 67.0)% 
90 °  (50.1, 300)%  −(−99.1, 57.4)% 

Weak dipoles with ( α 0 , β ) = ( 90 ° , 2 c ) disintegrate in the interval of 0 t 280 as observed in Ref. 31. Indeed, a decaying T-mode similar to that seen when α 0 = 60 ° elucidates the dipole dynamics for this parameter regime where the values of D and V rapidly increase over time (see the supplementary material). This destruction occurs much earlier than that achieved by a D-mode and is due to significant Rossby wave radiation. More specifically, the trajectory of the dipole exposes the vortex to interaction with the trailing Rossby waves, which enhances the rapid disintegration.

Finally, when α 0 45 ° for the weaker β = 3.5 c dipole, a D-mode does not appear over the considered time interval, and instead, we see components of decaying T-mode decelerate and drift apart, before the disintegration of the dipole. Since a D-mode appears when α 0 = 45 ° and β = 2 c, this suggests the existence of some bifurcation point in ( α 0 , β )-space, where there is a transition between D-mode and T-mode dominance. To confirm this, more research is needed. Finally, in the extreme case when α 0 = 90 °, weak dipoles with β = 3.5 c disintegrate much faster than for other initial tilts, as in Refs. 31 and 34.

To summarize, we obtained the following results:

  1. qA evolves as a decaying T-mode rather than a growing D-mode for weak dipoles when α 0 > 45 ° (and for α 0 = 45 ° when β = 3.5 c),

  2. Weak dipoles with β = 3.5 c rapidly disintegrate when α 0 45 °, with the dipole lifespan shortening as α0 increases,

  3. When β = 2 c and α 0 = 90 °, the dipole disintegrates much faster than with a growing D-mode, as a consequence of Rossby wave radiation,

  4. When β = 2 c and 45 ° < α 0 < 90 °, elliptic deformation of the separatrix and large enstrophy loss inhibit D-mode development and drive the adjustment to approximate steady propagation along the zonal axis.

Motivated by the ubiquity of isolated coherent vortices in geophysical (i.e., rotating and stratified) fluids, we investigated the dynamics of mesoscale vortex dipoles launched with different north-eastward tilts to the zonal direction. We considered an idealized equivalent-barotropic QG model in an oceanic configuration and numerically simulated the evolution of individual LRDs19 for a physically relevant range of parameters, focusing both on transient and long-time behaviors. While motivated by the important work of our forerunner (Ref. 31), we were able to extend their analysis using more advanced numerical schemes and techniques. These tools allowed us to reexamine their results and discover new dynamical behaviors, thus enriching our understanding of the long-term behavior of coherent vortices. In particular, we found dynamical sensitivities to both the initial-tilt angle, α0, and the initial intensity of the dipole, β / c.

Even strong dipoles (i.e., β = c; here, nondimensional planetary vorticity gradient and dipole speed, respectively) for moderate initial tilts (i.e., α 0 45 °) eventually developed a critical D-mode instability,30 resulting in essential deviations from the predictions of the kinematic strong-dipole theory in Ref. 35. This theory allows only for an alternating sign T-mode with the initial profile of qA, as depicted in Fig. 1(a), and decaying with time as in Figs. 7 and 8 with damping suggested in Ref. 36. For large tilt (i.e., α 0 > 45 °), our results are more aligned with the kinematic theory, while displaying meridional separation of dipole partners (illustrated by a splitting of the T-mode in Fig. 10) that accumulates over the adjustment period and prevents the above instability so that the D-mode did not develop for large tilt.

In the case of moderate initial tilt, weak dipoles experienced either a combination of significant partner separation with overall deceleration or complete destruction owing to fast developing of the D-mode. In the case of large initial tilts, more pronounced elliptical deformation of the dipole core—corresponding to meridional splitting of T-mode—was found to dominate over the D-mode, as the dipole underwent oscillations along the zonal axis until adjusting to steady-state eastward propagation. For extreme large tilts (e.g., α 0 = 90 °) and/or very weak dipoles ( β = 3.5 c), these structures disintegrated rapidly by radiating Rossby waves and proceeding to interact with them; leaving no time for the D-mode destruction mechanism.

Our numerical simulations employed very small implicit numerical viscosity in the background; however, the work presented in Ref. 41 showed that with explicit Newtonian viscosity (and β = 0), there exists an elliptical dipole structure that distinct circular dipole initializations converge toward as time progresses. These states are found to be unsteady, as characterized by decaying amplitude and increasing vortex size, which is broadly similar to our findings for large α0 and suggests that this adjustment would be stable in the inviscid limit (if implicit numerical viscosity was exactly zero).

In summary, we found and analyzed different dipole parameter regimes and reported evolution scenarios never previously discussed. All these cases may exist and co-exist in nature, but they are likely to be influenced by other unaccounted physical processes, such as realistic large-scale circulation, stratification, and topographic effects. Therefore, our results should act as a catalyst for further research.

See the supplementary material for animations of the dipole propagation and additional figures that support the results of this study.

J.D. acknowledges support from the Roth scholarship, Imperial College London. G.S. acknowledges support from the U.S. National Science Foundation (No. NSF-OCE-1828843). P.B. was supported by the NERC Grant No. NE/T002220/1, the Leverhulme Trust Grant No. RPG-2019-024, and by the Moscow Centre for Fundamental and Applied Mathematics (supported by Agreement No. 075-15-2019-1624 with the Ministry of Education and Science of the Russian Federation).

The authors have no conflicts to disclose.

Jack Davies: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Visualization (lead); Writing – original draft (lead). Georgi Sutyrin: Conceptualization (equal); Methodology (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal). Matthew Crowe: Data curation (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Pavel Berloff: Conceptualization (equal); Investigation (equal); Methodology (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal).

The data that support the findings of this work are available from the corresponding author upon reasonable request.

1.
V.
Kamenkovich
,
M.
Koshlyakov
, and
A.
Monin
,
Synoptic Eddies in the Ocean
(
D. Reidel Publishing Company
,
Dordrecht
,
1986
), Vol.
62
, p.
433
.
2.
J. C.
McWilliams
and
J. B.
Weiss
, “
Anisotropic geophysical vortices
,”
Chaos
4
,
305
311
(
1994
).
3.
J. C.
McWilliams
, “
The emergence of isolated coherent vortices in turbulent flow
,”
J. Fluid Mech.
146
,
21
43
(
1984
).
4.
A. P.
Ingersoll
, “
Atmospheric dynamics of the outer planets
,”
Science
248
,
308
315
(
1990
).
5.
L.
Bengtsson
and
J.
Lighthill
, “
Intense atmospheric vortices
,” in
Proceedings of the Joint Symposium (IUTAM/IUGG), Reading, UK, July 14–17, 1981
(
Springer Science & Business Media
,
2013
).
6.
T.
Dowling
and
E.
Spiegel
, “
Stellar and Jovian vortices
,”
Ann. N. Y. Acad. Sci.
617
,
190
(
1990
).
7.
Q.
Ni
,
X.
Zhai
,
G.
Wang
, and
C. W.
Hughes
, “
Widespread mesoscale dipoles in the global ocean
,”
J. Geophys. Res.
125
,
e2020JC016479
, https://doi.org/10.1029/2020JC016479 (
2020
).
8.
C. W.
Hughes
and
P. I.
Miller
, “
Rapid water transport by long-lasting modon eddy pairs in the southern midlatitude oceans
,”
Geophys. Res. Lett.
44
,
12-375
12-384
, https://doi.org/10.1002/2017GL075198 (
2017
).
9.
M.
Jalali
and
D.
Dritschel
, “
Stability and evolution of two opposite-signed quasi-geostrophic shallow-water vortex patches
,”
Geophys. Astrophys. Fluid Dyn.
114
,
561
587
(
2020
).
10.
M.
Rostami
and
V.
Zeitlin
, “
Eastward-moving equatorial modons in moist-convective shallow-water models
,”
Geophys. Astrophys. Fluid Dyn.
115
,
345
367
(
2021
).
11.
E. R.
Johnson
and
M. N.
Crowe
, “
Oceanic dipoles in a surface quasi-geostrophic model
,”
J. Fluid Mech.
958
,
R2
(
2023
).
12.
H.
Lamb
,
Hydrodynamics
, 2nd ed. (
Cambridge University Press
,
1895
).
13.
V.
Meleshko
and
G.
Van Heijst
, “
On Chaplygin's investigations of two-dimensional vortex structures in an inviscid fluid
,”
J. Fluid Mech.
272
,
157
182
(
1994
).
14.
R.
Khvoles
,
D.
Berson
, and
Z.
Kizner
, “
The structure and evolution of elliptical barotropic modons
,”
J. Fluid Mech.
530
,
1
30
(
2005
).
15.
P.
Luzzatto-Fegiz
and
C.
Williamson
, “
Determining the stability of steady two-dimensional flows through imperfect velocity-impulse diagrams
,”
J. Fluid Mech.
706
,
323
350
(
2012
).
16.
V.
Brion
,
D.
Sipp
, and
L.
Jacquin
, “
Linear dynamics of the Lamb-Chaplygin dipole in the two-dimensional limit
,”
Phys. Fluids
26
,
064103
(
2014
).
17.
M. E.
Stern
, “
Minimal properties of planetary eddies
,”
J. Mar. Res.
33
,
1
13
(
1975
).
18.
Z.
Kizner
and
D.
Berson
, “
Emergence of modons from collapsing vortex structures on the beta-plane
,”
J. Mar. Res.
58
,
375
403
(
2000
).
19.
V.
Larichev
and
G. M.
Reznik
, “
On two-dimensional solitary Rossby waves
,”
Dokl. Akad. Nauk
231
,
1077
1079
(
1976
).
20.
G.
Reznik
, “
Dynamics of localized vortices on the beta plane
,”
Izv., Atmos. Oceanic Phys.
46
,
784
797
(
2010
).
21.
E. R.
Johnson
and
M. N.
Crowe
, “
The decay of a dipolar vortex in a weakly dispersive environment
,”
J Fluid Mech.
917
,
A35
(
2021
).
22.
J. C.
McWilliams
, “
An application of equivalent modons to atmospheric blocking
,”
Dyn. Atmos. Oceans
5
,
43
66
(
1980
).
23.
N.
Butchart
,
K.
Haines
, and
J.
Marshall
, “
A theoretical and diagnostic study of solitary waves and atmospheric blocking
,”
J. Atmos. Sci.
46
,
2063
2078
(
1989
).
24.
J.
Nycander
, “
Refutation of stability proofs for dipole vortices
,”
Phys. Fluids A
4
,
467
476
(
1992
).
25.
G. E.
Swaters
and
G. R.
Flierl
, “
Ekman dissipation of a barotropic modon
,” in
Elsevier Oceanography Series
(
Elsevier
,
1989
), Vol.
50
, pp.
149
165
.
26.
J. C.
McWilliams
,
G. R.
Flierl
,
V. D.
Larichev
, and
G. M.
Reznik
, “
Numerical studies of barotropic modons
,”
Dyn. Atmos. Oceans
5
,
219
238
(
1981
).
27.
G. F.
Carnevale
,
G. K.
Vallis
,
R.
Purini
, and
M.
Briscolini
, “
Propagation of barotropic modons over topography
,”
Geophys. Astrophys. Fluid Dyn.
41
,
45
101
(
1988
).
28.
S.
Muzylev
and
G.
Reznik
, “
On proofs of stability of drift vortices in magnetized plasmas and rotating fluids
,”
Phys. Fluids B
4
,
2841
2844
(
1992
).
29.
G. E.
Swaters
, “
Spectral properties in modon stability theory
,”
Stud. Appl. Math.
112
,
235
258
(
2004
).
30.
J.
Davies
,
G. G.
Sutyrin
, and
P. S.
Berloff
, “
On the spontaneous symmetry breaking of eastward propagating dipoles
,”
Phys. Fluids
35
,
041707
(
2023
).
31.
J. S.
Hesthaven
,
J.-P.
Lynov
, and
J.
Nycander
, “
Dynamics of nonstationary dipole vortices
,”
Phys. Fluids A
5
,
622
629
(
1993
).
32.
J.
Pedlosky
et al,
Geophysical Fluid Dynamics
(
Springer
,
1987
), Vol.
710
.
33.
G. K.
Vallis
,
Atmospheric and Oceanic Fluid Dynamics
(
Cambridge University Press
,
2017
).
34.
G. G.
Sutyrin
,
J. S.
Hesthaven
,
J.-P.
Lynov
, and
J. J.
Rasmussen
, “
Dynamical properties of vortical structures on the beta-plane
,”
J. Fluid Mech.
268
,
103
131
(
1994
).
35.
J.
Nycander
and
M.
Isichenko
, “
Motion of dipole vortices in a weakly inhomogeneous medium and related convective transport
,”
Phys. Fluids B
2
,
2042
2047
(
1990
).
36.
J.
Nycander
, “
Dynamics of dipole vortices
,”
Verh. K. Akad. Wet., Afd. Natuurkd. Eerste Sect.
43
,
177
186
(
1994
).
37.
S.
Karabasov
,
P. S.
Berloff
, and
V.
Goloviznin
, “
Cabaret in the ocean gyres
,”
Ocean Modell.
30
,
155
168
(
2009
).
38.
I.
Shevchenko
and
P.
Berloff
, “
Multi-layer quasi-geostrophic ocean dynamics in eddy-resolving regimes
,”
Ocean Modell.
94
,
1
14
(
2015
).
39.
I.
Shevchenko
and
P.
Berloff
, “
On the roles of baroclinic modes in eddy-resolving midlatitude ocean dynamics
,”
Ocean Modell.
111
,
55
65
(
2017
).
40.
K. J.
Burns
,
G. M.
Vasil
,
J. S.
Oishi
,
D.
Lecoanet
, and
B. P.
Brown
, “
Dedalus: A flexible framework for numerical simulations with spectral methods
,”
Phys. Rev. Res.
2
,
023068
(
2020
).
41.
Z.
Kizner
,
R.
Khvoles
, and
D. A.
Kessler
, “
Viscous selection of an elliptical dipole
,”
J. Fluid Mech.
658
,
492
508
(
2010
).
Published open access through an agreement withJISC Collections

Supplementary Material