Jovian vortices and jets

We explore the conditions required for isolated vortices to exist in sheared zonal flows and the stability of the underlying zonal winds. This is done using the standard 2-layer quasigeostrophic model with the lower layer depth becoming infinite; however, this model differs from the usual layer model because the lower layer is not assumed to be motionless but has a steady configuration of alternating zonal flows [1]. Steady state vortices are obtained by a simulated annealing computational method introduced in [2], generalized and applied in [3] in fluid flow, and used in the context of magnetohydrodynamics in [4-6]. Various cases of vortices with a constant potential vorticity anomaly atop zonal winds and the stability of the underlying winds are considered using a mix of computational and analytical techniques.


Introduction
A great red spot on Jupiter has been observed for centuries, with its present manifestation dating back to telescopic observation in 1830. With the technological observational advancements and spacecraft measurements of modern times, additional Jovian vortical features have been discovered as well as such features in other planets. Many theoretical ideas have been proposed, yet it seems at minimum the red spot is a vortex enmeshed in zonal flow.
In this paper we investigate the conditions required for isolated vortices to exist in sheared zonal flows and the stability conditions for the maintenance of zonal winds. To this end we use the simple model of [7] to illustrate the basic concepts and to explore the sizes and shapes of the vortices. The model can be viewed as a version of the standard 2-layer quasigeostrophic (QG) model with the lower layer depth becoming infinite. Unlike the more common 1 1 2 -layer model, the lower layer is not assumed to be motionless but rather have a steady zonal flow, U (y). Since vortex stretching produced in the lower layer by vertical movement of the interface is negligible, this flow can remain unchanged. We will call this a 1 3 4 -layer model to distinguish it from either the motionless deep layer case or the one in which the deep layer evolves.
In Section 2 we will review how the 1 3 4 -layer model emerges from the 2-layer model. Also in this section we review the noncanonical Hamiltonian formalism [8][9][10] and briefly describe the Dirac bracket formalism, a Hamiltonian technique for the imposition of constraints. Both formalisms will be used later when we describe the simulated annealing (SA) procedure for obtaining steady states [2] and our generalization, the Dirac bracket simulated annealing (DBSA) procedure, introduced in [3].
Since the jets we consider can have β − U yy changing sign, by the Rayleigh criterion they may be subject to rapid instability and break-down in either the 1 1 2 -layer model or standard 2-layer model. We will re-examine the idea that steady deep flow can stabilize the jets [1] using the energy-Casimir method and relate it to the existence of isolated vortices in Section 3. Furthermore, in [11] it was shown that the "two-beta" model reduces the growth rates significantly and limits the unstable waves to small scales. In this system, the deep fluid has a strong reverse β effect because of large vertical extent of the flows in the direction parallel to the rotation vector, which is appropriate when the entropy gradients are small [12].
In Section 4 we describe steady states composed of localized vorticity anomalies, vortex patches, embedded in the layer zonal flows. From the equations for the dynamics relative to the jets we obtain integral conditions that are useful for identifying the allowed positions of the anomalies. Importantly, we obtain an artificial dynamics that is adaptable to a contour dynamics version of the simulated annealing technique [3], which we will subsequently use to explore vortices in sinusoidal shear flows.
Section 5 begins with a description of the Hamiltonian structure for contour dynamics in the usual incompressible two-dimensional Euler context. This structure, which applies to nonsingle-valued contours, was described briefly in [10], but is due to two of us (GRF and PJM) and has been used and described in talks over the past 25+ years. Next we describe the DBSA technique in the context of contour dynamics, then, as a warmup, show how to use it to construct the Kirchhoff ellipse and its generalization with a background shear. Lastly in Section 5, we describe constraints needed for application of DBSA to Jovian vortices.
In Section 6 a variety of Jovian vortices on jets are constructed by the contour dynamics DBSA technique with Dirac constraints chosen to be the linear momenta. The case with β = 0 and a depth independent background flow is first considered. The parametric dependence of vortices centered at y = 0, where the background flow reverses, is investigated. Using DBSA with β = 0, a vortex initially centered at y = 0 relaxes to a state consistent with the constraints. A triangular shape reminiscent of observations of Jupiter is seen to naturally emerge. Lastly in this section, a comparison of the contour dynamics solutions to DBSA applied to the full 1 3 4 -layer model is made. Finally, we conclude and summarize the paper in Section 7.

2.A. Layer model review
The standard 2-layer model of GFD [13] is composed of two coupled advection equations, where the potential vorticity in each layer is given by Here ∇ 2 ψ = ψ xx + ψ yy with subscript denoting partial derivative, [ψ, q ] = ψ x q y − ψ y q x is the Jacobian or ordinary bracket operator, and i = 1 corresponds to the upper layer. This formulation allows for two stretching terms, F i , and (non-standard) two values for β i . For convenience we let F 1 = F and F 2 = δF , whence δ = H 1 /H 2 is the ratio of the respective heights of the two layers.
The standard 1 1 2 -layer model is obtained by letting H 2 → ∞ keeping H 1 fixed, with the additional assumption that ψ 2 ≡ 0. However, in this limit, the bottom layer may remain dynamic with a vorticity q 2 = ∇ 2 ψ 2 + β 2 y that evolves independent of the upper layer. If the lower layer resides in a steady state, say ψ 2 (x), then we arrive at the 1 3 4 -layer model, with the upper layer dynamics in the presence of a steady deep layer flow governed by where without risk of confusion we have dropped the upper layer subscript 1. In Section 4 we will modify the dynamics of this 1 3 4 -layer model of (3) and (4) in such a way as to make it amenable to contour dynamics, which we will use subsequently in our analyses. Beforehand, let us now turn to the Hamiltonian description possessed by the model of (3) and (4).

2.B. Hamiltonian structure of 1 3 4 -layer model
The model is a Hamiltonian field theory with its Hamiltonian functional H naturally being the total energy, kinetic plus potential, written as a functional of the dynamical variable q: with dx = dxdy and the Green's function, G(x − x ), satisfying Because the dynamical variable q does not constitute a set of canonically conjugate field variables, the system takes noncanonical Hamiltonian form (see e.g. [8][9][10] for review) in terms of the following Poisson bracket: where J QG := −[q, · ], the usual Jacobian expression defined above, is the Poisson operator.
The bracket of (7) is a binary, antisymmetric operator, on functionals A, B, expressed in terms of their functional derivatives, δA/δq and δB/δq. Most importantly, the bracket of (7) also satisfies the Jacobi identity, With the Hamiltonian of (5) and the bracket (7), the time evolution of a functional A[q] is given by Using these in (8) and (7) yields the 1 3 4 -layer model of (3) and (4). The antisymmetry of the bracket ensures conservation of energy dH[q]/dt = 0, while invariants P µ associated with other Noether symmetries satisfy {P µ , H} = 0. In particular, if T is a function only of y, the zonal linear momentum, will be conserved: where the last equality follows given either periodic conditions or channel walls in the north and south. In the following we will mostly work in a reference frame translating at speed c (yet to be determined); this is equivalent to generating the motion with the Hamiltonian H c = H − cP rather than H since δH c /δq = −ψ − cy. Similarly, T may possess other symmetries giving rise to the possible class of invariants where φ µ ∈ {x, y, x 2 + y 2 }, corresponding to the momenta arising from two possible translational symmetries and L, the angular momentum arising from rotational symmetry, respectively.
Noncanonical Poisson brackets like (7) (7) and consequently the QG equations have Casimir invariants of the form with C(q) an arbitrary ordinary function; commonly used examples of such conserved functionals are the mean of the PV itself and the potential enstrophy

2.C. Dirac constraints and steady states
In our previous work [3] we used simulated annealing, a technique based on Hamiltonian structure, to obtain a variety of steady states for QG flows. Our work generalized a technique previously introduced in [2] by adding a smoothing metric and, importantly, Dirac constraints that allow for the relaxation to a larger class of steady states. Starting with an initial state, the technique produces a dynamical system that rearranges parcels of fluid conserving each bit's potential vorticity to achieve a maximum or minimum value of the Hamiltonian functional. This is done via a time-stepping of the PV by advecting with an artificial non-divergent velocity obtained from a so-called Dirac bracket, a generalized Poisson bracket that builds in arbitrary invariants, e.g., for two such invariants C 1,2 it has the form where {C j , B} D = 0 for any functional B. In Dirac's original work [14], the bracket {, } of (13) was the usual canonical Poisson bracket; however, his construction gives a bracket that satisfies the Jacobi identity for any bracket {, } which itself satisfies the Jacobi identity (see [15,16]).
Analogous to (7), associated with a bracket of the form of (13) is a Dirac Poisson operator, which we denote by J D . The relaxation dynamics is obtained by using a velocity field for advection obtained by essentially squaring J D . For PV-like dynamics we refer the reader to [3]. In Section 5 we tell this story more explicitly for the contour dynamics that we use in the present work. Before doing so, we first consider the stability of jets in the 1 3 4 -layer model setting, followed by a discussion of localized steady states.

3.A General form
Equilibrium states in some frame of reference satisfy δF/δq = 0, where with P µ defined by (11) being linear or angular momenta depending on the choice of the function φ µ (x), which correspond respectively to jets and vortices, C being a Casimir invariant yet to be chosen, and λ µ is a Lagrange multiplier for a chosen φ µ . We split up q = Q+q , with the associated ψ = Ψ + ψ , where Then, if Q is an equilibrium solution of (3) in a uniformly translating or rotating frame, corresponding to a jet or a vortex, it solves We these assumptions, (14) gives Equation (16) implies Ψ + λ µ φ µ and Q are functionally related -if we make the choice Thus, the quantity ∆F[q ] can serve as a Lyapunov functional for stability, by the energy-Casimir method. (See e.g. [9] for review and early references. See also [17,18] for detailed discussion.) The flow will be stable if q = Q is an extremum; i.e., if is either positive (except, of course, for q = 0) or sufficiently negative to overcome the first term in (17), which is positive.
Suppose C (q) ≥ D 0 ≥ 0; then we have the standard mean value theorem result: with s between 0 and s, and This is often referred to as Arnold's first theorem (A-1).
For the 1 3 4 layer model, we can also find cases with ∆F[q ] negative because of the existence of a bound on the energy and PV terms: In particular, where Z is defined by (12). This follows from where the first term of (18) is clearly positive definite. This means ∆F[q ] of (17), will be Therefore and the flow will also be stable; this is often referred to as Arnold's second theorem (A-2).
Linearized theory would end up in the same place with I just replaced by the Taylor-series

3.B Vortices and Jets with a linear PV-streamfunction relationship
A particularly simple choice, solutions of which we call "linear structures", has and Again, when we split into the background being a jet or a vortex, the term linear in q clearly vanishes, giving The expression of (20) is positive definite if α < 0, whence we infer stability. But the flow will also be stable if ∆F[q ] ≤ 0 with equality occurring only at q = 0. We use and the Fourier decomposition of q , denotedq , to obtain which will be negative if α < F . (Alternatively, the Poincaré inequality would yield this result.) Thus for linear structures we have the following two conditions for stability: For the specific case of zonal jets moving at speed c, we have φ = y and since T = F ψ 2 + βy. If Ψ has a term linear in y plus a periodic function, corresponding to a zonal flow that is periodic, ψ 2 will have the same character. This still leaves a great deal of freedom to choose the structure of the upper layer flows under the presumption that the deep flows are not well known. As noted before, c can be viewed as the rate of translation of the reference frame and will be the rate the vortices are moving when we look for timeindependent eddies atop zonal flows. For simplicity, however, we will take for the upper (setting both the length and velocity scales for nondimensionalization based on the upper layer jets). Then substituting (23) into (22) gives for the lower layer with Figure 1 shows the regions of stability for this solution, α < 0 or 0 < α < F , in terms of from available information from the Galileo probe [7] is that D = 1 which occurs for α = 1.
In that case, the flow will be stable if F > 1 meaning the length scale of the jets is larger than the deformation radius F −1/2 .
Intimately related to the energy-Casimir method is a Rayleigh criterion (see Refs. [17] and [19] for general discussion); taking However, so that we can make the flow stable if Q y does not change sign. Thus sufficiently strong β, i.e., greater than one in the case D = 1, will also stabilize the flow.
Finally, we note that the 1 1 2 layer model just has T = βy and ψ 2 = 0 so that (22) in the sinusoidal case implies α = 1 + F and (A.2) cannot hold. When β < 1, the flow is indeed unstable, and the jets break up.

4.A. Localized vortex in jets
The 1 3 4 -layer model written relative to the stationary jets in a frame moving at speed c takes the form, We seek steady states of this system that satisfy (26) with ∂q /∂t =0 . In particular, we are interested in localized vorticity anomalies, which from (26) must satisfy for an arbitrary function Q. If the vortex is decaying in all directions, the vanishing of q and ψ far from the vortex center implies and therefore from (27) outside closed streamlines. Far from the vortex center, ψ .
The conditions for stability C > 0 or C < −1/F ensure that ∇ 2 ψ is everywhere related to ψ by a positive coefficient. So we don't have far field wavelike behavior and can expect to find isolated vortex solutions.

4.B Linear case
For linear background flow, where C = −1/α and one has decaying modified Bessel function solutions for F > α. The conditions for stability are exactly those that permit isolated vortices to be embedded in the flow. The case of sinusoidal jets just has α = 1.
For comparison, in the 1 1 2 system, we would have which cannot have isolated solutions in all directions.

4.C. Integral conditions
From the x-moment of the equations for the PV anomalies where in proceeding from the third to the fourth equality (28) was used and the q ψ y -term vanishes by integration by parts. If we definẽ (which goes to zero in the far field), our moment equation tells us for steady propagation this leads to The constraint (31) places on the speed has to be considered along with the constraints implied by (22). In particular, for the linear, sinusoidal case, (22) implies α = −1/C (Q) = 1 and c = −β. A small vortex must then be centered where U = c = −β. Note that the zero in the denominator of matches with the zero in the numerator so that Q = −1 everywhere.

4.D Linear jets
Inserting Q = −α(Ψ + cy) into (26) gives where the final form that uses (30) is convenient for seeking steady states, as we shall see in Section 4.D. We will be looking for contour dynamics type solutions withq = q 0 within an area A and zero outside. We then want to solve The integral condition (31) gives while the far-field condition, from eliminating α from (25), has Given the parameters of the background flow, these two expressions for c imply that the north-south location of the vortex is determined, although it depends on the precise shape.
When finding steady states, we shall choose the centriod of the vortex, let the algorithm calculate c which will satisfy (31) and then use that to determine the correct β value. We can then use the β(y c ) values to find where the vortex will reside given instead the value of β.

4.C. Modified dynamics
From Section 4.D the problem of a vortex with a uniform PV anomaly under the linear jet assumption amounts to the following choice for Q(Z): with the characteristic function χ A (x) = 1 when x is in the patch area A and zero when it is outside. The boundary of the patch must be a contour of constant Z and the equation to solve is with Ψ + cy + ψ = const. on ∂A .
We will use contour dynamics to evaluate u given the boundary shape and the Dirac-bracket synthetic annealing of [3] for a modified dynamical system, adapted for contour dynamics in Section 5, to find the shape.
To put (32) into a form where the DBSA tools can be applied, we again writeq = q +αψ giving the dynamical equation However, if we are only interested in the steady states, we can just as well consider the modified dynamics of ∂ ∂tq because (35) and (36) possess the same steady states: if we construct a simulating annealing dynamics that relaxes to steady states of (36), then we obtain steady states of (35). However, it should be borne in mind that unlike (3), which preserves q on particles, (36) preservesq on particles, and our DBSA algorithm will do this as well. For the case of interest here, wherẽ q is a piecewise constant patch, in the modified dynamics the anomaly within the patch and the area of the patch are conserved, whereas that is not the case for the original equation (35). Therefore, the modified dynamics can be treated by the methods of Hamiltonian contour dynamics that we describe next.

5.A. Hamiltonian Structure of Contour Dynamics
The equations of contour dynamics [20,21] are an example of a reduction of the twodimensional Euler fluid equations. Consequently, contour dynamics inherits the Hamiltonian structure of vortex dynamics (see e.g. [8,9]), and can in fact be derived therefrom. The reduction is based on initial conditions where the dynamical variable is constant in a region bounded by a contour. Then for transport equations like that for two-dimensional vorticity it is known that this structure is preserved in time, with the dynamics restricted to be that of the moving bounding contour. The Hamiltonian structure of interest here is one with a noncanonical Poisson bracket, like that of Section 2, that does not require the contour bound a star-shaped region, i.e., the contour need not have a parameterization as a graph of an angle. Indeed the bounding contour is any plane curve (or curves) with an arbitrary parameterization. Here we will first describe the situation for the two-dimensional Euler fluid equations, where the scalar vorticity ω(x, t) = ∇ 2 ϕ is constant inside the contour, before showing how to apply this to the case of interest here.
The reduction to contour dynamics proceeds by replacing ω by a plane curve that bounds a vortex patch or patches X(σ) = (X(σ), Y (σ)). Here the curve parameter σ is not chosen to be arc length because arc length is not conserved by the dynamics of interest.
The noncanonical Poisson bracket for the contours is given, in its most symmetrical form, where we assume closed contours, although generalizations are possible. Observe that if the two functionals A and B are parameterization invariant, then so is their bracket {A, B}.
The bracket of (38) can be rewritten as where the noncanonical Poisson operator J for this bracket is the following skew-symmetric matrix operator: The dynamical equations for the contour are generated by inserting the following compact form for the Hamiltonian into the bracket of (38): whereτ andτ are the unit vectors tangent to the contour, with τ being that for the contour parameterized by σ , and φ(ρ) satisfies ∇ 2 φ(ρ) = G(ρ), where ρ = |x − x | and G is the twodimensional Green's function. Note, in (39) the argument of φ is |X−X |. This Hamiltonian can be obtained from that for the two-dimensional Euler equation by restricting to patchlike initial conditions and manipulating; furthermore, it can be shown to be parametrization invariant in accordance with our theory.
The area of a patch is evidently given as follows: Starting from an expression for the angular momentum L (physically minus the angular momentum) of (11) one can reduce to obtain its contour dynamics form which is a dynamical invariant following from Noether's first theorem, i.e., an invariant not due to bracket degeneracy but tied to the Hamiltonian of interest (39) via {L, H} = 0.
Again observe L is parametrization invariant, as expected on physical grounds.

5.B. Dirac Brackets and Simulated Annealing
Now we use the development of Section 5.A to construct a Dirac bracket analogous to that of (13) of Section 2.C, in order to apply our DBSA method. This is the contour dynamics version of the procedure in [3] with the Dirac bracket constructed using the contour dynamics bracket of (38). This will yield a system of the form where we have set X 1 = X and X 2 = Y and used repeated index notation over j, k = 1, 2, η H and η SA are numbers that weight the contributions of each of the terms of (43) to the dynamics, F[X] is a functional analogous to that of (14) but now for contour dynamics, G is a symmetric smoothing metric that we are free to choose and, most importantly, J CDD is the Poisson operator that arises from (13) upon using (38) with constraints C 1,2 , constraints that we will choose explicitly below. We found in [3] that in order to obtain a rich class of steady states it was necessary to use the Dirac constraints C 1,2 chosen judiciously for the desired state. Ordinary SA corresponds to the case where J CD is used in (43), while DBSA uses J CDD that enforces the Dirac constraints.
Relaxation proceeds under the dynamics of (43) in a manner analogous to the H-theorem relaxation of the Boltzmann equation to thermal equilibrium. The functional F generates relaxation dynamics to δF = 0 because dF/dt ≥ 0, which follows from (43). In Section 5.C, we will demonstrate how this works in the simpler context of plain contour dynamics, before using this simulated annealing technique to calculate Jovian vortices in Section 5.D.

5.C. The Kirchhoff ellipse with shear
As a first example we calculate the Kirchhoff ellipse, the well-known exact solution of the two-dimensional incompressible Euler equation. Because there are many steady states in rotating frames, e.g., the V-states of rigidly rotating vortex patches with m-fold symmetry [23], something more is need to select out the Kirchhoff ellipse. Extremization of δF = δ(H − ΩL) = 0, where L is given by (42), an expression for the angular momentum, and Ω a constant, is insufficient -it need not preserve L and simply goes to a circle. For this reason we employ Dirac constraints in our DBSA algorithm. To this end we choose L, already an invariant, to play a dual role as one of our Dirac constraints. The other Dirac constraint is chosen as in [3] to be the xy-moment that enforces 2-fold symmetry. For contour dynamics this is These are viable constraints because {L, K} = 0, a Dirac bracket requirement that ensures the denominator of (13) does not vanish.
Because the Dirac bracket is complicated we introduce the following shorthand notation: and obtain by direct calculation which give Similarly, we obtain which when inserted into (13) with F give the following vector field analagous to (40): Here the parameters η H and η SA have been set to zero and unity, respectively. added to the Hamiltonian functional (see [24]) and the relaxation to a variety of Moore and Saffman [25] ellipses -the steady version of Kida ellipses [26] -are obtained.

5.D. Application to Jovian vortices
Let us now describe how the formalism can be used to calculate Jovian vortices. For the modified dynamics of (36), vortex states will depend on the strength of the anomaly, q 0 , the initial radius, r 0 with A = πr 2 0 , and the initial center latitude y 0 . We use Ψ = cos(y) and let the procedure determine what the changes in c and U 1 need to be in order to have a steady equilibrium. We will use the linear momenta as Dirac constraints, which to within a sign are C 1 = xq and C 2 = yq. The DBSA routine then makes [cos(y) + ψ + a 1 x + a 2 y,q] = 0 .
For this case the integral conditions (cf. Section 4.C) are obtained by multiplying (51) by y and x and integrating, yielding respectively Consistent with symmetry, a 1 will be zero; comparing the second condition for a vortex patch of area A to the integral constraint of (33) of Section 4.D gives To match with the expressions from the full model, according to (25) this would need to be a 2 = −[β + F (U 1 − U 2 )]/α (or −β in the linear jet case with U 1 = U 2 and α = 1), but that will generally not be the case. The DBSA algorithm will give c(q 0 , A, y 0 ); we can adjust y 0 so that this agrees with the value stipulated by the far-field requirement. We shall simply look at the inverse problem: examining β(q 0 , A, y 0 ).

6.A. CD-DBSA
We first concentrate on the case in which the background flow does not vary with depth, U 1 = U 2 = 0, D = 1 so that α = 1 (cf. (23), (24), and (25)). When β = 0, the vortex center resides at y = 0 where the background flow changes sign, and it is symmetric both east-west and north-south as depicted in Figure 3. Observe in Figure 3(a) that it elongates as |q 0 | decreases. In Figure 3(b) where the size increases, the vortex becomes visibly different from an exact ellipse, since it feels the changes in the shear. It is known that the Moore and Saffman [25] ellipse in uniform shear (the steady version of Kida's solution [26]) has a smaller radius of curvature at the north and south for weaker shears relative to q 0 ; similarly, the solutions of Figure 3(b) have more curvature where the shear is small.
Although the solutions of Figure 3 were obtained using the full DBSA algorithm with Dirac constraints C 1,2 being the x and y moments, they could equally well be found by just applying SA; unlike the Kirchhoff ellipse [27] above which would, in the absence of the Dirac constraints, become circular. Essentially, the background shear locks in the position and orientation, with the amplitudes a i remaining zero throughout because of the symmetry.
However, when we consider β = 0 this is no longer the case and the constraints become critical.
For the β = 0 case, as mentioned in Section 5.D, we solve this problem in reverse by starting with a vortex centered at y = y 0 . Figure 4(a) shows that the standard synthetic annealing process moves this down to the y = 0 axis of symmetry and then reverts to the solution in Figure 3. In contrast, the constrained solution, Figure 4(b), remains centered at y 0 in the sense that the integral dx qy is conserved; this is because that is built into the Dirac bracket as C 2 . When the solution settles, we now have a finite value of a 2 , which depends on the offset y 0 and which then determines the compatible value of β = −αa 2 . Jupiter's Great Red Spot is, of course, in the southern hemisphere, so the sense of circulation for an anticyclone is reversed. So we can take the β = 0.6 solution in Figure 5 and reverse the direction of the shear flow and the vortex circulation; the result has a slightly triangular shape pointing towards the equator with relatively rapid westward flow around the north side. These features can be seen in the Red Spot.

6.B. Comparison with continuous case
For the continuous version (see [28]), with β = 0, DBSA is used directly on the original equations (3) and (4)  Because of the periodicity, the constraints are tapered to make them periodic at the boundaries, e.g., but otherwise the approach follows that in Ref. [3].
The case with the vortex centered is fairly straightforward, and gives solutions looking very similar to the CD solutions. Figure 6 shows an off-centered case with y 0 = 1. We have used a southern hemisphere situation, consequently, the signs of the cos(y) and q 0 terms are reversed. The estimated phase speed is c = −0.575.
When we include the beta-effect, the standard DBSA procedure works less well, perhaps because, at least initially, it tries to generate a net north-south flow which is problematic with the βv term. So we have taken a somewhat different tack: we solve the modified dynamics problem to find a steady state (α = 1) with the propagation speed c = − q 0 sin y/ q 0 . This can be rewritten by adding and subtracting − cos y + cy as  Thus, if we use as an initial condition the ψ 0 found from DBSA applied to the modified dynamics, the resulting structure should propagate at a speed c = −β.
The q 0 field is fairly similar, while the ψ fields, as seen in Figure 7 are virtually indistinguishable. It shows the rapid flow crossing north of the vortex, with some more northerly streamlines turning back and merging into the jet above. These features, along with the asymmetric bulge to the north, are noticeable in the movies of the flow near the Red Spot.

Summary and Conclusions
We have studied the 1 3 4 layer model with deep sinusoidal jets and a vortex in the upper layer. The criterion for the stability of the upper layer jets, which are free to move relative to the deep flow, is that the scale of the jets (in the sense of the inverse of the wavenumber) must be larger than the deformation radius. This is also the necessary condition for isolated vortices to exist within the upper layer. In the absence of deep flow, the jets are not stable and the condition for isolated vortices cannot hold. Thus, the model argues that the longlived vortices will be relatively shallow compared to the jets.
Using a Hamiltonian formulation of contour dynamics and synthetic annealing, we constructed vortex patch solutions in the sinusoidal flow in the absence of β. These are centered on the line where the shear flow has u = 0. Because of the deformation radius and the variations in shear, these are not precisely elliptical in contrast to the Moore and Saffman [25] case.
With β, the shape-preserving vortices will no longer be centered and will propagate. To be able to apply the synthetic-annealing procedure, we introduced a modified form of the dynamics which has the same steady solutions but evolves according to (36). It should be noted that this procedure no longer preserves the potential vorticity on each particle as they are rearranged. But the final solution is a valid steady state and shows that the vortex is now off-center and asymmetrical. It has a triangular component reminiscent of the Red Spot.
Physically, there are, of course, other processes acting in the Jovian atmosphere. We have used a quasi-geostrophic model, which can become inaccurate in regions with high Rossby number or, perhaps more appropriate for large baroclinic vortices, with changes in thickness between isopycnals which are not small compared to the mean thickness. The deep layer is not necessarily completely steady. This can lead to radiation of waves that drain energy from the spot; however, estimates of the rate from full two-layer calculations with the two-beta model [12] suggest it is very slow. Mergers with small spots may spin the large spots back up so that they can be maintained. The processes that maintain the deep jets