text:
# Emergence of intense jets and Jupiter Great Red Spot as maximum entropy structures.

## 1 Introduction

Atmospheric and oceanic flows are often organized into narrow jets. They can zonally flow around the planet like the jet streams in the earth stratosphere, or the eastward jet at 24 in the northern hemisphere of Jupiter ( Maxworthy 1984). Jets can alternatively organize into rings, forming vortices, like the rings shed by the meandering of the Gulf-Stream in the western Atlantic Ocean. The flow field in Jupiter most famous feature, the Great Red Spot, is an oval-shaped jet, rotating in the anticyclonic direction and surrounding an interior area with a weak mean flow ( Dowling and Ingersoll 1989), see figure 1(a). Robust cyclonic vortices have a similar jet structure ( Hatzes et al 1981 ), see figure 1(b).

Such jets and vortices are in a turbulent surrounding, and the persistence of their strength and concentration in the presence of eddy mixing is intriguing. The explanation proposed in this paper is based on a statistical mechanical approach: the narrow jet or vortex appears as the most probable state of the flow after a turbulent mixing of potential vorticity, taking into account constraints due to the quantities conserved by the dynamics, especially energy. Such a statistical theory has been first proposed for the two-dimensional Euler equations by Kuz’min (1982), Robert (1990) , Robert and Sommeria (1991), Miller (1990). See Brandt et al (1999) for a recent review and discussion. This theory predicts an organization of two-dimensional turbulence into a steady flow (with fine scale, ’microscopic’ vorticity fluctuations). Complete vorticity mixing is prevented by the conservation of the energy, which can be expressed as a constraint in the accessible vorticity fields. A similar, but quantitatively different, organization had been previously obtained with statistical mechanics of singular point vortices with the mean field approximation, instead of continuous vorticity fields (Onsager 1949, Joyce and Montgomery 1973).

Extension to the quasi-geostrophic (QG) model has been discussed by Sommeria et al (1991), Michel & Robert (1994), Kazantsev Sommeria and Verron (1998). This model describes a shallow water system with a weak vorticity in comparison with the planetary vorticity (small Rossby number), such that the pressure is in geostrophic balance, and the corresponding free surface deformation is supposed small in comparison with the layer thickness. For Jupiter the free surface would be rather at the bottom of the active atmospheric layer, floating on a denser fluid, as discussed by Dowling and Ingersoll (1989), see Dowling (1995) for a review. The gradient of planetary vorticity is accounted by a beta-effect. An additional beta-effect, depending on the latitude coordinate $`y`$, is introduced to represent the influence on the active atmospheric layer of a steady zonal flow in the deep interior, as discussed by Dowling and Ingersoll (1989).

The free surface deformability, representing the strength of the density stratification, is controlled by the Rossby radius of deformation $`R^{}`$. The two-dimensional Euler equation is recovered in the limit of very strong stratification for which $`R^{}\mathrm{}`$. We shall consider in this paper the opposite limit of weak stratification for which $`R^{}`$ is much smaller than the scale of the system $`L`$. This limit is appropriate for large scale oceanic currents, as the radius of deformation is typically 10-100 km. For Jupiter, $`R^{}`$ is estimated to be in the range 500-2500 km, while the Great Red Spot extends over 20,000 km in longitude, and 10,000 km in latitude, so the limit $`R^{}/L0`$ seems relevant. We show that in this limit the statistical equilibrium is made of quiescent zones with well mixed uniform potential vorticity, bounded by jets with thickness of order $`R^{}`$. This provides therefore a general justification of jet persistence. Some of the ideas used have been already sketched in Sommeria et al (1991), but we here provide a systematic derivation and thorough analysis.

In principle, the Quasi Geostrophic approximation breaks down for scales much larger than the radius of deformation, so that the limit $`R^{}/L0`$ seems inconsistent with the QG approximation. However the relevant scale is the jet width, which remains of order $`R^{}`$, so that the Quasi Geostrophic approximation remains valid in this limit. This point has been discussed by Marcus (1993) for the Great Red Spot, which he supposes to be a uniform Potential Vorticity ( PV ) spot surrounded by a uniform Potential Vorticity background ( we here justify this structure as the result of Potential Vorticity mixing with constraints on the conserved quantities ). Analyzing wind data in the Great Red Spot, Dowling & Ingersoll (1989) concluded that the QG approximation is good up within typically 30% error, which is reasonable to a first approximation. Statistical mechanics of the more general shallow water system (to be published), predicts a similar jet structure. The present Quasi Geostrophic results therefore provide a good description as a first approximation.

We first consider the case without beta-effect in section 2. We furthermore assume periodic boundary conditions (along both coordinates) in this section to avoid consideration of boundary effects. Starting from some initial condition with patches of uniform PV, we find that these patches mix with uniform density (probability) in two sub-domains, with strong density gradient at the interface, corresponding to a free jet. The coexistence of the two sub-domains can be interpreted as an equilibrium between two thermodynamic phases. We find that the interface has a free energy per unit of length, and its minimization leads to a minimum length at equilibrium. This results in a constant radius of curvature, in analogy with surface tension effects in thermodynamics, leading to spherical bubbles or droplets. The range of the vortex interaction is of the order $`R^{}`$, therefore becoming very small in the limit of small radius of deformation, so the statistical equilibrium indeed behaves like in usual thermodynamics with short range molecular interactions. This contrasts with the case of Euler equation, with long range vortex interactions, analogous to gravitational effects (Chavanis Sommeria and Robert 1996, Chavanis 1998).

Figure 7 summarizes the calculated equilibrium states, depending on the total energy and a parameter $`B`$ representing an asymmetry between the initial Potential Vorticity patch areas, before the mixing process. We obtain straight jets for a weak asymmetry and circular jets for higher asymmetry. Such circular a jet reduces to an axisymmetric vortex, with radius of order $`R^{}`$, in the limit of low energy.

We discuss the influence of the beta-effect or the deep zonal flow in section 3. The channel geometry, representing a zonal band periodic in the longitude $`x`$ is appropriate for that study. With the usual beta-effect $`\beta y`$, linear in the transverse coordinate $`y`$, statistical equilibrium is, depending on the initial parameters, a zonal flow, or a meandering eastward jet, or a uniform velocity $`v_m=R^2\beta `$ whose induced free surface slope cancels the beta-effect ( uniformization of Potential Vorticity ) on which circular vortices can coexist.

For more general beta-effects, due to the deep zonal flow, we find that the jet curvature depends on latitude $`y`$. In particular a quadratic beta effect $`ay^2`$ leads to oval-shape jets, similar to the Great Red Spot. Using the determination of the sublayer flow from Voyager data by Dowling and Ingersoll (1989), we show in section 4, that such a quadratic effective beta-effect is indeed a realistic model for Jupiter atmosphere in the latitude range of the Great Red Spot and the White Ovals, the other major coherent vortices on Jupiter. Using these data on beta-effect, as well as the shear in the zonal flow at the latitude of the Great Red Spot, the jet width and its maximum velocity, we deduce all the parameters of our model.

## 2 The case with periodic boundary conditions

### 2.1 The dynamical system

We start from the barotropic Quasi Geostrophic (QG) equation :

$$\frac{q}{t}+𝐯q=0$$
(1)
$$q=\mathrm{\Delta }\psi +\frac{\psi }{R^2}h(y)$$
(2)
$$𝐯=\widehat{𝐳}\psi $$
(3)
where $`q`$ is the potential vorticity ( PV ), advected by the non-divergent velocity $`𝐯`$, $`\psi `$ is the stream function, $`R`$ is the internal Rossby deformation radius between the layer of fluid under consideration and a deep layer unaffected by the dynamics. $`x`$ and $`y`$ are respectively the zonal and meridional coordinates ( $`x`$ is directed eastward and $`y`$ poleward ). The term $`h(y)`$ represents the combined effect of the planetary vorticity gradient and of a given stationary zonal flow in the deep layer, with stream function $`\psi _d(y)`$: $`h(y)=\beta y+\psi _d/R^2`$. This deep flow induces a constant deformation of the free surface, acting like a topography on the active layer. We shall therefore call $`h(y)`$ the ’topography’, and study its influence in section 3. Let us assume $`h(y)=0`$ in this section. We define the QG equations (1,2) in the non-dimensional square $`D=[\frac{1}{2},\frac{1}{2}]^2`$. $`R`$ is then the ratio of the internal Rossby deformation radius $`R^{}`$ to the physical scale of the domain $`L`$.

Let $`f_Dfd^2𝐫`$ be the average of $`f`$ on $`D`$ for any function $`f`$. Physically, as the stream function $`\psi `$ is related to the geostrophic pressure, $`\psi `$ is proportional to the mean height at the interface between the fluid layer and the bottom layer, and due to the mass conservation it must be constant (Pedlosky 1987). We make the choice

$$\psi =0$$
(4)
without loss of generality.

The total circulation is $`q=\mathrm{\Delta }\psi +\psi /R^2=\psi /R^2`$ due to the periodic boundary conditions. Therefore

$$q=0$$
(5)
We note that the Dirichlet problem (2) on $`D`$ with periodic boundary conditions has a unique solution $`\psi `$ for a given PV field.

Due to the periodic conditions for $`\psi `$, the linear momentum is also equal to 0,

$$𝐯=0$$
(6)
The energy

$$E=\frac{1}{2}_Dq\psi d^2𝐫=\frac{1}{2}_D[(\psi )^2+\frac{\psi ^2}{R^2}]d^2𝐫$$
(7)
is conserved ( we note that the first term in the right hand side of (7) is the kinetic energy whereas the second one is the gravitational available potential energy ).

The integrals

$$C_f(q)=_Df(q)d^2𝐫$$
(8)
for any continuous function $`f`$ are also conserved, in particular the different moments of the PV. In the case of an initial condition made of a finite number of PV levels, the areas initially occupied by each of these levels is conserved, and this is equivalent to the conservation of all the constants of motion (8).

### 2.2 The statistical mechanics on a two PV levels configuration.

#### 2.2.1 The macroscopic description.

The QG equations (1) (2) are known to develop very complex vorticity filaments. Because of the rapidly increasing amount of information it would require, as time goes on, a deterministic description of the flow for long time is both unrealistic and meaningless. The statistical theory adopts a probabilistic description for the vorticity field. The statistical equilibrium depends on the energy and of the global probability distribution of PV levels. Various previous studies (Sommeria & al 1991), (Kazantsev & al 1998) indicate that a model with only two PV levels provides a good approximation in many cases. The determination of the statistical equilibrium is then simplified as it depends only on the energy, on the two PV levels, denoted $`q=a_1`$ and $`q=a_1`$ and on their respective areas $`A`$ and $`(1A)`$ in $`D`$. The number of free parameters can be further reduced by appropriate scaling. Indeed a change in the time unit permits to define the PV levels up to a multiplicative constant, and we choose for the sake of simplicity :

$$\frac{a_1a_1}{2}=1$$
(9)
and define the non-dimensional parameter $`B`$ as :

$$B\frac{a_1+a_1}{2}$$
(10)
The condition (5) of zero mean PV imposes that $`a_1A+a_1(1A)=0`$. This means that $`a_1`$ and $`a_1`$ must be of opposite sign and, using (9) and (10), $`A=(1B)/2`$. The distribution of PV levels is therefore fully characterized by the single asymmetry parameter $`B`$, which takes values between -1 and +1. The symmetric case of two PV patches with equal area $`A=1/2`$ corresponds to $`B=0`$, while the case of a patch with small area (but high PV, such that $`q=0`$) corresponds to $`B1`$. Note that we can restrict the discussion to $`B1`$ as the QG system is symmetric by a change of sign of the PV.

The two PV levels mix due to turbulent effects, and the resulting state is locally described by the local probability (local area proportion) $`p(𝐫)`$ to find the first level at the location $`𝐫`$. The probability to find the complementary PV level $`a_1`$ is $`1p`$, and the locally averaged PV at each point is then

$$\overline{q}(𝐫)=a_1p(𝐫)+a_1(1p(𝐫))=2\left(p\frac{1}{2}\right)+B$$
(11)
where the second relation is obtained by using (9) and (10).

As we consider the evolution of two PV patches, the conservation of all invariants (8) is equivalent to the conservation of the area $`A`$ of the patch with PV value $`a_1`$ ( the area of the other PV level $`a_1`$ being $`1A`$ ). The integral of $`p`$ over the domain must be therefore equal to the initial area $`A`$ ( the patch with PV level $`a_1`$ is mixed but globally conserved ),

$$A\frac{1B}{2}=_Dp(𝐫)d^2𝐫$$
(12)
As the effect of local PV fluctuations is filtered out by integration, the stream function and the velocity field are fully determined by the locally averaged PV $`\overline{q}`$ as the solution of

$$\overline{q}=\mathrm{\Delta }\psi +\frac{\psi }{R^2};\psi periodic$$
(13)
$$and𝐯=\widehat{𝐳}\psi $$

Therefore the energy is also expressed in terms of the field $`\overline{q}`$ :

$$E=\frac{1}{2}_D\left[(\psi )^2+\frac{\psi ^2}{R^2}\right]d^2𝐫=\frac{1}{2}_D\overline{q}\psi d^2𝐫$$
(14)
Here the energy of the ’microscopic’ PV fluctuations has been neglected (replacing $`q`$ by $`\overline{q}`$), as justified in the case of Euler equation by Robert and Sommeria (1991). Indeed, considering a ’cutoff’ for the microscopic fluctuations much smaller than $`R`$, the small scale dynamics coincides with the Euler case.

The central result of the statistical mechanics of the QG equations (1,2) is that, under an ergodic hypothesis, we expect the long time dynamics to converge towards the Gibbs states defined by maximizing the mixing entropy

$$S=_D[p(𝐫)\mathrm{ln}p(𝐫)+(1p(𝐫))\mathrm{ln}(1p(𝐫))]d^2𝐫$$
(15)
under the constraints of the global PV distribution (12) and energy (14). It can be shown that the microscopic states satisfying the constraints given by the conservation laws are overwhelmly concentrated near the Gibbs state, which is therefore likely to be reached after a complex flow evolution. A good justification of this statement is obtained by the construction of converging sequences of approximations of the QG equation (1,2), in finite dimensional vector spaces, for which a Liouville theorem holds. This is a straightforward translation of the work of Robert (1999) for 2D Euler equations. The sequence of such Liouville measures has then the desired concentration properties as (1,2) enters in the context considered in Michel & Robert (1994).

#### 2.2.2 The Gibbs states

Following Robert & Sommeria (1991), we seek maxima of the entropy (15) under the constraints (12) and (14). To account for these constraints, we introduce two corresponding Lagrange multipliers, which we denote $`2\alpha `$ and $`C/R^2`$ for convenience in future calculations. Then the first variation of the functionals satisfies :

$$\delta S2\alpha \delta A+\frac{C}{R^2}\delta E=0$$

for all variations $`\delta p`$ of the probability field $`p`$. After straightforward differentiation we obtain:

$$\delta S=_D[\mathrm{ln}p\mathrm{ln}(1p)]\delta pd^2𝐫,\delta A=_D\delta pd^2𝐫$$

$$\delta E=_D\psi \delta \overline{q}d^2𝐫=_D2\psi \delta pd^2𝐫$$
(16)
where the expression of $`\delta E`$ has been obtained by integrating by part and expressing $`\overline{q}`$ by (11). Then we can write the first variation under the form $`_D[\mathrm{ln}p+\mathrm{ln}(1p)2\alpha +2C\psi /R^2]\delta pd^2𝐫`$ which must vanish for any small variation $`\delta p`$. This implies that the integrand must vanish, and yields the equation for the optimum state:

$$p=\frac{1\mathrm{tanh}(\alpha \frac{C\psi }{R^2})}{2},$$
(17)
and using (11) and (13), the partial differential equation

$$q=\mathrm{\Delta }\psi +\frac{\psi }{R^2}=B\mathrm{tanh}\left(\alpha \frac{C\psi }{R^2}\right)$$
(18)
determining the Gibbs states (statistical equilibrium). From now on we forget the $`q`$ over-line for the locally averaged PV and refer to it as the PV.

Therefore, we have shown that for any solution of the variational problem, two constants $`\alpha `$ and $`C`$ exist such that $`\psi `$ satisfies (18). Conversely it can be proved that for any such two constants, a solution to equation (18), in general not unique, always exists. Then $`p`$ associated with one of these solutions by (17) is a critical point of the ’free energy’ $`S(p)+2\alpha A(p)\frac{C}{R^2}E(p)`$ (i.e. its first variation vanishes). Then the Lagrange multipliers are not given but have to be calculated by prescribing the constraints (12) and (14) corresponding to the two parameters $`B`$ and $`E`$ respectively, given by the initial condition. Furthermore, among the states of given energy $`E`$ and asymmetry parameter $`B`$, we shall select the actual maxima.

Finally, let us find a lower bound for the parameter $`C`$ of the Gibbs states with non-zero energy (i.e. $`\psi `$ is not constant on $`D`$). Multiplying (18) by $`\mathrm{\Delta }\psi `$, integrating by part and defining $`f(C\psi )B\mathrm{tanh}(\alpha \frac{C\psi }{R^2})`$, we obtain :

$$C=\frac{_D\left((\mathrm{\Delta }\psi )^2+\frac{1}{R^2}(\psi )^2\right)d^2𝐫}{_Df^{}(C\psi )(\psi )^2d^2𝐫}$$

From which, using $`0<f^{}(C\psi )\frac{1}{R^2}`$ it follows that when $`\psi `$ is not constant :

$$C>1$$
(19)
### 2.3 The limit of small Rossby deformation radius

As suggested by oceanographic or Jovian parameters, we seek solutions for the Gibbs states equation in the limit of a small ratio between the Rossby deformation radius and the length scale of the domain : $`R<<1`$ with our non-dimensional coordinates 111 Modica (1987) considered the minimization of the functional $`E_ϵ(u)=_\mathrm{\Omega }[ϵ\left(u\right)^2+W_0(u)]𝑑𝐱`$ with the constraint $`_\mathrm{\Omega }u(𝐱)𝑑𝐱=m`$ in the limit $`ϵ0^+`$ where $`W_0`$ is a real function with two relative minima. He proved, in a mathematical framework, working with bounded variations functions, that if $`(u_ϵ)`$ are solutions of this variational problem, for any subsequence of $`(u_ϵ)`$ converging in $`L^1(\mathrm{\Omega })`$ as $`ϵ0`$, this subsequence converge to a function $`u_0`$ which takes only the values where $`W_0`$ reaches its minima ; with the interface between the corresponding subdomains having minimal area ( See Modica (1987) for a precise statement ). We note that the Euler equation of this variational problem may be the same as the Gibbs States equation (18) for a convenient choice of $`W_0`$. However as the variational problem itself is different this beautiful result cannot be used in our context.

#### 2.3.1 The uniform subdomains

Then we expect that the Laplacian in the Gibbs states equation (18) can be neglected with respect to $`\psi /R^2`$, except possibly in transition regions of small area. This transforms (18) into the algebraic equation :

$$q=\frac{\psi }{R^2}=B\mathrm{tanh}\left(\alpha \frac{C\psi }{R^2}\right)$$
(20)
Depending on the parameters, this equation has either one, two or three solutions, denoted $`\psi _1,\psi _0`$ and $`\psi _1`$ in increasing order (see figure 2 ). The case with a single solution would correspond to a uniform $`\psi `$, which should be equal to 0 due to the condition $`\psi =0`$. This is only possible for $`E=0`$. Otherwise, we have therefore two or three solutions, with different solutions occurring in subdomains. This condition of multiple solutions requires that the maximum slope for the right hand side of (20) must be greater than $`1/R^2`$ ; this is always realized due to the inequality (19). Furthermore $`\alpha `$ must be in an interval centered in $`CB`$ ( $`\alpha =CB`$ in the symmetric case of figure 2 ).

At the interface between two constant stream function subdomains, a strong gradient of $`\psi `$ necessarily occurs, corresponding to a jet along the interface. These jets give first order contributions to the entropy and energy, but let us first describe the zero order problem. Suppose that $`\psi `$ takes the value $`\psi _1`$ ( resp $`\psi _1`$ ) in subdomains of total area $`A_1`$ ( resp $`A_1`$ ). The reason why we do not select the value $`\psi _0`$ will soon become clear. Using (11) we conclude that the probability $`p`$ takes two constants values $`p_{\pm 1}`$ in their respective subdomains. The two areas $`A_{\pm 1}`$ ( measured from the middle of the jet ) are complementary such that :

$$A_1+A_1=1$$
(21)
Furthermore the constraint (4) of zero domain average for $`\psi `$ implies at zero order,

$$\psi _1A_1+\psi _1A_1=0,$$
(22)
or equivalently, using $`q_{\pm 1}=\psi _{\pm 1}/R^2`$, (11) and (21) :

$$2A_1\left(p_1\frac{1}{2}\right)+2(1A_1)\left(p_1\frac{1}{2}\right)=B$$
(23)
This can be obtained as well from the constraint on the (microscopic) PV patch area (12). The energy inside the subdomains reduces to the potential term $`\psi ^2/2R^2`$, since velocity vanishes. This area energy $`E_A`$ can be computed in terms of $`p_{\pm 1}`$ using $`q_{\pm 1}=\psi _{\pm 1}/R^2`$ and (11) :

$$\begin{array}{cc}\hfill E_A& \frac{\psi _1^2A_1+\psi _1^2A_1}{2R^2}=A_1e(p_1)+(1A_1)e(p_1)\hfill \\ \hfill \text{with}e(p)& R^2\left(2\left(p\frac{1}{2}\right)^2+2B\left(p\frac{1}{2}\right)+\frac{B^2}{2}\right)\hfill \end{array}$$
(24)
There is also an energy in the jet at the interface of subdomains, but it is small with respect to $`E_A`$. Indeed the velocity in the jet, of width $`R`$, is of order $`(\psi _1\psi _1)/R\psi /R`$, and the corresponding integrated kinetic energy is of order $`\psi ^2/R`$. This is small in comparison with the area energy $`E_A`$ (mostly potential) of order $`\psi ^2/R^2`$. A precise calculation will confirm this estimate in next sub-section.

We need to determine three unknown, the area $`A_1`$ and the probabilities $`p_{\pm 1}`$ of the PV level $`a_1`$ in each subdomain, while the constraints (23) and (24) provide two relations. An additional relation will be given by entropy maximization. As we neglect the jet area, the entropy reduces at order zero to the area entropy :

$$\begin{array}{cc}\hfill S_A& A_1s(p_1)+(1A_1)s(p_1)\hfill \\ \hfill \text{with}s(p)& p\mathrm{log}p(1p)\mathrm{log}(1p)\hfill \end{array}$$
(25)
Thus the zero order problem corresponds to the maximization of the area entropy (25) with respect to the 3 parameters $`p_{\pm 1}`$ and $`A_1`$, under the 2 constraints (23) and (24). A necessary condition for a solution of this variational problem is the existence of two Lagrange parameters $`\alpha _0`$ and $`C_0`$, associated respectively with the circulation constraint (23) and with the energy constraint (24), such that the first variations of the total free energy,

$$F_AS_A\frac{C_0}{R^2}E_A+\alpha _0\frac{\psi }{R^2},$$
(26)
vanish. Let us calculate $`F_A`$ using (23) and (24):

$$\begin{array}{cc}\hfill F_A& =A_1f(p_1)+(1A_1)f(p_1),\hfill \\ \hfill \text{with}f(p)& s(p)2C_0\left(p\frac{1}{2}\right)^22\left(C_0B\alpha _0\right)\left(p\frac{1}{2}\right)\frac{C_0B^2}{2}\alpha _0B.\hfill \end{array}$$
(27)
The vanishing of the variations with respect to $`p_1`$ and $`p_1`$ gives that $`f(p_{\pm 1})`$ are local minima of the free energy $`f(p)`$. It is easily proven that the function $`f`$ has two local minima and one local maximum ( for $`C_0>1`$ and $`(C_0B\alpha _0)`$ small enough ) ( see figure 3). The local maximum is achieved for $`p_0`$ corresponding to the value $`\psi _0`$. It is the reason why it has not been taken into account in this analysis. In addition, the vanishing of the first variations with respect to the area $`A_1`$ imposes the free energies $`f(p_{\pm 1})`$ in the two subdomains to be equal. This is like the condition of thermodynamic equilibrium for a chemical species shared by two coexisting phases.

In the expression (27) of $`f(p)`$, the entropy term $`s(p)`$ is symmetric with respect to $`p=\frac{1}{2}`$, as well as the quadratic term. Therefore if the linear term in $`\left(p\frac{1}{2}\right)`$ vanishes the two maxima are equal, with $`p_{\pm 1}`$ symmetric with respect to $`\frac{1}{2}`$. The addition of a linear term obviously breaks this condition of two equal maxima, so the coefficient of the linear term must vanish, thus :

$$\alpha _0=C_0B.$$
(28)
Since $`p_{\pm 1}`$ are symmetric with respect to $`\frac{1}{2}`$, we introduce the parameter $`u`$ by :

$$p_{\pm 1}=\frac{1}{2}(1\pm u).$$
(29)
Using (11),(23) we deduce:

$$\psi _{\pm 1}=R^2(B\pm u)$$
(30)
From (22) we state that the two constant stream function (30) have to be of opposite sign, so that $`u>|B|`$. Introducing (29) in the circulation constraint (23), and using (21), we get :

$$A_{\pm 1}=\frac{1}{2}\left(1\frac{B}{u}\right).$$
(31)
Using these results, the energy (24) becomes

$$EE_A=\frac{R^2}{2}(u^2B^2)$$
(32)
This relates the parameter $`u`$ to the given energy $`E`$ and asymmetry parameter $`B`$. Finally the condition that $`f(p_{\pm 1})`$ are maxima of $`f`$ leads to :

$$u=\mathrm{tanh}(C_0u),$$
(33)
which determines the ’temperature’ parameter $`C_0`$, as represented in figure 4. Therefore all the quantities are determined from the asymmetry parameter $`B`$ and from the parameter $`u`$, related to the energy by (32).

In the limit of low energy, $`u|B|`$, when for instance $`B>0`$, then $`A_1`$ goes to zero, so that $`\psi _1`$ tends to occupy the whole domain. This state is the most mixed one compatible with the constraint of a given value of $`B`$ (or equivalently a given initial patch area $`A=(1B)/2`$). In the opposite limit $`u1`$, we see from (30) that in the two subdomains $`q=\psi /R^2`$ tends to the two initial PV levels $`a_1=1+B`$ and $`a_1=1+B`$. Thus, this state is an unmixed state. It achieves the maximum possible energy $`E=\frac{R^2}{2}(1B^2)`$ under the constraint of a given patch area. We conclude that the parameter $`u`$, or the related ’temperature’ $`C_0`$, linked with the difference between the energy and the maximum accessible energy for the two given initial PV levels, characterizes the mixing of these two PV levels. We shall call $`u`$ the segregation parameter, as it quantifies the segregation of the PV level $`a_1`$ ( or its complementary $`a_1`$ ) between the two phases.

Let us now study the interface between the subdomains.

#### 2.3.2 Interior Jets

At the interface between two constant stream function subdomains, a strong gradient of $`\psi `$ necessarily occurs, corresponding to a jet along the interface. To study these jets, we come back to the Gibbs state equation (18). We expect the Lagrange parameters $`\alpha `$ and $`C`$ to be close to the zero order parameters $`\alpha _0`$ and $`C_0`$, computed in the previous sub-section, so we use $`\alpha =\alpha _0`$ and $`C=C_0`$ to calculate the jet structure. In such a jet, we cannot neglect the Laplacian term in (18), but a boundary layer type approximation can be used: we neglect the tangential derivative with respect to the derivative along the coordinate normal to the interface, $`\zeta `$. Accordingly, we neglect the inverse of the curvature radius of the jet with respect to $`1/R`$.

Thus, from the Gibbs states equation (18), using (28), we deduce the jet equation :

$$\frac{d^2\psi }{d\zeta ^2}+\frac{\psi }{R^2}=B\mathrm{tanh}\left(C_0\left(B\frac{\psi }{R^2}\right)\right)$$
(34)
As the stream function depends only on the normal coordinate $`\zeta `$, the velocity is tangent to the interface, forming a jet with a typical width scaling like $`R`$. We thus make the change of variables defined by :

$$\tau \frac{\zeta }{R};\varphi B+\frac{\psi }{R^2},$$
(35)
leading to the rescaled jet equation:

$$\frac{d^2\varphi }{d\tau ^2}=\mathrm{tanh}(C_0\varphi )+\varphi $$
(36)
The jet equation (36) is similar to a one dimensional equation of motion for a particle ( with the position $`\varphi `$ depending on a time $`\tau `$) under the action of a force $`dU/d\varphi `$ deriving from the potential,

$$U(\varphi )=\frac{\mathrm{ln}(\mathrm{cosh}(C_0\varphi ))}{C_0}\frac{\varphi ^2}{2},$$
(37)
represented in figure 2(b). In its trajectory the particle energy is conserved :

$$\frac{1}{2}\left(\frac{d\varphi }{d\tau }\right)^2+U(\varphi )=Cst$$
(38)
Let $`\varphi _i\psi _i/R^2B`$, $`i=1,0,1`$, corresponding to the solutions $`\psi _i`$ of the algebraic equation (20). From (30), we have $`\varphi _{\pm 1}=\pm u`$. Note that the stationary limit of (36), which must be reached for $`lim_{\tau \pm \mathrm{}}`$, yields again (33). Moreover, the particle energy conservation (38) imposes the integrability condition,

$$U(\varphi _1)=U(\varphi _1),$$
(39)
which is indeed satisfied due to the symmetry of the potential $`U`$. We note that the Lagrange parameter determination (28) and the symmetry of the probabilities (29) with respect to $`\frac{1}{2}`$ could have been deduced from this integrability condition (39) instead of minimizing the free energy (27) (we shall proceed in this way in section 3 to take into account the beta-effect).

The jet equation (36) has been numerically integrated. Figure 5 shows a typical stream function and velocity profile in the jets. Figure 6(a) shows how the jet width depends on the segregation parameter $`u`$. We note that the width of the jet is an increasing function of the mixing and therefore a decreasing function of the energy. Figure 6(b) shows the dependence in $`u`$ of the total non-dimensional energy $`e(u)=\frac{1}{2}_{\mathrm{}}^+\mathrm{}\left(d\varphi /d\tau \right)^2𝑑\tau `$ and of the maximum non dimensional jet velocity $`\left(d\varphi /d\tau \right)_{max}`$.

As the jet structure (36) does not depend on the coordinate tangent to the jet, we can define the jet entropy ( respectively energy, free energy ) per unit length $`S_{Jet}`$ ( respectively $`E_{Jet}`$, $`F_{Jet}`$ ). Multiplied by the jet length, these quantities are the first order corrections to the entropy ( respectively energy, free energy ). Using the change of variables (35), we calculate the jet entropy per unit length :

$$S_{Jet}=R_{\mathrm{}}^+\mathrm{}[s(p(\tau ))s(p_{\pm 1})]𝑑\tau $$

where $`s`$ is defined in (25), and $`p_{\pm 1}`$ are defined in (29). Using the probability equation (17) and (35) we obtain :

$$S_{Jet}=R_{\mathrm{}}^+\mathrm{}[\stackrel{~}{s}(\varphi ))\stackrel{~}{s}(\varphi _{\pm 1}))]d\tau $$
(40)
involving the function $`\stackrel{~}{s}(\varphi )\mathrm{ln}(\mathrm{cosh}(C_0\varphi ))C_0\varphi \mathrm{tanh}(C_0\varphi )`$. Similarly we straightforwardly calculate the potential and kinetic energy per unit length for the jet :

$$E_{Jet}^P=\frac{R^3}{2}_{\mathrm{}}^+\mathrm{}(\varphi ^2\varphi _1^2)𝑑\tau E_{Jet}^K=\frac{R^3}{2}_{\mathrm{}}^+\mathrm{}\left(\frac{d\varphi }{d\tau }\right)^2𝑑\tau $$

We use the integral (38) to calculate $`d\varphi /d\tau `$ and conclude :

$$E_{Jet}=R^3_{\mathrm{}}^+\mathrm{}[\stackrel{~}{g}(\varphi )\stackrel{~}{g}(\varphi _{\pm 1})]𝑑\tau ,$$
(41)
with $`\stackrel{~}{g}(\varphi )=\frac{\mathrm{ln}(\mathrm{cosh}(C_0\varphi )}{C_0}+\varphi ^2`$. Due to the symmetry of the jets, the jets provide no perturbation to the zero-order circulation, so there is no circulation term in the jet free energy expression : $`F_{Jet}=S_{Jet}C_0/R^2E_{Jet}`$. Then

$$F_{Jet}=C_0R_{\mathrm{}}^+\mathrm{}[\stackrel{~}{h}(\varphi )\stackrel{~}{h}(\varphi _{\pm 1})]𝑑\tau $$
(42)
with $`\stackrel{~}{h}(\varphi )=\varphi (\varphi \mathrm{tanh}(C_0\varphi ))`$.

Let us study the sign of $`F_{Jet}`$. As $`\varphi _1`$ verifies $`\varphi _1=\mathrm{tanh}(C_0\varphi _1)`$ and as $`\varphi (\tau )`$ is an increasing function of $`\tau `$ with $`lim_{\tau +\mathrm{}}\varphi (\tau )=\varphi _1`$ we conclude that $`\stackrel{~}{h}(\varphi )\stackrel{~}{h}(\varphi _1)>0`$ for any $`\tau >0`$. Thus $`F_{Jet}>0`$. Using the analogy with usual thermodynamics, the ‘surface tension’ is positive. This favors large ‘bubbles’ which minimize the interfacial length and therefore the corresponding free energy (27). Our initial hypothesis of well separated domains with uniform $`\psi `$ is thus supported, as discussed more precisely in next subsection.

#### 2.3.3 Selection of the sub-domain shape

The above analysis has permitted us to determine the areas of subdomains on which the stream function $`\psi `$ takes the constant values $`\psi _{\pm 1}`$, but the subdomains shape is still to be selected. There is an analogy with two phases coexisting in thermodynamic equilibrium, for instance a gas bubble in a liquid medium, for which a classical thermodynamic argument explains the spherical shape of the bubble by minimizing its free energy, proportional to the bubble area. Our system is isolated rather than in a thermal bath, but the jet energy is small (of order $`R`$) with respect to the total energy. Therefore the subdomain interior behaves like a thermal bath with respect to the jet, so the usual argument on free energy minimization applies. We shall now show this more precisely by directly maximizing the total entropy with constraints, taking into account the jet contribution.

The jet with length $`L`$ has an entropy $`S_{jet}L`$ and energy $`E_{jet}L`$. Since the total energy $`E=E_A(C)+LE_{Jet}`$ is given, the jet has also an indirect influence in the area energy $`E_A`$. This small energy change $`\delta E_A`$ results in a corresponding change in the area entropy $`\delta S_A=(C_0/R^2)\delta E_A`$, from the condition (2.2.2) of zero first variations. Note that there is no area change $`\delta A`$ since the jet is symmetric and has therefore no influence in the condition (12) of a given integral of $`p`$ (the difference in $`p`$ with the case of two uniform patches with boundary at the jet center has zero integral). Therefore, adding the direct and indirect contribution of the jet entropy leads to the total entropy

$$S=S_A(C_0)+\left(\frac{C_0}{R^2}E_{Jet}+S_{Jet}\right)L=S_A(C_0)F_{Jet}L$$
(43)
where $`S_A(C_0)`$ is the zero order, area entropy, obtained in the limit of vanishing jet width 222 This reasoning to obtain the first order entropy can be precised by evaluating explicitly the first order modification of the Lagrange parameter $`C`$ ( let say $`C_1CC_0`$ ) due to the jet energy, and the first order modification of the Lagrange parameter $`\alpha `$ ( let say $`\alpha _1\alpha \alpha _0`$ ) due to the jet curvature and computing the first order entropy from its definition (15). We have calculated the first order entropy in this way to actually obtain (43).

We deduce from (43) that the maximization of the entropy is achieved by minimizing the total free energy $`F_{Jet}L`$, which we have proved to be positive at the end of previous sub-section. Thus we conclude that the maximum entropy state minimizes the jet length, with a given area of the subdomains (31). The subdomains shape will therefore be a circle or a stripe. More precisely if $`A_1<1/\pi `$ the jet forms a circle enclosing the positive constant stream function domain ( the jet bounds a cyclonic vortex ), if $`1/\pi <A_1<11/\pi `$ two straight lines jets form a stripe and if $`A_1>11/\pi `$ the jet form a circle enclosing the negative constant stream function domain ( the jet bounds anti-cyclonic vortex ).

The different types of solutions can be summarized in a (E,B) diagram : figure 7. The outer parabola is the maximum energy achievable for a fixed B : $`E=R^2(1B^2)/2`$. The frontier lines between the straight jets and the circular jets corresponds to $`A_1=1/\pi `$ or $`A_1=1/\pi `$. It has been calculated using (31) and (32) : $`E=R^2B^2(2\pi 2)/(\pi 2)^2`$. Note that the maximum accessible energy is in $`R^2`$, but it has been scaled by the normalization condition (9) on PV levels, so the real energy is not bounded.

All this analysis assumes that the vortex size is much larger than the jet width $`\lambda `$, given by figure 6. In other words, the area $`A_1`$ or $`A_1`$ (31) must be larger than $`(2\lambda )^2`$. This is satisfied on the right side of the dashed line represented in figure 7. The dashed line itself corresponds to the equality, and the condition of large vortex is clearly not satisfied on its left side, for low energy. The position of the dashed line depends on the numerical value of $`R`$ (it has been here numerically computed for $`R=0.03`$), and it gets closer to the origin as $`R0`$. We shall now determine the statistical equilibrium in this case of low energy.

#### 2.3.4 Axisymmetric vortices

We have noted in subsection (2.3.1) that, when for instance $`B>0`$, in the limit of small energy ($`E0`$ or equivalently $`u|B|`$, for fixed $`B`$ and $`R`$), the area $`A_1`$ occupied by $`\psi _1`$ tends to 1, the whole domain. Therefore, in this limit, the complementary area $`A_1`$ tends to 0 and the vortex becomes smaller than the deformation radius, so we can no more neglect the curvature radius of the jet.

In this limit $`u|B|`$, as the vortex has a small area with respect to the total domain, it is not affected by the boundary conditions, so it can be supposed axisymmetric. From the general Gibbs states equation (18), we deduce the axisymmetric vortex equation :

$$\frac{d^2\psi }{d\zeta ^2}\frac{1}{\zeta }\frac{d\psi }{d\zeta }+\frac{\psi }{R^2}=B\mathrm{tanh}\left(\alpha \frac{C\psi }{R^2}\right)$$
(44)
As $`R`$ will be a typical scale of the vortex, we make the change of variable,

$$\zeta =Rr;\varphi =\frac{\alpha }{C}+\frac{\psi }{R^2},$$
(45)
leading to the rescaled axisymmetric vortex equation :

$$\frac{d^2\varphi }{dr^2}\frac{1}{r}\frac{d\varphi }{dr}+\varphi +\frac{\alpha }{C}=B+\mathrm{tanh}(C\varphi ))$$
(46)
From now on, we shall consider the case $`B>0`$ (the case $`B<0`$ is just the symmetric case of a negative vortex).

For this equation to describe a localized vortex, we impose $`lim_r\mathrm{}\varphi (r)=\varphi _1\alpha /C+\psi _1/R^2`$, where $`\psi _1`$ is the positive solution of the algebraic equation (20). Since nearly the whole fluid domain is covered by the asymptotic stream function $`\psi _1`$ outside the vortex, the condition of zero total circulation $`q=\psi /R^2=0`$ imposes that $`\psi _10`$ (it is of order $`R`$), so that $`\varphi _1=\alpha /C`$, and the algebraic equation (20) then leads to:

$$\alpha =\mathrm{arg}\mathrm{tanh}(B)$$
(47)
We can thus eliminate $`\alpha `$ in (46), leading to an equation depending on two parameters, $`B`$ and $`C`$,

$$\begin{array}{c}\hfill \frac{d^2\varphi }{dr^2}=\frac{1}{r}\frac{d\varphi }{dr}\mathrm{tanh}(C\varphi )+\varphi B+\frac{\mathrm{arg}\mathrm{tanh}B}{C}\\ \hfill \frac{d\varphi }{dr}(r=0)=0and\underset{r\mathrm{}}{lim}\varphi (r)=\frac{\mathrm{arg}\mathrm{tanh}B}{C}\end{array}$$
(48)
where the regularity condition at $`r=0`$ has been included. Let us consider, as in section (2.3.2), the analogy of equation (48) with a one particle motion with ’position’ $`\varphi `$ and ’time’ $`r`$.

The last four terms on the right-hand side of (46) can be written as the derivative $`dU/d\varphi `$ of the potential,

$$U(\varphi )=\frac{\mathrm{ln}(\mathrm{cosh}(C\varphi ))}{C}\frac{\varphi ^2}{2}+\left(B\frac{\mathrm{arg}\mathrm{tanh}B}{C}\right)\varphi ,$$
(49)
(represented in figure 8), while the first term can be interpreted as a friction effect. Indeed, an integration of (48) leads to:

$$U(\varphi _1)U(\varphi (r=0))=_0^+\mathrm{}\frac{1}{r}\left(\frac{d\varphi }{dr}\right)^2𝑑r<0$$
(50)
Thus, in figure 8(a), the hatched area on the right side must be greater than the one on the left (since $`(U(\varphi _1)U(\varphi _1))>U(\varphi (r=0)U(\varphi _1)>0`$). It is clear from figure 8 that this is possible only if $`\varphi _0<0`$ and $`\alpha /C<B`$, or, using (47), $`C>\alpha /B=\mathrm{arg}\mathrm{tanh}B/B`$. The value $`C=\alpha /B`$ corresponds to the integrability condition (39) when the effect of jet curvature is neglected. This effect is now taken into account by the departure of $`C`$ from this value, which we shall denote $`\mathrm{\Delta }CC\mathrm{arg}\mathrm{tanh}B/B`$. Then $`\mathrm{\Delta }C>0`$ and we expect to recover the results of section 2.3.2 in the limit $`\mathrm{\Delta }C0`$. Moreover, we must reach a uniform stream function at large distance, solution of the algebraic equation(20), so it must have three solutions. We see in figure 8 that the corresponding $`\mathrm{\Delta }C`$ must not exceed a maximal value, denoted $`\mathrm{\Delta }C_{max}`$ .

We can prove that for any $`B>0`$ and $`\mathrm{arg}\mathrm{tanh}B/B<C<\mathrm{arg}\mathrm{tanh}B/B+\mathrm{\Delta }C_{max}`$, equation (48) has a unique solution. Such solutions have been numerically obtained for $`B=0.75`$ and $`0<\mathrm{\Delta }C<\mathrm{\Delta }C_{max}`$. Corresponding stream function profiles are shown in figure 9.

As $`\mathrm{\Delta }C`$ is decreased from $`\mathrm{\Delta }C_{max}`$ to zero, two stages can be seen in figure 9. First the maximum value for the stream function is increased while the mean width of the vortex remains of the order of $`R`$. In a second stage, when $`\mathrm{\Delta }C`$ goes to zero, as we are closer to the integrability condition for big vortices (39), $`\varphi `$ remains longer in the vicinity of $`\varphi _1`$ so the vortex size increases. Note that the energy monotically increases as $`\mathrm{\Delta }C`$ is decreased, first by an increase in the vortex maximum stream function and then by an increase in size. Finally the case of a jet with negligible curvature studied in subsection 2.3.2 is reached when $`\mathrm{\Delta }C0`$.

In conclusion, we have shown that in the limit of small energy, with fixed $`B`$ and $`R`$, the Gibbs states are approximated by axisymmetric vortices, whose radial structure depends on the parameter $`\mathrm{\Delta }C`$, which monotonically decreases from $`\mathrm{\Delta }C_{max}`$ to 0 as energy is increased.

#### 2.3.5 The linear approximation for the Gibbs states

The previous discussion of axisymmetric vortices was concerned with the limit of small energy with fixed $`B`$ and (small) fixed $`R`$. We consider now the limit of small $`E`$ and $`B`$, i.e. the neighborhood of the origin in the phase diagram of figure 7. Then from (32), $`u|B|<<1`$. Figure 6 shows that for $`|u|<<1`$, the jet width diverges and therefore the jet tends to develop on the scale of the whole domain, so the approximation of a localized jet, or an isolated axisymmetric vortex, falls down.

In this limit of small $`E`$ and $`B`$, we can however linearize the Gibbs states equation (18), following the work of Chavanis and Sommeria (1996) for Euler equation. After linearization, solutions are expressed in terms of the eigenmodes of the Laplacian, and only the first eigenmodes can be entropy maxima.

These results are unchanged by the linear deformation term $`\psi /R^2`$, so the work of Chavanis and Sommeria (1996) directly applies here. With the periodic boundary conditions, the first eigenmode of the Laplacian, a sine function of one of the coordinates, for instance $`y`$, is thus selected. This corresponds to the low energy limit of the two jet configuration shown in figure 7. The next eigenmode, in $`\mathrm{sin}(\pi x)\mathrm{sin}(\pi y)`$, has the topology of the vortex states. A competition between these two modes is expected in the neighborhood of the origin for small $`E`$ and $`B`$. Note that the range of validity of the linear approximation is limited to a smaller range of parameters than in the Euler case, and this range of validity gets smaller and smaller as $`R0`$. The dominant solution with uniform subdomains and interfacial jets relies by contrast on the tanh like relation between PV and stream function, and it is genuinely non-linear.

## 3 The channel case

We now consider the channel geometry, which represents a zonal band around a given latitude. It is then natural to introduce a beta-effect, or a mean sublayer zonal flow ( topography ), with the term $`h(y)`$ in (1). We shall study the two cases of a linear $`h(y)`$ ( beta effect and/or uniform velocity for the sublayer flow ) or a quadratic $`h(y)`$. We shall follow the presentation for the periodic boundary conditions, stressing only the new features.

### 3.1 The dynamical system

Let us consider the barotropic QG equations (1, 2 and 3) in a channel $`D=[\frac{1}{2},\frac{1}{2}]^2`$ with the velocity $`𝐯`$ tangent to the boundary for $`y=\pm \frac{1}{2}`$ and 1-periodicity in the zonal direction. Thus we choose for the boundary conditions a constant $`\psi `$, denoted $`\psi _b`$, the same on the two boundaries $`y=\pm \frac{1}{2}`$. We note that, due to these conditions, the physical momentum (6) is equal to zero. It is always possible to satisfy this condition by a change of reference frame with a zonal velocity $`V`$ such that it moves with the center of mass of the fluid layer, and a corresponding change of the deep flow, resulting in an additional beta effect $`hh+\frac{Vy}{R^2}`$.

As in section (2), we need to specify the gage constant in the stream function $`\psi `$, and we generalize the integral condition (4) as,

$$\frac{\psi }{R^2}h(y)=0.$$
(51)
The total mass $`\psi `$ is then constant in time (but not the boundary value $`\psi _b`$ in general). With these conditions, the Dirichlet problem (2) has a unique solution $`\psi `$ for a given PV field $`q`$. We note that the scale unit is chosen such that the area of $`D`$ is equal to 1.

The integral of any functions of the potential vorticity (8) is still conserved. Let in particular $`\mathrm{\Gamma }`$ be the global PV, or circulation :

$$\mathrm{\Gamma }q=_D\mathrm{\Delta }\psi d^2𝐫=_D𝐯.d𝐥$$
(52)
By contrast with the doubly periodic boundary conditions, the circulation $`\mathrm{\Gamma }`$ is not necessarily equal to zero. The expression of the energy in terms of the PV (see equation (7)) is therefore modified (due to the boundary term in the integration by parts) :

$$E=\frac{1}{2}_D\left[(\psi )^2+\frac{\psi ^2}{R^2}\right]d^2𝐫=\frac{1}{2}_D(q+h(y))\psi d^2𝐫\frac{1}{2}\mathrm{\Gamma }\psi _b$$
(53)
Due to the invariance under zonal translation of the system, another conserved quantity exists :

$$M=_Dyqd^2𝐫$$
(54)
This constant moment fixes the ‘center of mass’ latitude for the PV field.

### 3.2 General form of the Gibbs states

et us consider the statistical mechanics on a two PV level configuration : the initial states is made of patches with two levels of potential vorticity, $`q=a_1`$ and $`q=a_1`$, occupying respectively the areas $`A`$ and $`(1A)`$ in $`D`$. We keep the normalization (9) and definition (10) for $`B`$. Now, since the circulation $`\mathrm{\Gamma }`$ is non-zero, the area $`A`$ is related to $`B`$ by $`A=(1B)/2+\mathrm{\Gamma }/2`$. The boundary term in the expression of the energy (53) leads to an obvious change in the energy variation (16). Let $`\gamma `$ be the Lagrange multiplier associated with the conservation (54) of the momentum $`M`$.

Adapting the periodic case computations, we then calculate the probability equation and the Gibbs state equation :

$$p=\frac{1\mathrm{tanh}\left(\alpha ^{}\frac{C\psi }{R^2}+\gamma y\right)}{2}$$
(55)
$$q=\mathrm{\Delta }\psi +\frac{\psi }{R^2}h(y)=B\mathrm{tanh}\left(\alpha ^{}\frac{C\psi }{R^2}+\gamma y\right)$$
(56)
with $`\alpha ^{}\alpha +C\psi _b/R^2`$. These results generalize (17) and (18) of section 2.

In the case of a Gibbs state depending on $`x`$, the Lagrange parameter $`\gamma `$ is related to a zonal propagation of the equilibrium structure. The statistical theory only predicts a set of equilibria shifted in $`x`$, but introducing the result back in the dynamical equation yields the propagation. Indeed the Gibbs state equation (56) is of the form $`q=f(\psi ,y)`$, which can be inverted (as it is monotonous in $`\psi `$, see Robert and Sommeria (1991)) to yield:

$$\psi =g(q)+R^2\frac{\gamma y}{C}$$
(57)
where $`g`$ is a function of the potential vorticity. From this relation we calculate the velocity using (3) : $`𝐯=R^2\frac{\gamma }{C}\widehat{𝐱}g^{}(q)\widehat{𝐳}q`$. As PV is advected ( equation (1) ) we obtain :

$$\frac{q}{t}+R^2\frac{\gamma }{C}\frac{q}{x}=0$$
(58)
Thus the PV field is invariant in a frame propagating with the zonal speed $`V_{sr}=R^2\frac{\gamma }{C}`$.

### 3.3 The limit of small Rossby deformation radius.

In this sub-section we propose to analyze the Gibbs state equation (56) in the limit of small deformation radius ( $`R<<1`$ ). The main difference with the periodic case resides in the latudinally depending topography $`h(y)`$, resulting in two effects. Firstly the subdomains of uniform PV are no more strictly uniform, and contain a weak zonal flow. Secondly the jet curvature is no longer constant in general, but depends on the local topography. Using the same boundary layer approximation as in the periodic case, the Laplacian term in the Gibbs state equation (56) will be neglected, except possibly in an interior jet and in the vicinity of the boundaries $`y=\pm \frac{1}{2}`$ ( boundary jets ).

Outside such jets, (56) reduces to an algebraic equation, like previously :

$$\frac{\psi }{R^2}h(y)=B\mathrm{tanh}\left(\alpha ^{}C\frac{\psi }{R^2}+\gamma y\right).$$
(59)
This is like (20), replacing the constants $`B`$ by $`B+h(y)`$ and $`\alpha `$ by $`\alpha ^{}+\gamma y`$. The three solutions can be still visualized by figure 8, but the position of the straight line with respect to the tanh curve now depends on $`y`$, due to the terms $`h(y)`$ and $`\gamma y`$. We assume that this dependence is linear in $`y`$ or varies on scales much larger than $`R`$ so that the Laplacian term remains indeed negligible. The zero order Lagrange parameters $`\alpha ^{}`$, $`C`$, $`\gamma `$, involved in this expression, can be obtained by directly maximizing the entropy by the same method as in section (2.3.1). A relation between the jet curvature and topography is then obtained at first order. This approach is developed in appendix (A).

However it is more simple to proceed differently : we start from the jet equation and show that its integrability condition provides the relation between the jet curvature and topography. To catch this effect, we take into account the radius of curvature of the jet, denoted $`r`$, like in section 2.3.4, but state that $`r`$ is constant across the jet, assumed thin. From the Gibbs state equation (56), using the boundary layer approximation, we thus obtain the jet equation in term of the transverse coordinate $`\zeta `$ :

$$\frac{d^2\psi }{d\zeta ^2}ϵ\frac{1}{r}\frac{d\psi }{d\zeta }+\frac{\psi }{R^2}h(y)=B\mathrm{tanh}\left(\alpha ^{}C\frac{\psi }{R^2}+\gamma y\right).$$
(60)
We have introduced $`ϵ=\pm 1`$ to account for the direction of curvature (keeping $`r>0`$). We define $`ϵ=1`$ (respectively -1) if the curvature of the jet is such that $`\varphi _1`$ ( respectively $`\varphi _1`$ ) is in the inner part of the jet. Note that, in the case of a vortex, as in our notations $`\psi `$ is proportional to the opposite of the pressure the case $`ϵ=1`$ ( resp $`ϵ=1`$ ) corresponds to a cyclone ( resp an anticyclone ).

The algebraic equation (59) depends on three Lagrange parameters, instead of two for the periodic case of previous section, but we have three additional constraints, the condition that $`\psi `$ has the same value at the two boundaries, the circulation constraint (52) and the momentum constraint (54). This will be achieved in general by boundary jets. Let us first study the interior jets.

To study the interior jet, we make the change of variables :

$$\tau \frac{\zeta }{R};\varphi \frac{\alpha ^{}+\gamma y}{C}+\frac{\psi }{R^2}$$
(61)
We assume the variations of $`y`$ in the jet width are negligible ( $`R<<`$ scale of variation of $`h(y)`$ ), so that $`y`$ is treated as a constant. Then we obtain the jet equation :

$$\frac{d^2\varphi }{d\tau ^2}ϵ\frac{R}{r}\frac{d\varphi }{d\tau }+\varphi +\frac{\alpha ^{}+\gamma y}{C}=B+h(y)+\mathrm{tanh}(C\varphi )$$
(62)
with $`\varphi \varphi _{\pm 1}`$ for $`\tau \pm \mathrm{}`$, where again $`\varphi _{\pm 1}`$ corresponds to the solutions of the algebraic equation (59), rescaled as

$$\varphi +\frac{\alpha ^{}}{C}B+\frac{\gamma y}{C}h(y)=\mathrm{tanh}(C\varphi ).$$
(63)
Let us consider, as in section (2.3.4), the analogy of the equation (62) with the motion equation of a particle in the potential :

$$U(\varphi )=\frac{\mathrm{ln}\mathrm{cosh}(C\varphi )}{C}\frac{\varphi ^2}{2}+\left(B+h(y)\frac{\alpha ^{}+\gamma y}{C}\right)\varphi .$$
(64)
Like in section 2.3.4., integration of (62) from $`\mathrm{}`$ to $`+\mathrm{}`$ imposes the integrability condition :

$$U(\varphi _1)U(\varphi _1)=ϵ\frac{R}{r}_{\mathrm{}}^+\mathrm{}\left(\frac{d\varphi }{d\tau }\right)^2𝑑\tau $$
(65)
The second term of the l.h.s. of equation (62) can be interpreted as a friction term: if $`ϵ=1`$, the ’particle’ starting from rest at $`\varphi _1`$ can reach a state of rest at $`\varphi _1`$ only if the difference of ’potential’ corresponds to the energy loss (65) by friction (if $`ϵ=1`$ the same is true in the reversed direction).

As in the periodic case we have made the thin jet assumption $`R<<r`$, so that the friction term ( rhs of (65)) is a correction of order $`R/r`$ : $`U(\varphi _1)U(\varphi _1)=O(R/r)`$. We first neglect it to get the zero order results, so we write $`U(\varphi )=U(\varphi _1)`$. Therefore the two hatched areas in figure 8 must be equal, like in figure 2. Due to the symmetry of the tanh function, this clearly implies that the central solution of the rescaled algebraic equation (63) must be $`\varphi _0=0`$, so that $`\alpha _0^{}/C_0B+\gamma y/C_0h(y)=0`$ (denoting the zero order Lagrange parameters by the index 0). This is possible at different latitudes $`y`$ only if $`\gamma y/C_0h(y)=0`$, or is of order $`R`$ (so that it can be neglected at zero order). Then the integrability condition becomes

$$\alpha _0^{}=C_0B.$$
(66)
Furthermore $`\varphi _{\pm 1}`$ are symmetric with respect to 0, of the form $`\varphi _{\pm 1}=\pm u`$, determined by equation (33), like in section 2.3. This parameter $`u`$ is again related to the energy by (32). Finally, the terms $`\gamma y/C_0h(y)=0`$ and the curvature term disappears in the jet equation (62), which therefore reduces to (36), discussed in section 2.3.2.

The first order solution outside the jet is obtained as a small correction $`\delta \varphi _{\pm 1}`$ to the zero order solutions $`\pm u`$, with also a small correction $`\alpha _1`$ and $`C_1`$ to the parameters $`\alpha _0^{}/C_0`$ and $`C_0`$,

$$\varphi _{\pm 1}=\pm u+\delta \varphi _{\pm 1}(y),\frac{\alpha ^{}}{C}=\frac{\alpha _0^{}}{C_0}+\alpha _1,C=C_0+C_1$$
(67)
From (65), we deduce that $`U(\varphi _1)U(\varphi _1)`$ has the sign of $`ϵ`$. As may be seen in figure 2, when $`\alpha _1`$ is positive, the line $`\varphi +\frac{\alpha ^{}}{C}`$ moves upward, so that $`U(\varphi _1)<U(\varphi _1)`$. Thus $`\alpha _1`$ has the sign opposite to the sign of $`U(\varphi _1)U(\varphi _1)`$ ; and we conclude that $`ϵ\alpha _1`$ is always positive. Introducing this expansion (67) in the algebraic equation (63), using the zero order results (33) and (66), we obtain :

$$C_0\delta \varphi _{\pm 1}(y)=\frac{\alpha _1C_0+C_1[B\pm C_0u(1u^2)]\gamma y+C_0h(y)}{1C_0(1u^2)}$$
(68)
Coming back to the stream function $`\psi `$, using (61), we deduce the corresponding velocity $`𝐯`$ by differentiation with respect to $`y`$. This velocity outside the jet is zonal, along the unit vector $`\widehat{𝐱}`$, and verifies:

$$𝐯=R^2\left(\frac{\frac{dh(y)}{dy}\gamma (1u^2)}{1C_0(1u^2)}\right)\widehat{𝐱}$$
(69)
It is therefore a constant plus a term proportional to the local beta-effect $`\frac{dh}{dy}`$. Notice that the corresponding shear $`\frac{d𝐯}{dy}`$ is stronger than the deep shear $`d_d^\psi /dy^2=R^2d^2h/dy^2`$ by the factor $`[1C_0(1u^2)]^1>1`$. The integrability condition (65) now provides the curvature of the jet. We can approximate the r.h.s. of this relation, of order $`R`$, by the zero order jet profile (36), denoting :

$$e(u)\frac{1}{2}_{\mathrm{}}^+\mathrm{}\left(\frac{d\varphi }{d\tau }\right)^2𝑑\tau .$$
(70)
( see figure 6(b) ). The l.h.s. can be expanded, using (67). We first expand the expression (64) of the potential, $`U(\varphi )=U_0(\varphi )+C_1/C_0\left[\varphi \mathrm{tanh}(C_0\varphi )\mathrm{log}\mathrm{cosh}(C_0\varphi )/C_0\right]+\left[h(y)\gamma y/C_0\alpha _1\right]\varphi `$. We can approximate $`\varphi \pm u`$ in the correction terms, and expand $`U_0(\varphi )=U_0(\pm u)+dU_0/d\varphi (\pm u)`$. The zero order equilibrium condition imposes that $`dU_0/d\varphi (\pm u)=0`$, so that (65) becomes

$$ϵu\left(h(y)\frac{\gamma y}{C_0}\alpha _1\right)=e(u)\frac{R}{r}$$
(71)
This equation (71), expresses the dependence in latitude $`y`$ of the curvature radius $`r`$ of the curve on which the jet is centered, thus defining the shape of the sub-domain interface as a function of the topography.

Without topography and for $`\gamma =0`$, we get a constant jet curvature. The same result was obtained in sub-section 2.3.3 by a different argument of free energy minimization. The parameter $`u`$, related to the energy by (32) and to $`C_0`$ by (33), quantifies the strength of the jet. By contrast, the vortex area is determined by the constraint on PV patch area (parameter $`B`$), but it is also related to the jet curvature, proportional by (71) to the small shift $`\alpha _1`$ in chemical potential and temperature. Likewise the equilibrium temperature at a liquid-gas interface slightly depends on the bubble curvature, due to capillary effects.

As explained in the end of section (3.2) the parameter $`\gamma `$ is linked to the zonal propagation speed of the structure. The term $`\gamma y`$ in (71), combined with a usual beta effect (linear topography term $`h(y)`$), leads to an oscillation with latitude $`y`$ of the jet curvature $`1/r`$, i.e. a meandering jet. Another possibility is an exact compensation of the beta-effect by the $`\gamma y`$ term, leading to a propagating circular vortex, and the selection between these two alternatives is discussed in next sub-section. An oval shaped, zonally elongated vortex, such as on Jupiter, is obtained when this compensation occurs, but with an additional quadratic topography $`h(y)`$. Indeed, to get a zonally elongated vortex, supposed latitudinally centered in zero, the radius of curvature of the jet must decrease for $`y>0`$ and increase for $`y<0`$. As a consequence, we deduce from (71) that the topography must be extremal at the latitude on which the vortex is centered (it actually admits a maximum in the cyclonic case and a minimum in the anticyclonic case ). Moreover, we deduce from (72) that the surrounding flow must have a zero velocity at the latitude on which the vortex is centered and that the shear is cyclonic when the vortex is a cyclone and anticyclonic when the vortex is an anticyclone. More generally, the curvature can be related to the zonal velocity outside the jet, eliminating the topography between (72) and (71),

$$𝐯=R^2\left(\frac{\gamma }{C_0}+\frac{ϵe(u)}{u\left(1C_0(1u^2)\right)}\frac{d}{dy}\left(\frac{R}{r}\right)\right)\widehat{𝐱}$$
(72)
Let us recall our approximations. Writing down the jet equation (60), when making the boundary layer approximation, we have assumed $`R<<r`$. We also assumed that in the jet width the topography can be considered as a constant. If $`1/\sqrt{a}`$ is a typical length scale for the topography variation this gives $`aR^2<<1`$. Moreover we have assumed that the effective topography effect $`h(y)\gamma y/C_0`$ remains small along the jet. If $`L_V`$ denotes the jet extension ( for example a vortex latitudinal size ) this approximation is valid as long as $`aL_V^2<<1`$.

#### 3.3.1 Beta effect or linear topography.

Let $`h(y)\beta y`$ in the following of this section. $`\beta `$ may mimic the beta-effect or a uniform velocity in the sublayer ( but we will refer it as the beta-effect ).

A first class of equilibrium states corresponds to a single solution $`\psi (y)`$ of the algebraic equation (59). This determines a smooth zonal flow, with possibly intense jets at the boundaries $`y=\pm \frac{1}{2}`$. The solution depends on the unknown parameters $`C`$, $`\alpha ^{}`$ and $`\gamma `$, which are indirectly determined by the energy $`E`$, momentum $`M`$, and the condition $`\psi =0`$. The limit of small energy corresponds to $`C\mathrm{}`$, for which we can neglect the term $`\psi /R^2`$ on the left hand side of (59), which then reduces to $`\psi =R^2/C[\mathrm{arg}\mathrm{tanh}(\beta yB)+\gamma y+\alpha ^{}]`$. This corresponds indeed to arbitrarily small values of $`\psi `$ (small energy) as $`C\mathrm{}`$.

When the particular energy value $`E=R^2\beta ^2/24`$ is reached, a uniform PV is possible, with $`\psi /R^2=\beta y`$. Then PV mixing is complete, which clearly maximizes the mixing entropy. In this case, $`\gamma =C\beta `$, so that $`\gamma y`$ cancels the term $`C\psi /R^2`$ in (59). Physically, the uniform westward zonal velocity,

$$v_m=R^2\beta $$
(73)
tilts the free surface with uniform slope by the geostrophic balance, and the corresponding topographic beta-effect exactly balances the imposed beta-effect.

For a still higher energy, a first possibility is that again

$$\gamma =C\beta $$
(74)
so that the beta effect exactly balanced by the $`\gamma y`$ term in the jet equation of previous subsection. This cancellation is directly obtained in the general Gibbs state equations (55) and (56). Indeed the modified stream function $`\psi ^{}=\psi +R^2\beta y`$ satisfies the same equations as in the doubly-periodic case. Therefore in the limit of small $`R`$, the Gibbs states are made of subdomains with uniform $`\psi ^{}`$ (uniform PV), separated by straight zonal jets or circular vortices.

However jets or a vortex persist in this sea of uniform PV due to the constraint of energy conservation. The vortex moves westward at the same velocity $`v_m`$, according to (58) , so they are just entrained by the background flow, without relative propagation (this can be physically understood by the cancellation of the beta-effect).

The selection of the subdomain areas and PV values is given like in the periodic case of section 2, just replacing $`\psi `$ by $`\psi ^{}=\psi \beta R^2y`$. Therefore we get again probabilities $`p_{\pm 1}=(1\pm u)`$ in the two subdomains with respective areas $`A_{\pm 1}`$ given by (31), and stream function,

$$\psi _{\pm 1}=R^2(B\pm u)R^2\beta y$$
(75)
From this relation, we can calculate the energy $`E=\frac{1}{2}(\psi _1^2A_1+\psi _1^2A_1)/R^2`$, so the energy condition (32) then becomes

$$E=\frac{R^2}{2}(u^2B^2)+\frac{R^2\beta ^2}{24}$$
(76)
Therefore these solutions with canceled beta effect can be obtained only beyond a minimum energy $`R^2\beta ^2/24`$, corresponding to the potential energy of the surface tilting associated with the drift velocity $`v_m`$. Then the excess energy will control the organization in two uniform PV areas.

The shape of these subdomains can be obtained again by minimization of the jet free energy. However, unlike in the periodic case, jets occur at the boundaries $`y=\pm \frac{1}{2}`$ as well as at subdomain interfaces. Indeed, such boundary jets are in general necessary to satisfy $`𝐯=\mathrm{𝟎}`$, or equivalently that the stream function $`\psi _b`$ must be equal at the two boundaries $`y=\pm \frac{1}{2}`$. In particular, the solutions (75) necessarily involve a stream function difference (or mass flux) $`R^2\beta `$ associated with the drift velocity $`v_m`$. This stream function difference must be compensated by boundary eastward jets with opposite total mass flux. We show in Appendix B that for two PV levels with similar initial areas a single eastward jet, separating two regions of uniform PV and weak westward drift, is the selected state (instead of two opposite jets in the periodic case). In the case of a strong PV level with a small initial area, the system organizes in a circular vortex like in the periodic case. In the limit $`u|B|`$, as one of the areas $`A_{\pm 1}`$ goes to zero, the jet approximation falls down. The corresponding analysis of axisymmetric vortices and of the linear approximation for the Gibbs states, as performed in section 2, is still valid here.

Up to now we have ignored the constraint of the momentum $`M`$ (54). This constraint imposes the latitude $`y_0`$ of the equilibrium structure ( a circular patch or a zonal band with uniform PV ). For instance in the case $`B>0`$, for which $`A_1>A_1`$ (as seen from (31)), we define $`y_0_{A_1}yd^2𝐫/A_1`$. Then

$$M_Dyqd^2𝐫=_{A_1}y(B+u)d^2𝐫+_{A_1}y(Bu)d^2𝐫=2uy_0A_1$$
(77)
We thus deduce the latitudinal position of the equilibrium structure :

$$y_0=\frac{M}{2uA_1}$$
(78)
In the case of a single eastward jet, the subdomain position has been already fixed by the area ($`y_0=A_1/2`$). Then the only possibility to satisfy a moment $`M`$ different from $`uA_1^2`$ is that the jet oscillates in latitude with some amplitude $`\mathrm{\Lambda }`$ (then $`MuA_1^2u\mathrm{\Lambda }^2`$ ). This is possible if $`\gamma \beta C`$ according to (71), which becomes

$$\frac{1}{r}=b(yy_0)$$
(79)
where $`b(u(\gamma +\beta C))/(2Ce(u)R)<0`$ and $`y_0\alpha _1/b`$. This equation clearly leads to a jet oscillating around the mean latitude $`y_0`$ ( as the curvature $`r`$ is positive for $`y<y_0`$ and negative for $`y>y_0`$ ; recall that the curvature is by definition positive when positive PV is in the inner part of the jet ). Note that this oscillation propagates eastward at speed $`R^2\frac{\gamma }{C}`$ ( given by (58) ). Since $`b<0`$, $`\frac{\gamma }{C}>\beta `$, this speed is eastward with respect to the background drift $`v_m`$ (73).

#### 3.3.2 Quadratic sublayer topography

As explained in section (3.3), in the limit of small Rossby deformation radius $`R`$, the Gibbs state equation has solutions consisting of a vortex bounded by a strong jet on the scale of $`R`$. This corresponds to the case of an initial patch with strong PV and small area ( the asymmetry parameter $`B`$ is sufficiently large ) with an energy sufficiently strong to get a structure of closed jet ( see figure 7). In the presence of a moderate topography $`h(y)`$, this internal jet is no more circular but its radius of curvature $`r<<R`$ depends on $`y`$ according to (71). We have seen in previous subsection that a linear topography $`h(y)=\beta y`$ leads to jets oscillating ( or to circular jets when $`\gamma =C\beta `$ ). We shall discuss here how a quadratic term in $`h(y)`$ modifies the shape of closed jets. We therefore assume a topography $`h(y)`$ of the form :

$$h(y)ay^2+by$$
(80)
This corresponds to a uniform deep zonal shear, with velocity $`v_d=R^2d(h\beta y)/dy=2aR^2y+b\beta `$. We focus our attention on vortex solutions, seeking close curves solutions of equation (71). The vortices will be typically oval shaped as the ones seen on Jupiter. We then study how this shape ( for instance the ratio of the great axis of the oval to the small one ) depends on the topography ( sublayer flow ) and on the jet parameters. Application of these results to Great Red Spot observations will be discussed in next section.

To make equation (71) more explicit, let $`s`$ be a curvilinear parameterization of our curve, $`𝐓(s)`$ the tangent unit vector to the curve and $`\theta (s)`$ the angular function of the curve defined by $`𝐓(s)=(\mathrm{cos}\theta (s),\mathrm{sin}\theta (s))`$ for any $`s`$. Then the radius of curvature $`r`$ of the curve is linked to $`\theta (s)`$ by $`1/r=d\theta /ds`$ and (71) yields the differential equations :

$`\{\begin{array}{cc}\frac{d\theta }{dS}\hfill & =dY^2+1\hfill \\ \frac{dY}{dS}\hfill & =\mathrm{sin}\theta (S)\hfill \end{array}`$ (81)
and $`{\displaystyle \frac{dX}{dS}}=\mathrm{cos}\theta (S)`$
with $`c^{}X=x`$, $`c^{}Y=yy_0`$, $`c^{}S=s`$ ; where $`d=(e^2(u)R^2a)/(ϵ\alpha _1^3u^2)`$, $`y_0e(u)R(C_0b\gamma )/(ϵ\alpha _1C_0u)`$. The space coordinates $`X`$, $`Y`$ and $`S`$ are here non dimensional and have been obtained by dividing the real coordinates by the scale $`c^{}e(u)R/(ϵ\alpha _1u)`$. Note that as explained in section (3.3) $`ϵ\alpha _1>0`$, so that $`c^{}>0`$. We further assume that $`a>0`$, so that $`d>0`$.

We first note that the two variables $`\theta `$ and $`Y`$ are independent of $`X`$. We will therefore consider the system formed by the two first differential equations (81). It is easily verified that this system is Hamiltonian, with $`\theta `$ and $`Y`$ the two conjugate variables and

$$H\mathrm{cos}\theta d\frac{Y^3}{3}+Y$$
(82)
the Hamiltonian. Thus $`H`$ is constant on the solution curves. We look for vortex solutions of our problem (81 and 81). Thus we require $`\theta `$ to be a monotonic function of $`S`$. Moreover the curves must close, that is $`X`$ and $`Y`$ must be periodic. For symmetry reasons, it is easily verified that the solutions of (81,81) with initial conditions $`\theta (0)=\frac{\pi }{2}`$, $`Y(0)=0`$ ( $`H=0`$ ) and some $`X(0)`$ are periodic. We prove in appendix (C) that these initial conditions are the only ones leading to closed curves. We also prove that the solutions of (81 and 81) when $`d>d_{max}\frac{4}{9}`$ does not define $`\theta `$ as a monotonic function of $`S`$. They contain double points and thus are not possible solution for our problem. Once given these initial conditions, we can easily prove that the structure has both a zonal symmetry axis and a latitudinal one passing through $`y_0`$.

To study the shape of the jets, we numerically solve equations (81,81) with initial conditions : $`\theta (0)=\frac{\pi }{2},y(0)=0`$. We obtain closed curves with oval shapes, as shown in figure 10. In figure 11 we have represented the width, the length and the aspect ratio of these vortices versus the parameter $`d`$. When $`d`$ tends to $`\frac{4}{9}`$, the vortex width tends to a maximum value : $`w_{max}=\frac{3}{2}`$ whereas the length diverges. In this limit, the vortices are thus very elongated.

## 4 Application to the Jupiter’s Great Red Spot and Oval BC

In previous sections, we have found maximum entropy states with the following properties :

* The fluid domain is partitioned in two subdomains with weak velocity, separated by jets whose width scales as the Rossby deformation radius. A strong initial PV level occupying a small area mixes in a subdomain with the form of a vortex bounded by an annular jet.
* In the presence of a parabolic topography $`h(y)`$ (due to a sublayer zonal flow with uniform shear), outside the jet, exists a zonal flow (69) with uniform shear. The velocity at the latitude of the vortex center vanishes in the reference frame of the vortex.
* The curvature of the jet is linked to the topography by (71). For the parabolic topography, solutions are oval shaped vortices, symmetric in latitude and longitude.

These properties of our solutions are the main qualitative properties of the Jovian vortices. Moreover, would this description be correct, it would predict that the topography has an extremum at the center of the vortex.

Dowling and Ingersoll (1989) derive the bottom topography, using the GRS and Oval BC velocity fields obtained from cloud motion. They analyze the results in the frame of a 1-1/2 shallow water model (SW), with an active shallow layer floating on a much deeper layer. This deep layer is in steady zonal motion which acts like a topography $`h_2`$. The SW topography $`h_2`$ is defined by $`fv_d=1/R_l(\lambda )(gh_2)/\lambda `$ where $`v_d`$ is the deep layer flow and $`R_l(\lambda )`$ is the latitudinal radius of curvature of Jupiter. Dowling and Ingersoll (1989) have deduced this SW topography $`h_2`$ by assuming the conservation of the shallow water potential vorticity $`(\omega +f)/h_1`$ ( $`h_1`$ is the upper layer thickness ) along the streamlines of the observed steady vortex flow. The vorticity $`\omega `$ is deduced from the measured velocity field, and the planetary vorticity $`f`$ is known, so that the variation of $`h`$ along each streamline is deduced. The pressure field is then obtained from the Bernoulli relation and the hydrostatic balance, leading to the field $`h_2`$. The result depends on the radius of deformation of Rossby ( $`R^{}=gh_0/f_0`$, where $`h_0`$ is the mean upper layer height and $`f_0`$ the mean Coriolis parameter ), a free parameter in this analysis. Three test values have been chosen, $`R_1^{}=1700,R_2^{}=2200`$ and $`R_3^{}=2600`$ km for the GRS and $`R_1^{}=1100,R_2^{}=1600`$ and $`R_3^{}=2000`$ km for the Oval BC ( we denote by a star superscript the physical parameters, to distinguish them from the non-dimensional quantities used earlier ). The height $`h_2`$ under each vortex has been found to depend only on latitude, and has been fitted as a quartic the planetographic latitude $`\lambda `$:

$$gh_2=A_0+A_1\lambda +A_2\lambda ^2+A_3\lambda ^3+A_4\lambda ^4.$$
(83)
The values obtained for the coefficients $`A_i`$ in the vortex reference frame, for each of the vortices and for each of the values $`R_1^{},R_2^{},R_3^{}`$ are reported in table 1 of Dowling and Ingersoll (1989).

Our model is the QG limit of this shallow water system. Starting from (SW) equations, we can derive the QG equations (1) by assuming the geostrophic balance and weak free surface deformation in comparison with the mean layer thickness $`h_0`$. The validity of this QG approximation has been discussed by Dowling and Ingersoll (1989) and was found reasonably good as a first approach, although not accurate. We furthermore use the beta-plane approximation, linearizing the planetary vorticity around a reference latitude $`\lambda _0`$ ( $`\lambda _0`$ is taken to be $`23^0`$ for the GRS and $`33.5^0`$ for the Oval BC). Therefore we write $`f=f_0+\beta y`$, with $`f_0=2\mathrm{\Omega }sin\lambda _0`$ and $`\beta 2\mathrm{\Omega }\mathrm{cos}\lambda _0/r_z(\lambda _0)`$ ( $`\mathrm{\Omega }`$ is the planetary angular speed of rotation : $`2\pi /\mathrm{\Omega }=9`$ h 55 mn 29.7 s and $`r_z(\lambda _0)`$ is the zonal planetary radius, which slightly depends on the latitude $`\lambda _0`$, due to the ellipsoidal planetary shape, see formula (4) of Dowling and Ingersoll (1989) ) . We then obtain the QG potential vorticity (2) with the QG topography $`h^{}(y^{})`$ linked to the SW topography (83) by:

$$h^{}(y^{})=\frac{gh_2}{f_0R^2}+\beta y^{}$$
(84)
We have computed the QG topography (84) using results of Dowling and Ingersoll (1989) for the SW topography (83) for the three values of the Rossby deformation radius $`R_1^{},R_2^{}`$ and $`R_3^{}`$, for the GRS and for the Oval BC. The result in figure 12 shows that, for both the GRS and the Oval BC, the QG topography has an extremum at a latitude which is nearly the center of the vortex. As far as we know, this fact has not been noticed in the literature. This result is in agreement with the predictions of our model. We note moreover that the two extrema of the topography are minima, thus our model predicts anticyclonic shears around the GRS and the Oval BC, as observed. Figure 84 shows a comparison of the QG topography derived from Dowling results with a quadratic approximation, in the case $`R^{}=R_2^{}=2200`$ km. This shows that the quadratic approximation $`h^{}(y^{})=a^{}y^2`$ is a good approximation on the latitudinal extension of the GRS. This also provides values of the parameter $`a^{}`$ (80) : $`a^{}=\mathrm{9.2\; 10}^{13}`$, $`a^{}=\mathrm{7.2\; 10}^{13}`$ and $`a^{}=\mathrm{6.4\; 10}^{13}`$ km-2 s-1 for $`R^{}=R_1^{},R_2^{}`$ and $`R_3^{}`$ respectively.

Let us deduce the corresponding non-dimensional parameters. First, the PV levels were normalized by (9), so our time unit $`T^{}`$ will depend on the real PV level difference : $`(a_1^{}a_1^{})/2\frac{1}{T^{}}`$. The other parameters are the Rossby deformation radius $`R^{}`$, the segregation parameter $`u`$ and the topography coefficient $`a^{}`$ (80).

We will consider $`R^{}`$ as a free parameter and use the following data from GRS observation :

* The jet width $`l^{}`$. Let us define the width of the jet $`l^{}`$ as the width on which the jet velocity is greater than one half of the maximum velocity. We use velocity measurement within the GRS of Mitchell et al (1981). They have used small clouds as tracers to measure velocities, and have observed that the velocity is nearly tangential to ellipses. Using a grid of concentric ellipses of constant eccentricity, the velocities have been plotted with respect to the semi-major axe $`𝒜`$ of the ellipse on which the measurement point lies. Results have been fitted by a quartic in $`𝒜`$ and may be seen in figure 13. Using these results, assuming that the velocity profile in the jet is symmetric with respect to its maximum, we choose as jet width $`l^{}=\mathrm{5.6\; 10}^3`$ km. In our model, the normalized jet width $`l=l^{}/R^{}`$ can be computed from the jet equation (34) as shown in figure 6(a). This determines the parameter $`u`$ from the parameter $`R^{}`$. The corresponding theoretical jet velocity can be compared to observations in figure 13, for $`R^{}=2500`$ km (the shape is not very sensitive to this parameter). The computed results for $`u`$ versus $`R^{}`$ is shown in figure 15.
* The maximum jet velocity $`v_{max}^{}`$. We will use the value $`v_{max}^{}=110`$ ms-1 ( Mitchell et al 1981). Using (35) and giving real dimension gives :

$$v_{max}^{}=\frac{R^{}}{T^{}}\frac{d\varphi }{d\tau }|_{max}(u).$$
(85)
$`d\varphi /d\tau |_{max}(u)`$ has been obtained by solving the non-dimensional jet equation (34), and is shown in figure 6(b). This now determine $`T^{}`$ from $`R^{}`$.
* The velocity shear surrounding the vortex. The ambient zonal shear measured at the latitude of the GRS from Limaye et al. (1986) is $`\sigma ^{}=1.5e5`$ s-1. Using (69) in its dimensional form, for a quadratic topography gives :

$$a^{}=\frac{\sigma ^{}}{2R^{}}\left(1\frac{C_0}{\mathrm{cosh}^2(C_0u)}\right)$$
(86)
This permits to compute $`a^{}`$ as a function of $`R^{}`$ (since $`u`$ has been determined, as well as $`C_0`$, related to $`u`$ by (33) ).

The computed results for $`a^{}`$ and $`T^{}`$ versus $`R^{}`$ are shown in figure 14. This shows that our determination of the topography is in agreement with the QG topography deduced from the shallow water model of Dowling and Ingersoll (1989) within a factor of two. The corresponding PV level difference $`a_1^{}a_1^{}`$ is comparable to the planetary vorticity $`f_0`$ at the latitude of the center of the GRS when $`R2400`$ km. For this value of $`R^{}`$, figure 15 shows that $`u`$ is very close to 1. Furthermore, as the GRS area is very small compared to the global area of a latitudinal band centered around the GRS, the non dimensional area occupied by the positive PV is very close to 1. Using this area expression (31) we conclude that $`B`$ is very close to 1. Using the definition of $`B`$ (10), we conclude that $`a_1^{}0`$ and $`a_1^{}f_0`$. As discussed below a forcing mechanism by convective plumes incoming from the sublayer is expected to yield this result.

The shape of the jet depends on the parameter $`d=\frac{e^2(u)}{ϵ\alpha _1^3u^2}\frac{a^{}R^^2}{a_1^{}a_1^{}}`$ and on the length scale $`c^{{}_{}{}^{}}\frac{e(u)R^{}}{ϵ\alpha _1u}`$. We can determine these two parameters from the previously determined values of $`a^{}`$,$`T^{}`$ and $`u`$, and from the observed half width of the Great Red Spot : $`y_{max}^{}=4900`$km. This permits to calculate $`c^{}`$, $`\alpha _1`$ and $`d`$ versus $`R^{}`$. Figure 16 show $`d`$ versus $`R^{}`$. The dot line represents the critical value $`d=\frac{4}{9}`$ beside which a vortex solution exists. The ratio of the length to the width of the GRS is approximatively $`2`$, which would correspond to $`d=0.441`$ ( figure 11 ) ; this is very close to the critical value $`\frac{4}{9}`$. From figure 16, our model predicts that the Rossby deformation radius is $`R^{}=1800`$ km. Figure 10 shows the actual shape of the vortex for $`d=0.441`$.

However, in the jet shape analyze, to obtain (81,81) we have supposed $`a^{}L_V^^2<<(a_1^{}a_1^{})/2`$ where $`L_V^{}`$ is the maximal latitudinal extension of the vortex ( the topography part of the PV remains negligible with respect to the PV ). For the value of $`R^{}`$ calculated above, we find $`a^{}y_{max}^{}{}_{}{}^{2}=3.910^5s^1`$ whereas $`a_1^{}a_1^{}=1.310^4s^1`$. We are thus at the limit of validity of our assumption.

Now we can reverse the procedure and propose a predictive model of the Great Red Spot. Assume a steady deep zonal flow with uniform shear, $`v_d=2a^{}R^2(y^{}y_I^{})`$, vanishing at an origin $`y_I`$ depending on our reference frame (which we shall choose in order to cancel the vortex drift). This flow is presumably generated by deep thermal convection but we are concerned here only with the dynamics of the upper layer, assumed stably stratified (due to cooling by radiative effects). We model this stratified upper layer as a shallow layer with radius of deformation $`R^{}2400`$ km. This layer is submitted to a total beta effect or ’topography’ $`h(y)=ay^2(\beta +2ay_I)y`$.

Assume that PV spots with value $`f_0`$, occupying small area proportion are randomly generated in this layer. This would be the result of intense incoming thermal plumes, as recently discussed by Ingersoll et al (2000), : conservation of the absolute angular momentum during the radial expansion leads to strong decrease of the local absolute vorticity, which comes close to zero. This means that in the planetary reference frame, a local vorticity patch with value $`f_0`$ is created. The opposite vorticity is globally created by the subducting flow, but it is close to 0 due the much larger area. This gives our time unit $`(2f_0)^1`$ and $`B=12A`$.

The outcome of random PV mixing with the constraint of the conservation laws is then a zonal velocity (69) in the observed upper layer and a vortex with area (31) with velocity profile shown by the dot curve in fig 14. The vortex moves with the upper velocity at $`y=0`$, it drifts with respect to the deep layer at velocity (73) so that the beta effect is suppressed. The shape is an oval symmetric in $`x`$ and $`y`$, with aspect ratio computed from the shape parameter $`d`$ ( figure 11 ).

Note that a slightly larger area , or stronger energy could lead to very elongated vortices. Then argument of free energy minimization show that this would lead to a single eastward zonal jet. This may explain the jet observed in the Northern hemisphere of Jupiter at the same latitude as the GRS.

For smaller Jovian vortices such as the Oval BC, the size of the vortices is comparable with the Rossby deformation radius. Thus such vortices may be described as done for axisymmetric vortices ( section 2.3.4 ). This explain why such vortices do not have a quiescent core as the GRS.

Let us describe dark brown cyclonic spots ( ’Barges’ ) at $`14`$ N on Jupiter. Their first interest is to stress that cyclonic vortices embedded in cyclonic shear exist on Jupiter. Let us go further. The greater of these barges is studied from Voyager observations by Hatzes at al (1981). The meridional velocities measured at the latitude of the center of the barge ( Hatzes et al (1981), figure 7 ) show a boundary jet organization around the perimeter of the barge ( $`v_{max}=25`$ ms-1 ), see figure 1(b). The surrounding shear is such that the shear velocity at the maximum latitude of the barge is the same as the maximum jet velocity. Thus our approximation $`aL_V^2<<1`$ is not good. We can however explain the elongated shape, similar to figure 10(b) obtained for $`d`$ very close to $`d_{max}`$.

We conclude that the Gibbs state equation (56) derived from maximization of entropy of the QG model (1) is in the limit of small Rossby deformation radius a model that explains the main qualitative features of Jovian vortices. The statistical mechanics itself explains the organization of a turbulent flow in coherent structures.

## 5 Conclusion

Our first result is to provide a general explanation for the emergence and robustness of intense jets in atmospheric or oceanic turbulent flows. In the absence of topography or beta-effect turbulence mixes potential vorticity in subdomains, and such jets occur at the interface of these subdomains, with a width of the order of the deformation radius. From a thermodynamic point of view, this is like coexistence of two phases. Indeed the vortex interaction becomes short ranged in the limit of small deformation radius, and statistical mechanics leads to a thermodynamic equilibrium between two 'phases', with different concentrations of the 2 Potential Vorticity levels. Another approach leading to the same result is to consider the general partial differential equation (18) characterizing the equilibrium states. This equation reduces to the algebraic equation (20) in the limit of small deformation radius. The two uniform subdomains correspond to two solutions $`\psi _1`$ and $`\psi _1`$ of this equation. At the interface of these subdomains, the general pde reduces to the equation (36), whose solution determines the jet profile. In addition, a solvability condition of this equation confirms the relation of equilibrium between the two 'phases', which was obtained in the thermodynamic approach.

All our results have been obtained for a 2 potential vorticity level case, but cases with more levels would lead to qualitatively similar results, although the quantitative analysis would be more involved due to the additional parameters.

In the presence of beta effect or topography, for low energy, purely zonal flows with gentle variation in latitude are obtained. A critical energy is the energy of the state where the zonal flow just compensates the beta-effect. For this state the PV is strictly uniform in the whole domain. For greater energy, two well mixed domain separated by jets appears, as in the without topography case. However the PV is no more strictly uniform in the well mixed subdomain : a zonal flow exists. Also the jet curvature depends on latitude. With ordinary beta effect this yields an intense eastward jet, purely zonal or wavy depending on the constraint on the momentum $`M`$. With the quadratic beta-effect generated by a deep shear, this can produce an oval shaped vortex. This vortex then drift so as to compensate the beta-effect. In other words, in the vortex reference frame the equivalent topography $`h(y)`$ admits an extremum, and this is in agreement with the data of Dowling and Ingersoll (1989).

Our quasi analytical approach therefore explains most of the basic features of the Great Red Spot and other Jovian vortices. It can be developed into a more accurate predictive model along the following lines. First the approximation $`R<<r`$ of thin jet is convenient for a qualitative understanding but is only marginally satisfied. However this limitation can be overcomed by numerical determination of the equilibrium state equation (18) by methods like used by Turkington and Whitaker (1995) or using relaxation equations toward the maximum entropy state as described by Robert et Sommeria (1992). Furthermore, extension to the more general shallow water model is desirable, as the Rossby number ( $`0.36`$ where it a maximal ) is not very small. This can be formally achieved ( in preparation ).

Finally the results rely on an assumption of ergodicity, or complete potential vorticity mixing consistent with the constraints on the conservation laws. Various numerical and laboratory experiments in the case of Euler equations ( see e.g. Brands et al. 1989 ) indicate that mixing may not be global but be restricted to active regions. Organization into local vortices, rather that at the scale of the whole domain, is more likely with a small radius of deformation, as vortex interactions leading to coalescence are then screened. This is observed for instance in the numerical computations of Kukharkin Orszag and Yakhot (1995). By contrast, the zonal shear here promote vortex encounters ( as observed in Voyager data ) and we expect a much better relaxation toward the global statistical equilibrium, which involves always a single vortex in a given shear zone ( as it minimizes the interfacial free energy ).

## 6 Acknowledgments

The authors thank R. Robert for collaboration on statistical mechanics approach and for useful comments on the present work.

## Appendix A Determination of the Gibbs state by direct entropy maximization, in the presence of a topography ( beta-effect ) .

In section (3.3), we studied the limit of small Rossby deformation radius in the Gibbs state equation (56) by considering the jet equation (60) and its integrability condition (65). We deduced that the Gibbs states are composed of subdomains in which $`\psi `$ verifies the algebraic equation (59) separated by an interfacial jet whose curvature verifies (71). The aim of this annex is to prove that these results can be obtained by directly maximizing the entropy, adapting the method used in section (2.3.1).
Let us make the following assumptions :

1. In the limit of small Rossby deformation radius, the probability $`p`$ of finding the PV level $`a_1`$ takes two values $`p_{\pm 1}(y)`$, depending only on $`y`$. We are looking for vortex solutions. The vortex shape is described by the length $`l(y)`$ on which the probability $`p`$ takes the value $`p_1(y)`$ ( See figure 17 ).
2. The two subdomains where $`p`$ take the two values $`p_{\pm 1}(y)`$ are separated by a jet. The probabilities $`p_{\pm 1}(y)`$ are supposed to be close to the their values without topography $`p_{\pm 1}=\pm u`$, such that the free energy per unit length of the jet is well approximated by the one calculated without topography (42). If $`L_V`$ denotes the vortex size, $`1/\sqrt{a}`$ a typical length on which topography varies, we will show this approximation to be valid as soon as $`aL_V^2<<1`$.
3. The boundary conditions can be relaxed ; that is no boundary term appears in the variation of the free energy at the order considered here. ( See discussion concerning boundary jets in section (3.3).

Given these hypothesis, the Gibbs states is described by the 3 functions $`p_{\pm 1}(y)`$ and $`l(y)`$. We will determine them by maximizing the entropy $`S`$ (15) under the 3 constraints : energy (53), mass conservation (51) and momentum (54). A necessary condition for a solution to this variational problem is the existence of 3 Lagrange parameters $`C`$,$`\alpha `$ and $`\gamma `$ such that the first variations of the free energy :

$$FS\frac{C_0}{R^2}E+\alpha _0\frac{\psi }{R^2}+\gamma M$$
(87)
vanish. Using (15),(53),(51) and (54) ; (2) and (11) and the above hypothesis ( see figure 17 ), a direct calculation shows that the free energy (87) is up to a constant :

$$\begin{array}{cc}\hfill F& =_{y_{min}}^{y_{max}}\left[f(p_1(y),y)(1l(y))+f(p_1(y),y)l(y)\right]𝑑y+_{\frac{1}{2}}^{y_{min}}f(p_1(y),y)𝑑y+\hfill \\ & _{y_{max}}^{\frac{1}{2}}f(p_1(y),y)𝑑y+LF_{Jet}(u)\hfill \\ \hfill \text{with}f(p,y)& (p\mathrm{log}p+(1p)\mathrm{log}(1p))2C_0\left(p\frac{1}{2}\right)^22\left(C_0B\alpha _0\right)\left(p\frac{1}{2}\right)\hfill \\ & 2\left(C_0h(y)\gamma y\right)\left(p\frac{1}{2}\right)\hfill \end{array}$$
(88)
and where $`L`$ is the jet length, $`F_{Jet}(u)`$ is the jet free energy per unit length (42), calculated without topography.
Considering first variations of the free energy (88) under variations of $`p_1(y)`$ ( resp $`p_1(y)`$ ) proves that $`f/p(p_1(y),y)=0`$ and that $`f/p(p_1(y),y)=0`$. A direct calculation shows that :

$$2\left(p_{\pm 1}\frac{1}{2}\right)=\mathrm{tanh}\left(2C_0\left(p_{\pm 1}\frac{1}{2}\right)+C_0h(y)\gamma y+C_0B\alpha \right)$$
(89)
Using (11) and (2) ; recalling that we neglect the Laplacian term, a straight calculation shows that (89) is equivalent to the algebraic equation (59).

Let us consider now first variations of the free energy (88) under small variations $`\delta l(y)`$ of $`l(y)`$. Using that the length of the jet is given by $`L=2_{y_{min}}^{y_{max}}\sqrt{1+\frac{1}{4}\left(dl/dy\right)^2}𝑑y`$, a straightforward calculation shows that $`\delta L=_{y_{min}}^{y_{max}}\delta l(y)/r𝑑y`$ where $`r`$ is the radius of curvature of the jet. We thus deduce from first variations of the free energy (88) :

$$\frac{F_{Jet}(u)}{r}=f(p_1(y),y)f(p_1(y),y)$$
(90)
Hypothesis (2) : $`aL_V^2<<1`$ permitted us to consider $`F_{Jet}(u)`$ as independent of $`y`$. In accordance with this hypothesis we evaluate $`f(p_1(y),y)f(p_1(y),y)`$ at order zero, with $`p_{\pm 1}=\frac{1}{2}(1\pm u)`$ at this order. We obtain : $`f(p_1(y),y)f(p_1(y),y)=2u(\alpha 1+C_0h(y)\gamma y)`$. Moreover using the free energy per unit length expression (42), (37) and (38), one can show that $`F_{Jet}(u)=2e(u)C_0R`$ where $`e(u)`$ is defined by (70). These two last results show that (90) is equivalent to (71) the expression for the radius of curvature $`r`$ found by the integrability condition for the jet.

## Appendix B Boundary jets in the channel case. Beta-effect or linear sublayer topography.

Let us derive the boundary jet properties in the case of a beta-effect (or linear sublayer topography). For the sake of simplicity, we will treat the case of zero circulation $`\mathrm{\Gamma }=0`$. Let $`\varphi _b^+`$ and $`\varphi _b^{}`$ be the values of $`\varphi `$ on the $`y=\frac{1}{2}`$ and $`y=\frac{1}{2}`$ boundary respectively. The boundary jet satisfies the jet equation (36), but with boundary conditions $`\varphi (\tau =0)=\varphi _b^\pm `$ and $`\varphi (\tau +\mathrm{})=\varphi _1or\varphi _1`$. Thus, using (38), we deduce that :

$$\frac{1}{2}\left(\frac{d\varphi }{d\tau }\right)^2\left(y=\pm \frac{1}{2}\right)+U(\varphi _b^\pm )=U(\varphi _1)=U(\varphi _1)$$
(91)
(where the last equality comes from the integrability condition for the interfacial jet).

We can relate $`d\varphi /d\tau `$ to the derivative $`d\psi /d\zeta `$ normal to the boundary ( denoting the normal coordinate $`\zeta =\pm \frac{1}{2}y`$ ), using (61),

$$\frac{d\varphi }{d\tau }\left(y=\pm \frac{1}{2}\right)=\frac{1}{R}\frac{d\psi }{d\zeta }\left(y=\pm \frac{1}{2}\right)R\frac{\gamma }{C}$$
(92)
Furthermore the condition $`\mathrm{\Gamma }=0`$ imposes

$$\mathrm{\Gamma }=_{y=\frac{1}{2}}\frac{d\psi }{dy}𝑑x+_{y=\frac{1}{2}}\frac{d\psi }{dy}𝑑x=\frac{d\psi }{d\zeta }\left(y=\frac{1}{2}\right)+\frac{d\psi }{d\zeta }\left(y=\frac{1}{2}\right)=0,$$

so that $`\frac{d\psi }{d\zeta }(y=\frac{1}{2})=\frac{d\psi }{d\zeta }(y=\frac{1}{2})`$. Then (91) becomes :

$$U(\varphi _b^\pm )=U(\varphi _1)\frac{1}{2}\left(\frac{1}{R}\frac{d\psi }{d\zeta }\left(y=\frac{1}{2}\right)R\frac{\gamma }{C}\right)^2$$
(93)
so that

$$U(\varphi _b^+)=U(\varphi _b^{})$$
(94)
Note that the two values $`\varphi _b^\pm `$ cannot be equal. They must indeed satisfy, from (61) and the condition of zero mass flux ($`\psi ^+=\psi `$),

$$\varphi _b^+\varphi _b^{}=\gamma /C$$
(95)
To solve these two equations (94) and (95), let us have a look at the potential $`U`$ in figure 19. We have to compare $`\varphi _b^+\varphi _b^{}`$ and $`\varphi _1\varphi _1`$. We thus distinguish two cases. Using (95) and $`\varphi _{\pm 1}=\pm u`$, we calculate : $`(\varphi _b^+\varphi _b^{})/(\varphi _1\varphi _1)=\gamma /(2Cu)`$. We recall that $`u`$ is of order unity.

* Low beta-effect case : $`\gamma /C<2u`$.
For an x-independent statistical equilibrium, one possibility is a zonal band with $`\varphi _1`$ inside the band and $`\varphi _1`$ near both boundaries, with $`\varphi _{b1}^+`$ and $`\varphi _{b1}^{}`$ at the boundaries. The symmetric solution, with $`\varphi _{b1}^+`$ and $`\varphi _{b1}^{}`$ at the boundaries and $`\varphi _1`$ near the boundaries, has the same free energy ( due to the symmetry of the potential $`U`$ ). However the solution which maximizes the total free energy of the jets corresponds to $`\varphi _1`$ in the lower part of the domain, and $`\varphi _1`$ in the upper part, with $`\varphi _{b0}^+`$ at the upper boundary ( $`y=\frac{1}{2}`$ ) and $`\varphi _{b0}^{}`$ at the lower boundary ( $`y=\frac{1}{2}`$ ), and a single eastward interior jet.
* High beta-effect case : $`\beta >2u`$
In this case, we note that $`\varphi _b^+\varphi _b^{}>\varphi _1\varphi _1`$. Then the two equations (95) and (94) determine a unique solution for $`(\varphi _b^+,\varphi _b^{})`$ ( see figure 19 ). We are then necessarily in the case with $`\varphi =\varphi _1`$ near the boundaries $`y=\frac{1}{2}`$ and $`\varphi =\varphi _1`$ near the other one. It again involves a single eastward internal jet. This jet can oscillate in latitude due to the momentum constraint ( according to (79) ).

Once $`\varphi _b^\pm `$ are fixed, we use equation (61) on the two boundaries $`y=\pm \frac{1}{2}`$ to conclude that $`\alpha =\frac{1}{2}(\varphi _b^++\varphi _b^{})`$. Then using the integrability condition (66) we calculate :

$$\psi _b=R^2\left(B+\frac{\varphi _b^++\varphi _b^{}}{2}\right)$$
(96)
This closes the determination of our parameters.

## Appendix C Boundary condition for the equation (81,81) of the oval-shaped vortices boundary

We want to see whether equations (81,81) with $`d>0`$, defining the curve formed by the jet, admit periodic solutions both in $`x`$ and $`y`$ corresponding to vortices, or not. As stressed in section (3.3.2), equations (81) derive from the Hamiltonian (82) ; $`y`$ and $`\theta `$ are the two conjugated variables. Let us study the phase portrait of $`H`$. For $`\theta `$ in $`[0,2\pi [`$, there are 4 critical points : $`P_1=(0,1/d),P_2=(0,1/d),P_3=(\pi ,1/d),P_4=(\pi ,1/d)`$. By linearization around these fixed points, one easily prove that $`P_1`$ and $`P_3`$ are stable fixed points whereas $`P_2`$ and $`P_4`$ are hyperbolic fixed points. This permits to draw the phase portrait : figure 18. Using $`H`$ (82), we obtain that the unstable manifolds are given by $`12/(3\sqrt{d})=H`$ and $`1+2/(3\sqrt{d})=H`$ respectively. The parameter $`d`$ governs a transition of the phase space structure. This transition occurs when the two unstable manifolds merge ; this permits to compute the transition value for $`d`$ : $`d=\frac{4}{9}`$. We are looking for vortices solutions of (81,81). We recall that $`\theta `$ is the angle with the $`x`$ axe. We thus impose $`\theta `$ to be a monotonic function of $`s`$ on a trajectory. Thus areas $`c`$ of the phase portrait (figure 18 ) are forbidden ; they would correspond to $`\theta `$ varying in a finite interval. We won't address their analyze there, but oscillating jet solutions can be found in these areas. Area $`b2`$ : it can be shown that, as on such trajectories, as $`\theta `$ is not strictly increasing, double points exists, thus forbidding area $`b2`$. Conversely areas $`a`$ and $`b1`$ are admissible trajectories, giving $`y`$ as periodic functions of $`\theta `$.We have to had the condition that $`x`$ (81) must also be a periodic function of $`\theta `$ for the curve to close. Let us denote $`\mathrm{\Delta }x`$ the $`x`$ variation when $`\theta `$ runs in $`[0,2\pi ]`$. We thus impose the condition $`\mathrm{\Delta }x=0`$. Using (81) and (81) we calculate :

$$\mathrm{\Delta }x=_0^L\mathrm{cos}\theta ds=_0^{2\pi }\frac{\mathrm{cos}\theta d\theta }{dy^2(\theta )+1}=_{\frac{\pi }{2}}^{\frac{\pi }{2}}\mathrm{cos}\theta \frac{d[y^2(\theta )y^2(\theta \pi )]}{(dy^2(\theta )+1)(dy^2(\theta \pi )+1)}d\theta =0$$
(97)
( The last expression is obtained rewriting the integral as a sum on $`[\frac{\pi }{2},\frac{\pi }{2}]`$ plus a sum on $`[\frac{\pi }{2},\frac{3\pi }{2}]`$ ; and performing a variable change ). Let us study the sign of $`y^2(\theta )y^2(\theta \pi )`$. Using (82) we deduce that $`d\frac{y^3(\theta )}{3}+y(\theta )=H\mathrm{cos}\theta `$ and $`d\frac{y^3(\theta \pi )}{3}+y(\theta \pi )=H+\mathrm{cos}\theta `$. From these two relations we conclude that $`y(\theta \pi )=y(\theta )`$ implies $`\mathrm{cos}(\theta )=0`$ and that $`y(\theta \pi )=y(\theta )`$ implies $`H=0`$. Thus if $`H0`$ ; $`y^2(\theta )y^2(\theta \pi )`$ does not change sign on $`[\frac{\pi }{2},\frac{\pi }{2}]`$. Moreover, on the areas a) and b1) $`(dy^2(\theta )+1)`$ does not change sign. Thus if $`H0`$, the argument of the last integral of (97) does not change sign and $`\mathrm{\Delta }x`$ can not be zero. We thus conclude that the only solution where $`x`$ is a periodic function of $`\theta `$ is the solution corresponding to $`H=0`$. This solution is the one obtained from (81,81) with initial conditions $`y(0)=0`$ and $`\theta (0)=\frac{\pi }{2}`$. As we have previously excluded the area $`b2`$ ( when $`d>\frac{4}{9}`$ ) we conclude that no vortex solution exists when $`d>\frac{4}{9}`$. We conclude that equations (81,81) with $`d>0`$, defining the curve formed by the jet, admit one periodic solution both in $`x`$ and $`y`$ corresponding to vortices, only when $`d<d_{max}=\frac{4}{9}`$. This solution corresponds to $`H=0`$ (82). The vortex then admits a latitudinal and a zonal axes of symmetry.

Figure Captions.

Figure 1 : Annular jets observed in the atmosphere of Jupiter. a) Velocity field in the Great Red Spot of Jupiter (200 South), from Dowling and Ingersoll (1989). b) Velocity field in the cyclonic Barge of Jupiter (140 North) from Hatzes et al (1981)

Figure 2 : (a) Graphical representation of the algebraic equation (20), with the rescaled variable $`\varphi \alpha /C+\psi /R^2`$. The three solutions are at the intersection of the curve (left-hand side) and straight line (right-hand side). Here the integrability condition $`\alpha =C_0B`$ for the differential equation (36)

Figure 3 : The free energy density $`f(p)`$ (27) versus the probability $`p`$. For $`C_0>1`$ and $`(C_0B\alpha _0)`$ small enough $`f(p)`$ has two local minima and one local maximum, allowing to obtain two values $`p_{\pm 1}`$ in the maximization of entropy under constraints.

Figure 4 : The parameter $`u`$ versus the Lagrange parameter $`C_0`$, as the solution of (33).

Figure 5 : Typical stream function profile in a jet ($`u`$ = 0.75 ) versus the transverse coordinate $`\tau =\zeta /R`$ ( Top ) and corresponding velocity profile ( Bottom ).

Figure 6 : Jet properties versus the segregation parameter $`u`$. a) Jet width, defined as the width of the region with velocity greater than half the maximum jet velocity. b) Maximum velocity $`\left(d\varphi /d\tau \right)_{max}`$ and jet kinetic energy $`e(u)`$ ( dotted line).

Figure 7 : Phase diagram of the Gibbs states versus the energy $`E`$ and the asymmetry parameter B. The outer line is the maximum energy achievable for a fixed B : $`E=\frac{R^2}{2}(1B^2)`$. The frontiers line between the straight jets and the circular jets corresponds to $`A_1=1/\pi `$ or $`A_1=1/\pi `$. it as been calculated using (31) and (32) : $`E=R^2B^2(2\pi 2)/(\pi 2)^2`$. The dot line represents the frontiers between axisymmetric vortices and the circular jets. We define it as the energy value for which the circular vortex area $`A_1`$ or $`A_1`$ (31) is equal to $`(2l)^2`$, where $`l`$ is the typical jets width ( figure 6 ). Such a line depends on the numerical value of R the ratio of the Rossby deformation radius to the domain scale. It has been here numerically calculated for R = 0.03.

Figure 8 : (a) Graphical representation of the algebraic equation (20), with the rescaled variable $`\varphi \alpha /C+\psi /R^2`$, like in figure 2, but in the case of a Gibbs state with an axisymmetric vortex ($`\mathrm{\Delta }C>0`$). Then the rhs hatched area is greater than the lhs one.
(b) The corresponding potential $`U(\varphi )`$, given by (49), is asymmetric to compensate for the friction term in equation (50).

Figure 9 : Various axisymmetric stream-function profiles for decreasing $`\mathrm{\Delta }C`$ ( $`\mathrm{\Delta }C=[\mathrm{0.90.60.30.10.050.030.01}]`$ and $`B=0.75`$.

Figure 10 : a ) Typical sub-domain shape with a topography $`h(y)=ay^2`$. The parameter $`d`$ has been chosen such that the ratio of the length on the width be 2 ; as on Jupiter's GRS. b ) Typical sub-domain shape with a topography $`h(y)=ay^2`$ when the parameter $`d`$ id very close to its maximum value $`d=\frac{4}{9}`$. The shape is then very elongated, with latitudinal boundaries quasi parallel, as for instance the Jovian cyclonic vortices ( 'Barges' ) described by ( Hatzes et al 1981 ).

Figure 11 : a ) Sub-domain non dimensional length and width versus the parameter $`d`$ ( topography $`h(y)=ay^2`$ ). b ) Sub-domain aspect ratio versus the parameter $`d`$.

Figure 12 : QG topography ( units $`s^1`$ ) versus latitude computed from data of Dowling and Ingersoll (1989) : a) under the GRS ; b) under the Oval BC.

Figure 13 : Velocity profile within the GRS from ( Mitchell et al 1981 ). They have observed that the velocity is nearly tangential to ellipses. Using a grid of concentric ellipses of constant eccentricity, the velocities have been plotted with respect to the semi-major axis $`𝒜`$ of the ellipse on which the measurement point lies. Results have been fitted by a quartic in $`𝒜`$. $`l^{}`$ is the jet width, defined as the width on which the jet velocity is greater than half the maximum velocity.

Figure 14 : A) Coefficient $`a`$ for a quadratic topography $`h(y)=ay^2`$ versus $`R^{}`$ computed from our QG model (86). The three cross show the coefficient $`a`$ computed from QG topography deduced from Dowling SW observed results. B) Difference of PV levels $`a_1a_1`$ versus $`R^{}`$ computed from our QG model (85). The dot line represents the planetary vorticity $`f_0`$ at the latitude of the center of the GRS. This show that one of the PV levels may be interpreted as vorticity generated by convection-plumes from the sublayer.

Figure 15 : The segregation parameter $`u`$ versus the Rossby deformation Radius $`R^{}`$ for the GRS. $`u`$ has been computed using the actual jet maximum velocity and width ( see section (4) ).

Figure 16 : The non dimensional parameter giving the shape of the curve : $`d`$ ( see (81) ) with respect to $`R^{}`$ in our model of the GRS. The dot line represents the critical value $`d=\frac{4}{9}`$ below which a vortex solution exists. The ratio of the length to the width of the GRS is approximately $`2`$. From figure 11 we conclude that this correspond to $`d`$ very close to the critical value $`\frac{4}{9}`$. From this figure, our model predicts that the Rossby deformation radius is $`R^{}=1800`$ km ( see section (4) for comments ).

Figure 17 : Definition of $`l(y)`$.

Figure 18 : Phase portraits of the Hamiltonian H (82) for $`y_0=0`$, governing the jet shape via differential equations (81) ( two periods in $`\theta `$ ). For vortices, we are looking for periodic solutions in $`y`$. Thus only trajectories of areas a) and b) are under interest. Conversely trajectories of area c) could correspond to oscillating jets. The parameters $`d`$ governs a transition between two type of phase portraits. A) For $`d<\frac{4}{9}`$ ( here $`d=0.075`$ ), trajectories of area a) can define $`y`$ as a function of $`\theta `$ corresponding to convex vortices. B ) For $`d>\frac{4}{9}`$ ( here $`d=0.075`$ ), for trajectories of area a), the curve $`y(\theta )`$ admits double points. Thus they can not define vortex boundaries.

Figure 19 : Resolution of the equations (93 and 95). The long-doted and the doted lines represent the two cases discussed in appendix (B).

.