Mathematical Modeling of Vortex Interaction Using a Three-Layer Quasigeostrophic Model. Part 1: Point-Vortex Approach

: The theory of point vortices is used to explain the interaction of a surface vortex with subsurface vortices in the framework of a three-layer quasigeostrophic model. Theory and numerical experiments are used to calculate the interaction between one surface and one subsurface vortex. Then, the conﬁguration with one surface vortex and two subsurface vortices of equal and opposite vorticities (a subsurface vortex dipole) is considered. Numerical experiments show that the self-propelling dipole can either be captured by the surface vortex, move in its vicinity, or ﬁnally be completely ejected on an unbounded trajectory. Asymmetric dipoles make loop-like motions and remain in the vicinity of the surface vortex. This model can help interpret the motions of Lagrangian ﬂoats at various depths in the ocean.

In this paper, we use a quasigeostrophic model to study the features of the interaction between one vortex of the upper layer and one/two vortices of the middle layer of a three-layer rotating fluid. The vortex of the upper layer is a prototype of a surface ocean vortex (many such vortices occupy the upper 500-600 m of 4000 m deep oceans). The vortices of the middle layer represent intrathermocline lens vortices, observed in all oceans, but especially common in the Northeastern Atlantic at depths of 600-1600 m [64,65]. Note that for the first time the idea of modeling intrathermocline lenses with point vortices of the intermediate layer was proposed by Hogg and Stommel in [37]. The lower layer will contain no vortex and will correspond to the deep ocean.
In [66,67], the SEMANE and MEDTOP cruise data were analyzed to study the interactions of intrathermocline vortices with a cyclonic surface vortex; however, due to intermittent data collection at sea, only a few snapshots of these interactions were obtained. This circumstance is the motivation for N l , N n l, n = 1 n l i f m = i 3 j=2 q ij s jm κ n m κ l i ln r nl mi + K 0 γ j−1 r nl mi (1) with canonical variables g n m = x n m and p n m = κ n m y n m /2. Here, K 0 is a modified Bessel function of zero order; γ 1 , γ 2 are the parameters inversely proportional to the first and second deformation radii, respectively; r nl mi = In addition to the Hamiltonian (1), the system has integral invariants: i.e., the linear impulses P x , P y and the angular momentum M. In Cartesian coordinates, the equations of motion of the vortices have the form

Numerical Modeling of Vortex Interaction
In all numerical experiments, which solve Equations (2) and (3), we will assume that one cyclonic vortex is located in the upper layer (N 1 = 1, κ 1 1 > 0) and either one anticyclonic lens can be located in the middle layer (N 2 = 1, κ 1 2 < 0), or two vortices of opposite signs N 2 = 2, κ 1 2 < 0, κ 2 2 >0 . We also introduce the notation for the effective intensity of a point vortex Λ

Cyclonic Surfer Vortex and Anticyclonic Intrathermocline Lenses
Let us consider now Λ 1 1 = −Λ 1 2 ; this is an analogue of a heton (a two-layer dipolar structure) [33]. We assume that the upper and middle layer vortices lie initially on the x axis and are separated by a distance L. We obtain from (3) .
i.e., (1) at L = 0, the structure will remain in place; (2) at L > 0, the translational velocity of this two-layer vortex is a non-monotonic function of L, and its maximum is achieved at L ≈ 1.75·Rd 1 (green line in Figure 1). If Λ 1 1 −Λ 1 2 , then the vortices always rotate with angular velocity relative to the center of vorticity with coordinates Figure 1 demonstrates the numerically obtained function V(L) for pointwise vortices in both the upper and middle layers (4) of the three-layer model.
We assume that the upper and middle layer vortices lie initially on the x axis and are separated by a distance L. We obtain from (3) i.e., (1) at = 0, the structure will remain in place; (2) at > 0, the translational velocity of this two-layer vortex is a non-monotonic function of L, and its maximum is achieved at ≈ 1.75 • 1 (green line in Figure 1). If 1 1 ≠ − 2 1 , then the vortices always rotate with angular velocity relative to the center of vorticity with coordinates        (2) and (3) on the same time interval. Obviously, the end points of the straight segments repeat the profile of the green line in Figure 1. Note, that the non-monotonic character of V(L) expresses the fundamental difference between baroclinic and barotropic dipoles, which have a singularity in velocity at L = 0. Panels (b,c) show the initial trajectories of vortices at the same initial location but with a stronger cyclonic vortex in the upper layer. When the ratio Λ 1 2 /Λ 1 1 is changed, the position of the center of vorticity changes, according to (6).

Collinear Initial Configuration
Let us consider the simplest case when three vortices (one in the upper layer and two in the middle layer) initially form a symmetric collinear structure. The cyclonic vortex with effective intensity Λ 1 1 > 0 is still located in the upper layer; in the middle layer there are two vortices of opposite signs on opposite sides of the surface vortex, with Λ 1 2 = −Λ 2 2 > 0. In such a situation two main interactions can govern the motion of the vortices of the middle layer: (a) an intralayer interaction, leading to uniform and rectilinear motion; (b) an interlayer interaction, leading to counterclockwise rotation under the influence of the surface cyclone. Obviously, the upper layer vortex does not undergo intralayer interaction; the pair of middle layer vortices advects it in the direction of its self-propagation. All these mechanisms act simultaneously and contribute to the observed motion of the vortex structure. Figure 3 shows examples of trajectories of such a three-vortex structure. Here, as before, L is the horizontal distance between the upper layer vortex and one of the middle layer vortices. Thus, Mathematics 2020, 8, 1228 5 of 13 the initial distance between the vortices of the middle layer is 2L. The figure shows a change in regimes as the distance L is increased, from a predominant intralayer interaction to a dominance of interlayer interaction.
interaction, leading to uniform and rectilinear motion; (b) an interlayer interaction, leading to counterclockwise rotation under the influence of the surface cyclone. Obviously, the upper layer vortex does not undergo intralayer interaction; the pair of middle layer vortices advects it in the direction of its self-propagation. All these mechanisms act simultaneously and contribute to the observed motion of the vortex structure. Figure 3 shows examples of trajectories of such a three-vortex structure. Here, as before, is the horizontal distance between the upper layer vortex and one of the middle layer vortices. Thus, the initial distance between the vortices of the middle layer is 2 . The figure shows a change in regimes as the distance is increased, from a predominant intralayer interaction to a dominance of interlayer interaction.  First, we consider the case of intrathermocline dipole, a particular case where the middle layer vortices form a dipolar structure, i.e., Λ 1 2 /Λ 2 2 = −1. Let us start with the simplest case when the vortices of the middle layer form a pair, symmetric with respect to the y axis; we consider different distances (2L) between the vortices in the pair.
In this case, for any value of L, the middle layer vortex motion deviates to the right of a straight trajectory due to the action of the upper layer cyclone. The cyclone itself, initially stationary, begins to move counterclockwise locally due to the action of the middle layer vortex pair. Numerical simulations show that for L ≤ 0.6, the dipole rotates around the surface cyclone and escapes to the left, while the cyclone slows down and stops over time without having performed a complete revolution. For L ≥ 0.7, the surface vortex always traps the middle layer dipole, and the entire vortex structure rotates counterclockwise near the center of the domain for a finite duration. The middle layer cyclone trajectory lies inside the trajectory of the surface cyclone, and the anticyclonic lens rotates along an outer orbit. After this duration, the middle layer dipole moves away, and the surface cyclone stops.  Figure 4 shows a gallery of the vortex trajectories at different time intervals (in the different panels) until the dipole structure leaves the central region (which is defined by R ≤ 20 around the origin). During the stage of bounded motion, the global vortex structure is composed of a two-layer cyclone with a "tilted axis" and a peripheral anticyclonic lens, all performing a counterclockwise rotation.  Figure 5a shows that the period of rotation of the vortex structure inside the bounded region monotonically increases with , while the number of full revolutions is an irregular function of . Note that this effect, known as "chattering" [9], is observed in other problems of vortex dynamics (for example, [68]). Two examples with surprisingly large values of are shown in panels (c) and (f) of Figure 4.   Figure 5a shows that the period of rotation T of the vortex structure inside the bounded region monotonically increases with L, while the number of full revolutions N is an irregular function of L. Note that this effect, known as "chattering" [9], is observed in other problems of vortex dynamics (for example, [68]). Two examples with surprisingly large values of N are shown in panels (c) and (f) of Figure 4.
The transition between bounded and unbounded trajectories occurs when two velocity contributions become similar in magnitude: the first is the vortex velocity due to intralayer interaction, proportional to the singular function 1/L, and the second is the velocity of vortices due to interlayer interaction, proportional to the regular non-monotonic function (1/L − K 1 (L)). During the evolution of the vortex structure, the ratio between these contributions constantly changes; the prevalence of the first one leads to a transition towards unbounded trajectories, whereas the prevalence of the second one leads to bounded trajectories. The behavior of the vortices during such evolutions is illustrated in Figure 4b, which shows the time changes in the y-coordinates of the cyclones of the middle and upper layers at L = 4.8, when N = 8: initially the unbounded motion becomes bounded when y 1 2 decreases below y 1 1 ; then the amplitude of y 1 1 in the first half of the finite cycle initially increases to a maximum value, and y 1 2 decreases to a minimum; in the second half of the cycle, the changes in these functions are reversed; finally, y 1 2 reaches a limit value at which the subsurface vortex breaks out of the closed region, and the motion becomes unbounded again.  Figure 5a shows that the period of rotation of the vortex structure inside the bounded region monotonically increases with , while the number of full revolutions is an irregular function of . Note that this effect, known as "chattering" [9], is observed in other problems of vortex dynamics (for example, [68]   The transition between bounded and unbounded trajectories occurs when two velocity contributions become similar in magnitude: the first is the vortex velocity due to intralayer interaction, proportional to the singular function 1⁄ , and the second is the velocity of vortices due For large values of N (for example, when L = 4.3 and N = 241, Figure 4f), there are several such subcycles, but at the end of each of them (except of the last one), the amplitude of y 1 2 has not yet reached its limit value.
Another interesting property in the case of large N is that before reaching the unbounded regime, the annular regions filled by the trajectories of the surface cyclone and middle layer anticyclone expand significantly (as shown by the red and blue rings in Figure 4c,f).
Thus, when an intrathermocline pair runs into a surface cyclone, starting from a certain value of L (here, at B = −16, this value is L ≈ 1.7) a temporary trapping occurs, followed by an expulsion of the dipole. Note that stationary localized regimes, observed for initially collinear vortices, as in Figure 3c,d, do not occur when a middle layer dipole runs into a surface cyclone.
Next, we consider the case of asymmetric middle layer vortices. Up to now, it has been assumed everywhere that the two middle layer vortices have the same effective intensity. However, observations [66,67] show that cyclones at intermediate depths are generally weaker than anticyclones. Such an asymmetry in the distribution of potential vorticity is now studied: we assume that The interaction with the surface cyclone changes dramatically since the weak subsurface cyclone now rotates around the stronger anticyclone in the middle layer. In the general case, this weak subsurface cyclone follows a cycloidal trajectory both around the surface cyclone and around its anticyclonic partner. An example of this is shown in Figure 6. The trajectory of the middle layer cyclone (in green) is much wider than that of the anticyclonic vortex (in blue). that 2 1 2 2 ⁄ = −1 2 ⁄ . The interaction with the surface cyclone changes dramatically since the weak subsurface cyclone now rotates around the stronger anticyclone in the middle layer. In the general case, this weak subsurface cyclone follows a cycloidal trajectory both around the surface cyclone and around its anticyclonic partner. An example of this is shown in Figure 6. The trajectory of the middle layer cyclone (in green) is much wider than that of the anticyclonic vortex (in blue). Depending on the value of parameter , there is a countable number of periodic motions of the three vortices along closed trajectories, called "absolute choreographies" [69]. Examples of the first six choreographies are presented in Figure 7. Depending on the value of parameter L, there is a countable number of periodic motions of the three vortices along closed trajectories, called "absolute choreographies" [69]. Examples of the first six choreographies are presented in Figure 7.  Thus, we have a discrete set of values L n , decreasing with n, which corresponds to a family of n-fold symmetric configurations.    Finally, Figure 9 shows how the transition between one-mode to two-mode choreography occurs when the distance changes with a constant step ∆ = ( 1 − 2 ) 5 ⁄ . In order not to confuse the picture when each trajectory fills its annular region, panels (b-e) show only the first four loops (in time), and the number indicates the order of their formation. For → 2 , the loops thicken in the vicinity of the stationary position so that they coincide completely when the limit is reached.  Figure 7 (blue markers), i.e., at Λ 1 2 /Λ 1 1 = 1/4, and the same at Λ 1 2 /Λ 1 1 = 1/8 (brown markers) and Λ 1 2 /Λ 1 1 = 1/2 (green markers). Figure 8 shows that values L n decrease with increasing n and with decreasing of intensity of subsurface vortices.
Finally, Figure 9 shows how the transition between one-mode to two-mode choreography occurs when the distance L changes with a constant step ∆ = (L 1 − L 2 )/5. In order not to confuse the picture when each trajectory fills its annular region, panels (b-e) show only the first four loops (in time), and the number indicates the order of their formation. For L → L 2 , the loops thicken in the vicinity of the stationary position so that they coincide completely when the limit is reached.

Discussion and Conclusions
We recapitulate the main results of this work and provide their oceanographic interpretation.


If the cyclone of the upper layer and the anticyclonic lens of the middle layer are separated by some distance, then such a two-layer vortex can either move forward (when its total effective vorticity is zero) or rotate relative to the center of vorticity (when its total effective vorticity is nonzero). In any case, both vortices can move far enough from the original location (Section 3.1).  Figure 9 additionally shows that, for any L L N , for large times, all three trajectories fill the annular regions, as in Figure 6. Only a countable number of stationary periodic solutions (choreographies) exist at L = L N and Λ 1 2 /Λ 2 2 −1.

Discussion and Conclusions
We recapitulate the main results of this work and provide their oceanographic interpretation.
• If the cyclone of the upper layer and the anticyclonic lens of the middle layer are separated by some distance, then such a two-layer vortex can either move forward (when its total effective vorticity is zero) or rotate relative to the center of vorticity (when its total effective vorticity is nonzero). In any case, both vortices can move far enough from the original location (Section 3.1). Note that loop-like motions of SOFAR floats seeded in lenses are often observed in the ocean (for example, see [70][71][72]). This loop-like float motion can be explained by the position of the float on the periphery of the lens. It can also be explained by the lens describing loops when interacting with a weak cyclonic partner; this is difficult to determine experimentally.
We believe that the results obtained here can be useful in analyzing experimental measurements of surface and subsurface vortices in the ocean.

Conflicts of Interest:
The authors declare that they have no competing interests.