Models of interacting pairs of thin, quasi-geostrophic vortices: steady-state solutions and nonlinear stability

Abstract We study pairwise interactions of elliptical quasi-geostrophic (QG) vortices as the limiting case of vanishingly thin uniform potential vorticity ellipsoids. In this limit, the product of the vertical extent of the ellipsoid and the potential vorticity within it is held fixed to a finite non-zero constant. Such elliptical “lenses” inherit the property that, in isolation, they steadily rotate without changing shape. Here, we use this property to extend both standard moment models and Hamiltonian ellipsoidal models to approximate the dynamical interaction of such elliptical lenses. By neglecting non-elliptical deformations, the simplified models reduce the dynamics to just four degrees of freedom per vortex. For simplicity, we focus on pairwise interactions between identical elliptical vortices initially separated in both the horizontal and vertical directions. The dynamics of the simplified models are compared with the full QG dynamics of the system, and show good agreement as expected for sufficiently distant lenses. The results reveal the existence of families of steadily rotating equilibria in the initial horizontal and vertical separation parameter space. For sufficiently large vertical separations, equilibria with varying shape exist for all horizontal separations. Below a critical vertical separation (stretched by the constant ratio of buoyancy to Coriolis frequencies N / f), comparable to the mean radius of either vortex, a gap opens in horizontal separation where no equilibria are possible. Solutions near the edge of this gap are unstable. In the full QG system, equilibria at the edge of the gap exhibit corners (infinite curvature) along their boundaries. Comparisons of the model results with the full nonlinear QG evolution show that the early stages of the instability are captured by the Hamiltonian elliptical model but not by the moment model that inaccurately estimates shorter-range interactions.


Introduction
Observations of the world's oceans indicate complex eddying motions which may account for a transport of heat, salt, momentum and trace constituents comparable to that of the large-scale thermohaline circulation (see Olson 1991, Stammer 1998, Garrett 2000, Carton 2001, Siegel et al. 2001, Chelton et al. 2007, Carton 2010, Chelton, Gaube et al. 2011, Souza et al. 2011, Beron-Vera et al. 2013, Dong et al. 2014, Zhang et al. 2014, for a sample of the vast literature). Indeed, the perception of the role of ocean vortices has transformed dramatically in just 30 years following major improvements in ocean observations, particularly at the sea surface (Ebbesmeyer et al. 1986, Wunsch and Stammer 1996, Samelson et al. 2014. Those observations have revealed a plethora of vortices drifting across ocean basins, interacting with each other, with the major ocean currents and jets, and with bottom topography and marine boundaries. Collectively they induce a turbulent, highly-irregular flow. There is growing evidence that the deep ocean, too, exhibits a ubiquitous field of vortices (e.g. the "meddies" spinning out of the outflow of the Mediterranean at roughly 1000 m depth) (McDowell and Rossby 1978, Hebert et al. 1990, Prater and Sanford 1994, Serra et al. 2005, Casal et al. 2006, Ambar et al. 2008, L'Hégaret et al. 2014. Vortices are important because they exhibit highly anomalous transport: unlike other fluid motions, vortices typically carry material over large distances with little change. This is particularly true of passive tracers, but also for a key dynamical quantity, potential vorticity, which is therefore a convenient marker for ocean vortices (Zhang et al. 2014). The extent to which vortices contribute to the overall ocean circulation is difficult to quantify, but there is no doubt vortices are a common dynamical feature.
The ubiquity of vortices in the oceans, in the Earth's atmosphere and in other planetary atmospheres has motivated an extensive body of research aimed at understanding vortex characteristics, such as their structure, stability, the impact of the surrounding flow, and interactions with other vortices, including merger. Of relevance to the present study, where we consider flows with a non-trivial vertical dependence, is the early work of Polvani (1991), who examined vortex interactions in distinct vertical layers. Steadily-rotating equilibrium states, consisting of two identical vortices of (equal or opposite) uniform potential vorticity in adjacent equal-thickness layers, were found. Such solutions are unstable when the horizontal distance between their centers is sufficiently small (depending on the Burger number, Bu, a measure of the relative importance of rotation to stratification). This instability often leads to vortex "alignment", where the vortex centers come together at the expense of shedding filaments and small vortices in the periphery to conserve angular momentum.
Subsequent research has considered more realistic three dimensional structure and the effects of vortex shape, i.e. the height-to-width aspect scaled by the relevant Coriolis to buoyancy frequency ratio (equal to √ Bu or to the ratio of Rossby and Froude numbers, both considered small). In general, strong vortex interactions only occur when the vortex centers are sufficiently close (depending on parameters), irrespective of the sign, shape and intensity of the vortices (see, e.g. Viera 1995, Sutyrin et al. 1998, Dritschel 2002, Reinaud and Dritschel 2005, Martinsen-Burrell et al. 2006, Bambrey et al. 2007, Ozugurlu et al. 2008, Reinaud and Carton 2009, Reinaud and Dritschel 2009. Here, as in previous studies, we simplify the physical system by making the Quasi-Geostrophic (QG) approximation, tantamount to imposing both hydrostatic and geostrophic balance at leading-order in the equations of motion. Ignoring generally weak dissipative and diabatic effects, QG flow is governed by the following set of equations (Vallis 2008) ∂q ∂t + u · ∇q = 0, q = ψ , u = − ∂ψ ∂y , ∂ψ ∂x , 0 , (1a-c) where q is the potential vorticity, ψ the streamfunction, u the horizontal velocity, and is the three-dimensional Laplacian in which z has been stretched by N/f , the ratio of buoyancy and Coriolis frequencies, both constant. The QG model is essentially the simplest model that includes the effects of both stratification and planetary rotation.
The QG approximation filters out (typically weak) higher frequency inertia-gravity waves and leaves only low-frequency nonlinear potential vorticity (PV) advection. The motion is layer-wise two-dimensional but the flow field itself is determined by the threedimensional distribution of the QG PV. This distribution in general includes boundary contributions associated with buoyancy anomalies (or surface potential temperature anomalies in the atmospheric context, (e.g. Dritschel and Saravanan 1994)). In an adiabatic (frictionless and unforced) flow, both the interior PV and the surface buoyancy fields are materially conserved, making them natural markers for vortices.
The QG model is formally identical to the familiar two-dimensional Euler equations for the vorticity, with the crucial difference that ψ is the full three-dimensional Laplacian of the streamfunction ψ (x, y, z, t), as opposed to the two-dimensional Laplacian of ψ(x, y, t) that appears in the Euler equations. A natural question is whether the QG equations possess any solutions similar to the classical Kirchoff-type solutions of the two-dimensional Euler equations (Kida 1981), which are characterized by a single elliptical region of uniform, non-zero vorticity. This question was answered in the affirmative by Meacham Meacham (1992) and by Zhmur and Shchepetkin (1991) who discovered a class of three dimensional solutions to the QG equations in which the PV takes a uniform value inside an ellipsoidal domain. The discovery of these exact single-ellipsoid solutions originated a large literature concerned with using these solutions to study two or more interacting vortices (see Miyazaki et al. 2001, Dritschel et al. 2004, McKiver 2015. A particularly important research topic has been to determine the conditions under which ellipsoidal vortices can merge, as discussed in Reinaud and Dritschel (2005).
Interestingly, it appears that the QG equations do not possess any flat (verticallylocalized) elliptical solutions with uniform PV. However, Dritschel (2011) has recently discovered a class of non-uniform QG elliptical vortex solutions. These solutions are obtained as the limit of a sequence of ellipsoidal vortices in standard position, each with uniform PV Q and z-semi-axis length c, letting c → 0 while keeping Qc fixed. Dritschel studied this type of limit in the context of the "surface QG model," in which the elliptical vortex is located on the surface of a semi-infinite, three-dimensional "ocean." In this case, one needs to use the method of images and consider the limit of two ellipsoidal vortices placed symmetrically above and below the surface (Dritschel 2011). However, the limiting procedure applies equally well to any single ellipsoidal interior QG vortex and leads to a singular (two-dimensional) elliptical sheet of non-uniform PV.
A natural question to ask is whether Dritschel's exact solutions can be used in order to study the interactions among several vortices, using suitable approximation methods. In general, systems of flat, singular vortices contain fewer free parameters than systems containing the same number of three-dimensional ellipsoidal vortices, because one need not consider the height-to width ratio of each vortex nor its vertical tilt. This greatly simplifies, for instance, the study of vortex interactions leading to merger or vertical alignment. The fact that the PV is no longer uniform in the ellipses poses no inherent difficulties. Physically, the situation is akin to that posed in multi-layer QG models (e.g. Polvani 1991, Sokolovskiy andFilyushkin 2015), where instead of examining finite thickness, uniform vorticity patches in piece-wise continuous density profiles, we study vanishingly thin vorticity distributions in constant stratification.
We consider one of the simplest possible situations: the interaction between a pair of elliptical two-dimensional vortices that are initially lying on parallel planes orthogonal to the z-axis. We develop a moment approximation to the dynamics of the two vortices, along the lines of well-known moment models for the two-dimensional Euler equations (see Melander et al. 1986, for instance). Thus, our study is related to that of Sutyrin et al. (1998), who studied the alignment problem for uniform elliptical thin disks of vorticity. In their Ansatz, since the vorticity density is taken to be uniform, each individual vortex does not correspond to an exact QG solution. Here, by contrast, by considering non-uniform initial conditions, we obtain long-lived interacting vortices that for certain separations settle in equilibrium states. The calculation of these equilibria is arguably our most interesting new result.
In addition to deriving the truncated moment model, we also adapt the Hamiltonian ELlipsoidal Model (HELM) originated by Dritschel et al. (2004) to the evolution of QG ellipsoidal vortices of vanishing vertical thickness, i.e. to elliptical lenses. To assess the validity of the approximations employed, we compare predictions from both simplified models to numerical solutions of the full QG equations for the same two-vortex initial conditions.
In the next section, we review the single elliptical solution Dritschel (2011) modified to the present problem. In section 3 we provide a brief derivation of both the moment and Hamiltonian elliptical models. The reduced model results are then compared with the full QG dynamics in section 4.1 for two identical vortices on different horizontal planes. These results reveal the existence of steadily rotating equilibria, which are the subject of section 4.2. Equilibria at any horizontal separation exist for sufficiently large vertical separation. Below a critical vertical separation (stretched by the constant factor N/f ) that is comparable to the mean radius of either vortex, a gap opens up in horizontal separation where no equilibria are possible. The significance of this for stability is explored in section 4.3, which also examines the nonlinear QG evolution of unstable equilibria. The early stages of this instability are captured by the elliptical model but not by the moment model, due to its inaccurate estimate of shorter-range interactions. A few conclusions and ideas for future work are provided in section 5.

A single, non-uniform QG vortex
Any ellipsoid with uniform PV density, Q, together with its self-induced velocity field constitutes a steadily-rotating solution to (1), as originally discovered by Meacham (1992) and by Zhmur and Shchepetkin (1991). Suppose the centroid of the ellipsoid is at the origin and denote the semi-axis lengths by a, b and c (with a ≥ b ≥ c). Assume that the shortest semi-axis (of length c) is aligned with the z axis, and let the uniform potential vorticity inside the ellipsoid have a given value Q. Dritschel's "elliptical" limit (Dritschel 2011) results from letting c → 0 while simultaneously letting Q → ∞, in such a way that the product Qc remains fixed. The streamfunction is given by where D is the ellipsoidal domain. Defining we can re-express the integrals as where D 0 is the elliptical intersection of D with the plane z = 0. Rescaling z by c and using the notation Finally, taking the Dritschel limit c → 0, holding Qc fixed at ω m , and carrying out the z integration we find This is exactly the same streamfunction found by Dritschel (2011), except that in his case ω m takes twice the value because he considers two ellipsoids of uniform PV that collapse into an elliptical patch of buoyancy sitting at the surface of a semi-infinite domain. Here, instead, the elliptical patch is immersed in an infinite fluid. We can identify ω m with the maximum magnitude of the PV density at the center of a "flat" elliptical lens with non-uniform PV density which gives rise to a well-defined, regular flow field everywhere in three-dimensional space, including at its edge. The rotation rate (the analogue to the Kirchhoff rotation rate for the two-dimensional case) is obtained from the limiting value of the known ellipsoidal rotation rate and equals where λ = b/a ≤ 1 and R D is the elliptic integral of the second kind (in Carlson's symmetric form) Explicit formulas for the integrated PV, the enstrophy, the angular impulse J and the energy E can be found in Dritschel (2011), again with the proviso that the value of ω m in that paper is double that studied here. The same reference contains also detailed expressions for the flow field associated with the ellipse, which inherits from the ellipsoid the property of being linear in the interior of the vortex.

QG elliptical vortex models
We now consider a system of N flat QG lenses of the kind obtained in the previous section. We denote by X k ≡ (X k , Y k , Z k ) the centroid of the kth elliptical vortex, by (a k , b k ) its semi-axes (with λ k ≡ b k /a k ) and by φ k the angle from the x-axis to the major semi-axis. Then the initial PV density takes the form with ω k ≡ ω (k) m and on each elliptical domain D k and zero outside. Since in the QG model the velocity field has no z-component, the motion of each vortex will be confined to its respective initial plane z = Z k . However, in general, the lateral motion of the vortices cannot be determined exactly without solving the full QG equations numerically. We now introduce two loworder approximate models which avoid this and preserve the simple elliptical form of each vortex for all time.

Moment model
In this section we introduce a moment model for the interaction between any number of QG vortices. For the sake of simplicity we work out the equations only in the case N = 2. It is straightforward to generalize these equations to arbitrary N. Moment models have a long history in the literature on the 2D Euler equations, starting with the work by Melander et al. (1986). They are based not only on the assumption that each vortex remains elliptical at all times, but also that the major axis of each ellipse be much smaller than the minimum distance between the two vortex centroids. The analysis in the QG case follows closely the one in Melander et al. (1986): (i) The streamfunction for the "far" field produced by each ellipse is Taylor-expanded with respect to the ratio between the major axis of the ellipse and the distance between the centroids. The only difference from the two-dimensional case is that the Green's function used for the streamfunction is three dimensional, and thus leads to different coefficients in the Taylor series. (ii) The "near" field is given at leading order by Dritschel's single-vortex solution. Thus, the PV (density) inside each ellipse is non-uniform and the self-induced angular velocity of each ellipse does not obey Kirchoff's well-known formula for the 2D Euler equations, but rather (3).
Let us consider again the PV density in (4), now with N = 2. By linearity, the corresponding streamfunction is the sum of two terms of the form in (2): where g k is defined in (5). For each domain D k , it is convenient to consider a local Cartesian coordinate system (ξ , η) with origin at the centroid X k . Then, we can rewrite the previous expression as where E k is simply D k translated by X k . Let ψ i denote the streamfunction in the plane of the ith ellipse, and let (Z i − z) 2 = (Z 2 − Z 1 ) 2 ≡ h 2 be the (fixed) squared vertical separation between the two vortices. We may write the two terms in the sum in (7) as follows Here ψ self i is the contribution to the streamfunction from the term with k = i, that is, by the ellipse itself; while ψ far i is the contribution from the other ellipse (k = i). Since we are assuming that the vortices remain elliptical, ψ self i must be given by Dritschel's single vortex solution discussed in the previous section. Since we are also assuming that the two vortices are well separated, we can approximate the effect of one vortex on the other by Taylor-expanding ψ far i under the assumption when x belongs to the ith ellipse and the index j refers to the other ellipse. This gives where are the (weighted) geometric moments of order (m, n) in the local system of coordinates (ξ , η). Notice that the first-order moments J vanish because the integrand is odd. This leaves us with The moments can be calculated using elliptic coordinates in terms of λ i and φ i . A summary of the calculation is given in appendix A; it leads to i is the angular momentum or impulse: The difference of the same two moments is Using these results, we can now derive the model equations for two interacting elliptical lenses. These equations form a set of 6 ODEs for the aspect ratios λ 1 , λ 2 , the tilting angles φ 1 , φ 2 , the inter-centroid distance R and the angle of rotation θ.
In order to define the centroid of each vortex, we define its "weighted area" by The last two equalities reflect the fact that this quantity coincides with the volumeintegrated PV of the "flattened ellipsoid" divided by ω i (see Dritschel 2011, equation (10)), and therefore is conserved. It is then natural to define each centroid by and consider its derivative with respect to time. Observe that for x fixed, the total time derivative of this integral is zero because the function g i (x − X i , y − Y i ) moves together with D i . Therefore, Since the self-induced centroid velocity of each vortex is zero (by symmetry), we need only consider the effect of the other vortex on the centroid motion. Thuṡ where U far i is the velocity field generated by the "far-field" streamfunction in (11). Taylor expanding U far i around the centroid X i , we obtaiṅ Next we truncate the series at the second-order moments, and substitute ψ far i from (11). Observe that, to the desired order of approximation, the second derivatives inside the square brackets in (22) need to be applied only to the term proportional to J (0,0) j in (11). Switching to relative polar coordinates and using the symmetry θ 1,2 = θ 2,1 + π ≡ θ, the centroid evolution equations becomė Finally, these equations can be easily combined into two equations for θ and R: where we recall that T = √ R 2 + h 2 . We next derive the evolution equations for the local geometrical moments. Starting from their definition, we havė where (u, v) is the velocity field inside the ith ellipse. By linearity, it can be written as the sum of two terms. The first is U far i , the field produced by the other, "distant" vortex, to be obtained again from the streamfunction ψ far i . The second is the velocity field produced by the ith ellipse itself. For this we take Dritschel's exact solution for a single elliptical vortex. Thus, we writė whereJ (m,n) i,self is the moment's rate of change due to the uniform self-induced vortex rotation.
For the lower moments, we obtain the following set of ODEs: which are formally identical to the equations for the two-dimensional Euler case found in Melander et al. (1986), except that they cannot be simplified using (∂ 2 x + ∂ 2 y )ψ = 0 because in the QG model ψ = 0 only in the three-dimensional sense. Adding and subtracting equations (28) and (29) and exploiting formulas (15) and (16), we obtaiṅ The only time dependence in Dritschel's solution comes from the steady rotationφ i = Ω i , whereas λ i is constant. Therefore, from equations (15) and (16), it follows thatJ i,self = 0 andḊ i,self = −4J (1,1) i Ω i . We can also combine (12)- (14) in order to find J (1,1) i without solving (30): Substituting (11), (15) and (16) into (31) and canceling the factors 2A 2 i λ i /15π, we find that Therefore, the aspect ratios λ i , i = 1, 2, evolve according tȯ Similarly, from (32) we find that which simplifies, using (35), tȯ where Ω i is the steady angular velocity obtained by applying (3) to the ith vortex. For reference, the complete, non-dimensional set of equations defining the moment model is given in appendix B.
The final model is similar to that presented in Sutyrin et al. (1998) for uniform elliptical vortices, the main difference being that here all the integrals over the elliptical patches include the PV density function g k (x, y). At leading order in the expansion the resulting ODEs are identical. The far-field approximation cannot distinguish between uniform and non-uniform vorticity distributions. However, for our non-uniform vortices, the self-induced angular velocity Ω i of each vortex (assigned by (3)) is exact, whereas it only approximates the (not necessarily constant) self-induced rotation frequency for thin patches of constant vorticity.

The Dritschel-Reinaud-McKiver HELM model
We now introduce a more sophisticated model for the evolution of QG elliptical lenses along the lines of the HELM developed by Dritschel et al. (2004) for the evolution of a system of QG ellipsoidal vortices. Once again, the vortices are assumed to remain ellipsoidal at all times, but there is no separation assumption like the one we introduced for the moment model. Therefore, the Dritschel-Reinaud-McKiver model has more general validity; the trade-off is that the resulting ODE's are more complicated.
At any given time the boundary of the ith ellipse is determined by a time-dependent 2 × 2 symmetric matrix B i such that Here only, we abuse the bold-faced notation for 3-vectors employed throughout, by using x − X i to indicate the 2-vector in the z = Z i plane obtained by removing the third component of x − X i , which is identically equal to zero when x is a point on the ellipse. We remark that B i has eigenvalues a i , b i , the lengths of the semi-axes of the ellipse, whose direction is assigned by the corresponding unit eigenvectorsâ i ,b i . We recall that the HELM model also starts from a representation of the form in (38), but for an ellipsoid; in that case B i is a 3 × 3 symmetric matrix and x − X i is a "true" 3-vector pointing from the centroid of the ellipsoid to any point on its surface. Then, the HELM model consists of 9 ODEs per vortex, 3 describing the evolution of the centroid X i and 6 describing the distinct entries of B i . By contrast, for the limiting case of an elliptical lens, there are only 6 unknowns, the 3 components of X i and the 3 distinct entries of B i , which are denoted as Nevertheless, one can directly verify that the Hamiltonian formalism still applies: B i and X i satisfy the evolution equations In these formulas, κ i = ω (i) m a i b i /3 and L is the symplectic matrix The Hamiltonian H has the form H = H V + H I where H V is the sum of the selfinteraction energies of the individual vortices, and H I is the interaction energy between all possible pairs of vortices. Following Dritschel et al. (2004) where R F is the elliptic integral of the first kind (in Carlson's form) When substituted in (41) this term produces a rotational motion of each vortex with the same constant angular velocity as if it where a single-vortex solution of the kind described in the previous section; the calculation is the same as in Dritschel et al. (2004) (but with c = 0). The interaction Hamiltonian is given by and can be obtained formally from the interaction Hamiltonian for a system of ellipsoids just by changing the order of the semi-axes (so that a ≥ b ≥ c) and then taking the limit c → 0. Whereas (40) together with (41), (43) and (44) do describe the evolution our system of elliptical vortices, they are not practical because the integrals in (44) cannot be calculated exactly. The same difficulty also arises in the original HELM formulation, where it is overcome by a further approximation: the integrals over the ellipsoids are replaced by discrete sums akin to Gaussian quadratures. This is done via a two-step process: first, one shows that the potential generated by a three-dimensional ellipsoid of uniform PV is identical to the potential due to a certain two-dimensional "focal ellipse" of PV in the middle-major axis plane of the ellipsoid (parallel to the y-z coordinate plane). Next, the potential of this elliptical sheet is approximated in terms of a suitable numerical quadrature formula, based on the principle of matching as many spatial moments of the focal ellipse as possible with the chosen number of quadrature points.
In our context the first step is unnecessary: our elliptical vortices are already twodimensional, and coincide with the limit as c → 0 of the focal ellipses parallel to the x-y plane that would generate the same potential as the ellipsoids with c > 0. Therefore, we only need to use the HELM quadrature scheme in order to write the Hamiltonian in discrete form, except that now the points have to lie in the x-y plane. This gives where n is the number of quadrature points, and σ l and x i,l are the strength and location of the lth point in the ith elliptical vortex. Let wherex l = a l ρ l cos ϑ l ,ỹ l = b l ρ l sin ϑ l .
The values of σ l , ρ l , ϑ l for l = 1, . . . , n, and various values of n are found in Dritschel et al. (2004) (with the proviso that each of theirỹ l -z l pairs becomes ax l -ỹ l pair here). The derivatives of H I have the same form as in Dritschel et al. (2004), namely where ∂ x i,l /∂B i can be calculated following the same exact procedure as in Dritschel et al. (2004). We obtain where we have introduced the matrices Finally, by substituting (48) back into (41) and then into (40) we obtain the desired set of 6N ODEs for the time evolution of our system of elliptical vortices. The resulting equations are relatively cumbersome, but we expect them to be more accurate and more generally applicable than the moment model equations obtained in the previous subsection. This is because the present derivation does not rely on any assumption on the separation between the vortices. This does not mean that the model is accurate for any such separation, only that an assumption is not necessary. The present elliptical model only assumes the vortices remain elliptical, and computes all interaction terms which ensure this. In the remainder of this paper we will focus on the case N = 2 and solve these equations numerically. Following Dritschel et al. (2004), who compared different choices of quadrature points, n, for modeling ellipsoid interactions and found very little sensitivity in their results for n ≥ 7, we use n = 7 discretization points in what follows.

QG contour dynamics simulations
The multi-layer QG contour surgery algorithm detailed in appendix A of Dritschel (2002) was applied to study the interaction of lenses of general shape confined to two layers. Contour surgery is a regularized form of contour dynamics wherein sufficiently thin filaments are removed automatically, to control the flow complexity and to enable topological reconnections between regions having the same field isolevel (Dritschel 1998). In the QG algorithm, each layer is given a finite thickness δz, and the vertically-integrated PV over this thickness is used to compute the velocity on a finite set of contour nodes located at the middle of each layer. In the limit δz → 0, the vertically-integrated PV reduces to the PV density defined above. For practical reasons, δz cannot be made identically zero (due to the necessity to choose a correspondingly small time step). Here, we take δz = 0.05 in all results reported, including in the related numerical code to find the generally non-elliptical equilibria (see below). A further necessary approximation is to replace the continuous PV distribution with a piecewise uniform one. Then, the internal distribution of each vortex is modelled by a finite set of contours. Here, 10 contours are used in each lens initially, each with an identical jump in PV. All simulations reported use the standard dimensionless node spacing parameter μ = 0.2 and large-scale length L = √ 2/4 (Dritschel 1998). The weak effect of these approximations can be gauged by examining how well the algorithm reproduces the steady rotation of an exact elliptical lens (Dritschel 2011). Starting with a vortex having an aspect ratio λ = 0.5 and area π, after one period we find a phase error of −0.039 radians, corresponding to an error in the rotation rate of approximately −0.63%. During this time, the r.m.s. variation in the aspect ratio is 1.3 × 10 −4 . This gives an indication of how well the numerical model approximates the exact dynamics of an infinitely thin lens.
The aspect ratio and phase are computed by matching the second-order moments J (2,0) , J (0,2) and J (1,1) of the numerical solution with those of an elliptical lens of the same area and center. Similarly, the vortex center can be obtained from the first-order moments J (1,0) and J (0,1) , though these are zero in the comparison above. Below, for two lenses, all of these moments are used to compare the moment and elliptical model solutions with the full QG solutions.

Results
We focus on the simplest configuration consisting of a pair of identical elliptic vortex lenses, initially centered at ( ± x c , 0, ± z/2) and with major axes aligned with the x axis (φ i (0) = 0), see figure 1. Without loss of generality, the lens areas are set to A = πab = π and ω m = 1.

Comparison of model dynamics
The evolution of the vortex aspect ratio λ 1 (t) = λ 2 (t) ≡ λ(t) for the two models and the full QG dynamics are shown for x c (0) = (1, 1.25, 1.5) in figure 2 for moderately elongated lenses, with λ(0) = 2/3, and for a vertical separation z = 0.4. At this aspect ratio, the scaled rotation period of an isolated ellipse is 5.33. Consistent with the asymptotics used in their derivation, both reduced models show close agreement with the full QG solution for  large enough horizontal separation distances, here x c (0) ≥ 1.5. The lenses initially repel each other and become more circular before returning, periodically, to their initial shape. Compared to the moment model, the HELM solution provides a better approximation to the extreme λ values observed in the full simulation and shows a slightly slower phase drift over the ∼6 oscillation periods plotted.
The comparisons deteriorate significantly when x c (0) is reduced to 1.25. The full QG solution indicates strong deformation and "merger" of the vortex centers (vertical alignment of the vortex cores) for these initial conditions, whereas both reduced models continue to show periodic repulsion and compression, albeit with significant differences between the two in predicted amplitudes and frequencies. For even smaller values of the initial horizontal separation, x c (0) = 1, all three solutions indicate extreme and rapid elongation of the lenses, resulting in the merger of the vortex centers in the full QG simulation and its elliptically-constrained analogue in the reduced models. Notably, even at small x c (0), the initial behavior of the full solution is closely matched by HELM. In the moment model, despite the severe truncation of the interaction dynamics, quantitatively similar behavior is observed at early times.
A more complete comparison of the model behavior in the (x c (0), z) parameter space is provided in figure 3, where the evolution of the aspect ratio λ(t) and the horizontal separation R h (t) = 2x c (t) are shown. For each z, similar changes in behavior are observed for decreasing values of x c (0). As expected, both HELM and moment models show excellent agreement with the nearly periodic QG solutions for large initial horizontal vortex separations where the lenses repel. Consistent with the asymptotic expansions employed, at smaller values of x c (0) the moment model shows better agreement with both HELM and QG when z is larger.
In each model, there exists a critical separation distance below which the lenses attract each other and experience significant stretching. This distance decreases with increasing vertical separation. For z = 0.2 and 0.4, vortex merger (measured by a rapid decrease in R h (t)) occurs at larger initial separations in the full QG solution than in either HELM or the moment model. Differences between the reduced models and the full QG simulations are greatest at these intermediate separations. However, and perhaps unexpectedly, for initial separations well below the critical value, HELM solutions accurately shadow the dynamics of the full equations during their initial stretching phase where, presumably, the vortices remain nearly elliptic.

Co-rotating steady state solutions
As shown in figures 2 and 3, in all models the lenses repel each other (and contract) for large initial separations and attract each other (and elongate) for small enough initial separations. Therefore, for each initial configuration (φ i (0) = 0, z) there exists a critical fixed value of the vortex centroid x c (0) = x c where both R and λ remain constant for all times. Physically, such steady states correspond to co-rotation of the two vortex system whereθ =φ 1 =φ 2 . In the simplified moment model (B.2)-(B.7),λ i =μ = 0 only occurs when φ i = θ ± mπ. While these solutions, for any model, are only strictly "steady" in the proper rotating frame, we will refer to them as steady states, or 'co-rotating equilibria'.
In any model, for pairs of identical lenses with fixed initial phases, φ i (0), steady state solutions are described by a function of the form f (λ, z, x c ) = 0. Root-finding then provides a full picture of the dynamics in the three-dimensional initial condition space. In the moment model, f (λ, z, x c ) =θ −φ 1 is explicitly defined by (B.4) and (B.6). In HELM and in the contour-based QG model, steady states can only be determined numerically by iterative continuation procedures. Such procedures must be started close to an equilibrium state, e.g. at large separation when the vortices are approximately circular. Then, by decreasing either z or x c , entire families of solutions can be determined, together with their rotation rates Ω. When z is held fixed at a small to moderate value, when decreasing x c a turning point is reached beyond which there are no further equilibria with smaller x c . In this case, special care must be taken to continue the solution branch now for increasing x c to find additional equilibria. Thus, there may be (at least) two equilibria for the same values of x c and z. A similar situation arises when starting near aligned states with x c ≈ 0. Such states need not be nearly circular, but can be found from nearly circular states with large z.
The QG contour equilibria are found using an algorithm closely analogous to the one presented in the appendix of Reinaud and Dritschel (2002). The principal difference is that here we have used a variety of control parameters: in place of varying z or x c , we have used the minimum or maximum x coordinate for the vortex centered at x c , the total angular impulse, J = J 1 + J 2 , of the vortex system, the aspect ratio λ, and the width of one of the vortices. This was necessary to pass through the turning point in x c over the wide parameter space investigated. The steady states obtained were verified to be stationary in the appropriate rotating frame.   For z = 0.9, the moment model captures the location and elliptical shape of the steady states for x c 1.1 with similar accuracy. Despite the breakdown in the validity of the expansion for small separations, the moment model also predicts a steady state branch as x c → 0, albeit with the wrong shape. As found above, the largest differences occur at intermediate horizontal separations where the moment model grossly over-predicts the elongation of the lenses.
As the vertical separation z decreases, the comparison between the HELM and QG results deteriorates, especially for x c < 1. Differences are more extreme in the moment model and are not shown here. For the smallest z, no solutions are found for x c < 0.8 approximately. In this case, a turning point occurs which is barely visible at the end of the QG branch. As discussed below, the solutions near the turning point are always found to be unstable in the full QG dynamics.
To accurately describe the steady states near turning points or near the termination of solution branches, we turn next to the full QG dynamics. An overview of the range of steady state solutions found is presented in figure 5. First of all, the top row shows the effect of keeping z fixed at a moderately large value and varying x c . For both small and large x c , the vortices are nearly circular, as exhibited by their aspect ratio λ in figure 4. At intermediate x c , the vortices are most deformed, and here exhibit an e.g. shape in the full QG dynamics. The higher curvature at the outward tips of the vortices indicates the presence of stagnation points a little further out along the x axis (this has been confirmed by examining the streamfunction field). These stagnation points can collide with the vortex boundaries, as in figure 5(i) and (nearly) in figure 5(e), and explain why there are no further (singly-connected) equilibria for smaller or larger x c , respectively, in these cases.
For all values of z, circular solutions exist for x c = 0. As we show below, for a limited range, (approximately 0.57 < z < 0.75), non-circular co-rotating solutions can be found for small horizontal separations and again for large separations. For these vertical displacements, there exists a range of x c values wherein no steady-solutions are found. Nearly circular steady-states exist for all z for large horizontal separations. As x c decreases, a turning point occurs as the vortices elongate and deform into a non-elliptical shapes. Examples of steady states at the turning point in x c are shown in figures 5(f)-(i) for decreasing z. The limiting states change from outward to inward pointing corners as the vertical separation decreases. All states near and beyond the turning point have been found to be unstable (see below).
A complete picture of the steady states obtained using the full QG dynamics in the (x c , z) parameter space is provided in figure 6. In figure 6(a), we show where steady state solutions exist and indicate their deformation by color, while in figure 6(b), the stability of a number of these states is indicated (see next subsection for details). The horizontal lines in figure 6(a) shows families with constant z while the predominantly vertical lines show families with constant total angular impulse J. There is a conspicuous gap for small to moderate x c which first appears for z < 0.795. At right hand edge of this gap for For comparison, moment model predictions for co-rotating states are also shown in figure 6(a). The curves correspond to solutions, in (x c , z), ofφ 1 =φ 2 =θ for fixed values of λ, here 0.5 and 0.75. While mimicking the general shape of the steady-state diagram, unlike the full QG model, the moment model predicts the existence of co-rotating states for small x c for all values of z. In addition, on the outer, large x c boundary, the moment model produces steady-states at significantly smaller values of separation than the full dynamics.

Stability of co-rotating solutions
The stability of solutions is determined by tracking the nonlinear dynamics of numerically determined co-rotating states. Such states are subject to numerical perturbations which either grow or decay. The time evolution of a few selected unstable steady states is illustrated in figure 7 for the full QG equations. The images are arranged in order of increasing z from top to bottom. For the smallest vertical separation z = 0.2 in figure 7(a), we start with a limiting state with an inward pointing corner on each vortex. The tips subsequently elongate and partially wrap around the other vortex in its own horizontal plane. This causes each vortex to become significantly more circular, albeit with surprising complexity near each vortex edge. This complexity is characteristic of surface QG flows (Harvey andAmbaum 2010, Scott 2010) and is due to relatively strong shortrange interactions compared to 2D Euler flows. Despite this, the vortices settle into more circular forms with a fragment of the other vortex above or below it. Turning to z = 0.4 in figure 7(b), we see a similar (but much slower) evolution toward instability, and much less complexity in the final state. Again, a fragment of each vortex is found either above or below the other at late times. For z = 0.5 in figure 7(c), we see the evolution of an unstable state near a turning point in x c (0). In this case, each vortex is split into two substantial but unequal parts. The larger parts align near the center and the smaller parts are left orbiting around in a near tripole configuration. Many small, weak structures are found as well, but these have little impact on the final state which persists for at least 20 rotation periods. For z = 0.6 in figure 7(d), we show the evolution of another unstable state near a turning point in x c (0). The instability here is highly asymmetric leading to an aligned structure off-center and a large satellite in one layer only. This state also persists for long times. For z = 0.7 in figure 7(e), we start with a near limiting steady state with outward pointing corners. The instability here is relatively weak, leading first to a filament being shed from the outer tip of the blue vortex (and rolling up into many small-scale vortices), followed by the same event on the red vortex. The main vortices settle down to a slightly more aligned and less deformed state orbited by small-scale debris. Finally, for z = 0.8 in figure 7(f), a similar evolution is observed, albeit the instability is now stronger, with more substantial tongues of PV shed from the tips of each vortex. Again, the vortices settle down to a more aligned and less deformed state, despite the significant amount of debris found orbiting the vortices.

Discussion
We have studied the behavior of thin lens-like vortices with their potential vorticity concentrated on sheets in separate horizontal layers. This idealized situation reduces the three-dimensional dynamics of interacting vortices to effectively a two-dimensional problem. We have gone further here and developed approximate models for interacting lenses, taking each lens to be an elliptical distribution of PV density. Each lens must have non-uniform PV density to remain exactly elliptical in isolation (Dritschel 2011), but this poses no difficulty in deriving reduced models of interacting lenses. The reduced models, here a "moment" model and an "elliptical" model, both approximate the interaction between each pair of lenses. This interaction generally does not preserve the elliptical shape of each lens, but the error is small for sufficiently large vertical and/or horizontal separations, as verified by direct comparison with the full QG dynamics. In the moment model, both the vertical and the horizontal separation must be large compared with the semi-major axis length of either vortex. In the elliptical model only the three-dimensional distance between the vortex centers must be large. The elliptical model does not explicitly carry out an asymptotic expansion, but rather extracts the part of the interaction which preserves the elliptical shape of each lens. As a result, the elliptical model is substantially more accurate than the moment model, particularly for shorter-range interactions.
The simple problem of two identical lenses in two different layers was explored in depth. The reduced models were shown to accurately capture the vortex centers, aspect ratios and phases diagnosed from the full QG dynamics for sufficiently large separations. The models also indicated that, for particular parameter choices, steadily-rotating solutions ("steady states") may exist. This was confirmed through explicit calculation and comparison with the full QG dynamics, whose predictions closely match those of the reduced models for large separation distances. A complete survey of the possible steady states was carried out in the QG model, revealing a large gap in parameter space where no solutions are possible. Essentially, when the vertical separation falls below the mean radius of either vortex (approximately), states with moderate to small horizontal separations no longer exist. Either a turning point occurs in horizontal separation, or the vortices reach a limiting state characterized by an inward or an outward pointing corner on their boundaries.
The stability and nonlinear evolution of these states has also been addressed. We have shown that all solutions near the edge of the gap in parameter space are unstable. These instabilities manifest themselves in a variety of ways, from shedding thin tongues of PV and becoming more aligned and circular, to vortex splitting and subsequent partial alignment of the vortex parts. Vortex splitting can be either symmetric, resulting in a tripole state, or asymmetric, resulting in an aligned off-centered structure and a satellite in one layer only. In all cases, the late time states persist for long times.
In the future, it would be interesting to consider asymmetric configurations having vortices of different sizes and intensities. Additionally, a combined elliptical-ellipsoidal model should be developed to explore the interaction of a surface buoyancy anomaly (a lens with singular PV) with a three-dimensional interior vortex.