Stable three-dimensional vortex families consistent with Jovian observations including the Great Red Spot

Abstract Detailed observations of the velocities of Jovian vortices exist at only one height in the atmosphere, so their vertical structures are poorly understood. This motivates this study that computes stable three-dimensional, long-lived planetary vortices that satisfy the equations of motion. We solve the anelastic equations with a high-resolution pseudo-spectral method using the observed Jovian atmospheric temperatures and zonal flow. We examine several families of vortices and find that constant-vorticity vortices, which have nearly uniform vorticity as a function of height, and horizontal areas that go to zero at their tops and bottoms, converge to stable vortices that look like the Great Red Spot (GRS) and other Jovian anticyclones. In contrast, the constant-area vortices proposed in previous studies, which have nearly uniform areas as a function of height, and vertical vorticities that go to zero at their tops and bottoms, are far from equilibrium, break apart, and converge to constant-vorticity vortices. Our late-time vortices show unexpected properties. Vortices that are initially non-hollow become hollow (i.e. have local minima of vertical vorticity at their centres), which is a feature of the GRS that cannot be explained with two-dimensional simulations. The central axes of the final vortices align with the planetary spin axis even if they align initially with the local direction of gravity. We present scaling laws for how vortex properties change with the Rossby number and other non-dimensional parameters. We prove analytically that the horizontal mid-plane of a stable vortex must lie at a height above the top of the convective zone.


Introduction
Jupiter and Saturn have robust, long-lived vortices at their mid-latitudes.The Great Red Spot (GRS) of Jupiter, with current mean diameter ∼13 000 km, has survived for more than 9600 of its own vortex turn-around times since the time it became monitored continuously (starting with Dawes 1857), and, if it is the same vortex observed by Hooke (1665) and Cassini (1666), then it has survived at least 21 000 turn-around times.The GRS is not unique in its longevity.Three Jovian vortices known as the White Ovals, with diameters ∼9000 km, each survived for more than 6000 turn-around times before merging to form the current Red Oval or Little Red Spot.In contrast, the Atlantic meddies can survive for up to ∼150 vortex turn-around times (Richardson, Bower & Zenk 2000), while terrestrial tropical cyclones over open water typically survive only ∼3 vortex turn-around times.More than 200 persistent mid-latitude vortices with diameters greater than 1000 km have been identified on Jupiter, and more than 50 on Saturn (Trammell et al. 2014).These vortices are governed by the familiar laws of fluid dynamics, so their ubiquity and persistence are of interest to and should be explainable by fluid dynamicists, yet there remain many unanswered questions.Why are these vortices so prevalent in some planetary atmospheres but not others?What determines their longevity?Why are observed mid-latitude planetary vortices and Atlantic meddies predominately anticyclones, whereas the polar vortices on Jupiter and Saturn, and tropical storms on Earth, are cyclonic?For Jovian and Saturnian mid-latitude vortices where data are available, the vortices are all embedded in east-west zonal flows that have the same sign vorticity as the vortices, but are the zonal flows necessary for their formation or persistence?Must the vortices be forced continuously by underlying thermal convection to persist?
Although we now have high-resolution velocity maps (Simon et al. 2014;Wong et al. 2021) for many of the vortices due to advances in correlation imagery techniques (Asay-Davis, Shetty & Marcus 2006;Shetty, Asay-Davis & Marcus 2006;Asay-Davis et al. 2008, 2009), the maps provide only the horizontal (north-south and east-west) velocities, and are given at only one height in the atmosphere, the top of the visible clouds (which, we have argued, is near 1 bar in the atmosphere; Marcus et al. 2019), but others believe to be higher in the atmosphere (Fletcher et al. 2010(Fletcher et al. , 2016(Fletcher et al. , 2021;;Wong et al. 2021).We believe that to answer the questions posed above, we need to understand the three-dimensional structures of these vortices, which to date remain unknown and controversial.For example, it is uncertain whether the GRS lies entirely at heights above the Jovian convective zone or whether its bottom lies deep within it.Despite recent advances in observations that are designed specifically to probe multiple depths of the Jovian atmosphere, such as radio measurements with the Very Large Array radio telescope (de Pater et al. 2016), gravity field anomaly (Parisi et al. 2021) and microwave (Bolton et al. 2021) measurements by the Juno spacecraft, and infrared measurements with the James Webb Space Telescope (de Pater et al. 2022), there is no consistent three-dimensional picture of the vortices.We believe that the best way to obtain an understanding of the three-dimensional vortex dynamics is with numerical simulations such that the simulated vortices are consistent with the horizontal velocities at the one available height where they can be measured in high-resolution and with the three-dimensional equations of motion.
The few partial answers that have been found to the questions posed above were generally provided by two-dimensional numerical simulations using reduced equations of motion (shallow-water, quasi-geostrophic, etc.) (Dowling & Ingersoll 1988, 1989;Marcus 1993;Showman 2007;Li et al. 2020).However, these studies beg the question of whether vortices in two dimensions behave as they do in three dimensions (and we show in this paper that in some important and relevant ways they do not).There have been few numerical calculations of three-dimensional, Jovian-like vortices, and most used the Boussinesq approximation (Williams 1997;Hassanzadeh, Marcus & Le Gal 2012;Lemasquerier et al. 2020), which is appropriate for a fluid without large vertical density variations, but not for the Jovian atmosphere.Even if we use one of the smallest published estimates of the vertical thickness of the GRS, 190 km, the mass density within the GRS varies by a factor of 100 from its top to its bottom.(Fletcher et al. (2010) argue that the top of the GRS is at a height above 140 mbar, while Parisi et al. (2021) estimate that the bottom of the GRS is between 150 km and 375 km, which corresponds to 40 and 340 bar.Using 40 bar for the bottom, we obtain that the thickness of the GRS is 190 km, and then obtain the difference of the densities at the top and bottom from figure 2 below.)Therefore, it is not surprising that Boussinesq simulations have velocities that differ qualitatively from the observed GRS velocities at the cloud tops, and the results of these studies may not be applicable to Jovian vortices.Liu & Schneider (2010) used a global circulation model based on the hydrostatic primitive equations for a compressible ideal gas atmosphere, with 30 levels of vertical resolution, to compute zonal flows on Jupiter.Their focus was on the east-west jets, not the vortices.During the simulations, vortices formed spontaneously, but they were predominately cyclones, rather than anticyclones.Some numerical three-dimensional studies used the anelastic equations (or a combination of Boussinesq and anelastic equations) to simulate the Jovian convective zone and the stable-stratified layer at heights above it containing the visible cloud tops, and they were able to create convection-forced three-dimensional vortices at the cloud-top levels.However, these vortices survived for fewer than 5 turn-around times (Heimpel, Gastine & Wicht 2016;Yadav & Bloxham 2020).Cho et al. (2001) solved the three-dimensional quasi-geostrophic equation using contour dynamics and 8 vertical levels.Their focus was on correlating the horizontal contours of their simulation with the cloud patterns of the GRS, and they did not report the vertical structure or the longevity of their vortex.Their calculation was initialized with a large vortex, rather than creating one via convective forcing.The vortex was allowed to evolve freely in time.Morales-Juberías & Dowling (2013) followed this line of reasoning to examine a freely evolving vortex to see if the vertical vorticity (as a proxy for the Jovian clouds) would form the observed patterns of the GRS.They solved the primitive isentropic coordinate equations (with the EPIC code) initialized with a three-dimensional ellipsoidal vortex.Their vortex was forced by neither convection nor any diabatic heating from radiative transfer of latent heating or cooling from cloud formation, sublimation or evaporation.They successfully produced long-lived vortices.However, they carried out their calculations such that the lower parts of their vortices were in an ambient atmosphere that was much more stably stratified with respect to instability than the Jovian atmosphere, and they did not report on the vertical structures of their vortices (possibly because they used only 14 levels in the vertical direction of their calculation, four of which were sponge layers at the top of the calculation used to prevent numerical instabilities).Palotai, Dowling & Fletcher (2014) followed up the study of Morales-Juberías & Dowling (2013), and modelled how ammonia clouds interact with these vortices.They found ammonia clouds forming inside anticyclones, consistent with Jovian observations.
In this paper, we follow the approach taken by Morales-Juberías & Dowling (2013), and evolve initial vortices, without forcing, diabatic heating, heat transfer or any form of thermal or viscous dissipation with the equations of motion, to see if they survive, qualitatively change, or are destroyed.The goal of this study is to produce non-decaying vortices so that in follow-up studies we can use these vortices to examine their longevity and slow evolution when realistic forcing and dissipation are included.Here, our goal is to examine the vertical structures and other properties of the unforced/undissipated vortices.We use the anelastic equations, which allow for global and local convection (the latter of which we find is important for reproducing several features of Jovian vortices).The calculations are well resolved in the horizontal and vertical directions, with 385 vertical collocation points.Rather than comparing our simulated vortices to Jovian cloud patterns that they might produce, we compare their horizontal velocities and temperature fields at the cloud-top height to observations.The simulated vortices mimic several prominent features of the GRS, including the shielded vorticity, the quiescent interiors with warm, high-velocity annular rings at their peripheries, a cool top, first observed by Voyager (Flasar et al. 1981), and a warm bottom, first proposed by Marcus et al. (2013a).
Our calculations are carried out using the observed atmospheric temperature and pressure (de Pater et al. 2019;Moeckel, de Pater & DeBoer 2023), which consists of a highly stratified region at heights above a convection zone (see § 3 for details).The atmosphere is highly stratified at the top of our computational domain and has N/f 1, where N is the Brunt-Väisälä frequency of the atmosphere, and f is the local Coriolis frequency.The convection zone at the bottom of our computational domain is near adiabatic ( N/f 0).We include the convection zone in our domain because many researchers expect that the vortex bottom extends nearly to its top, which we set to be 10 bar in our calculations.Some observers believe that the bottom of the large Jovian vortices lies within the convection zone (Bolton et al. 2021).We want to allow simulations of such scenarios.In addition to the choice of hydrostatic background, our method differs from previous studies because of the numerical method that we use (see the Appendix for details).Most numerical methods used for computing atmospheric flows require small time steps when N/f is large, making the computation of Jovian vortices numerically challenging.However, our time-stepping algorithm allows us to use a relatively large time step, which makes the simulations feasible with our available computational resources.Using this method, we are able to have a large Courant number compared to previous calculations, without sacrificing accuracy or stability.
Our philosophy for computing stable vortices that look like Jovian vortices incorporates the same set of ideas that we used for computing vortices in non-forced and non-dissipated rotating, stratified Boussinesq flows (Hassanzadeh & Marcus 2013;Mahdinia et al. 2017).Those Boussinesq flows have multiple families of equilibrium vortex solutions, and within each family there is a continuum of vortex solutions.The fact that each family is a continuum not only allows the vortices to have different Rossby and Burger numbers, but also allows the vortices to have a continuous range of sizes and a continuous range of locations at different heights within the flow.To compute stable vortices in stratified, rotating Boussinesq flows, we used an initial-value code that was initialized with a vortex that was 'close to' a stable equilibrium, Typically, but not always, if the initial vortex was sufficiently 'close to' a stable equilibrium, then the flows settled into a late-time vortex whose appearance looked similar to the initial condition.The approach of the initial condition to equilibrium was facilitated by the fact that the boundaries of the computational domain were designed to absorb, rather than reflect, the Poincaré waves that were generated and radiated by the initial vortex as it adjusted rapidly to come into hydrostatic and geostrophic or cyclostrophic equilibrium.Because the interior of the flow's computational domain was dissipationless, once the rapid adjustment and shedding of Poincaré waves ended, the late-time equilibrium vortices stopped dissipating and lasted indefinitely.If the initial condition of the initial-value code was a vortex that was not sufficiently 'close to' a stable equilibrium, and if the flow remained a vortex at late times, then the late-time vortex did not look like the initial vortex.Frequently, the initial vortex split into multiple vortices, some of which survived, but none of which looked like the 984 A61-4 initial vortex.Here, in this study of Jovian-like vortices in a rotating, stratified Jovian-like atmosphere governed by the anelastic equations of motion, we take the same approach by using an initial-value code and choosing an initial condition that 'looks like' a Jovian vortex and is sufficiently 'close to' a stable equilibrium', where 'looks like' and 'close to' are defined below.
One of the goals of this paper is to show how to create an initial three-dimensional vortex that is sufficiently 'close to' a stable equilibrium of the anelastic equations of motion so that at late times it converges to a vortex similar in appearance to the initial condition and also 'looks like' a Jovian vortex because its cloud-top velocity is similar in appearance to that of an observed Jovian vortex.
The remainder of this paper is as follows.In § 2, we review the anelastic equations.In § § 3 and 4, we review the Jovian observations that we use to determine the initial atmosphere's vertical thermal structure and the structure of the east-west zonal system of winds in which we embed the initial vortex.In § 5, we show how to create initial vortices that are approximate equilibria so that they are sufficiently 'close to' an equilibrium that they evolve to a stable vortex that looks like the initial condition that itself 'looks like' a Jovian vortex.In particular, we consider four families of anticyclonic vortices, one of which was a family proposed by Parisi et al. (2021) as a model of the GRS to interpret the Juno gravitometer data.In § 6, we present the results of running the initial-value code with the families of initial anticyclones.Not all families survive.In this section, we also compare the velocity and temperature fields of the final, quasi-steady vortices to Jovian vortices.In § 7, we determine how various properties of the final vortices scale with their Rossby numbers and with their vertical aspect ratios (the vertical thickness divided by the mean horizontal diameter of the vortex).We also present a semi-analytic formula that explains the shapes of the vortices in their vertical-east/west planes.In § 8, we take a closer look at the vertical structures of the vortices, and show how local convective instabilities constrain the locations of the tops and bottoms of the vortices.A summary and discussion of our results, along with a plan for future work, are in § 9.

Useful coordinate systems
A coordinate system used commonly for studying planetary scale flows is spherical coordinates (r, θ, φ), where r is the radius, θ is the latitude, and φ is the longitude.However, planetary vortices are local structures in planetary atmospheres.The characteristic horizontal length L of a typical planetary vortex is small compared to the planet's radius.The largest planetary vortex, the GRS of Jupiter, spans ∼11 • (longitude) by 7.4 • (latitude), or 13 000 km × 9200 km (Simon et al. 2018;Wong et al. 2021), so L 13 000 km.Because L is small compared to the mean Jovian radius r 0 72 000 at the 1 bar level (near the tops of the visible clouds), we use a local Cartesian coordinate system (cf.Pedlosky 1982, chapter 6).Let (r 0 , θ 0 , φ 0 ) be the centre of the vortex of interest.Then the local Cartesian coordinates in figure 1 (2.1) Note that ẑ is a radially outward unit vector and is along the vertical direction in which the flow is stratified; x is eastwards; and ŷ = ẑ × x is northwards.When computing the The locations of the south pole and the equator are marked.The origin of both Cartesian coordinates (black dot) is at the centre of the vortex (r 0 , θ 0 , φ 0 ).The angle shown is −θ 0 because the GRS is in the southern hemisphere.For the Cartesian coordinates (x, y, z) in red, x is eastwards, ẑ is radially outwards, and ŷ = ẑ × x is along the local north direction.For the Cartesian coordinates (x c , y c , z c ) in green, xc is eastwards, ẑc is along the direction of the planetary rotational axis, and ŷc = ẑc × xc .This illustration is chosen with θ 0 < 0, which would be appropriate for the GRS of Jupiter, which lies in the southern hemisphere.
GRS at θ 0 −23 • , we can ignore Jupiter's small oblateness and neglect the curvature correction terms because (L/r 0 ) 2 0.06 in (2.1)-(2.3).To exploit further approximations in the following sections, we will also need to consider the horizontal length scale L, the characteristic distance associated with the horizontal derivative of the vortex.For example, the operator (∂ 2 /∂x 2 + ∂ 2 /∂y 2 ) acting on the velocity or pressure of a vortex will be of order 1/L 2 .For some vortices, L can be considerably different from L. For example, the GRS has a relatively quiet interior surrounded by a high-velocity collar of thickness ∼2000 km, so although L 13 000 km, we have L 2000 km.(The width of the high-velocity collar might be related to the Rossby deformation radius L R .Two-dimensional, quasi-geostrophic studies (Marcus 1988;Dowling & Ingersoll 1989;Shetty & Marcus 2010) have used a 'depth-averaged' value of L R ∼ 2000 km to try to simulate qualitatively the velocity field of the GRS at cloud-top level.However, L R ≡ H P N/f , where N is the Brunt-Väisälä frequency, H P is the local vertical pressure scale height, and f is the local value of the Coriolis parameter.Because N varies by three orders of magnitude over the height of the GRS, so does L R .See § 3.2.) Local Cartesian coordinates are useful because much of the dynamics is controlled by gravity in the z direction, which strongly stratifies the flow in that direction.However, Coriolis forces due to the planetary spin are also important, especially for low Rossby number flows.In the absence of baroclinicity and stratification, the vortices would obey the Taylor-Proudman theorem and be invariant along the direction parallel to the planet's spin axis.Therefore, it is convenient to introduce another local Cartesian coordinate system (x c , y c , z c ), which has the same origin as the (x, y, z) coordinates but is aligned with the spin of the planet (see figure 1), such that ẑc is in the direction of the planet's rotation vector, x c ≡ x, and ŷc is orthogonal to xc and ẑc such that the right-hand rule is satisfied.
To convert between the two local Cartesian coordinate systems: (2.4) y = y c sin θ 0 + z c cos θ 0 , (2.5) z = −y c cos θ 0 + z c sin θ 0 , (2.6) (2.7) (2.9) Note that the (x, y, z) coordinates of the straight line parallel to the spin axis of the planet that passes through the point (x , y , z ) obey x = x , (2.10) (2.11) In other words, along the line described by (2.10)-(2.11),x c and y c are fixed, while z c varies.We will need these relations to construct our initial conditions of the vortices.

Anelastic approximation
We denote the solutions of the steady hydrostatic (i.e. the velocity v = 0) equations with overbars.The Euler equation and ideal gas equation become d P(z) dz = − ρ(z) g, (2.12) where P is pressure, ρ is mass density, T is temperature, g is the acceleration of gravity pointing in the −ẑ direction, and R is the ideal gas constant of the hydrogen-helium mixture of the planetary atmosphere of interest, which for Jupiter at 1 bar is R = 3.60 × 10 3 m 2 s −2 K −1 (de Pater et al. 2019).For the anelastic approximation, it is useful to introduce the potential temperature Θ, with where P ref is a reference temperature, and C p is the heat capacity at constant pressure.
To define the hydrostatic solution completely, an energy equation is needed.However, as will be discussed in § 3, we replace that equation with the observed value of T(z).Using the observed T(z), all of the other barred thermodynamic quantities can computed with (2.12)-(2.15).

Anelastic equation of motion
There are different versions of the anelastic equations (Ogura & Phillips 1962;Gough 1969;Adrian 1984;Bannon 1996;Brown, Vasil & Zweibel 2012), and we chose the one developed by Bannon (1996) because it conserves energy and we have found it to be robust in astrophysical calculations (Barranco & Marcus 2005;Marcus et al. 2015Marcus et al. , 2016;;Barranco, Pei & Marcus 2018).(This equation set is very similar to the Lantz-Braginsky-Roberts (LBR) equations (Lantz 1992;Lantz & Fan 1999;Braginsky & Roberts 1995).The only difference is that our thermal computational variable is potential temperature, but the one in the LBR equations is entropy.Brown et al. (2012) show that the LBR equations capture dynamics correctly in subadiabatic atmosphere because of energy conservation, whereas some other versions of anelastic equations do not.)In the Cartesian coordinates (x, y, z) shown in figure 1, these equations are 0 = ∇ • ( ρv), (2.20) where N is the Brunt-Väisälä frequency, with and f = f ẑc is the full Coriolis vector, which is along the direction of ẑc , or f = In writing the gravitational acceleration as −gẑ, rather than a vector in the spherical radial direction, we have ignored the small curvature correction terms.We also ignored the higher-order curvature terms in the Coriolis vector and the centrifugal force (Pedlosky 1982).In (2.21), the viscous term is dropped.In (2.22), the thermal conduction and radiative transfer terms are not considered.As the system of interest has a high Reynolds number, the dissipation length scale is much smaller than the numerical resolution, and we use hyperviscosity and hyperdiffusivity to perform numerical dissipation at the smallest resolvable scales.Dropping the viscous, thermal conduction and radiative transfer terms means that the convective velocities below the top of the underlying convection zone are not computed.However, the near-adiabatic temperatures of the underlying convection zone are used in these calculations, and the Brunt-Väisälä frequency goes smoothly from the observed large value at the top of the computational domain down to zero at the top of the convective zone.(See § 3.1 and figure 2(a).)Moreover, although we have omitted the convective velocities within the convection zone itself, we do compute accurately the intermittent convective velocities above the top of the convection zone when and where the flow becomes locally superadiabatic and unstable to local convection -see § 8.We do not believe that our calculations are harmed by the lack of convective velocities below the top of the convection zone because the bottoms of our computed vortices lie well above the convection zone and never penetrate down to the top of the convective zone.The details of our numerical solver and boundary condition are discussed in the Appendix.
To compute our numerical solutions using an initial-value code in § 6, we solve numerically (2.20)-(2.24).However, in § 5, where we compute initial conditions of the vortices that are 'close to' equilibrium flows, it is more useful to replace the vector momentum equation (2.21) with an equivalent set of three scalar equations: the vertical component of the vorticity equation, the horizontal divergence of the horizontal component of the momentum equation, and the vertical component of the momentum equation.These equations are where ω is the vorticity vector, v ⊥ is the horizontal velocity, and ω ⊥ is the horizontal vorticity.Also, is the horizontal Laplacian operator, and is the horizontal divergence of any given vector V .

Our non-use of the hydrostatic energy equation
The reason why we use the observed T(z), rather than solving the energy equation to find it, is because the former is relatively easy to obtain, and the latter is difficult to solve.The equation governing the energy or temperature of the flow is where k is the coefficient of heat conduction, and C v is the heat capacity at constant volume of the fluid.Because the flow is optically thin at heights above the visible cloud tops, the radiative heat transfer responsible for heating and cooling is complex (especially in regions with different cloud decks where there are phase changes).
A second complication of this equation is that, like (2.22), it will have convective velocities in any region where N(x, y, z, t) 2 < 0. To truly simulate Jovian flows, we would need to use (3.1) with correct parameter values of k(x, y, z), the cloud compositions, and the radiation physics such that the top of the convection zone would be at 10 bar, with large convective velocities in the convection zone.Those velocities would mix the potential temperature, making the time-averaged temperature within the convection zone approximately adiabatic with N 0. (Local, intermittent convection has been inferred in Jupiter at heights above 10 bar, but it is not generally believed that the flow is fully convective at heights above 10 bar.)While experience has shown us that we can compute convection with the anelastic equations (Marcus 1978), andYadav &Bloxham (2020) and other authors have done so, we believe that the dynamics of three-dimensional vortices is sufficiently complex that this first study should be carried out with the vortices in a stable rather than turbulently convective ambient atmosphere.Also, like Morales-Juberías & Dowling (2013), we wish to explore freely evolving Jovian vortices, rather than those forced by convection.For these reasons, we artificially suppress global convection (but allow local convection) in this study by choosing a T(z) that agrees with observations and is weakly subadiabatic (i.e. with N = 0 + ) beneath the top of the convective zone at 10 bar.Thus, T(z) is close to its actual value in the convection zone, but the fluid deeper than 10 bar is not filled with turbulent convection.(This replacement is what was done in the early days of calculating stellar structure (Schwarzschild 2015).Typically, one solved the steady-state energy equation along with the other hydrostatic equations, and in any region in which the vertical gradient of the temperature made the flow unstable to convection, the temperature was replaced with an adiabat such that N = 0.) Our underlying philosophy is that some of the long-lived vortices computed with this approximation might be destroyed or modified by convection.However, our belief is that there are no families of unforced vortices that could be computed with a true convective zone that would not have a counterpart using this approximation.We suspect that vortices computed with this approximation that extend into the true convective zone would be destroyed if convection were actually present, or at least, have their bottom parts that extend into the convective zone truncated.(See our discussion in § 9.2.)

Constructing hydrostatic fields based on observations
The thermal structure of the Jovian atmosphere is crucial for the existence of Jovian vortices.However, directly simulating the radiative transfer and convective processes that determine the Jovian thermal structure is complex, as described in § 3.1.Therefore, for simplicity, we use the observed zonal-averaged values of T( P) at the GRS latitude θ 0 = 23 • S, as specified by de Pater et al. (2019) and Moeckel et al. (2023), along with (2.12) and (2.13) to obtain first guesses of P(z), P( T) and P( ρ) as shown by the blue dots in figure 2. Although the temperature data appear to be smooth, when we differentiate to

Choice of zonal flow
The Jovian east-west zonal flow v zonal ( y, z) x far from the large vortices is, by definition, averaged over longitude, and is therefore independent of x.At the height of the visible cloud tops, it has been measured repeatedly (Limaye 1986;Porco et al. 2003;Asay-Davis et al. 2011;Tollefson et al. 2017) and observed to be nearly independent of time.The observations are shown as a function of latitude θ or y in figure 3.For the calculations here, we need to know v zonal for all values of z that contain a Jovian vortex.The only direct measurement of v zonal at elevations other than the cloud tops was from the Galileo probe (Atkinson, Pollack & Seiff 1998), which measured the velocity at only one location, and that result is controversial because the location was a 'hot spot' with anomalous properties (Marcus et al. 2019).Furthermore, because the probe entered at latitude 7.3 • N jovigraphic, it measured velocities along a descent path more aligned with the y c axis than the z c axis.For low Rossby number flows, the thermal wind equation (Pedlosky 1982), or the equivalent density wind equation (Kaspi et al. 2020), can be used to determine the vertical derivative of the zonal flow.However, the use of those equations requires accurate thermal T or density ρ data as a function of ( y, z).The temperature and density measurements are sufficiently uncertain that the results are ambiguous (Fletcher et al. 2016(Fletcher et al. , 2021)).Observations based on Juno gravity measurements (Kaspi et al. 2020) suggest that the zonal flow beneath the cloud tops depends primarily on y c and has an exponential decay in the z direction, with e-folding length 1471 km, which is much greater than the vertical size ∼200 km of the Jovian vortices studied here.
Determining the exact vertical structure of the zonal flow is not the focus of this study.Therefore, for simplicity, we choose a zonal flow that is independent of z c , satisfies the steady anelastic equations, and agrees with the observations at the cloud tops.This simple form of zonal flow is consistent with the Juno gravity measurement (Kaspi et al. 2020), which shows that the e-folding depth of Jovian zonal flow is ∼1471 km, which is much bigger than the vertical size of the domain in this study.It is also consistent with the recent James Webb Space Telescope zonal wind measurement (Hueso et al. 2023), which shows that the zonal wind near the top of our vortex (∼140 mbar) is similar to the cloud-top measurement at the latitude of our interest (23 • S).
Any zonal flow v = v zonal ( y, z) x is a solution to the steady anelastic equations subject to the requirements that Therefore, the only constraint on our choice of zonal flow is that where F( y sin θ 0 − z cos θ 0 ) matches the Jovian observations in figure 3 at z = 0 (the value of z near the visible cloud tops).Two different zonal flows v zonal ( y c ) are used in this study.One corresponds to the blue curve in figure 3, which is the observed Jovian zonal velocity, and the other is an approximation (orange curve) that uses a zonal velocity with constant shear.(The observed zonal velocity is smoothed to avoid numerical instabilities.Both choices of v zonal ( y c ) are artificially made periodic near the y direction boundaries.See the Appendix.)We include a study of vortices with the constant-shear zonal flow (orange line in figure 3) because we want to determine whether the local maximum and minimum of the zonal flow have a strong effect on the ability of a vortex to be hollow (as proposed by Shetty, Asay-Davis & Marcus 2007) and whether locations where the Rossby Mach number becomes supersonic control the flow's stability (as proposed by Dowling 2014Dowling , 2019Dowling , 2020;;Afanasyev & Dowling 2022).We note that we carried out numerous simulations in which the initial condition consisted of a zonal flow, finite-amplitude 'noise', and no vortices.In no case did the noise destabilize the zonal flow, so we argue that the zonal flows used in these computations and these boundary conditions are linearly stable and also stable to a variety of finite-amplitude perturbations.We note that if the zonal flows of Jupiter extend into the convective zone, as many researchers believe, then the stability of the actual Jovian zonal flows cannot be determined by our calculations.

Initial condition set-up
5.1.Strategy The unforced and undissipated equations that govern the flow do not have a unique solution.For example, P = ρ = Θ = T = v = 0 is a solution with no zones and no vortices.Purely zonal flows v zonal ( y, z) ≡ v zonal ( y, z) x with no vortices also exist with the zonal vorticity equal to ω zonal ≡ ∇ × v zonal , and with zonal pressure Pzonal , zonal temperature Tzonal and zonal potential temperature Θzonal determined by (4.1) and (4.2) with v = v zonal ( y, z).In addition, there are families of solutions that have vortices embedded in the zonal flow.For these flows, it is useful to decompose the velocity into contributions due to the vortex as v vortex ≡ v − v zonal , with similar decompositions for pressure, temperature and potential temperature.For the vertical vorticity, the contribution due to the vortex is ω z,vortex ≡ (ω − ω zonal ) • ẑ, where ω zonal • ẑ ≡ −∂v zonal /∂y.We also define ω z,vortex,0 (x, y) ≡ ω z,vortex (x, y, z = 0).Throughout the remainder of this paper, we define z = 0 to be 1 bar, which we take to be the nominal height of the visible Jovian cloud tops where the horizontal velocity is observed.(Some observers believe that the cloud tops can be as high as 500 mbar.)

Stacked models
Here, we show how to construct an initial vortex that is 'close to' an equilibrium solution of the governing equations of motion and that 'looks like' a Jovian vortex because its ω z,vortex,0 (x, y) is in agreement with Jovian observations at the cloud-top level.To find these near equilibria, we note that the vertical velocities v z of planetary vortices are small compared to their horizontal velocities, and difficult to measure.Many studies assume v z = 0 and the quasi-geostrophic equation ( (5.1) (5.5) These equations, along with (2.23)-(2.24),are equivalent to the governing equations of motion (2.20)-(2.24)for a steady flow with v z = 0. To construct an approximate numerical solution, we begin by solving (5.1) and (5.2) for v ⊥ (x, y, z = z ) for each collocation point z .We can do this because (5.1) and ( 5.2) at one value of z are decoupled from those at other values of z.Equations (5.1) and (5.2) are equivalent to the steady Euler equations for a two-dimensional, incompressible velocity for a constant density fluid, and there are many ways to solve them.For example, contour dynamics (Shetty et al. 2006) allows the computation of steady equilibrium flows even if they are unstable.The most popular method is to solve the two-dimensional Euler equation with an initial-value solver, and let the flow come to a steady state.In general, if the zonal flow v zonal ( y, z) x is given, then a multi-parameter family of steady or slowly drifting vortices, with ω vortex embedded in the zonal flow, can be found.Once a solution is found at every collocation point in z, the two-dimensional solutions can be stacked (as described in the next subsection) on top of each other to create a steady three-dimensional solution to (5.1) and (5.2).(It does not appear to be a problem if the initial two-dimensional vortices at different heights z drift at different speeds.For all cases that we have examined so far, the initial vortex -which is not designed to be a solution to the full equations, but only to be close to an attracting solution -comes to a statistically steady vortex when observed in some Galilean frame moving in the x direction.)After finding v ⊥ (x, y, z), we substitute it into the right-hand side of (5.3) and solve for P(x, y, z).We then use (5.4) with P(x, y, z) to find Θ(x, y, z), and use Θ(x, y, z) in (2.23)-(2.24) to find ρ(x, y, z) and T(x, y, z).The only governing equation that is unsatisfied is (2.22), which is why this is only an approximate solution.In § 6, we show numerically that with some judicious choices, this approximate solution is close to an equilibrium solution, and explain why.
For the case of a vortex embedded in a uniform zonal shear flow and β ≡ 0, Moore & Saffman (1971) found a family of analytic, stable two-dimensional solutions to (5.1)-( 5.3) consisting of an elliptical vortex with uniform vertical vorticity ωz = ω 0 inside the vortex.For v zonal = −σ yx, they found where is the horizontal aspect ratio (length of of the ellipse's major diameter in x to its minor diameter in y).As σ increases, the zonal flow shear stretches the vortex in the x direction, and increases.Note that Moore-Saffman vortices, like all steady, inviscid, constant-density, two-dimensional solutions to Euler's equation, have streamlines that coincide with their vorticity isocontours, so v • ∇ω z = 0.In the limit of ω 0 → 0, we have → ∞, which is the case where velocity induced by the vortex is negligible.The flow streamlines are strictly zonal.One of the properties of Moore-Saffman vortices that we will exploit later in this section is that for constant σ and ω 0 , the ellipticity is determined by (5.6), but the diameter of the vortex is a free parameter.Therefore, regardless of the size of the vortex in (5.6), it is a solution to (5.1)-( 5.3) with β = 0.
Calculations of Moore-Saffman vortices in two dimensions show that they are robust, so they are therefore potentially useful building blocks of initial three-dimensional vortices even though they have β ≡ 0, and even though their vertical vorticities ω z,vortex are uniform and do not look like the GRS, which has a pronounced local minimum of ω z,vortex at or near its centre (i.e. the GRS is 'hollow').However, we do make one change to the Moore-Saffman vortices when we use them to build a non-equilibrium initial vortex.We smooth the outer edge of the vortex, and include an annular region of opposite-signed ω z,vortex at the outer edge, to create shielded vortices.Observations at the cloud-top level show that the Jovian vortices are shielded (Choi et al. 2007;Shetty & Marcus 2010;Grassi et al. 2018;Wong et al. 2021;Scarica et al. 2022) and that the horizontal thickness of the opposite-signed annular shield of the GRS is ∼2000 km.To create our Moore-Saffman-like vortices at the cloud-top level, we define the following elliptical coordinates for convenience: where (x 0 , y 0 ) is the centre of the vortex.We define the velocity field of the vortex as (5.9) (5.10) (5.11)where R v is the elliptical radius of the vortex (the locus where the outer annular, opposite-signed shield begins), and L r is the characteristic thickness of the annulus.We use L r = 2000 km throughout this study.The value of ω 0 is determined by via (5.6), and ω z,vortex (x 0 , y 0 ) = ω 0 .In the limit L r → ∞, the vortex given by (5.9)-(5.11) is an exact Moore-Saffman vortex with vorticity ω 0 .Of course, the three-dimensional initial vortices constructed from neither Moore-Saffman vortices nor the vortices in (5.9)-(5.11)are exact equilibria of (5.1)-(5.3),but we will show that they work well enough.

Projection methods
Two-dimensional vortices that satisfy (5.1) and (5.2) can be stacked along any arbitrary direction and satisfy the three-dimensional governing equations (5.1)-(5.4)and (2.23)-(2.24)(but not (5.5)).Two evident directions along which to have the centres of each two-dimensional vortex lie are the z axis as in figure 4(a) or the z c axis, parallel to the planet's rotation axis, as in figure 4(b).One would choose the former direction if the Rossby number were large, so that the stratification was more important than the Coriolis forces (like a tornado), and choose the latter orientation for vice versa.For a low-Rossby-number barotropic flow, the Taylor-Proudman theorem would make the flow independent of z c , and the vortex would be aligned along the z c axis.In either case, the strong stratification of the Jovian atmosphere makes the flow baroclinic, so that regardless of the orientation of the stacking, to satisfy (5.1)-(5.4)and (2.23)-(2.24),all of the two-dimensional vortices lie in x-y planes rather than x c -y c planes, and the flows have v z = 0, rather than v z c = 0. See figure 4.
Although we will show in the next section that a three-dimensional vortex is well approximated by a stack of two-dimensional vortices with v z = 0, constructing three-dimensional vortices from observations is challenging because currently the two-dimensional velocity can be measured only at the level of the visible cloud tops.The velocity elsewhere has little observational constraint.For this reason, at height z = 0, we define a reference horizontal velocity v vortex,0 (x, y) and vertical vorticity ω z,vortex,0 (x, y).
In future studies, we will set ω z,vortex,0 (x, y) equal to the observed values of the Jovian vortex that we wish to study.Here, we forego observational velocities, and instead use the modified Moore-Saffman velocity in (5.9)-(5.11)because in this study we are more concerned with finding solutions to the full three-dimensional equations of motion where a known and controllable reference ω z,vortex,0 (x, y) is given.We are particularly interested in whether the vortex's final quasi-equilibrium v vortex (x, y, z = 0) is approximately the same as its initial reference v vortex,0 (x, y).In using (5.9)-(5.11), the value of ω 0 is slaved to the value of by (5.6) and the value of σ .In this paper, we use σ = −9.1 × 10 −6 s −1 when calculating ω 0 for all of the simulations, regardless of whether they use the constant-shear zonal velocity or the Jovian zonal velocity shown in figure 3.
To create the two-dimensional velocities needed for the stacked vortex at the collocation points z / = 0, rather than using observations (which do not exist) or other two-dimensional solutions to (5.1) and (5.2) (for which we have no guidance on what solutions to compute), we create new two-dimensional velocities and vertical vorticities at the other collocation points, by projecting the reference vertical vorticity ω z,vortex,0 (x, y) ≡ ω z,vortex (x, y, z = 0) to other heights z.(The two-dimensional velocity at any z is determined uniquely by ω z,vortex (x, y, z) because ∇ ⊥ • v vortex = 0 at all z.)For example, we could create a columnar stacked vortex aligned with the z axis by defining ω z,vortex (x, y, z) = ω z,vortex,0 (x, y) for all z.Of course, the GRS is not an infinite column, and we require that it has a well-defined top and bottom.One way to do that is by defining a family of vortices aligned with the z axis (which we call the CA-Rad family): CA-Rad family ω z,vortex (x, y, z) = ω z,vortex,0 (x, y) F(z) F( 0) , (5.12) with which projects ω z,vortex,0 (x, y) along the z axis.Here, F(z) is the projection function, which in this case determines how the magnitude of the vorticity changes as a function of z.The choice of this functional form of F(z) is arbitrary, but note that an assumed Gaussian vertical dependence of the pressure, vertical vorticity or angular velocity of the ocean and planetary vortices is quite common (cf.Morales-Juberías & Dowling 2013; Yim, Billant & Ménesguen 2016;Mahdinia et al. 2017).We do not expect the ω z,vortex of the final quasi-equilibrium vortex to have an exact Gaussian vertical dependence.However, by carefully choosing F(z) and the other properties of the initial vortex, we are trying to 'nudge' the final solution of our initial-value code to a hoped-for attracting basin that looks like a Jovian vortex.We use F(z) and not F(z c ) in all projections.Because z c = z sin θ 0 , the two choices are equivalent when H top and H bot are rescaled.We have labelled the family of initial conditions represented by (5.12) as CA-Rad because it projects ω z,vortex along the planet's spherical radial axis z (as in figure 4a), while preserving the horizontal area of the vortex and decreasing the magnitude of ω z,vortex (x, y) as |z − z 0 | increases.
To construct an initial condition stacked as in figure 4(b) to create a family of vortices in the CA-Rot family, we set where we have used the fact that (2.10) and (2.11) define the locus of a point along a line parallel to the planet's spin axis z c .Rather than requiring that the tops and bottoms of vortices be defined as the locations where their ω z,vortex go to zero, we can define them by having their horizontal areas go to zero.To do that, we map the x-y plane at the cloud-top level to other depths via a hybrid coordinate (x , y ).For example, CV-Rad Family ω z,vortex (x, y, z) = ω z,vortex,0 (x , y ), (5.15) x, (5.16) produces a family of vortices, CV-Rad, oriented as in figure 4(a), where the magnitude of ω z,vortex (x, y, z) is constant along z but the horizontal area of the vortex decreases in a self-similar way from z 0 to the tops and bottoms of the vortices.A similar family of initial vortices, CV-Rot, but oriented as in figure 4(b), can be created with CV-Rot family ω z,vortex (x, y, z) = ω z,vortex,0 (x , y ), (5.18) (5.20) If we had used Moore-Saffman vortices, rather the modified velocity in (5.9)-(5.11),then the two-dimensional vortices at every height z within the CV-Rad and CV-Rot families would be exact solutions to (5.1)-( 5.3) with β = 0, rather than just good approximations, so we might expect that the initial vortices in these two families are close to the equilibrium solutions to the full equations.However, the vortices in the CA-Rad and CA-Rot families are far from equilibrium and not approximate solutions to (5.1)-(5.2).The shear σ of the zonal flow at the tops and x (km) (×10 4 ) -5 0 5 x (km) bottoms of these vortices is the same as it is at z = 0, but ω z,vortex at the tops and bottoms is much smaller than it is at z = 0. Thus we expect that as these vortices evolve in time, their tops and bottoms become stretched and elongated in the x direction, and possible destroyed.Our motivation for examining initial vortices of the CA-Rot family in the next section is that Galanti et al. (2019) and Parisi et al. (2021) used them, or vortices similar to them, as GRS models (they did not solve the governing equations of motion) to interpret Juno data.
Examples of ω z,vortex (x, y = 0, z) for vortices from the CV-Rot and CA-Rot families are shown in figure 5.
For all four families of vortices, z = z 0 is the height of the x-y plane where the initial anticyclones have their maximum values of Pvortex / ρ and |v vortex |.It is also the x-y plane where the initial CA vortices have their maximum values of | ωz |, and the x-y plane where the initial CV vortices have their maximum values of area.We define the mid-plane z mid (t) of the vortex as the height z where the anticyclone has its maximum Pvortex / ρ, and note that z mid (t) is a function of time, with z mid (t = 0) = z 0 .Typically, the x-y plane at z mid (t) is where a vortex in the CA family will have its biggest vertical vorticity, and where a vortex in the CV family will have its biggest horizontal area, so it is useful to track the flow fields in the mid-plane.

Numerical results with β ≡ 0
We examined several initial anticyclonic vortices from each of the four vortex families discussed in § 5.3, with parameter values for F(z) (in (5.13)) in the ranges 120 km ≤ H bot ≤ 180 km, 50 km ≤ H top ≤ 52.5 km and 0 ≤ z 0 ≤ 15 km (equivalently, between 500 mbar and 1 bar), and with parameter values for ω z,vortex,0 (x, y) (in (5.9)-(5.11)) in the range −0.3 ≤ Ro ≤ −0.1, or equivalently, 1.2 ≤ ≤ 1.5.We chose R v = 4600 km and L r = 2000 km.Here, the initial Rossby number is defined as Ro ≡ ω 0 /(2f sin θ 0 ).Note that ω 0 is slaved to the value of by (5.6) and the value of σ .In this paper, we use σ = −9.1 × 10 −6 s −1 when calculating ω 0 for all of the simulations, regardless of whether they use the constant-shear zonal velocity or the Jovian zonal velocity shown in figure 3. We set θ 0 = 23 • S and f = 3.5 × 10 −4 s −1 .Throughout the paper, we report time in the unit of Earth days.Note that the figures of vertical slices are stretched vertically for graphical purposes.The simulated vortices are pancake vortices with very small vertical depths.

Orientation: aligned with the z or z c axis?
The planetary vortices of interest in this study have low Rossby number (Ro 0.3), and all the vortices that we examined have their final axes aligned approximately with the planet's spin axis (z c axis), even if they were initially aligned with the local gravity (z axis).This is consistent with our previous findings (Zhang & Marcus 2022).For example, figure 6 shows how the central axis of a vortex initially in the CV-Rad family changes from the beginning to the end of a simulation.For all heights z, we define the y location of the vortex centre as We use Pvortex / ρ as the weighting factor because it is almost never negative and is a good proxy for the stream function of the horizontal velocity.If dz/dy ori = ±∞, then the central axis is aligned with the z axis, and (using (2.11)) if dz/dy ori = tan θ 0 , then it is aligned with the z c axis.Initially, the vortex in figure 6 is aligned with the z axis, but by the end of the simulation, for most values of z, the slope is dz/dy ori tan θ 0 , and the vortex is aligned with the z c axis.However, it must be noted that because the vertical thickness of the vortex is so much smaller than its diameter (∼0.01), the change in the location of the central axis during the simulation is extremely small.The orange curve defining y ori (z) in figure 6 does not align perfectly with the z c axis.A small portion of y ori (z) (approximately −20 km < z < 0 km) is unaligned.Our simulations show that whether the central axis of a vortex aligns with the direction of the local gravity (z axis) or with the z c axis depends on more than the value of the local Rossby number.However, we postpone that discussion for a future paper.All of our numerical simulations with initial vortices in the CV-Rad and CA-Rad families ended with their central axes aligned approximately with the z c axis.For this reason, we confine the remainder of this section to the evolution of vortices that are initially in the CV-Rot and CA-Rot families.6.2.CV-Rot Using (5.18), the maximum vertical vorticity ω z,vortex of the initial vortices is the same at each height, but the area of the vortex changes.We examined several initial CV-Rot vortices, and in all cases, the initial vortex converges to a large, quasi-steady vortex that drifts longitudinally.The vertical size of the vortex remains similar to its initial condition.An initially deeper CV-Rot vortex evolves to a deeper final state unless there is convectional instability locally (see § 8 for details).
To illustrate how this family of vortices evolves, § 6.2.1 shows an example of an initial vortex embedded in a constant-shear zonal velocity (case CV1).In § 6.2.2, we examine the same initial vortex embedded in the observed Jovian zonal velocity (case CV2).Both cases have H bot = 120 km, H top = 50 km, = 1.5 (equivalently, Ro = −0.1)and z 0 = 8.2 km (equivalently, 700 mbar; note that z = 0, the top of the visible clouds, corresponds to P 1 bar).One reason why we chose this value z 0 was so that there are no initial, local convective instabilities (see § 8).

Vortex embedded within a zonal velocity with constant shear (case CV1)
Figure 7 shows ω z,vortex /| f sin θ 0 | in the x-y plane at the mid-plane at four different times.Because the initial condition is not in exact equilibrium, parts of it are torn apart by the zonal flow and roll up into small vortices.However, most of the initial vortex resists the zonal shear and remains intact.By the end of the simulation, all of the small vortices have merged back into the main vortex or become sheared out in the x direction modifying the zonal flow, and there is only one large longitudinally drifting vortex.Figure 8 shows ω z,vortex /| f sin θ 0 | in the x-z plane at y = 0, which contains the central axis of the vortex.Throughout the simulation, the heights of the top and bottom of the vortex remain approximately constant.The hollowness, or local minimum, of ω z,vortex at the central rotation axis of the vortex is visible at most values of z.
The truly surprising result is that the vortex becomes hollow.The initial ω z,vortex in figures 7(a) and 8(a) is non-hollow, but in the centre and along the central spin axis of the vortices in figures 7(b-d) and 8(b-d), ω z,vortex is smaller than it is in the outer parts of the vortices.In our previous two-dimensional, quasi-geostrophic studies of vortex dynamics, we never found an initially non-hollow vortex spontaneously becoming hollow.Moreover, the only way that we found to prevent an initially hollow vortex from becoming non-hollow was by adding an artificial 'bottom topography' to the governing quasi-geostrophic equation of motion (Youssef 2000;Shetty et al. 2006).However, Barranco & Marcus (2005) found stable hollow vortices in three-dimensional protoplanetary disk simulations, which indicates that the hollow vortex structure is indeed a three-dimensional effect.
Figure 9 shows Θvortex / Θ in the x-z plane passing through the centre of the vortex.The cool top and warm bottom of the vortex occur because the major force balance in the z component of the momentum equation is hydrostatic balance:  (Marcus 1993).In the first 40 days, the vortex remains mostly intact, but casts off many filaments that roll up into small vortices that are similar to those in the CV-Rot family.After 120 days, all of the small vortices have merged with the large vortex.The centre of the vortex is prominently hollow with a local minimum of ω z,vortex at its centre.The vortex remains shielded, but most of the negative ω z,vortex has moved to the vortex's northern and southern edges.We illustrate this case because it is qualitatively similar to all the simulations that we carried out with CV-Rot vortices embedded in constant zonal shear.There is little qualitative change of the vortex after 120 days.
An anticyclone has a high pressure centre, so Θvortex < 0 at heights above the mid-plane, and Θvortex > 0 below the mid-plane to maintain the balance.The cooler temperature of the top of the GRS and other Jovian anticyclones has been observed many times (Flasar et al. 1981;Fletcher et al. 2010Fletcher et al. , 2016)).
Figure 10 shows Θvortex (x, y)/ Θ in the mid-plane of the vortex, corresponding to the horizontal white region in the middle of figure 9.It is not known whether the mid-planes of the Jovian vortices are at, below or above the visible cloud-top level at z = 0 (i.e. at ∼1 bar).Initially, the mid-planes of our vortices are at z 0 , and although the height of the mid-planes changes in time, they remain near z 0 .Figure 10 shows a thin annulus of high Θvortex / Θ at the outer edge of the vortex, approximately coincident with the shield of negative ω z .This is an unexpected result, and the figure shows that this annulus was not present in the initial conditions.Because there is no corresponding anomaly in Pvortex / P at the location of this annulus (and because Tvortex / T = Θvortex / Θ + (R/C p ) Pvortex / P), there is a high-temperature annulus of Tvortex / T at the outer edge of the vortex where the characteristic temperature of the annulus is approximately 1 K warmer than the surrounding fluid.Because this annulus was not present initially, and because N 2 > 0, the energy equation (2.22) shows that the annulus of high Θvortex was created by a downward (1) that the bright rings were created by down-welling velocities in the annuli that could not be observed directly; (2) that the adiabatic heating caused by the unseen down-welling locally dried out the atmosphere, making the annular regions free of clouds; and (3) that the 5 µm radiation from the underlying, high-temperature atmosphere, which would normally be blocked by the clouds, would be observed.Figure 10 supports these hypotheses.
Figure 11 shows the maximum values of |v vortex |, ω z,vortex , Pvortex / ρ and v z,vortex in each x-y plane as functions of z for the quasi-steady final vortex at day 120.The figure also shows ( N2 − N 2 ) as a function of z along the central axis of the vortex.Figures 11(a,c) show that the vertical profiles of |v vortex | and Pvortex / ρ are similar, as would be expected from geostrophic balance in (2.21).The mid-plane of the final vortex at 120 days is not very different from its initial location at z = 0. Figure 11(b) shows that the vertical vorticity does not have a sharp maximum at the vortex mid-plane, but is nearly constant from z = 0 down to z = −100 km.This is evidence that the final vortex is a CV-Rot vortex.Figure 11(d) shows that v z is 0.2 % of the full velocity at most in the three-dimensional simulation.Figure 11(e) shows N2 − N 2 = −(g/ Θ)(∂ Θ/∂z), which is the stratification due to the net velocity (zonal flow plus vortex) along the central axis of the vortex.The value is non-dimensionalized by ( f sin θ 0 ) 2 .When this quantity is greater than zero, it means that the vortex has locally mixed the fluid in a way to de-stratify the vortex; the mixing has increased the potential temperature at the bottom of the vortex at the expense of the potential temperature at the top of the vortex.(In other words, mixing has increased the potential mass density at the top of the vortex at the expense of the mass density at the bottom of the vortex, making it more prone to local convective instabilities.)In regions where Ñ2 vortex /g < 0, mixing has super-stratified the fluid, making it more stable to convection.

Vortex embedded in the Jovian zonal velocity (case CV2)
Note that case CV1 is conducted in a constant-shear zonal flow without the locations of maximum and minimum zonal velocity.On the other hand, Shetty et al. (2007) suggested that the latter is important to make hollow vortices in stable equilibrium in a quasi-geostrophic system with bottom topography.We repeated the numerical calculation with the same initial vortex that we used in the CV1 case, but embedded in the Jovian zonal velocity (blue curve in figure 3), rather than with the straight orange line (which has no local maximum or minimum velocities).The purpose of this simulation is to see if including the locations of maximum and minimum zonal velocity can significantly affect the simulation result.Figures 12-14 show the vertical vorticity and potential temperature of this new vortex labelled CV2.To our surprise, the two final vortices are qualitatively similar, other than the fact that the horizontal area of ω z of the vortex is smaller than that of the CV1 vortex at day 120.They are both shielded and both developed warm, thin annular rings at their outer edges like the large Jovian anticyclones.In addition, they both developed hollow interiors like the GRS.Our result indicates that the locations (or existence) of the maximum and minimum zonal velocities are not important for the vortex features that we capture in both simulations (hollowness, warm ring, etc.).Although we did not intend to recover the velocity and temperature field of any specific Jovian vortex in this study, the velocities of the CV1 and CV2 quasi-steady vortices are qualitatively similar to the GRS because of their hollow interiors.This can be seen in figure 15, which shows the ω z,vortex of the GRS and of the CV2 vortex along their north-south minor principal axes as functions of y, at the cloud-top plane.The interiors of the vortices are the regions in y where ω z,vortex > 0; the shields are the regions just outside interiors where ω z,vortex < 0. The ω z,vortex of a non-hollow shielded vortex would decrease monotonically as a function of |y| from its maximum value at the latitude of the vortex centre at y = 0 until it became negative in the shield, and then it would increase to zero at the outer edge of the shield.Both the GRS and CV2 have pronounced local minima of ω z,vortex at y = 0 surrounded in the north and south by large 'shoulders' of positive ω z,vortex .A measure of the hollowness of a vortex is the ratio of the largest value of ω z,vortex in the shoulders divided by the local minimum value of ω z,vortex at the vortex centre.Both the GRS and CV2 have similar ratios, ∼2.However, the GRS is larger than CV2, and its region of hollowness is also larger.6.3.CA-Rot Using (5.14) to create an initial CA-Rot vortex makes the horizontal area of the vertical vorticity ω z,vortex of the initial vortex the same at each height, but the magnitude of ω z,vortex decreases away from the mid-plane, vanishing at the tops and bottoms of the vortex.None of our simulations that began with CA-Rot vortices produced a quasi-steady CA-Rot shows that near the mid-plane, the vortex is de-stratified, making it more prone to convective instabilities; whereas at the vortex's top and bottom, the fluid is superstratified and more resistant to convection.Note that the mid-plane location is z = 0.58 km (975 mbar) at day 120, slightly lower than the height of its initial location at z 0 = 8.2 km (700 mbar).See text for details.vortex at the end of the simulation.The initial vortex either quickly filamented into smaller vortices due to the shear of the zonal flow being much greater than ω z,vortex at the tops and bottoms of the vortices, or broke apart due to a local convective instability (as discussed in § 8 because the top of the initial vortex is too close to the mid-plane).In all cases, we found that the flow creates one or more quasi-steady CV-Rot vortices at the end of the simulation.This is shown in figures 16 and 17 for a vortex embedded in a zonal flow with constant shear (case CA1), and in figures 18 and 19 for a vortex embedded in the Jovian zonal flow shown in figure 3 (case CA2).For both cases, H top = 120 km, H bot = 50 km, = 1.5 (equivalently, Ro = −0.1)and z 0 = 8.2 km (equivalently, P = 700 mbar).Initially, the CA-Rot vortices are column-like in the x-z plane, with the area of the vortex on each x-y plane identical.For both cases, the initial vortices fragment within the first 40 days into filaments that roll up into smaller vortices.However, during that same time, the initial vortex and the smaller vortices that it spawns evolve into CV-Rot vortices.By day 80, many of the smaller vortices have merged with the large vortex, and by day 120, only a few vortices remain unmerged.We truncated the calculations after day 120 because the purpose of the calculations is to show that the CA-Rot vortices are far from equilibrium and that they quickly evolve into CV-Rot vortices.Based on two-dimensional vortex dynamics in zonal flows, we expect that if we continue calculation, the vortices whose central latitudes are closer together than the semi-minor radius of the largest vortex will eventually merge together.

How the energy equation comes to equilibrium
In § 5.1, we showed that our initial conditions of the vortices were such that they satisfy the continuity equation (2.20), all three components of the steady momentum equation (2.21), and the ideal gas equations (2.23) and (2.24).However, because v z ≡ 0, none of our initial conditions satisfied the energy equation (2.22).Numerically, we found that with the exception of (2.22), the full equations of motion (2.20)-(2.24)were nearly satisfied because the initial characteristic times for v, ρ and P to change were long compared to the turn-around time of the vortex, but the initial characteristic time for Θ to change was much faster.Of the two terms on the far right-hand side of (2.22),only (v ⊥ • ∇ ⊥ ) Θ can initially change ∂ Θ/∂t because v z is initially zero.When the vortex reaches its final steady state, the magnitude of the (v ⊥ • ∇ ⊥ ) Θ term has decreased approximately by a factor of 3 from its initial magnitude, and is balanced primarily by the N2 v z Θ/g term.Thus even though v z is small compared to v ⊥ (see (7.1) below), v z cannot be zero and plays an essential role in creating the final steady vortex.(If v z = 0, then the only way that the steady version of (2.22) could be satisfied would be if the isocontours of Θ coincide with the two-dimensional horizontal streamlines.)7. Scaling analyses and the shape of vortices

Scaling analysis of the vertical vorticity equation
In § 5.2, we showed that a steady vortex with v z ≡ 0 has (v ⊥ • ∇ ⊥ )(ω z + βy) = 0.However, none of our computed vortices has v z ≡ 0, and it is important to understand how small (v ⊥ • ∇ ⊥ )(ω z + βy) is, how it scales with the small dimensionless parameters of the flow, such as the Rossby number Ro ≡ ω z ∞ /(2f | sin θ 0 |), and what terms balance it in the vertical component of the anelastic vorticity (2.26).To understand how the various terms in (2.26) scale, we must first understand how v z / v ⊥ scales, where angle brackets around a quantity are defined to mean the 'characteristic value' of that quantity.We also define the characteristic horizontal length of the vortex L ≡ v ⊥ 2 / ω z 2 evaluated at the vortex mid-plane.We also define the characteristic vertical size of a vortex H ≡ v ⊥ 2 / ω ⊥ 2 evaluated at the vortex mid-plane where the L 2 norms are taken over the horizontal area of the vortex.For the types of vortices considered here, which are not vertically confined by walls, the vertical and horizontal components of the velocity have the same vertical scale H, and horizontal scale L. In our calculations, we found that H is the approximate vertical scale height of ρ and P, smaller than the prescribed vertical depth H top or H bot .Our result shows clearly that the vertical characteristic length scale of Jovian vortices can be different from their depth, which is similar to the findings that the Rossby deformation radius is less than the semi-minor radius of the GRS (Marcus 1988;Dowling & Ingersoll 1989;Shetty & Marcus 2010).
Note that these scaling relationships differ from those used for quasi-geostrophic homogeneous flows, where it is assumed that the vertical scale height of the horizontal component of the velocity is much greater than that of the vertical component (see Chapter 6 of Pedlosky 1982).Our numerical calculations suggest the following scalings: Note that the scaling equations (7.4)-(7.10)are for the seven terms on the right-hand side of (2.26).This tells us that when the flow is steady, although Thus the two-dimensional streamlines and the isocontours of (ω z + βy) are nearly coincident, as argued in § 5.2.We remind the reader that this needed coincidence was the basis on which we argued that the CA family of vortices was far from equilibrium, whereas the CV family was near equilibrium.Scaling equations (7.4)-(7.10)also tell us that the dominant balance in (2.26) is between 20b) in the mid-plane.These two terms do not exactly cancel each other near the centre because they both become small there (due to the hollowness of the vortex), and the other terms in (2.26) become comparable to them.Note that for most Jovian vortices and our numerical simulations here, Ro ∼ 0.1.For the simulations presented here, H/L ∼ 0.01.The value of H/L for Jovian vortices is not well constrained by observations.Also note that L /r 0 0.03 for the largest Jovian vortices.

Shape of the vortex
The shapes of ω z,vortex of the steady vortices in the x-z plane all have a distinctive 'ice cream cone' shape (e.g.figures 8(d), 13(d), 17(d) and 19(d)).We can explain these by recalling that .11)and using the hydrostatic approximation (6.2), To eliminate ( P/ ρ) in (7.13), we use the small-Rossby number limit (quasi-geostrophic approximation) of ( 5.3) with β = 0: where A(z) has dimensions of length squared.Our numerical calculations show that √ A(z) is approximately equal to the mean cross-sectional radius of the vortex.(We first found this relationship empirically in Boussinesq vortices; Hassanzadeh et al. 2012.)Because all the final, statistically steady vortices computed here, regardless of their initial conditions, are CV-Rot vortices, it is reasonable to approximate ω z (z)| averaged as independent of z, and assign it a value ω A .Thus (7.13)-(7.15)become Given d ln Θ/dz, N2 (z), N(z)| averaged and ω A , we can solve (7.16) for √ A(z), the mean radius of the vortex, using the boundary conditions that A = 0 at the top and bottom of the vortex.Figure 21 shows the results of solving (7.16) for √ A(z) for the CV1 and CV2 vortices, and that the equation works well in describing the vortex shapes.is approximately zero.We define z cross (x, y) to be the location near the mid-plane where Θ(x, y, z) = 0 (where the red and blue lines cross in figure 22a).The requirement that at all locations in the fluid along with the hydrostatic pressure equation puts bounds on z top (x, y) and z bottom (x, y).Two integral relationships of hydrostatic equilibrium can be determined by integrating the approximate anelastic hydrostatic equation (6.2) from the bottom of the vortex at z bottom (x, y) to z cross (x, y), and from z cross  8.3) show that the areas of the two lobes (between the vertical dotted red line and the solid blue curve in figure 23c) above z cross and below z cross are the same.Let us modify the vortex in figure 23(a) for z > z cross in order to minimize z top , while keeping the vortex at z < z cross unchanged, and therefore the area of the lower lobe in figure 23(c) unchanged.The modified vortex is shown in figures 23(b,d).To maintain hydrostatic equilibrium, the areas of the upper and lower lobes in the modified vortex in figure 23(d) must also be the same.The modified vortex in figure 23(b) must have ∂Θ/∂z ≥ 0 to be stable, but to minimize z top , it needs to have ∂Θ/∂z as small as possible.Thus the vortex with the minimum z top has ∂Θ/∂z = 0 at z ≥ z cross .Our numerical simulations with initial vortices with a z top smaller than the minimum value implied by figure 23 show that the vortices either break apart or increase the heights of their tops as they evolve in time.An argument similar to that illustrated in figure 23 shows that there is maximum upper bound to z bottom , and we have validated that argument numerically.
Figure 24 is a schematic like figure 23(a).It shows that z cross must be greater than z convection for all x and y.If z cross < z convection as pictured in figure 24, then at z convection , we have Θ < Θ.This means that ∂Θ/∂z < 0 for some values of z in the range z cross < z < z convection , making the flow locally unstable to convection.Thus z cross > z convection .
Because our simulations show that z cross is approximately the height of the mid-plane of the vortex (see figure 9), this suggests that the mid-plane is always at a height above the convection zone.
9. Conclusions, discussion and future work 9.1.Summary Although the main objective of this study was not to replicate any specific planetary vortex, our simulations were carried out in the anticipation that the stable, three-dimensional vortices would qualitatively resemble Jovian vortices, including the Great Red Spot (GRS).Our focus is on understanding the vertical structure of planetary vortices, because their horizontal structures are often well constrained by observations (but typically at only one height in the atmosphere), whereas their vertical structures remain unobserved and poorly understood.Our main results are as follows.
(i) This is the first numerical calculation of Jovian anticyclonic vortices with temperatures and horizontal velocities similar to the GRS that are solutions to the three-dimensional, anelastic equations of motion.The vortices are computed with the observed atmospheric vertical thermal stratification or Brunt-Väisälä frequency N(z) so that the atmosphere is strongly stratified at its top and nearly adiabatic at its bottom.The calculations use the observed (cloud-top level) east-west zonal shear flow.The vortices are shielded (i.e.no net circulation of the vertical vorticity), and the heights of their tops and bottoms are within observational bounds.The bottoms of the stable vortices that we computed were at approximately 30 bar, or 150 km beneath the top of the visible clouds.The horizontal velocity field at the cloud-top level (the only height at which we have direct observations) agrees qualitatively with observations.The calculations were computed with no dissipation other than hyperviscosity and hyperdiffusivity (i.e.no viscosity, radiative transfer terms, thermal conductivity or Newton cooling), and after 120 days, the vortices do not appear to decay with time.(ii) One of the reasons why we are able to compute these vortices is that the initial-value code used to integrate the governing equations can take large time steps with no loss (xi) Although we do not compute whether the bottom of a vortex can survive when it is within a fully turbulent convection zone, we prove that the mid-plane of the vortex must lie at a height above the top of the convective zone.(xii) The cross-section of the vortex in the vertical-horizontal plane has an 'ice cream cone' shape that can be explained by a simple expression.(xiii) The magnitudes of the Rossby numbers in parts of the Jovian-like vortices that we have simulated are approximately from −0.3 to −0.1; so if quasi-geostrophic approximations are used to model these vortices, then the approximations would be expected to be accurate only to within 30 %.In particular, a comparison of the thermal wind equation with the curl of the momentum equation, from which it is derived, shows that if the thermal wind equation is used to approximate the relationship between the horizontal temperature gradients within the vortices and the vertical gradients of the horizontal velocities, then those approximations are accurate only to within 30 %. 9.2.Discussion and future work We recognize that our study is not exhaustive and that there may be other stable families of Jovian-like vortices.Our calculations do not include β effects or a fully turbulent convective zone, and we will address those problems in future calculations because we do not anticipate that they will be computationally challenging.Micro-scale physics, such as the formation and sublimation of clouds, with associated diabatic heating, etc., is not included in this study, which might be important.We plan to examine some of these effects by incorporating a subgrid-scale cloud model in our future calculations.
The calculations presented here use a zonal east-west flow that does not vary in the vertical direction.It is set to the observed zonal flow at ∼1 bar.In future calculations, we propose to make it a function of height, but currently there are no consistent observational data to guide our choices (Fletcher et al. 2021).However, we believe that future observations using the James Webb Space Telescope (de Pater et al. 2022) will provide zonal velocities at another height, between 50 and 700 mbar in the atmosphere, to guide us and also to provide horizontal velocity fields of the GRS and other Jovian vortices that can be used to test our numerical calculations.
The focus of our future calculations will be to determine under what circumstances and for what parameter values initial vortices in the CV-Rot family evolve into long-lived, quasi-steady vortices.For those vortices that do become quasi-steady, what properties (e.g.thickness of the annular shield, location of the top, horizontal area, height of the mid-plane) of the initial vortex are inherited by the final quasi-steady vortex, which change, and why?Can we initialize the flow with a large-area CV-Rot vortex that evolves into a long-lived vortex as large as the GRS, and with a horizontal velocity field at ∼1 bar that matches it quantitatively?Can we start with small-area CV-Rot vortices that evolve into long-lived vortices with sizes and horizontal velocity fields at ∼1 bar that match quantitatively the current Red Oval, or any of the three previous White Ovals, which are not hollow?We also plan to add radiative transfer terms to the energy equation to mimic the ∼4.5 year thermal dissipation time of the atmosphere near 1 bar to determine whether the numerically computed Jovian vortices are still long-lived.
The second reason why the computational domain in y is larger than the physically relevant domain shown in figure 7 and elsewhere is that the calculations within the computational y domain are periodic in y, but the zonal flow v zonal ( y) is not periodic.Therefore, in the regions outside the physically relevant domain, we artificially modify v zonal ( y) to make it periodic in y.This is a well-known and very versatile computational trick, which does not seem to affect the results (see Marcus 1993).
Figure 1.Two sets of Cartesian coordinates.The thin quarter-circle indicates the cloud-top level of the planet.The locations of the south pole and the equator are marked.The origin of both Cartesian coordinates (black dot) is at the centre of the vortex (r 0 , θ 0 , φ 0 ).The angle shown is −θ 0 because the GRS is in the southern hemisphere.For the Cartesian coordinates (x, y, z) in red, x is eastwards, ẑ is radially outwards, and ŷ = ẑ × x is along the local north direction.For the Cartesian coordinates (x c , y c , z c ) in green, xc is eastwards, ẑc is along the direction of the planetary rotational axis, and ŷc = ẑc × xc .This illustration is chosen with θ 0 < 0, which would be appropriate for the GRS of Jupiter, which lies in the southern hemisphere.

Figure 2 .
Figure 2. Comparison of the vertical profiles of the background thermodynamic functions based on the observed T( P) (blue dots) and the smoothed functions (red curves).The black dashed line is the top of the convection zone at 10 bar.We define z = 0 as the plane where P = 1 bar, near the top of the visible clouds.
terms due to radiative transfer and viscous heating, (3.1)

Figure 3 .
Figure 3.The observed Jovian zonal velocity (blue) at a radius of the planet R corresponding to a pressure of approximately 1 bar versus a constant-shear zonal velocity (orange) with v zonal = −σ y, where σ = −9.1 × 10 −6 s −1 .The velocity is shown as a function of planetographic latitude θ and of the north-south length y, with y = 0 chosen to be the approximate centre of the GRS.

Figure 4 .
Figure 4. Two examples of three-dimensional vortices made from stacked two-dimensional vortices (red ellipses).The local Cartesian axes with z oriented in the direction opposite to gravity are in black.The planetary rotational axis z c is in green, and θ 0 is the latitude at the centre of the vortex.(a) A vortex stacked along radial axis z.(b) A vortex stacked along rotational axis z c.Note that both images show that the two-dimensional vortices lie in x-y planes, not x c -y c planes.The vorticity of each two-dimensional vortex is parallel to the z axis, not the z c axis.For Jovian vortices, where the ratios of their vertical thicknesses to horizontal diameters are ∼0.01, the figure is highly exaggerated.It would be difficult to observe directly the differences between the vortices in the two images if plotted with the true horizontal-to-vertical aspect ratio.

Figure 5 .
Figure 5. Examples of ω z,vortex (x, y = 0, z)/( f sin θ 0 ) of two initial vortices in the x-z plane at y = 0.The horizontal axis is x, and the vertical axis is z = z c sin θ 0 .Both vortices have H top = 50 km, H bot = 120 km, ω 0 = 3.0 × 10 −5 s −1 , = 1.5, R v = 4600 km and L r = 2000 km.(a) A constant-vorticity, area-varying CV-Rot vortex.(b) A constant-area, vorticity-varying CA-Rot vortex.Note that the planetary vortices are pancake vortices.For example, our vertical domain size is ∼300 km, while the horizontal domain size is ∼100 000 km.The images are stretched along the vertical direction for graphical purposes.The vertical stretching applies to all figures of vertical slices presented.

Figure 6 .
Figure 6.The change in the orientation of the central axis, y ori (z), of a vortex from initial condition to final state.The black dashed line has slope tan θ 0 and is parallel to the planet's spin, or z c axis.The blue vertical line shows the orientation of the initial vortex: y ori ≡ 0 for all z, which is parallel to the direction of gravity or z axis.The continuous orange curve shows z as a function of y ori for the final vortex.Its central axis is aligned with the spin axis of the planet for almost all z, with the exception of the mid-plane near z = −20 km.The initial vortex is part of the CV-Rad family, is in near equilibrium, and is constructed using (5.15) with H bot = 120 km, H top = 50 km and z 0 = 8.2 km.For constructing ω z,vortex,0 (x, y), we chose Ro = −0.1, or equivalently, = 1.5.

Figure 7 .
Figure 7. Plots of ω z,vortex (x, y)/| f sin θ 0 | at the mid-plane of the vortex (case CV1).See supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.132for more details.At all times, the vortex mid-plane is defined as the height z where the anticyclone has its maximum Pvortex / ρ.The time is in units of Earth days: (a) t = 0, (b) t ∼ 40, (c) t ∼ 80, (d) t ∼ 120.The colour bar shows the value of ω z,vortex (x, y)/| f sin θ 0 |.The turn-around time of this vortex is ∼4 days; that of the GRS is ∼6 days(Marcus 1993).In the first 40 days, the vortex remains mostly intact, but casts off many filaments that roll up into small vortices that are similar to those in the CV-Rot family.After 120 days, all of the small vortices have merged with the large vortex.The centre of the vortex is prominently hollow with a local minimum of ω z,vortex at its centre.The vortex remains shielded, but most of the negative ω z,vortex has moved to the vortex's northern and southern edges.We illustrate this case because it is qualitatively similar to all the simulations that we carried out with CV-Rot vortices embedded in constant zonal shear.There is little qualitative change of the vortex after 120 days.

Figure 8 .
Figure 8. Plots of ω z,vortex (x, y = 0, z)/| f sin θ 0 | in the x-z plane of the vortex (case CV1).See supplementary movie 2 available at https://doi.org/10.1017/jfm.2024.132for more details.The colour bar is as in figure 7.In this plane at day 40, the shed filaments have mostly negative values of ω z , and the hollowness along the central axis of the vortex has developed and extends from near the top of the vortex to near the bottom of the vortex.At day 120, some of the negative ω z in the initial shield from other x-z planes has become concentrated in this plane at y = 0, and fills much of the east-west domain.Here, (a) t = 0 days, (b) t ∼ 40 days, (c) t ∼ 80 days, (d) t ∼ 120 days.

Figure 9 .
Figure 9. Plots of Θvortex (x, y = 0, z)/ Θ (with values as shown in the colour bar) in the x-z plane passing through the central axis of the vortex (case CV1).The mid-plane of the vortex lies in the horizontal band of white where Θvortex (x, y, z) 0. The high-Θvortex annular ring at the edge of the vortex can be seen in (b-d) by noting that the white horizontal band has high 'shoulders' at the two edges of the vortex.The shoulders indicate that they contain higher Θvortex / Θ than the surrounding fluid at the same z.The shoulders are not present in (a) at t = 0.The top of the vortex in (d) is at z = 60 km, which corresponds to atmospheric pressure 35 mbar.The tops of the large Jovian vortices are thought to be at a height above 140 mbar (Fletcher et al. 2010) or heights above z = 37 km in this figure.Here, (a) t = 0 days, (b) t ∼ 40 days, (c) t ∼ 80 days, (d) t ∼ 120 days.

Figure 10 .
Figure 10.Plots of Θvortex (x, y)/ Θ on the mid-plane of the vortex (case CV1).This plane corresponds the white horizontal band in the middle of the vortex in figure 9.The characteristic difference between the temperature in the high-Θvortex annular ring at the outer edge of the vortex and the ambient flow is approximately 1 K.The high-Θvortex annulus in (b-d) is not present in (a) at t = 0.The high-Θvortex annular ring at the edge of the vortex corresponds to the high 'shoulders' in figures 9(b-d).Here, (a) t = 0 days, (b) t ∼ 40 days, (c) t ∼ 80 days, (d) t ∼ 120 days.

Figure 11 .
Figure 11.Vertical profiles of the vortex in case CV1.(a) Maximum |v vortex | at each x-y plane as a function of z.(b) Maximum ω z,vortex at each x-y plane as a function of z.(c) Maximum Pvortex / ρ at each x-y plane as a function of z.(d) Maximum |v z,vortex | at each x-y plane as a function of z: v z,vortex ∞ / v ⊥,vortex ∞ × 100.(e) Plot of ( N2 − N 2 )/( f sin θ 0 ) 2 along the z axis of the vortex.Plots (a-c) are normalized by their maximum values; (d) is normalized by the maximum |v vortex | of the vortex.The green dashed line in (e) indicates where N2 − N 2 = 0. Plot (e)shows that near the mid-plane, the vortex is de-stratified, making it more prone to convective instabilities; whereas at the vortex's top and bottom, the fluid is superstratified and more resistant to convection.Note that the mid-plane location is z = 0.58 km (975 mbar) at day 120, slightly lower than the height of its initial location at z 0 = 8.2 km (700 mbar).See text for details.

Figure 12 .Figure 13 .Figure 14 .
Figure 12.Plots of ω z,vortex (x, y)/| f sin θ 0 | of the CV2 vortex plotted as we did for the CV1 vortex in figure 7. The two horizontal arrows at the right-hand side indicate the latitudes where the zonal flow in figure 3 has local maximum and minimum values.The arrows point in the direction of the local zonal flow.Here, (a) t = 0 days, (b) t ∼ 40 days, (c) t ∼ 80 days, (d) t ∼ 120 days.

Figure 15 .
Figure 15.Comparison of the ω z,vortex of the CV2 vortex at the cloud-top plane and the observed GRS along the principle, minor, north-south or y axis: (a) CV2, (b) the observed GRS profile from Wong et al. (2021).The black arrows indicate the locations of the extrema of westward/eastward zonal velocity.The hollowness ratio (defined as the maximum value of ω z,vortex in the 'shoulders' divided by its local minimum value at the vortex centre -see text for details) of the GRS and CV2 are both ∼2.The main difference between the GRS and CV2 is that the GRS is larger.
.14)For each z within the vortex, averaging (7.14) over the horizontal area within the vortex, we define A(z) with ( f sin θ 0 ) ω z (z) Figure 21.Comparing 2 √ A(z) (dashed line) to the vertical vorticity ω z,vortex (x, y = 0, z)/| f sin θ 0 | of the (a) CV1 and (b) CV2 vortices in the x-z plane going through the centre of each vortex.The factor of 2 is to make the boundary close to the visual boundary of the vortex.
Figure23(a)  shows schematically the potential temperature Θ of an anticyclone similar to the one in figure22(a).Figure23(c)shows the corresponding Θ/ Θ. Equations (8.2) and (8.3) show that the areas of the two lobes (between the vertical dotted red line and the solid blue curve in figure23c) above z cross and below z cross are the same.Let us modify the vortex in figure23(a) for z > z cross in order to minimize z top , while keeping the vortex at z < z cross unchanged, and therefore the area of the lower lobe in figure23(c) unchanged.The modified vortex is shown in figures 23(b,d).To maintain hydrostatic equilibrium, the areas of the upper and lower lobes in the modified vortex in figure23(d) must also be the same.The modified vortex in figure 23(b) must have ∂Θ/∂z ≥ 0 to be stable, but to minimize z top , it needs to have ∂Θ/∂z as small as possible.Thus the vortex with the minimum z top has ∂Θ/∂z = 0 at z ≥ z cross .Our numerical simulations with initial vortices with a z top smaller than the minimum value implied by figure23show that the vortices either break apart or increase the heights of their tops as they evolve in time.An argument similar to that illustrated in figure23shows that there is maximum upper bound to z bottom , and we have validated that argument numerically.Figure24is a schematic like figure23(a).It shows that z cross must be greater than z convection for all x and y.If z cross < z convection as pictured in figure24, then at z convection ,

Figure 24 .
Figure24.Schematic with curves for Θ(z) and Θ(z), and lines for z cross and z convection , plotted as in figure23(a).If z cross < z convection , then Θ(z convection ) < Θ(z cross ), so ∂Θ/∂z < 0 somewhere in the range z cross < z < z convection , and the vortex is locally unstable to convection.