Mathematical Modeling of Vortex Interaction Using a Three-Layer Quasigeostrophic Model. Part 2: Finite-Core-Vortex Approach and Oceanographic Application

: The three-layer version of the contour dynamics / surgery method is used to study the interaction mechanisms of a large-scale surface vortex with a smaller vortex / vortices of the middle layer (prototypes of intrathermocline vortices in the ocean) belonging to the middle layer of a three-layer rotating ﬂuid. The lower layer is assumed to be dynamically passive. The piecewise constant vertical density distribution approximates the average long-term proﬁle for the North Atlantic, where intrathermocline eddies are observed most often at depths of 300–1600 m. Numerical experiments were carried out with di ﬀ erent initial conﬁgurations of vortices, to evaluate several e ﬀ ects. Firstly, the stability of the vortex compound was evaluated. Most often, it remains compact, but when unstable, it can break as vertically coupled vortex dipoles (called hetons). Secondly, we studied the interaction between a vertically tilted cyclone and lenses. Then, the lenses ﬁrst undergo anticlockwise rotation determined by the surface cyclone. The lenses can induce alignment or coupling with cyclonic vorticity above them. Only very weak lenses are destroyed by the shear stress exerted by the surface cyclone. Thirdly, under the inﬂuence of lens dipoles, the surface cyclone can be torn apart. In particular, the shedding of rapidly moving vortex pairs at the surface reﬂects the presence of lens dipoles below. More slowly moving small eddies can also be torn away from the main surface cyclone. In this case, they do not appear to be coupled with middle layer vortices. They are the result of large shear-induced deformation. Common and di ﬀ ering features of the vortex interaction, modeled in the framework of the theory of point and ﬁnite-core vortices, are noted.


Introduction
Vortices are ubiquitous features of many fluid flows (geophysical, industrial, biological flows). They are often long living and they carry substantial energy or chemical tracers across the fluid. They are also the key to the description of developed turbulence [1,2]. Though the flow around a singular circular vortex is laminar and predictable, more than three vortices evolve chaotically, and the flow around them is intermittent. By interacting, deforming and splitting, a few initial vortices can produce many other flow features (filaments, smaller or larger vortices). Thus, a continuous spectrum of energized scales develops in space and time. The interaction between vortices, of different sizes, is an These theorems are central to vortex dynamics, a specific study of which follows. They allow the full description of the flow by the potential vorticity distribution which is advected by the velocity field that it itself generates. After presenting the model equations and numerical method (Section 2), we will study the interaction of a surface vortex with a deeper and smaller vortex (Section 3). Due to vorticity symmetry, we chose the surface vortex to have positive potential vorticity. In Section 4, we study the interaction of the same surface vortex, with one or two pairs of deep and small vortices. Our aim in these various experiments and cases is to determine the fate of the vortices: do they remain aligned vertically, or coherent (in close vicinity), or do they irreversibly separate? The latter evolution indicates that under given conditions, tall columns will not spontaneously form. Furthermore, do these vertical vortex interactions destroy the horizontal coherence of the initial vortices. This too would prevent the growth of intense and large vortices and send energy downscale. In the final Section 5, a discussion and conclusion are provided.

Model
As in the first part [19], we used a three-layer quasi-geostrophic model with the following parameters characteristic of the North Atlantic Ocean basin, and water stratification: the total depth was 4 km, the thicknesses of the upper, middle and lower layers were H 1 = 600 m, H 2 = 1000 m and H 3 = 2400 m, respectively (see Figure 1). The first and second deformation radii calculated from the average density stratification correspond to the horizontal scale where the Coriolis and buoyancy effects are similar in magnitude. These deformation radii are equal to: Rd 1 = 32 km and Rd 2 = 15 km [20]. The time scale T * is chosen as the rotational period of an isolated circular vortex patch. Assuming that the characteristic horizontal scale R = Rd 1 and that the maximum velocity of a vortex of unit radius is 40 cm/s, then we get T * ≈ 9 days. This allows us to render our model dimensionless. These theorems are central to vortex dynamics, a specific study of which follows. They allow the full description of the flow by the potential vorticity distribution which is advected by the velocity field that it itself generates. After presenting the model equations and numerical method (Section 2), we will study the interaction of a surface vortex with a deeper and smaller vortex (Section 3). Due to vorticity symmetry, we chose the surface vortex to have positive potential vorticity. In Section 4, we study the interaction of the same surface vortex, with one or two pairs of deep and small vortices. Our aim in these various experiments and cases is to determine the fate of the vortices: do they remain aligned vertically, or coherent (in close vicinity), or do they irreversibly separate? The latter evolution indicates that under given conditions, tall columns will not spontaneously form. Furthermore, do these vertical vortex interactions destroy the horizontal coherence of the initial vortices. This too would prevent the growth of intense and large vortices and send energy downscale. In the final Section 5, a discussion and conclusion are provided.

Model
As in the first part [19], we used a three-layer quasi-geostrophic model with the following parameters characteristic of the North Atlantic Ocean basin, and water stratification: the total depth was 4 km, the thicknesses of the upper, middle and lower layers were H1 = 600 m, H2 = 1000 m and H3 = 2400 m, respectively (see Figure 1). The first and second deformation radii calculated from the average density stratification correspond to the horizontal scale where the Coriolis and buoyancy effects are similar in magnitude. These deformation radii are equal to: Rd1 = 32 km and Rd2 = 15 km [20]. The time scale T * is chosen as the rotational period of an isolated circular vortex patch. Assuming that the characteristic horizontal scale R = Rd1 and that the maximum velocity of a vortex of unit radius is 40 cm/s, then we get T * ≈ 9 days. This allows us to render our model dimensionless.
In this model, we consider one or several intra-thermocline vortex/vortices, idealized as a patch in the middle layer, with uniform positive or negative potential vorticity (PV). The surface vortex will be represented as a patch of uniform, positive, potential vorticity in the upper layer. Note that the three-layer quasi-geostrophic model does not completely reproduce real ocean vortices, but this model has been shown to satisfactorily explain vortex interactions such as the merging of vortices of the same sign belonging to one layer [16], the interaction of vortices located at different layers [22], and the effect of the bottom topography on the dynamics of intrathermocline lenses [20]. In this model, we consider one or several intra-thermocline vortex/vortices, idealized as a patch in the middle layer, with uniform positive or negative potential vorticity (PV). The surface vortex will be represented as a patch of uniform, positive, potential vorticity in the upper layer.
Note that the three-layer quasi-geostrophic model does not completely reproduce real ocean vortices, but this model has been shown to satisfactorily explain vortex interactions such as the merging of vortices of the same sign belonging to one layer [16], the interaction of vortices located at different layers [22], and the effect of the bottom topography on the dynamics of intrathermocline lenses [20].

Mathematical Model
The quasi-geostrophic model assumes the predominance of horizontal motion for ocean dynamics (see, e.g., [23]). Two dimensionless numbers characterize quasi-geostrophic dynamics: the Rossby number Ro = U/fR, where U is the magnitude of the horizontal velocity, f is the Coriolis parameter and R is the horizontal scale of motion; and the Burger number Bu = (N 2 H 2 /f 2 R 2 ), where N is the buoyancy frequency (static stability frequency) and H is the vertical scale of motion. In the quasi-geostrophic model, Ro is small, and Bu is of order unity. Here, we will also assume that R is much smaller than the scale of planetary variations of f, so that f is constant (f -plane approximation). Therefore, the fluid is hydrostatic, and the motion is governed at first order by the balance of the Coriolis acceleration and of the horizontal pressure gradient, in the horizontal momentum equations. A stream-function is then obtained by dividing the dynamic pressure by density and by the Coriolis parameter. The curl of the Euler momentum equations on a rotating planet provide a vorticity equation. In a reasonable approximation of the ocean, we assume that the model stratification is layer wise (layers of constant density). Then, the horizontal velocity is uniform vertically in each layer. Finally, we obtain that the quasi-geostrophic model is governed by the conservation of potential vorticity in each layer: where the indices i = 1, 2, 3 denote the upper, middle and bottom layers, respectively, Π i , ψ i and h i represent the potential vorticity (PV), stream-function and dimensionless thickness for the i-th layer, respectively, ∇ 2 = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 and d i /dt = ∂/∂t + u i ∂/∂x + v i ∂/∂yare two-dimensional Laplace and total time operators, respectively. Potential vorticity is thus composed of relative vorticity (the curl of horizontal vorticity, the Laplacian of the stream-function), and of vortex stretching (the vorticity induced by the stretching or squeezing of fluid columns via the Coriolis acceleration; the matrix A times the stream-function vector Ψ). We also used the formal vector Π with 3 components Π i . Using the method of diagonalization described in [24] reduced (2) to system for the components of auxiliary vectors Z and Φ: where E is a unit matrix and matrices Q and S have the forms: The Froude Numbers F 1 and F 2 here and in (2) are expressed in terms of the eigenvalues λ i : of the spectral problem Aq + λq = 0; with eigenvectors q (1) , q (2) , q (3) . For the non-zero eigenvalues: Setting, as before, H 1 = 600 m, H 2 = 1000 m and H 3 = 2400 m; and R = Rd 1 = 32 km, Rd 2 = 15 km we obtain from (5) and (6) F 1 = 0.14 and F 2 = 0.7378 with h 1 = 0.15, h 2 = 0.25, h 3 = 0.6.
We can determine the components ϕ j of the auxiliary vector Φ from (3) using the Green's function G i : where K 0 is the modified Bessel function of zero order (below the notations K 1 will be used for the modified Bessel functions of the first order), r = (x − x ) 2 + (y − y ) 2 (here, (x, y) and (x , y ) are the coordinates of the point of observation and the point of integration, correspondingly), γ 1,2 = −λ 2,3 . The streamfunctions in the layers are convenient to represent using the regular functions G 1 = G 2 − G 1 , G 2 = G 3 − G 1 , and G 1 which has an integrated singularity at r → 0 : where q il and s il (l = 1, 2, 3) are the elements of the matrices Q and S introduced above. Equation (7) allows us to calculate the stream functions for each of the layers if the distributions of Π i (x, y) are known.

Numerical Models
The most convenient and effective method to solve the equations for our problem is the contour dynamics method which assumes that the potential vorticity in each layer is zero, everywhere except for a finite set of areas in which it takes constant values. In this case, the improper integral in (7) is replaced by a set of integrals over finite areas. Let us, then consider that in each layer, PV is the where N i is the number of vortex patches in the ith layer and Π j i are the potential vorticities of vortex patches, which are nonzero constants inside compact areas S j i (i = 1, 2, 3; j = 1, . . . , N i ). Thus, instead of (7), we can write: Then, using the Stokes theorem for a transition the double integrals to contour ones we obtain: where C n i are the contours of areas S n i described with the continuously varying parameter ν n i in the positive direction (counterclockwise) along each contour C n i , and: where the dot superscript represents differentiation with respect to ν n m . Then, we select a set of M n ij reference points (Lagrangian nodes) on each contour C n i , and can write the equations of motion: where r n ij is the radius vector of jth Lagrangian node of contour C n i bounding the nth vortex patch in the layer i-th layer. The velocities in the right and sides of Equation (10) are determined from (9) using relations of the form: They are related by the means of contour integrals, which allows us to use the version of the contour dynamics method [25] generalized to the case of three-layer rotating fluid in [26].
The numerical model [26] that we use here was added by the "contour surgery" procedure (CS) [27].
Note that the equations of motion of point vortices in [19] are obtained from (9) and (11) in the limiting case: Then, using the Stokes theorem for a transition the double integrals to contour ones we obtain: where are the contours of areas described with the continuously varying parameter in the positive direction (counterclockwise) along each contour , and: where the dot superscript represents differentiation with respect to .
Then, we select a set of reference points (Lagrangian nodes) on each contour , and can write the equations of motion: where is the radius vector of jth Lagrangian node of contour bounding the nth vortex patch in the layer i-th layer. The velocities in the right and sides of Equation (10) are determined from (9) using relations of the form: They are related by the means of contour integrals, which allows us to use the version of the contour dynamics method [25] generalized to the case of three-layer rotating fluid in [26].
The numerical model [26] that we use here was added by the "contour surgery" procedure (CS) [27].
Note that the equations of motion of point vortices in [19] are obtained from (9) and (11) in the limiting case: where is an undisturbed circular area and: where S l i is an undisturbed circular area and: where δ is Dirac delta function and κ l i is the intensity of a point vertex.

Physical Model of the Present Problem
This model is composed of immiscible layers and the interaction between vortex patches, situated in different layers, occurs via deformed interfaces under/above them. The conservation of PV leads to the inevitable formation of nonzero relative vorticity in the adjacent layer. The behavior of interfaces is described by the equations: Due to the geostrophic and hydrostatic balances, the vortex patch of the upper layer with a positive/negative PV will initiate a local curvature of the interface between the upper and the middle layers η 1 , with the convexity directed upward/downward. The cyclonic/anticyclonic vortex of the middle layer will generate a configuration of η 1 and η 2 in the form of double concave/double convex lens. As an example, Figure 2 shows a vertical section of these surfaces' shapes, calculated by the Equation (14) for y = 0 and t = 0 for the case N 1 = 1, Π 1 1 > 0, N 2 = 2, Π 1 2 = −Π 2 2 > 0 inside the initially circular regions; and N 3 = 0. Note that in the figure, the vertical scale is significantly increased for clarity. The dimensionless radii of the vortices in the figure are as follows: R 1 1 = 4 and R 1 2 = R 2 2 = 1. middle layer will generate a configuration of η1 and η2 in the form of double concave/double convex lens. As an example, Figure 2 shows a vertical section of these surfaces' shapes, calculated by the Equation (14) for y = 0 and t = 0 for the case 1 = 1, 1 1 > 0, 2 = 2, 2 1 = − 2 2 > 0 inside the initially circular regions; and 3 = 0. Note that in the figure, the vertical scale is significantly increased for clarity. The dimensionless radii of the vortices in the figure are as follows: 1 1 = 4 and 2 1 = 2 2 = 1. Note that in the interaction of vortex patches, the key role is played not by the local potential vorticity , but by the "effective" PV: P = ℎ ( -is the area of the horizontal section of the corresponding column with a PV equal to ). Hereafter, we will vary the locations of the middle layer, or intrathermocline lenses (ITLs) and dipole structures relative to the position of the upper layer cyclone. Using the contour dynamics method, we studied the features of their joint evolution. The various numerical experiments aim at determining: 1. If the surface cyclone over a large intrathermocline vortex forms a stable structure.  Note that in the interaction of vortex patches, the key role is played not by the local potential vorticity Π Hereafter, we will vary the locations of the middle layer, or intrathermocline lenses (ITLs) and dipole structures relative to the position of the upper layer cyclone. Using the contour dynamics method, we studied the features of their joint evolution. The various numerical experiments aim at determining: 1.
If the surface cyclone over a large intrathermocline vortex forms a stable structure.

2.
If the surface cyclone remains coherent and/or if it sheds filaments, small eddies, or small eddy pairs. 3.
The evolution of the intrathermocline vortices and their ability to deform the upper cyclone with a characteristic surface signature.

4.
If the intra-layer or inter-layer interactions dominate and in particular, if horizontal vortex pairs or vertical vortex pairs (hetons) or the vertical alignment of like-signed vorticity patches, are the prevalent mechanism.

Vertically Aligned Vortices
Following [28,29], in the simplest case, when the centers of the surface vortex and the lens lie on the same vertical axis, we obtain an axisymmetric vortex. We can study its stability with respect to the perturbations of vortex patch contours. Figure 3 shows the neutral stability curve of a perturbation with azimuthal mode m = 2, in part of the plane of the variables r 2 ≡ r 1 2 and Π 2 ≡ Π 1 2 . The analysis [28,29] showed that perturbations with higher modes are stable, within the selected intervals of the parameters r 2 and Π 2 . Following [28,29], in the simplest case, when the centers of the surface vortex and the lens lie on the same vertical axis, we obtain an axisymmetric vortex. We can study its stability with respect to the perturbations of vortex patch contours. Figure 3 shows the neutral stability curve of a perturbation with azimuthal mode m = 2, in part of the plane of the variables 2 ≡ 2 1 and 2 ≡ 2 1 .
The analysis [28,29] showed that perturbations with higher modes are stable, within the selected intervals of the parameters 2 and 2 . The stability diagram shows that vortex structures, having a wide and weak PV component in the middle layer, can be unstable.
The nonlinear evolution of the linearly unstable vortex structure with 1 = −0.4444 • 2 is shown in Figure 4. The numerical calculation shows that both contours retain their circular shape for a rather long time. Then, from a weak numerical noise, the second azimuthal mode of deformation amplifies, and forms two peripheral vortices in the upper layer. The velocity shear that the upper and lower vortices exert on each other lead to their (dipolar) breaking. They form a heton (a two-layer vortex dipole); the filaments that they shed roll up via horizontal shear instability into smaller, peripheral vortices. It is noteworthy to mention that the remnant of the upper cyclone breaks up into two vortex patches, vertically aligned with the two fragments of the middle layer vortex (called here ITL). This illustrates the dynamical influence of this ITL in the upper layer (reflecting the ratio of integrated vorticities). Furthermore, the substantial fragmentation of the initial vortex patches can be associated with an energy transfer towards smaller scales (a direct energy cascade) and with the substantial stirring of water masses. Mixing then occurs when very small vorticity patches reach the threshold size for contour removal, an "ersatz" of real viscosity.   This example clearly shows that the presence of a sufficiently large ITL, located under the surface vortex, can lead to the destruction of the latter into smaller components due to instability. According to Figure 3, an axisymmetric two-layer structure, containing an anticyclonic lens of unit radius, for example, will always be stable, keeping its shape.  The stability diagram shows that vortex structures, having a wide and weak PV component in the middle layer, can be unstable.

Vertically Shifted Vortices
The nonlinear evolution of the linearly unstable vortex structure with P 1 = −0.4444·P 2 is shown in Figure 4. The numerical calculation shows that both contours retain their circular shape for a rather long time. Then, from a weak numerical noise, the second azimuthal mode of deformation amplifies, and forms two peripheral vortices in the upper layer. The velocity shear that the upper and lower vortices exert on each other lead to their (dipolar) breaking. They form a heton (a two-layer vortex dipole); the filaments that they shed roll up via horizontal shear instability into smaller, peripheral vortices. It is noteworthy to mention that the remnant of the upper cyclone breaks up into two vortex patches, vertically aligned with the two fragments of the middle layer vortex (called here ITL). This illustrates the dynamical influence of this ITL in the upper layer (reflecting the ratio of integrated vorticities). Furthermore, the substantial fragmentation of the initial vortex patches can be associated with an energy transfer towards smaller scales (a direct energy cascade) and with the substantial stirring of water masses. Mixing then occurs when very small vorticity patches reach the threshold size for contour removal, an "ersatz" of real viscosity.
This example clearly shows that the presence of a sufficiently large ITL, located under the surface vortex, can lead to the destruction of the latter into smaller components due to instability. According to Figure 3, an axisymmetric two-layer structure, containing an anticyclonic lens of unit radius, for example, will always be stable, keeping its shape.  This example clearly shows that the presence of a sufficiently large ITL, located under the surface vortex, can lead to the destruction of the latter into smaller components due to instability. According to Figure 3, an axisymmetric two-layer structure, containing an anticyclonic lens of unit radius, for example, will always be stable, keeping its shape.  At first glance, it seems that there is a qualitative similarity of the translational velocity of the point and finite-core hetons: both profiles are nonmonotonic, and in the second case, the maximum velocity is reached near the boundary of the surface vortex, which is quite obvious. Now, let us consider in more detail the case of finite-core vortices. At very small and very large distances L (compared with the vortex radius), the centers of the vortices move qualitatively in the same way as  [19] and finite-area vortex patches in the same layers (black-red line) at r 1 = 4; P 1 = 0.6 and r 2 = 1; P 2 = −0.6 on the distance L between vortices in the first case and their centers in the second case. The part of the curve in which the surface vortex loses its compactness and partially collapses is marked in red. Markers on the black (L = 0.71, L = 7.00) and red (L = 1.00, L = 6.50) lines correspond to the experiments presented in Figure 6. the point vortices, and the contours of the vortex patches during their motion keep a shape close to circular. The situation changes in the vicinity of the red part of the curve in Figure 5 and, even more, in the red part itself. Figures 6 and 7 show the particular motions of such hetons (two-layer vortices). Obviously, the total length of the trajectories in Figure 6 is proportional to the translational velocity of the heton (compare with Figure 2a in [19]). Panels 6a and 6d show the contours of the vortex patches when the values of correspond to points on the black line in Figure 5 in the vicinity of the red part of the curve. During the motion of the two-layer pair, the ITL's contour retains its circular shape, while the contour of the surface vortex undergoes deformations, but this vortex remains compact. Here, as well as in Figures 7-11, the initial configurations of the vortex spots are shaded. This allows us to determine the direction of the motion of all the vortices.  At first glance, it seems that there is a qualitative similarity of the translational velocity of the point and finite-core hetons: both profiles are nonmonotonic, and in the second case, the maximum velocity is reached near the boundary of the surface vortex, which is quite obvious. Now, let us consider in more detail the case of finite-core vortices. At very small and very large distances L (compared with the vortex radius), the centers of the vortices move qualitatively in the same way as the point vortices, and the contours of the vortex patches during their motion keep a shape close to circular. The situation changes in the vicinity of the red part of the curve in Figure 5 and, even more, in the red part itself. Figures 6 and 7 show the particular motions of such hetons (two-layer vortices). Obviously, the total length of the trajectories in Figure 6 is proportional to the translational velocity of the heton (compare with Figure 2a in [19]). Panels 6a and 6d show the contours of the vortex patches when the values of L correspond to points on the black line in Figure 5 in the vicinity of the red part of the curve. During the motion of the two-layer pair, the ITL's contour retains its circular shape, while the contour of the surface vortex undergoes deformations, but this vortex remains compact. Here, as well as in Figures 7-11, the initial configurations of the vortex spots are shaded. This allows us to determine the direction of the motion of all the vortices.   Motion of a two-layer heton (vertical dipolar vortex), with the same parameters r 1 , r 2 , P 1 , P 2 , as in Figure 6, but with L = 4 at the indicated dimensionless times. On the contrary, in panels (b) and (c) of Figure 6, the surface vortex loses its compactness over time, and a part of its mass is transferred into filaments and into small-scale vortices. This is especially evident in Figure 6b, where the deviation of the vortex motion from rectilinear is due to the fragmentation of the vortex patch. We can conclude that the formation of small-scale structures (filaments in particular) and the vorticity stirring are maximal when the lens (middle layer vortex) lies close to the surface vortex rim. Indeed, this is the location of maximal velocity shear.

Vertically Shifted Vortices
In Figure 7, we examine in more details the intense interaction of the upper cyclone and ITL when the translational velocity of the pair is maximal at = 4, and initially, when the center of the lens is located exactly below the rim of the surface cyclone. In this case, the lens immediately significantly affects the surface vortex contour; a part of this contour begins to "wrap" around the lens, and later breaks up into a few small vortices. One of the small vortices (about 1.5% of the initial volume) remains in the vicinity of the lens, while the largest part of the surface cyclone (about 70% of the initial volume) remains compact and is advected by the "strong" lens in a collective clockwise rotation.
Note that in the experiment presented in Figure 7, the surface vortex and the lens have an equal contribution to the global motion of the two-vortex structure, whereas in the last experiment, the lens eventually plays a dominant role as the surface vortex weakens due to filamentation.
Obviously, when 2 decreases, the dominant role of the lens decreases, and calculations show that a change in the nonlinear regime takes place at 2 ≈ 0.575. A series of experiments in Figure 8 shows the behavior of lens when the surface vortex is dynamically dominant. In this Figure, the potential vorticity of the lens is smaller by a factor of two (panel 8a), of four (panel 8b) and of eight (panel 8c), compared with the case of Figure 7.
For the twice weaker lens, a global rotation of the two-layer vortex occurs in the opposite direction to the case of Figure 7. The collective motion of the two-vortex structure is then largely determined by the surface cyclone, though, under the shearing influence of the lens, this cyclone partially collapses and loses about 17% of its volume. In the case of the four times weaker lens, the areal vorticity loss by the cyclone via filamentation is already just over 5%, and in the case of the eight-time attenuation, the cyclone remains compact, but as before, the cyclone does not remain stationary under the influence of the lens. The area filled with superimposed cyclone contours in Figure 8c significantly exceeds the original, and the motion of vortex patches is then qualitatively close to the motion of point vortices (compare with Figure 2b,c in [19]). and 32 km). To show the generality of this evolution, we consider the case when this radius is 1.5 times larger, i.e., r1 = 6 (i.e., 192 km and 32 km). Figure 9b shows that no qualitative changes compared with Figure 9a were observed. In this case, 91.1% of the original volume remains in the central core. Note, that in Figure 9b, for convenience, the dimensions of the surface vortices are shown on the same scale as in Figure 9a, but their actual increase by 1.5 times is visible by the corresponding decrease in the size of the intrathermocline vortices. A comparison of Figure 9a,b with Figure 3a,b in [19] suggests that the theory of point vortices satisfactorily explains the main qualitative features of the infinite mode of distributed vortices. However, there are important differences: (1) in the middle layer, finite-core dipoles which "escape" from under a large-scale surface cyclone recede not along a rectilinear (ballistic) trajectory, but along a divergent spiral; (2) due to the vertical alignment mechanism, these dipoles acquire a two-layer structure with a component in the upper layer. Though vertical alignment seems commonplace for vortex interactions, it has a really unexpected consequence in the ocean: a small subsurface dipole which escapes the influence of a large surface vortex can tear a fragment from it and carry this fragment far away.

Collinear Initial Configuration
We add that bounded solutions, as in Figure 3c,d in [19], for finite-core lenses located under the surface vortex (at least within the framework of the considered geometric parameters) fail. Two examples in Figure 10 will suffice to illustrate the similarities and differences between the two representations of vortices at = 1 and = 6 (strong and weak pairs). In Figure 10a, a strong deformation of the surface vortex was generated by the shear stress exerted by the middle layer dipole. External fluid enters the cyclone, while the dipole is deflected west of the cyclone. At > 6, the dipole moves away, tearing and advecting away two fragments of this cyclone. Most of the surface cyclone is left behind, but with substantial filamentation and intrusion of external fluid. Qualitatively, the results of the interaction here and in Figure 10b (with the same value of 1 ) are similar.

Impact of the Two External Intrathermocline Vortices on the Surface Cyclone
In Figure 10b, in the case of a weakly intrathermocline dipole, takes place the capture of a cyclonic vortex of middle layer by the surface cyclone, to form a two-layer cyclone; the anticyclonic lens rotates on its periphery.

Case of Asymmetric Middle Layer Vortices
To date, it has been assumed everywhere that the two middle layer vortices have the same potential vorticity. However, observations [15,17] 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 2 1 2 2 ⁄ = −1 2 ⁄ . Two experiments at different values of show the main evolutions and their differences with the discrete model. Figure 11a, for = 2, shows that for relatively close vortices of ITL structure, the intralayer interaction is dominant: both middle layer vortices rotate clockwise mutually and counterclockwise around the surface cyclone as point vortices in Figure 6 in [19]. Their slowing down during their We note two circumstances. First, in Figure 11a, all vortex patches retain their compactness and behave almost like point vortices, but in the Figure 11b, noticeable deformations and even the partial destruction of the surface vortex contour are observed. A small fragment, separated from the cyclone, moves along a peripheral trajectory; it follows the motion of the middle layer anticyclone with about half a period lag. Second, Figure 11b shows that the modal evolutions are not characteristic of point On the contrary, in panels (b) and (c) of Figure 6, the surface vortex loses its compactness over time, and a part of its mass is transferred into filaments and into small-scale vortices. This is especially evident in Figure 6b, where the deviation of the vortex motion from rectilinear is due to the fragmentation of the vortex patch. We can conclude that the formation of small-scale structures (filaments in particular) and the vorticity stirring are maximal when the lens (middle layer vortex) lies close to the surface vortex rim. Indeed, this is the location of maximal velocity shear.
In Figure 7, we examine in more details the intense interaction of the upper cyclone and ITL when the translational velocity of the pair is maximal at L = 4, and initially, when the center of the lens is located exactly below the rim of the surface cyclone. In this case, the lens immediately significantly affects the surface vortex contour; a part of this contour begins to "wrap" around the lens, and later breaks up into a few small vortices. One of the small vortices (about 1.5% of the initial volume) remains in the vicinity of the lens, while the largest part of the surface cyclone (about 70% of the initial volume) remains compact and is advected by the "strong" lens in a collective clockwise rotation.
Note that in the experiment presented in Figure 7, the surface vortex and the lens have an equal contribution to the global motion of the two-vortex structure, whereas in the last experiment, the lens eventually plays a dominant role as the surface vortex weakens due to filamentation.
Obviously, when P 2 decreases, the dominant role of the lens decreases, and calculations show that a change in the nonlinear regime takes place at P 2 ≈ 0.575. A series of experiments in Figure 8 shows the behavior of lens when the surface vortex is dynamically dominant. In this Figure, the potential vorticity of the lens is smaller by a factor of two (panel 8a), of four (panel 8b) and of eight (panel 8c), compared with the case of Figure 7.
For the twice weaker lens, a global rotation of the two-layer vortex occurs in the opposite direction to the case of Figure 7. The collective motion of the two-vortex structure is then largely determined by the surface cyclone, though, under the shearing influence of the lens, this cyclone partially collapses and loses about 17% of its volume.
In the case of the four times weaker lens, the areal vorticity loss by the cyclone via filamentation is already just over 5%, and in the case of the eight-time attenuation, the cyclone remains compact, but as before, the cyclone does not remain stationary under the influence of the lens. The area filled with superimposed cyclone contours in Figure 8c significantly exceeds the original, and the motion of vortex patches is then qualitatively close to the motion of point vortices (compare with Figure 2b,c in [19]).

Collinear Initial Configuration
Here, after a short and oceanographic justification, we study two cases of unbounded motions to evaluate the influence of the finite area of vortices (mainly, the cyclone of the upper layer).
In the Gulf of Cadiz, meddies and their cyclonic partners are often generated near the Portimao Canyon, on the northern side of the gulf. Then, they move southward under the cyclonic gyre of the gulf. This is why we initialize a middle layer dipole under a surface cyclone.
The formation of a new type of vortex structure occurs when the intrathermocline dipole is initially located below the surface vortex core. Figure 9a shows that a sufficiently "strong" middle layer pair tries to self-propel outwards of the surface vortex. However, due to the vertical alignment and hetonic [30] mechanisms, the lenses tear the surface cyclone apart: most of its core is left behind (72.7% of the initial volume) and two fragments of this surface cyclone pair up with the lenses and drift away. The two-layer structure, including the lenses, has a complex form: one part consists of two cyclonic vortex patches in the upper and middle layers, and the second part has a hetonic nature. In all the cases of finite-core vortices studied above, the dimensionless radii of the surface cyclone r 1 and intrathermocline vortices r 2 are 4 and 1, correspondingly, (at the selected scale, they are 128 km and 32 km). To show the generality of this evolution, we consider the case when this radius is 1.5 times larger, i.e., r 1 = 6 (i.e., 192 km and 32 km). Figure 9b shows that no qualitative changes compared with Figure 9a were observed. In this case, 91.1% of the original volume remains in the central core. Note, that in Figure 9b, for convenience, the dimensions of the surface vortices are shown on the same scale as in Figure 9a, but their actual increase by 1.5 times is visible by the corresponding decrease in the size of the intrathermocline vortices.
A comparison of Figure 9a,b with Figure 3a,b in [19] suggests that the theory of point vortices satisfactorily explains the main qualitative features of the infinite mode of distributed vortices. However, there are important differences: (1) in the middle layer, finite-core dipoles which "escape" from under a large-scale surface cyclone recede not along a rectilinear (ballistic) trajectory, but along a divergent spiral; (2) due to the vertical alignment mechanism, these dipoles acquire a two-layer structure with a component in the upper layer. Though vertical alignment seems commonplace for vortex interactions, it has a really unexpected consequence in the ocean: a small subsurface dipole which escapes the influence of a large surface vortex can tear a fragment from it and carry this fragment far away.
We add that bounded solutions, as in Figure 3c,d in [19], for finite-core lenses located under the surface vortex (at least within the framework of the considered geometric parameters) fail.

Impact of the Two External Intrathermocline Vortices on the Surface Cyclone
Now, let two ITLs of opposite signs initially be located relatively far from the surface vortex (in all the examples considered below, y 1 2 = y 2 2 = −16 and x 2 2 = −x 1 2 = L are assumed) and move towards it.

Case of Intrathermocline Dipole
Let us start with the simplest case when the vortices of the middle layer form a pair i.e., P 1 2 /P 2 2 = −1, symmetric with respect to the y axis; we consider different distances (2L) between the vortices in the pair.
Two examples in Figure 10 will suffice to illustrate the similarities and differences between the two representations of vortices at L = 1 and L = 6 (strong and weak pairs).
In Figure 10a, a strong deformation of the surface vortex was generated by the shear stress exerted by the middle layer dipole. External fluid enters the cyclone, while the dipole is deflected west of the cyclone. At t > 6, the dipole moves away, tearing and advecting away two fragments of this cyclone. Most of the surface cyclone is left behind, but with substantial filamentation and intrusion of external fluid. Qualitatively, the results of the interaction here and in Figure 10b (with the same value of r 1 ) are similar.
In Figure 10b, in the case of a weakly intrathermocline dipole, takes place the capture of a cyclonic vortex of middle layer by the surface cyclone, to form a two-layer cyclone; the anticyclonic lens rotates on its periphery.

Case of Asymmetric Middle Layer Vortices
To date, it has been assumed everywhere that the two middle layer vortices have the same potential vorticity. However, observations [15,17] 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 P 1 2 /P 2 2 = −1/2. Two experiments at different values of L show the main evolutions and their differences with the discrete model. Figure 11a, for L = 2, shows that for relatively close vortices of ITL structure, the intralayer interaction is dominant: both middle layer vortices rotate clockwise mutually and counterclockwise around the surface cyclone as point vortices in Figure 6 in [19]. Their slowing down during their loops was noted in [15,17]. Another dynamical regime takes place at L = 4 (Figure 11b): the vortex motion closely resembles the mode 1 structure of Figure 7a in [19].
We note two circumstances. First, in Figure 11a, all vortex patches retain their compactness and behave almost like point vortices, but in the Figure 11b, noticeable deformations and even the partial destruction of the surface vortex contour are observed. A small fragment, separated from the cyclone, moves along a peripheral trajectory; it follows the motion of the middle layer anticyclone with about half a period lag. Second, Figure 11b shows that the modal evolutions are not characteristic of point vortices only.

Discussion and Conclusions
By analogy with [19], we note the main results of the work from the point of view of their oceanographic interpretation: • A two-layer vortex with a vertical axis, consisting of a large-scale cyclone of the upper layer and a sufficiently large anticyclone of the middle layer, can be unstable with respect to small perturbations of the shape of the vortex patches and decompose into smaller two-layer vortex structures. In the case when the intratermocline anticyclone is of the order of Rd 1 , the two-layer vortex behaves stably (Section 3.1). Note that such a vortex structure consisting of point vortices is also always stable.

•
If the centers of the cyclone of the upper layer and the anticyclonic lens of the middle layer are spaced, 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). These properties are also characteristic of point vortices. For finite-core vortices, in addition, there is a certain interval of distances between their centers, in which the motion of the vortex patches is accompanied by a partial destruction of the surface cyclone and a deviation from the trajectory predicted by the theory of point vortices occurs (Section 3.2).

•
If intrathermocline dipole is initially located below the surface vortex core, the formation of a such new type of vortex structure occurs when most of the surface cyclone core is left behind and two of its fragments pair up with the lenses and drift away. The two-layer structure including the lenses has a complex form: one part consists of two cyclonic vortex patches in the upper and middle layers, and the second part has a hetonic nature (Section 4.1).

•
If the intrathermocline dipole is running onto the surface vortex, then (like in the point vortex case [19]) we have two regimes: (a) the pair passes under the surface vortex, changing its direction in its vicinity and capturing part of its mass; (b) the cyclonic vortex of the intrathermocline dipole is captured by the surface cyclone, and the anticyclonic lens rotates at the periphery of the newly formed bilayer cyclonic structure; and, unlike the case of point vortices, the dipole of the middle layer cannot leave the vicinity of the surface cyclone. (Section 4.2.1).

•
If an asymmetry in the distribution of potential vorticity is observed in a vortex structure incident on a surface cyclone, then it makes loop-like movements and, at the same time, rotates around the cyclone of the upper layer. It is striking that in the case of finite-core vortices (depending on the initial distance between them), states can be formed that resemble the modal structures obtained for point vortices (19) as quasi-stationary solutions (Section 4.2.2).
These findings from vortex dynamics deserve further measurements at sea, both from intensive surveys and with Lagrangian floats. It is noteworthy that since cyclonic eddies, as a rule, are weaker than anticyclonic ones at intermediate depths in the ocean [15,17], loop-shaped trajectories enveloping a surface cyclone can be most probable in nature. Further work will investigate the databases of surface drifter trajectories and of deep Argo float trajectories as an oceanographic application of this paper.