Instability of Lenticular Vortices: Results from Laboratory Experiments, Linear Stability Analysis and Numerical Simulations

: The instability of surface lenticular vortices is investigated using a comprehensive suite of laboratory experiments combined with numerical linear stability analysis as well as nonlinear numerical simulations in a two-layer Rotating Shallow Water model. The development of instabilities is discussed and compared between the different methods. The linear stability analysis allows for a clear description of the origin of the instability observed in both the laboratory experiments and numerical simulations. While global qualitative agreement is found, some discrepancies are observed and discussed. Our study highlights that the sensitivity of the instability outcome is related to the initial condition and the lower-layer ﬂow. The inhibition or even suppression of some unstable modes may be explained in terms of the lower-layer potential vorticity proﬁle.


Introduction
Mesoscale vortices are a major component of the global oceanic circulation. They are particularly common in oceanic regions of high mesoscale activity and can have long lifetimes, up to the order of a year (e.g., [1][2][3][4]). These coherent structures play an active role in the energetics of the ocean general circulation, in the transport of biological species, heat and salt anomalies, and in air-sea interaction (e.g., [5,6]). For example, eddies that detach from currents separating oceanic regions with different properties may propagate from one side to the other, transporting into a new region water with anomalous properties. These warm-and cold-core vortices are associated with the intersection of the isopycnals with the surface of the ocean, forming a front (e.g., [7]), but can also be sub-surface intensified [8] or form lenticular vortices at depth, a well known example of such vortices being the meddies [3]. Typical eddy size in the ocean are a few Rossby deformation radii (between 1 and 5 for eddies described in details in the literature), with Rossby numbers smaller than unity but finite (typically below 0.3 in absolute value) (e.g., [4,[9][10][11]). Smaller eddies can also be found, down to the submesoscale and with higher Rossby numbers (e.g., [12,13]). Warm-core surface anticyclonic vortices are the subject of this study and will be referred as lenticular vortices.
Early modelling studies of the dynamics of such vortices date back to the 1970s, starting with experimental studies [14,15]. Subsequently Griffiths and Linden [16] performed laboratory experiments in a rotating tank, using two different methods for generating the vortex. The first method (referred to as constant-flux) consists of slowly injecting a lighter fluid at the surface of the tank filled with heavier fluid until the desired volume for the vortex is obtained (this technique was already used by Gill et al. [15]). The second technique, the so-called constant-volume method previously used by Saunders [14], consists of releasing a given volume of lighter fluid initially contained in a suspended cylinder.
Under the action of buoyancy and rotation, the patch of lighter fluid adjusts towards an axisymmetric vortex which eventually destabilizes and breaks up.
The issue of the stability of oceanic surface vortices has also been addressed in theoretical and numerical studies. The (multi-layer) Rotating Shallow Water (RSW) model has been widely used to investigate their dynamics, representing the eddies as patches of lighter fluid overlying a deep ocean. Such model is appealing because it allows for an explicit representation of the density anomaly (compared to the exterior region) in the eddy, as well as the outcropping of the isothermal or isopycnal (e.g., [7,17]). This outcropping can be associated with specific dynamical features (including instabilities) associated with the frontal configuration [18][19][20]. In the two-layer configuration, Ripa [21] derived a semi-analytical criterion for the instability of vortices in solid-body rotation. He classified the instability regimes in terms of wave-wave resonances, showing that lenticular vortices with constant (zero in their case) Potential Vorticity (PV) over a quiescent layer of fluid are prone to instability through resonance between a Poincaré-like wave in the upper layer and a Rossby wave in the lower layer, for realistic oceanic eddy size and intensity. Paldor and Nof [22] investigated numerically the linear stability of zero-PV lenticular vortices and found unstable modes with azimuthal wavenumber m ≥ 2. More recently, Cohen et al. [23] extended this analysis to vortices with constant non-zero PV.
The fact that eddies with very long lifetimes are observed seems to conflict with results from these idealized studies, which found instability for typical oceanic regimes. Part of the explanation could be related to the flow surrounding the vortex, and in particular to the circulation below it. Dewar and Killworth [24] showed (although not in the context of lenticular vortices) that a co-rotating flow under an eddy decreases the growth rate of the instability, and Benilov [25] showed that a constant-PV lower-layer flow suppresses baroclinic instability in a QG model, the explanation being that Rossby wave propagation is precluded in the lower layer in such a configuration. In a two-layer RSW model, Cohen et al. [26] also reported no finding of unstable mode for lenticular vortices over a constant PV layer.
While the nonlinear development of the instability of isolated but non-outcropping vortices (i.e., vortices surrounded horizontally by a ring of opposite-signed vorticity) has been investigated numerically in numerous papers, either in layered or continuously stratified models (e.g., [27][28][29][30][31]), few simulations of the evolution of lenticular vortices have been reported to our knowledge. Verzicco et al. [32] reproduced the constant-volume experiments of Griffiths and Linden [16] and the subsequent evolution of the vortex for a few cases, using numerical simulation of the Boussinesq equations. Numerical simulations of the two-layer shallow water equations were performed by Thivolle-Cazat et al. [33]. They used the MICOM code in conjunction with data assimilation of laboratory experiments measurements (on the Coriolis platform in Grenoble, France). They did not focus on lenticular vortices, which represented only 5 runs out of 20 experiments. This rather scarce number of numerical studies is partly due to the difficulty of having a reliable numerical scheme for the shallow water equations that deals with the outcropping of one of the layers. On the other hand, fully three-dimensional numerical simulations still require large computational resources, making a large number of runs with different parameters at high resolution impractical.
In this paper, we report results on the instability of constant-PV surface lenticular vortices in a two-layer configuration. These results are based on laboratory experiments covering a wide range of vortex parameters (aspect ratio and PV), supplemented by extensive nonlinear numerical simulations of the two-layer RSW equations. These are compared to results from a linear stability analysis using the latter model. Besides allowing to capture outcropping and related frontal dynamics, the two-layer model is very close to the experimental setup used in this study, thus allowing for a more accurate comparison. A description of the model equations and methods, including the laboratory experiments setup, the nonlinear numerical simulation code and the linear stability analysis, are given in Section 2. In Section 3, we discuss the linear stage of the developing instability observed in both the nonlinear numerical simulations and laboratory experiments, and compare to results from the linear stability analysis. Then Section 4 addresses the nonlinear saturation of the instability, mostly in the numerical simulations, and contains a detailed analysis of one run illustrating the dynamics. A discussion on the impact of a lower layer flow is provided in Section 5, and conclusions are presented in Section 6. The two-layer RSW equations in polar coordinates (r, θ) are, e.g., [34] Dv i Dt

Problem Formulation and Methods
where i is the index of the layers (1 and 2 for the upper and lower layer, respectively), D/Dt = ∂ t + u∂ r + (v/r)∂ θ the Lagrangian derivative, v i = (u i , v i ) the horizontal velocity, p i the pressure and h i the layer thickness. This system of equations is supplemented by a dynamical condition at the interface: In the following, we will use either the rigid lid assumption, which implies h 1 + h 2 = H 0 with H 0 a constant corresponding to the nondimensional total thickness of fluid at rest (see below for the nondimensionalization); or a free-surface upper boundary condition. The rigid lid assumption allows for a simpler discussion of the formulation of the linear stability analysis and the spectra of unstable modes, while capturing the essential mechanisms at play in both the laboratory experiments and the nonlinear numerical simulation.Relaxing this approximation yields additional unstable modes that mostly take place for non-realistic dynamical parameters, as will be discussed separately. We shall consider here the stability of surface lenticular vortices and nondimensionalize the equations accordingly. We use the size of the eddy L e (i.e., the radius of outcropping of the interface) as the typical horizontal length scale and the inverse Coriolis parameter f −1 as a typical time scale. The velocity is non-dimensionalized by L e f (i.e., the nondimensional velocity is the Rossby number Ro = U/ f L). The thickness h i is non-dimensionalized with the value H sc associated with a baroclinic (or rather upper-layer) deformation radius equal to the size of the vortex: g H sc / f = L e , where g = λ 2 g is the reduced gravity. That is, the depth scale used for non-dimensionalization is H sc = f 2 L 2 e /g 2 and, correspondingly, the pressure scale is ρ i L 2 e f 2 . The barotropic Burger number of the vortices investigated is given by gH 0 / f 2 L 2 e and is large (greater than 10 for most experiments, the minimal value being 3), justifying the use of the rigid lid approximation for the linear stability analysis.

Vortex Profiles
Lenticular vortices are steady axisymmetric solutions of the equations with zero radial velocity and azimuthal velocity obeying cyclogeostrophic balance. They satisfy in each layer the following system of (nondimensional) equations: which is valid both for the free-surface and the rigid lid upper boundary conditions (the difference arising in the relationship between P i and H i ). In the notation used in this paper, upper-case variables stand for background state whereas lower case are for perturbations (and for some parameters). The family of upper-layer profiles we consider is defined by a constant value of potential vorticity Q 1 . This aims at reproducing the profiles obtained in a constant-volume experiment, where some fluid initially contained in a bottomless cylinder is released and adjusts through the combined action of Coriolis and buoyancy forces. The PV of the upper layer, initially contained in the cylinder, is constant with the initial thickness of the upper layer fluid contained in the cylinder) and nearly remains so throughout the adjustment and subsequent evolution. This assumption of constant upper PV is also relevant for typical oceanic profiles [7,9]. In the lower layer, we will primarily consider two different situations: V 2 = 0 (quiescent lower layer, with Q 2 = 1/H 2 ), which is the most trivial solution, and Q 2 piecewise constant. The latter configuration aims at capturing more accurately the profile obtained in constantvolume experiments. This method consists of releasing a volume of buoyant fluid initially contained in a suspended cylinder at the center of a rotating tank. The fluid initially below the cylinder in the lower layer has a positive PV anomaly, so that PV has a piecewiseconstant distribution. In the course of the adjustment, this volume of fluid slightly shrinks inward and acquires a weak cyclonic circulation. The lower layer PV is then given by where r * c is the location of fluid parcels that were initially at the radius of the cylinder before adjustment. With the rigid lid approximation, the momentum equations for the background state (4) in each layer reduce to: Equations (5) and (7) come with the boundary conditions H 1 (r = 1) = 0 and V i (r = 0) = 0, i = 1, 2. In the lower layer V 2 = 0 at r = 1, as mass conservation implies that fluid parcels in the lower layer at r = L e are unaffected during the adjustment. The resulting boundary-value problem is solved numerically with a standard boundary value problem solver, using a linear expansion for V i near the singularity at r = 0 of the form V i (r → 0) ∝ r. Examples of vortex profiles with different values of Q i are given in Figure 1.

Range of Parameters Investigated and Estimates of Realistic Values
Typical dynamical parameters for major ocean eddies, such as the Gulf Stream or Agulhas rings, were reported in (e.g., Olson [2], Goni et al. [11]), sometimes with an interpretation based on a two-layer formulation [10,35]. In the latter study, the authors reported typical values for two Agulhas rings of 400 m for the vertical displacement of the isopycnal, L e in the range 110 km to 130 km and maximum absolute azimuthal velocity in the range 0.6 m/s to 0.9 m/s. The deformation radius is of the order of 30 km in this area [36] (computing g h 0 f , with g = 0.018 m/s 2 as reported by the authors, gives the same value). The resulting Burger number ranges between 0.05 and 0.07, while the Rossby number is between around 0.04 and 0.07. These typical dynamical parameters were confirmed by observations from satellites in ths region [11], who reported typical eddy size ranging from 50 km to 170 km.
For a Gulf stream warm core ring (82B), the reported values were around 500 m for the isothermal displacement, V max around 0.5 m/s and L e ∼ 60 km. The deformation radius is around 20 km in this region, so the typical Burger number is around 0.11, as well as the typical Rossby number.
In our scaling, the potential vorticity (at the center of the vortex, and assuming it is constant) and is thus equal to the inverse Burger number in the limit of small Rossby number. Therefore, the nondimensional PV associated with the above-mentionned eddies ranges between 9 and 20. More generally, sizes of ocean rings range between 1 and 5 deformation radius [2], so that Bu ∼ [0.04, 1] and Q 1 ∼ [1,25]. Typical realistic values of δ do not exceed 0.3. Of course, smaller eddies with larger Rossby numbers are known to exist (e.g., [4,13]), and these would have smaller values of Q 1 (as well as δ).  , lower-layer PV Q 2 (lower left panel) and upper and lower layer azimuthal velocity V 1 and V 2 (upper and lower right panels, respectively). The vortex parameters are Q 1 = 10, δ = 0.2, and the solution with a quiescent lower layer (black) is compared to the solution with a piecewise-constant lower-layer PV (blue). The analytical solution for Q 1 = 0, V 2 = 0 is plotted as a grey dashed line for reference.

Setup of the Laboratory Experiments
Lenticular vortices were produced through the adjustment of a buoyant fluid initially confined in a bottomless cylinder in a rotating tank, the so-called constant-volume experimental setup [16,37]. The mechanical lift-up of the inner cylinder, which is co-aligned with the rotation axis of the tank, was triggered after the water in the tank had reached solid-body rotation. A schematic of the experimental setup is given in Figure 2.
The rotating tank was a square with side 99 cm and depth 51 cm. The rotation rate could be varied between 0.5 rad s −1 and 1.35 rad s −1 within an accuracy of 0.01%. We used two different cylinders with radius R c 5.5 cm and 20 cm. Dyed fresh water was used to visualize the evolution of the vortex using a CCD camera placed above the tank. The lower layer was salt water with concentration in the range 3.6 to 80 g L −1 , thus setting the value of the reduced gravity. The total fluid depth was varied from 8 to 40 cm while the thickness of lighter fluid in the cylinder was 3.5 ≤ h 0 ≤ 7 cm, covering a range of initial depth ratio 0.088 ≤ δ 0 ≤ 0.67. In total, 42 different experiments were performed, with initial Burger numbers θ 0 = g h 0 /( f R 2 c ) varying between 0.015 and 0.88. A table summarizing the parameters of the laboratory experiments is given in Appendix D. We relate these values to the nondimensional values of the potential vorticity Q 1 and aspect ratio δ. This is done using the semi-empirical fact that the vortex spreads out by approximately one baroclinic deformation radius, Rd = g h 0 / f , during the initial adjustment, as was noticed by Griffiths and Linden [16], so that L e ≈ R c (1 + √ θ 0 ). Thus, the nondimensional upper-layer PV becomes: For the depth ratio we use the qualitative formulas given by Griffiths and Linden [16]: Formulas (9) and (10) are used to infer the vortex parameters in the laboratory experiments, allowing for a comparison with the nonlinear numerical simulations and results from the linear stability analysis.

Linear Stability Analysis
For comparison with observations from laboratory experiments and nonlinear numerical simulations, we perform a linear stability analysis of the base flow described above. For simplicity, and since it does not affect the main results, we use the rigid-lid approximation here -extension to the free surface case is discussed in Appendix C. The nondimensional equations are linearized around the basic state using the standard normal mode decomposition (u i , v i , p i ) ∝ e i(mθ−ωt) , introducing an additional change of variable u → iu for convenience. This gives the linear eigenproblem for the eigenvalue ω defined for r ≤ 1: Boundary conditions for this problem are regularity of the solution at r = 0 and r = 1.
In the exterior domain (r > 1), the rigid lid condition enforces a non-divergent flowcf. (13)-with p 1 − p 2 → h 2 = 0. The flow in the exterior region can be obtained analytically, using a condition of decay as r → ∞. This leads to boundary conditions at r = 1 for the interior solution (see Appendix A for details): The resulting eigenproblem is solved using a collocation method with Chebyshev polynomials. Regularity of the solution is ensured by the form of the interpolants. The azimuthal momentum and mass conservation equations are multiplied by r for the numerical implementation. Numerical convergence in the case of an eddy over a quiescent lower layer is reached for resolutions as low as 10 points.

Nonlinear Simulations: Numerics and Initial Conditions
For the nonlinear simulations, we used a high-resolution finite-volume numerical code [38]. It integrates the two-layer RSW Equations (1) and (2) with a free surface in conservative form, using a well-balanced entropy-satisfying scheme. The latter ensures physically appropriate behavior in cases of shock formation (namely a decrease in energy) without requiring explicit dissipation. In this regard, the numerical code belongs to the family of implicit Large Eddy Simulation (iLES). In particular, it copes well with outcropping (vanishing layer thickness), which is required for the simulation of lenticular vortices. This code was previously used in the context of frontal dynamics by, e.g., Gula et al. [18], Ribstein and Zeitlin [19], Gula and Zeitlin [39]. The numerical simulations are initialized using the solutions for the lenticular vortices considered in Section 2.1.2, with a quiescent lower layer, adapted to the free surface condition and suitably perturbed to trigger the instability. The perturbation is obtain by applying a θ-dependent map to the radial variable when computing the solution H i (r), V 1 (r), as described in Appendix B. A weak radial velocity field is then computed so that the divergence of the velocity vanishes, in order to limit the generation of ageostrophic motions and especially of inertia-gravity waves. The azimuthal wavenumbers of the perturbation range from 2 to 10 (unless otherwise stated), and the perturbation's relative magnitude is of the order of a few percent. A series of numerical experiments were carried out to span the δ, Q 1 plane. Values used in the initial conditions are given in Table 1. The size of the domain is four time the diameter of the vortices, and the numerical grid has 400 2 points.

Early Stage of the Instability
In this section, we report the behavior of the vortices at early stages of the instability from the laboratory experiments and nonlinear simulation and compare it with the results from the linear stability analysis.

Preliminary Considerations from Linear Stability Analysis
We first briefly discuss the results from the linear stability analysis, neglecting the background current (V 2 = 0) for simplicity. A diagram of the wavenumber m and growth rate σ = (ω) of the most unstable mode in the (δ, Q 1 ) plane is given in Figure 3. An extensive discussion of the linear stability of these kind of vortices can be found in (Cohen et al. [23], their "regular" regime).
Instability is found for large values of the aspect ratio δ and/or the potential vorticity Q 1 , with the wavenumber of the most unstable mode increasing with both of these parameters, while the growth rate depends predominantly on δ. The latter dependence is in agreement with what was found for non-outcropping eddies in two-layer models (e.g., [24,40]). From the wave resonance point of view, the "standard" unstable modes (in the upper and rightmost part of the diagram) result from the interaction between a frontal-type wave in the upper layer, sometimes referred to as a mixed Rossby-Poincaré or Yanai mode, and a Rossby wave in the lower layer [18,21,23,39]. The latter is supported by the presence of the PV gradient associated with the slope of the interface since the lower-layer PV is given by Q 2 = 1/H 2 . The growth rate thus depends on the efficiency of the mutual amplification between the upper-and lower-layer waves, which is primarily driven by their spatial structure (the shape of the modes) and their phase speed difference, as phase-locking is a necessary condition for mutual amplification. A review of two-layer and frontal instabilities is given in Zeitlin [34].
There is a zone of relative stability in the lower-left corner of the diagram, delimited by the m = 2 contour. In this region, however, small zones of instability also exist with rather high (4 < m < 10) wavenumbers and growth rates not exceeding 2 × 10 −2 f −1 . These instabilities involve higher radial wavenumber Poincaré modes in the upper layer. The latter have different phase speed, therefore experiencing phase locking for different (higher) values of the azimuthal wavenumber compared to the unstable modes described above. Finally, one sees a weak region of very high growth rate (up to 0.4 f −1 ) with wavenumber m = 10 for very large depth ratio (δ > 0.85) and weak Q 1 , which we identify as a frontal-frontal (i.e., implying a Rossby-Poincaré wave in the lower layer instead of a Rossby wave) instability [19].  . Growth rate (colours) and wavenumber (black contours) of the most unstable mode, as found from the linear stability analysis, in the (δ, Q 1 ) parameter space, for lenticular vortices with a quiescent lower layer (V 2 = 0). Following the discussion of realistic parameters (see Section 2.1.3), typical values for large ocean rings are log(Q 1 ) ∈ [0, 1.5] and log(δ) < −0.5.

Observations in Laboratory Experiments
Just after the withdrawal of the cylinder, the freshwater spreads radially and the flow adjusts to an axisymmetric lenticular vortex. The detailed dynamics of this initial adjustment has been carefully discussed in several papers [16,33,37] and is not the subject of the current paper. One should keep in mind that dissipative processes and mixing occur, especially close to the edge of the vortices, diminishing the azimuthal velocity and reducing the shear between the two layers there [37]. This is almost certainly the main source of differences between the experimental adjusted lenticular vortex profile and the inviscid solution. After this period, the flow eventually starts to develop non-axisymmetric features that consist of m arms developing out of the initial vortex and rolling up in a cyclonic sense, forming dipolar structures that can eventually detach and propagate away.
Snapshots of different stages of the initial adjustment and the development of an instability with m = 6 are shown in Figure 4. The parameters of this eddy (see table in Appendix D) are δ ≈ 0.17, Q 1 ≈ 58, Bu ≈ 0.02. That is a very large eddy (around 7 deformation radii) with a Rossby number of order 0.1. Here and below, when restoring dimensional values, we will use Rd ∼ 25 km and f ∼ 10 −4 s −1 , which is fairly typical of mid-latitudes. Using these, the maximum velocity is 1.75 m/s, which is large (but one should keep mind that all these parameters are estimated within at least a factor of 2 uncertainty). In the laboratory experiments, the number of arms is often well defined and their development rates are of the same order. Sometimes, at high-wavenumber, one or two arms may have an inhibited growth compared to the others (see Figure 4d). The wavenumber of the developing unstable mode, as well as the timescale of the growth of the instability, are different between the various experiments (2 ≤ m ≤ 7 for the investigated range of parameters). A comparison of the prediction from the linear stability analysis with the experimental observations is given in Figure 5 (left panel) and shows good agreement. Possible reasons for the observed discrepancies are the following: • The parameters δ, Q 1 in the laboratory experiments are deduced from qualitative considerations and semi-empirical relations, neglecting the details of the initial adjustment; • Friction effects decreases the azimuthal velocity, especially during the initial adjustment and near the edge of the eddy [37]. • The initial adjustment leaves out high-wavenumber perturbations to the vortex before the instability starts developing (see Figure 4a). The vortex may be unstable with respect to several wavenumbers, which are unevenly excited by this perturbation, so that a mode that is not the most unstable may initially gain more energy than the most unstable one and develop faster. Another possible source of differences is that the lower-layer flow is neglected in the linear stability analysis. However, as will be shown in Section 5, this is very unlikely to be the dominant cause of discrepancies.

Results from Nonlinear Numerical Simulations
Our numerical simulation reproduces the observations reported previously from the laboratory experiments, as far as the development of the instability is concerned, at least qualitatively. We compute a growth rate associated with the development of the instability based on the perturbation of the upper layer thickness with respect to the angular mean at the initial time. A clear period of linear growth is observed in every simulation. Figure 5 (right panel) shows a comparison of the growth rate calculated in the numerical simulations and that predicted by the linear stability analysis here adapted to the free surface for a more rigorous comparison. Note that, for the range of (δ, Q 1 ) investigated, the impact of the free surface compared to a rigid lid approximation is negligible (see Appendix C for an extended discussion).
Rather good agreement is found in most cases. In particular, no instability occurs for the two runs outside the linear instability zone and the growth rate of the most unstable mode increases with the depth ratio and to some extent with the upper layer PV. Nonetheless, the growth rate measured in numerical simulations is systematically smaller that the linear prediction, with a relative error ranging between 10% and 40% and roughly increasing with Q 1 . Several explanations are possible for this discrepancy. It is most likely because of numerical dissipation near the edge of the vortex: although the numerical code does not use explicit viscosity, the entropy condition used in the Riemann solver is associated with a decrease of the energy across shocks. Hence, energy dissipation occurs at the edge of the lens. In addition, there is a non-negligible time before the unstable mode starts developing, because the initial amplitude of the perturbation is weak and the shape of the perturbation is different from that of the unstable modes. Short-time non-normal modes are very likely to be excited by such a perturbation, before the long-time unstable modes start developing and dominate the flow (e.g., [41]), so that the angular mean profile changes before the unstable mode starts developing.
The development of the instability is investigated in a numerical experiment where we used the linear unstable modes as a perturbation to the initial state. This procedure has been routinely used (e.g., [30,39]) and enables a clear investigation of the dynamics associated with the unstable mode in isolation from other phenomena. The vortex parameters are δ = 0.2 and Q 1 = 12, corresponding to Ro ∼ 0.25 and Bu ∼ 0.1, i.e., it has a typical size of roughly 3 deformation radius. The dimensional maximum velocity reaches 1.8 m/s, which is roughly twice the observed typical maximum velocity in ocean rings. The pressure and velocity perturbations associated with the developing unstable mode for this vortex initialized with an m = 3 perturbation in the upper layer are shown at t = 40 f −1 , roughly after twice the growth timescale, in Figure 6. (Note that here and below, for clarity, fields in the upper layer are only represented where the layer thickness is greater than H 0 /10 3 ). The standard picture of a frontal mode in the upper layer and Rossby mode in the lower layer (with the velocity along the pressure isolines) clearly emerges, and matches the structure of the linear unstable mode (not shown, but one can compare with Gula and Zeitlin [39] in the rotating annulus): the m = 3 growing perturbation manifests itself by the non-axisymmetric shape of the vortex in the left panel.

Sensitivity to the Initial Conditions
When initially perturbed with several wavenumbers, several unstable modes with differ'ent wavenumbers may develop, as might be expected from the linear stability analysis. This is visible in Figure 7 (left panel), where the evolution of the modal amplitude of the perturbation in one numerical run is plotted and exhibits linear growth for two wavenumbers. It is different from what is observed in the laboratory experiments, where a single wavenumber most often emerges. We suspected this could be due to the difference in the initial perturbation. We thus consider the destabilization of the eddy with δ = 0.2 and Q 1 = 12 discussed previously, in five nonlinear simulations using 3 different types of initial perturbation. The first consists of a random noise with different azimuthal wavenumbers; the second consists of noise with the same radial structure as before but only one wavenumber (we consider m = 2 and m = 3); the last one consists of the linearly unstable mode, as shown previously (two runs). The basic vortex profile is linearly unstable with respect to wavenumbers 2 and 3, and the associated growth rates are σ 2 = 4.3 × 10 −2 f and σ 3 = 3.4 × 10 −2 f , respectively. In the nonlinear simulations, we find growth rates of 3.8 × 10 −2 and 2.9 × 10 −2 for m = 2 and m = 3, respectively, both in the simulation initialized with several wavenumbers and with a single wavenumber, as shown in Figure 7. The development of this instability is associated with a transfer from potential to kinetic energy. The evolution of the relative energy anomaly, with the potential energy on the same domain with the lens subtracted, is shown in Figure 7 (right panel) for the three different runs with different initial perturbation (multiple wavenumbers, m = 2 and m = 3). The energy transfer starts at the initial instant, as seen from the slight but constant decrease in potential energy. However, there is no visible increase in the kinetic energy at this time, probably due to dissipation. Nonlinear saturation, as seen in the decrease of the linear growth rate (cf. Figure 7), occurs around t = 100 f −1 here and is associated with an enhancement of the energy transfer, which stops later on (at about t = 200 f −1 depending on the simulation). During the same time, the dipoles form and start moving away from the initial vortex position. For the m = 2 instability, t ≈ 100 f −1 corresponds to when the vortex becomes pinched near its center, while it is associated with the arms beginning to roll up for the m = 3 instability. Rather surprisingly, the development of the instability, as seen in the evolution of the energy, is actually slower for the simulation perturbed with several wavenumbers compared to the m = 2 perturbation alone associated with the strongest transfer. To our knowledge, this kind of sensitivity to the initial perturbations has not been reported before.

Nonlinear Saturation of the Instability
We now turn to the mature stage and nonlinear saturation of the instability, relying mostly on nonlinear numerical simulations and comparing with laboratory experiments when appropriate. The range of parameters investigated in the former is the same as in the latter (see Table 1). In addition, we performed a set of simulations with Q 1 = 12 and δ = 0.2 already discussed in the previous section. The latter serves as a reference for the discussion of the nonlinear saturation of the instability, and is presented in detail. The outcome of the nonlinear saturation of the instability is not unique: it may be divided into 3 different scenarios, depending on whether the dominant unstable mode consists of (1) one wavenumber m = 2, (2) one wavenumber m > 2, or (3) a superposition of several wavenumbers. For the same initial vortex, one of the above situation can be picked by choosing the proper initial perturbation, provided the vortex is unstable with respect to several wavenumbers with nearby growth rates. This is particularly true (see previous section) for the vortex with δ = 0.2, Q 1 = 12 described below.

Nonlinear Saturation for m = 2 Dominant
The destabilization of the vortex in this case is similar to the standard "dipolar vortex breaking" observed as the outcome of barotropic instability (e.g., [29,30], and references therein). Developing arms start rolling up in a cyclonic sense (although the circulation within the arms is anticyclonic) and form a dipole with a lower-layer cyclonic circulation associated with the positive PV anomaly shredding from the initial core. The destabilized vortex is shown in Figure 8 (first row) at two different stages ,and the corresponding lower-layer PV is shown in Figure 9 below. The lower layer is nearly balanced (in the geostrophic sense) and circulation may thus be inferred from the pressure gradient. The self-mutual interaction of the poles is baroclinic, as each pole is in a different layer, and is more heton-like [29,42]. As the dipoles radiate away, they detach and eventually leave a monopole at the location of the initial vortex, which looks like an elliptical "cat's eye" vortex. In some case, weaker satellites result from the filaments between the radiating dipole and the monopole (on the right above the monopole and below it on the left in Figure 8, upper-right panel). The magnitude of this monopole is larger for larger values of δ, and is almost negligible for small values (within the parameter range investigated in our numerical simulations and leading to a wavenumber 2 unstable mode). It is worth mentioning that the radial momentum of the dipoles almost vanishes once they detach from central monopole. Almost all of the lower-layer initial patch of PV splits off and moves along with the baroclinic dipoles, or is stirred during of the detachment of the latter, so that the PV anomaly under the remaining elliptical monopole is very weak, enhancing its stability (see Section 5.2 for a discussion of the impact of a vanishing lower-layer PV gradient). This concludes the nonlinear saturation of the instability, as a marginally stable state is thus reached. The situation shown at t = 260 f −1 in Figure 8

Nonlinear Saturation with Wavenumber(s) m > 2
Cyclonic arms also develop for vortices that are unstable with respect to one or several mode number(s) larger than 2. However, they hardly ever detach from the initial vortex nor propagate away from it. Rather, as can be seen in the second and last row of Figure 8, the baroclinic dipoles formed rotate on their own and stay in the same neighbourhood. At later times, they adjust toward monopoles (as in Figure 8), or eventually remerge with the main anticyclone.
It is interesting to note the differences between the development and outcome of the destabilization when several wavenumbers develop simultaneously. In particular, the arms often detach at different times. In the case given in Figure 8, at t = 130 f −1 , only one dipole in the process of forming is visible, whereas there are as many as the wavenumber of the initial perturbation in the two upper rows. At t = 260 f −1 , one can see that another dipole has detached (to the right of the main anticyclone) and additional, much weaker patches are visible, e.g., to the left of the main anticyclone. Again, this is representative of many cases.
The saturation of the instability in the numerics is quantitatively different from what is observed in the laboratory experiments. In the latter, the initial vortex is always completely destroyed (meaning that all the lighter fluid initially in the eddy went into the dipoles), while we do observe a remaining (somewhat weak) monopole at the initial location of the eddy in the numerics. This sheds light on the uncertainties associated with both numerical and laboratory experiment, and especially on the initialization of the flow, as discussed throughout the paper.
Although a detailed comparison of the different models is beyond the scope of this manuscript, the scenario of the breakup of the vortex is qualitatively similar to nonoutcropping QG or RSW vortices with in either a gaussian shape or a piecewise-constant PV distribution (e.g., [24,31,40]).

Discussion: Impact of a Lower-Layer Flow
To date, in our nonlinear numerical simulations and linear stability analysis, we have only considered basic states with the lower layer at rest. In this case, the lower-layer PV is given by Q 2 = 1/H 2 , and its radial variation allows for the propagation of a Rossby wave that can resonate a upper-layer wave to give rise to an instability. This is highlighted in Figure 9 which shows the evolution of the lower PV (Note that, due to discretisation errors in the vicinity of the outcropping region, PV is not exactly conserved and develops small sharp perturbations, including in the lower layer. For visualization purpose, we smoothed the computed PV field using a Gaussian kernel of width equal to 4 grid cells.) during the development of the instability for the vortex with Q 1 = 12, δ = 0.2 previously shown in Figure 8. The role of the lower-layer PV in the instability is obvious from its phase quadrature with the upper layer at early stages of the instability (best visible in simulations initialized with well-defined azimuthal mode number-see the left and center upper panels), typical of baroclinic instability.
However, the actual inviscid solutions associated with the experimental configuration we used are different. In particular, the lower-layer PV has a piecewise-constant radial profile, as shown in Figure 10. It is also different in the other standard experimental method for producing lenticular vortices, the "constant-flux" experiment [16]. In the latter configuration, where fluid is gradually injected from a small source at the surface, conservation of the lower-layer PV implies that its value is constant and the same as the planetary PV outside the eddy. While fluid is being injected, fluid parcels in the lower layer are pushed radially outward and thereby acquire negative relative velocity. As a result, the lower layer is co-rotating with the upper layer (cf. Figure 10, lower right panel, orange line). (An approximate analytic solution for this configuration, neglecting the centrifugal term in the lower layer, was derived by Griffiths and Linden [16].) On the contrary, as discussed in Section 2.1, a weak counter-rotating circulation forms in the constant-volume experiment.
A solution with a constant lower-layer PV is a plausible configuration for realistic ocean eddies. Indeed, deep ocean flows are weaker than surface currents and thereby associated with very weak departure of the PV from the background planetary vorticity. Thus, as an eddy propagates above this quiescent deep layer, it will slightly moves the isopycnal down which will implies the development of a weak co-rotating flow in the deep layer, through PV conservation. Although very qualitative, this is somewhat in agreement with the deep extension of eddies reported from various observations (e.g., [2,8,9]). Based on this remark, one can argue that constant-PV lower layer is at least as plausible as a lower-layer at rest.
One would expect that the vanishing of the PV gradient and the associated suppression of Rossby waves would remove the instability, as was shown for nonisolated eddies [25,43] and more recently for warm-core eddies in the RSW model [26]. We thus consider lenticular vortices with profiles corresponding to the two situations mentioned, which are shown in Figure 10.

Piecewise-Constant Lower-Layer PV
The jump in Q 2 at r = r * c yields a cusp in the velocity profile, as follows from Equation (5) (see Figure 10 around r = 0.6). As spectral methods do not cope well with discontinuities, we smooth the profile of Q 2 (and, correspondingly, the cusp in V 2 ) using a third-order polynomial that matches both the values of Q 2 and its derivative (which is zero) at r = r * c ± 0.05 for the numerical resolution of the linear eigenproblem. The linear stability calculation is performed over a restricted range of values for (δ, Q 1 ): 0.01 < δ < 0.9 and 1 < Q 1 < 80.
The instability regime we found with this configuration is marginally changed by the lower layer flow compared to the quiescent lower layer. Therefore, it is not shown here, and we refer the reader to Figure 3. The most striking difference is that the growth rate of the instability vanishes for Q 1 < 5 (log(Q 1 ) < 0.7). For Q 1 > 5, the growth rate is very similar for low values of δ, which is not surprising as the lower-layer flow remains weak. For larger values of the depth ratio (δ > 0.1), the growth rate for piecewise-constant Q 2 is lower than that for V 2 = 0, the more so as δ increases. We found weakly unstable modes. (σ < 2 × 10 −2 for δ 0 < 0.9) for log(δ) > −0.5, which have critical levels associated with the lower layer flow. However, their physical relevance is unclear and the basic flow in this region tends to be unrealistic, as the lower layer velocity peaks at values similar (in absolute value) to the upper layer flow and the initial value of the depth ratio (before the adjustment) is very high: in particular, neglecting viscous effects in the lower layer is not appropriate in this configuration.
We also tried running four nonlinear simulation initialized with this modified vortex profile for δ = 0.2 and 0.4, and Q 1 = 12 and 25. These simulations did not exhibit significant difference compared to their counterpart with V 2 = 0, confirming the previous results from the linear stability analysis.

Constant Lower-Layer PV
The equations for the background state are unchanged, except that the lower layer PV is now Q 2 = 1/H 0 , and the boundary condition for the lower-layer velocity is different: The only solution that reproduces the inviscid "constant flux" experiment adjusted state has Q 1 = 0, which is associated with solid-body rotation V 1 = −r/2 (e.g., [16]). We investigated, however, a broader range of values for Q 1 , considering that eddies produced in oceanic flows, e.g., by detaching from boundary currents, can have non-zero PV values.
From the linear stability analysis, we found no unstable mode in the range of parameters investigated: 5 × 10 −2 < δ < 0.8 and 0 < Q 1 < 40, except for the highest values of δ and low values of Q 1 where the frontal-frontal instability occurs (not shown, but see Figure 3, lower right corner). This result is in agreement with the fact that the absence of gradient of PV in the lower layer prevents the propagation of a Rossby mode, thus suppressing the instability resulting from its resonance with the upper-layer modes. This also confirms some results previously reported in the literature [24][25][26]. It is worth mentioning that we still found ageostrophic instability for very high values of δ and wavenumbers m ≥ 10, which was not reported in Cohen et al. [26] as they did not investigate the corresponding range of parameters (note that δ > 0.8 is unrealistic in oceanic flows, except over shallow regions such as continental shelves).
We verified the inhibition of the instability by the constant lower PV was verified by means of nonlinear simulations for two sets of parameters: δ = 0.2, Q 1 = 12 and δ = 0.4, Q 1 = 20 (not shown). The typical dimensional parameters of the first eddy were already given above, and the second case corresponds to a fairly large (Bu ∼ 0.03) and intense (Ro ∼ 0.3) vortex. These simulations shown that the vortices eventually destabilized, but only after a significant part of the lower layer velocity sustaining the constant PV dissipated, thereby restoring a PV anomaly in the lower layer. Once a significant lowerlayer PV was reached (around half the magnitude of the lower PV in the quiescent lower layer configuration), the eddy started destabilizing, with a typical timescale that was clearly extended (by at least a factor 2) compared to the simulations with V 2 = 0. These results tend to confirm the stabilization of the vortex by the constant lower-layer PV. To obtain more insight about the impact of the lower-layer flow, we investigate the linear stability of profiles that lie in between the two extreme configurations for the lower layer: (1) no flow (V 2 = 0, Q 2 = 1/H 2 ); and (2) constant PV (V 2 > 0, Q 2 = 1/H 0 ). We consider a class of solutions with Q 1 = 0 with the lower layer velocity proportional to the constant-PV solution with a factor α (as in, [43]). The layer thickness are recomputed afterwards by solving (7), with δ held constant. Note that the value of H 0 , and hence the dimensional size of the vortex L e , vary with α. (Other "paths" for the evolution of the vortex profile can be used, e.g., by keeping the dimensional volume of the lens of fluid constant, which would probably be a better approximation to experimental flows. However, we do not aim at modelling precisely the evolution of the eddy under bottom friction.) The evolution of the growth rate as a function of α is given in Figure 11 for two different values of δ. One sees that the growth rate of all the unstable modes vanishes continuously as the lower layer flow evolves towards a zero PV-gradient solution.

Conclusions
In this study, we investigated the instability of constant-PV warm-core eddies in a two-layer Rotating Shallow Water model, by means of a linear stability analysis, laboratory experiments and nonlinear numerical simulations. We now summarize the results of this work.
A linear stability analysis was first carried and a stability diagram was constructed. It shows the dependence of the growth rate and wavenumber of the most unstable mode with respect to the depth ratio δ of the eddy and the value Q 1 of the upper Potential Vorticity (PV), for an eddy over a quiescent lower layer and with the rigid-lid approximation. We confirmed the following results, previously reported in the litterature ( [19,23], and references therein): • The instability domain is dominated by the hybrid instability, associated with resonance between a frontal mode and a lower-layer Rossby wave. • The growth rate of the most unstable mode increases with the depth ratio and slightly decreases (increases) with the value of the PV of the eddy for small (large) values of the depth ratio. The wavenumber of the most unstable mode increases with both depth ratio and the eddy PV. Several unstable modes with close growth rates and wavenumbers co-exist for m > 2. We also reported some results that were not discussed in former studies, to our knowledge: • Some unstable modes may be found in the zone of stability in the lower-left corner of the δ, Q 1 parameter space. • For very large aspect ratio, we found a very unstable ageostrophic unstable mode with high wavenumber (m ≈ 10), associated with the resonance between Poincarelike waves.
Moving to a linear stability analysis, we then investigated the impact of the lower flow: • If the lower-layer PV is constant (relevant for oceanic application and "constant-flux" experiments, and associated with a weak co-rotating lower-layer circulation), the hybrid instability vanishes as no lower-layer PV gradient supporting the Rossby wave propagation is present. As a result, the vortex is found to be stable, as was previously suggested (e.g., [24,25]) and confirmed by a linear stability analysis in the same model by Cohen et al. [26] with a different numerical method. We further confirmed this inhibition of the instability by means of nonlinear numerical simulations. The ageostrophic instability for very high depth ratio persists (from the linear stability analysis), which was not previously reported.

•
If the lower layer is piecewise-constant (corresponding to the inviscid adjusted state in "constant-volume" laboratory experiments and associated with a weak counterrotating lower-layer circulation), the stability properties are hardly changed for Q 1 > 5, as shown by both the linear stability analysis and the nonlinear simulations. The growth rate is slightly lower for large depth ratio. For Q 1 < 5, no significant unstable mode has been found, but weakly unstable mode may exist for high values of the depth ratio.
The destabilization of the vortex was studied by means of "constant-volume" laboratory experiments (release of lighter fluid initially contained in a cylinder at the center of a rotating tank filled with denser fluid) and numerically using a high-resolution finitevolume scheme of the two-layer RSW equations. For a wide range of vortex parameters, the wavenumber of the unstable mode observed in the laboratory experiments, as well as the growth rate measured in the numerical simulations, showed good agreement with the linear theory prediction. However, from a more quantitative point of view, some discrepancies were observed: • The wavenumber observed in laboratory experiments was sometimes different than the linear prediction, especially for large values of the depth ratio; • The growth rate in the numerical experiments is lower that expected from the linear stability analysis (between 10% and 40%, increasing with Q 1 ).
These differences could be caused by friction effects: the Ekman damping in the lower layer (laboratory experiments) becomes important for large values of δ, while friction near the edge of the eddy can have an impact, especially for large Q 1 where the velocity profile peaks at the edge. More importantly, we observed that the wavenumber of the developing unstable mode is sensitive to the initial perturbation. In the numerical simulations, while single-valued high azimuthal wavenumber perturbation is not enough to trigger instability, a superposition of several wavenumbers leads to the concurrent development of several unstable modes (if present). On the other hand, instability seems to be fed by highwavenumber perturbation in the laboratory experiments although, in most cases, it exhibits a single well-defined wavenumber at mature stage.
The nonlinear saturation of the instability, afterwards the initial development of (typically m) arms rolling-up in a cyclonic sense, can be classified as follows: • for m = 2, two baroclinic dipoles form and propagate away from the initial eddy, leaving eventually a weak monopole at the center; • for m > 2, several dipoles propagate with a strongly curved direction. As a result, they sometimes stay in the neighbourhood of the initial eddy, interacting in a complex manner with each other. Again, a weak monopole can be left at the location of the initial location of the eddy, depending on the parameters.
Here again, differences are observed between the numerical simulations and the laboratory experiments. In particular, the destabilization of the vortex in the numerical simulations appears more chaotic when several unstable modes (with different wavenumbers) develop, as they nonlinearly interact with each other. Moreover, it is very rare to observe a leftover monopole at the center of the domain in the laboratory experiments, as almost all the initial lighter fluid moved away inside the baroclinic dipoles. Although it is difficult to identify the origin of these discrepancies, one should bear in mind that dynamical parameters of the eddies produced in the laboratory experiment are estimated from simple heuristic arguments. As a result, for instance, the spectra of unstable modes associated with the actual eddy produced could be different resulting in one mode being dominant over the others, which would lead to a "simpler"development of the instability as observed.
The results presented in this paper strengthen the need to investigate the impact of the environment of eddies to better understand their instability properties. In particular, we found that the lower-layer flow may act to suppress the instability. In this direction, investigating the role of smaller-scale turbulence (e.g., random fluctuations in the lower layer) could provide very interesting information regarding the stability of ocean eddies. Hence, upon insertion of this form in (A4), the streamfunction satisfies a Laplace equation with the condition ψ → 0 as r → ∞. This gives the following relation: which is used as a boundary condition at r = 1 to ensure continuity of the eigenfunction across r = 1. Finally, (A1) and (A2) can be combined to express u and v as functions of p and its first derivative: Using u = v, these equations yield a relation between p and its first radial derivative which provides the inhomogeneous linear boundary condition (14), again by imposing continuity of p and dp/dr at r = 1.
(δ, Q 1 domain was shown). The relative difference between the free-surface and rigid-lid, for a value of the density ratio of 0.99, is less that 10% for Q 1 > 5 and of the order of 15% for lower Q 1 . This error associated with the rigid lid approximation increases with δ and decreases with Q 1 . The only remarkable emerging features are frontal-frontal modes (e.g., [19]) that manifest themselves emerging from the bottom of the figure for very deep eddies (around log(δ) = −0.2). Table A1. Parameters of the laboratory experiments: experimental configuration are in the first 7 columns, estimated nondimensional parameters of the adjusted vortex in the next 4 columns, and observed most unstable mode number in last column. Symbols are explained in the text. Bold values are for the experiment shown in Figure 4.