Deformation and destruction of north-eastward drifting dipoles

We systematically study the evolution of Larichev-Reznik dipoles in an equivalent-barotropic quasigeostrophic 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; (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.


I. INTRODUCTION
5][6] Within this broad class of phenomena, we focus on vortex dipoles (i.e., pairs of opposite-sign vortices) and their long-term evolution.][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. 16ilding from this, a 2D stationary solution (so-called modon) on the b-plane (taking into account planetary rotation and sphericity) was obtained because Stern 17 was found to be unstable. 18Furthermore, a more general solution of the classical quasigeostrophic (QG) equivalent-barotropic model (given by rotating top-layer shallowwater 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). 19hen 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 speed 20 to prevent the vortex losing energy to generated waves. 21In 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. 27Proof of their stability, however, remains evasive. 24,28,29ecently, the phenomenon of spontaneous symmetry breaking was discovered in eastward weak-intensity LRDs leading to their ultimate destruction. 30The corresponding unstable growing perturbation has even symmetry about the zonal axis (similar to the unstable modes obtained for the evolving viscous LCD 16 ) 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 Amode 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.

II. MODEL DESCRIPTION A. Equivalent-barotropic framework
To connect to previous studies, we adopt the equivalentbarotropic quasigeostrophic (QG) model on the beta-plane with flat bottom relief, 32,33 D P where b represents the meridional gradient of the Coriolis parameter; P ¼ Pðx; ŷ; 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 with Rd denoting the internal Rossby radius of deformation and ŵ ¼ ŵðx; ŷ; tÞ denoting the velocity streamfunction.For compactness, we make use of the differential operators, 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 þ xc ð tÞ and ŷ ¼ Ŷ 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 Û , respectively, where the Jacobian operator is defined as for functions A and B, and the absence of hat notation indicates the following nondimensional quantities: Focusing on steady vortices with purely zonal drift yields the following governing equation: where 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 ; where q A denotes the A-component, which has even PVA relative to the zonal axis, and q S 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 w 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 R, where dR ¼ dXdY; and the summations 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 Re À1 r 4 w, where Re ¼ ð L Û Þ À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.

B. LRD initialization
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, where ðr; #Þ are standard polar coordinates.Here, r ¼ 1 defines a circular separatrix ( L is chosen to be the vortex radius); p 2 ¼ b=c þ S > 0; a 0 ¼ 0 ; 180 (eastward and westward propagation, respectively); and k is a positive constant satisfying the nonlinear equation, where J l and K l are the order-l 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 betaplane 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., w !0, as r ! 1 (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: where FðrÞ ¼ ðp=kÞ 2 ðJ 1 ðkrÞ=rJ 1 ðkÞ À 1Þ À 1; r 1 ; ÀK 1 ðprÞ=rK 1 ðpÞ; r > 1; ( and the corresponding PVA components are as follows: where In this work, we investigate tilted dipoles for S ¼ 1, i.e., L ¼ R d , launched at angles 0 < a 0 90 to the zonal axis.As a consequence of this, q A 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 Scomponents 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).
In the (eastward propagating) case a 0 ¼ 0 and for dipoles with b=c ! 2, the component q A (i.e., D-mode) was initially zero but appeared due to spontaneous symmetry breaking and shown to oscillate and grow exponentially over time. 30This 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 a 0 6 ¼ 0 .

C. Kinematics of tilted dipoles
For dipoles, let us consider the coordinates of the positive and negative PVA extrema (X 1 , Y 1 ) and (X 2 , Y 2 ), respectively.The dipole tilt is, then, characterized by where D is the distance between extrema, and a is the evolving angle, such that að0Þ ¼ a 0 .The dipole center is defined as (X c , Y c ), with while the zonal and meridional drift components of the dipole center, respectively, give the propagation velocity.These components allow us to obtain the corresponding dipole speed as For very strong dipoles (c ) b), 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 ¼ D 0 ) and in the dipole center velocity V ¼ c.This theory tells us that the evolution of the dipole angle, aðtÞ, is given by the following equation: where M d ¼ Ð QðrÞr 3 dr ' 3V is the magnitude of the LRD moment, which is described by Eqs. ( 13)-( 15) with b ¼ 0 and normalized by the area of fluid trapped inside the separatrix.Combining Eq. ( 21) with 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 Y m is given in Ref. 35 as 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 In this case, the dipole evolution calculated from Eqs. ( 22)-( 24) displays decaying oscillations, which depend on the damping parameter k.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 highresolution numerical simulations.

III. NUMERICAL EXPERIMENTS A. Methodology
We solve Eq. ( 1) on a doubly periodic computational domain of size ð4L Y ; L Y Þ ¼ ð60; 15Þ with 4N Â N ¼ 8192 Â 2048 nodes.Motivated by mesoscale observations, we consider dipoles on the scale of the internal Rossby radius of deformation, that is, L ¼ Rd (consequently, S ¼ 1).Furthermore, to draw comparisons with Ref. 31, we nondimensionalize so that c ¼ 0.1 (so that Û ¼ 10ĉ, where ĉ is the dimensional equivalent of c) and explore the following parameter space: ða 0 ; bÞ ¼ 5 ; 30 ; 45 ; 60 ; 90 ½ Â 1; 2; 3:5; 4 where a 0 45 is referred to as small tilt, and a 0 > 45 is referred to as large tilt.We also use the terms weak and strong to correspond with c < b 4c and b c, respectively.Benchmark solution animations for this parameter space are shown in Fig. 3 (Multimedia view).
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. 40This 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.

B. Overview of results
The dipole evolution is analyzed by considering both the trajectory of the dipole center (x c , y c ) 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 (X 1 , Y 1 ) and (X 2 , Y 2 ) are evaluated using local 2D-parabolic interpolation.
Comparing with the kinematic theory, we found new behavior in the evolution of dipoles with b=c ¼ 1 (strong dipoles) and with b=c ! 2 (weak dipoles).These modifications to the dipole dynamics are mostly related to the growing D-mode on top of decaying Tmode, 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. 30Here, 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.
Weak dipoles for a 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.
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Þ þ jq A ðX 2 ; Y 2 ; tÞj 2 max X;Y q A ðX; Y; 0Þ ; describes the time evolution of the A-component amplitude in the dipole normalized by its initial intensity.

C. Strong dipole evolution
Figure 6 shows the evolution of the meridional position of the dipole center, y c , 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.
In the case of strong dipoles with a 0 ¼ 5 , the modeled dipole trajectory is compared with the kinematic theory trajectory (Fig. 7), where we used k ¼ 0:2 and normalized DX $ D sin a and y c with D 0 sin a and Y m [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.
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 a 0 ¼ 5 ; 30 ; 60 over the time of integration 0 t 400 and seven intersections for a 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 a, is expected to decay at subsequent crossings together with the angle aðt ' Þ, similar to jDXðt ' Þj=D sin a 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 a 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 aðt ' Þ for the earlier crossings.However, a significant increase in the value of q m can be seen in Fig. 8 at ' ¼ 7; 8 (t > 300) for a 0 ¼ 5 and a 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.Moreover, the growth rate between the final zonal axis crossings is for a 0 ¼ 5 , which is smaller than D-mode growth for weak eastward LRDs. 30Note that the increase in q m between the last crossings is weaker, r $ 0:003, for a 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 a 0 ¼ 45 [Fig.4(b)].
For an initial launch of a 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 a 0 ¼ 45 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.
Similar T-mode persistence with alternating sign value is observed in the dynamics when a 0 ¼ 90 (see the supplementary material).However, despite approximately linear depreciation in q m 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 b ¼ 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 < a 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 < a 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.FIG. 7. ðDX; y c Þ-spiral for ða 0 ; bÞ ¼ ð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.

D. Slow destruction of weak dipoles
Our numerical simulations indicate that weak dipoles released at an angle a 2 ½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 a 0 (corresponding to a 0 ¼ 90 when b ¼ 0:2 and a 0 ¼ 60 ; 90 when b ¼ 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 ða 0 ; bÞ ¼ ð½5 ; 30 ; 45 ; 2cÞ, we find a phase transition in q A at intersections with the zonal axis (y c ¼ 0).In particular, when a 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. 30bserved 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 (b=c ¼ 2), launched at a 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 a 0 ¼ 30; 45 (see Fig. though such opposite-sign couples remain self-advecting eastward until t ¼ 400, the amplifying D-mode will inevitably disintegrate the dipole longer simulations. For even weaker dipoles with b ¼ 0:35, a growing D-mode emerged for a 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 < a 0 < 45 (and for a 0 ¼ 45 when b ¼ 2cÞ.
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. b ¼ 3:5c dipoles completely disintegrate within the time interval of 0 t 400, whereas b ¼ 2c dipoles begin to propagate westward much later and would completely disintegrate if simulated for longer time.

E. Adjustment or fast destruction of weak dipoles
Weak dipoles with ða 0 ; bÞ ¼ ð60 ; 2cÞ have q A 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 b ¼ c.This deformation is similar to the meridional separatrix stretching captured in simulations for dipoles with initially circular separatrix in Ref. 14 meridional separation at early times, this elliptical deformation does not change significantly in the range of t 400 (see Table II), and these small alterations are likely due to implicit numerical viscosity. 41Furthermore, the much lower deceleration experienced compared to cases with smaller values of a 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.
Weak dipoles with ða 0 ; bÞ ¼ ð90 ; 2cÞ disintegrate in the interval of 0 t 280 as observed in Ref. 31.Indeed, a decaying T-mode similar to that seen when a 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 a 0 !45 for the weaker b ¼ 3:5c 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 a 0 ¼ 45 and b ¼ 2c, this suggests the existence of some bifurcation point in ða 0 ; bÞ-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 a 0 ¼ 90 , weak dipoles with b ¼ 3:5c disintegrate much faster than for other initial tilts, as in Refs.31 and 34.
To summarize, we obtained the following results: 1. q A evolves as a decaying T-mode rather than a growing D-mode for weak dipoles when a 0 > 45 (and for a 0

IV. CONCLUSIONS AND DISCUSSION
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 equivalentbarotropic QG model in an oceanic configuration and numerically simulated the evolution of individual LRDs 19 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, a 0 , and the initial intensity of the dipole, b=c.
Even strong dipoles (i.e., b ¼ c; here, nondimensional planetary vorticity gradient and dipole speed, respectively) for moderate initial tilts (i.e., a 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 q A , as depicted in Fig. 1(a), and decaying with time as in 7 and 8 with damping suggested in Ref. 36.For large tilt (i.e., a 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 Tmode-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., a 0 ¼ 90 ) and/or very weak dipoles (b ¼ 3:5c), 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 b ¼ 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 a 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.

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

FIG. 2 .
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. 3 .FIG. 4 .
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 a 0 increasing from top to bottom.The strongest LRD with c=b ¼ 1 is shown in the left column, while other columns show progressively weaker LRDs, with the right-most column corresponding to c=b ¼ 0:25.Multimedia available online.
¼ 45 when b ¼ 3:5c), 2. Weak dipoles with b ¼ 3:5c rapidly disintegrate when a 0 !45 , with the dipole lifespan shortening as a 0 increases, 3. When b ¼ 2c and a 0 ¼ 90 , the dipole disintegrates much faster than with a growing D-mode, as a consequence of Rossby wave radiation, 4. When b ¼ 2c and 45 < a 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.

TABLE I .
Percentage change in partner separation and dipole deceleration for b ¼ c, with first entries corresponding to t ¼ 200 and second entries corresponding to t ¼ 400.FIG. 5. ða 0 ; bÞ-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-modeinduced destruction with t s > 100 (where t s 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.

TABLE II .
Percentage change in partner separation and dipole deceleration for b ¼ 2c, with first entries corresponding to t ¼ 200 and second entries corresponding to t ¼ 400.When a 0 ¼ 90 , the dipole disintegrates into the background completely at t $ 280.