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.

1. Introduction

Jupiter and Saturn have robust, long-lived vortices at their mid-latitudes. The Great Red Spot (GRS) of Jupiter, with current mean diameter ${\sim }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 $\sim$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 $\sim$150 vortex turn-around times (Richardson, Bower & Zenk 2000), while terrestrial tropical cyclones over open water typically survive only $\sim$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, 2016, 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 $\bar {N}/f \gg 1$, where $\bar {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 ($\bar {N}/f \simeq 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 $\bar {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 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.

2. Coordinate systems, hydrostatic equilibrium and anelastic equations of motion

2.1. Useful coordinate systems

A coordinate system used commonly for studying planetary scale flows is spherical coordinates $(r,\theta,\phi )$, where $r$ is the radius, $\theta$ is the latitude, and $\phi$ is the longitude. However, planetary vortices are local structures in planetary atmospheres. The characteristic horizontal length ${\mathcal {L}}$ of a typical planetary vortex is small compared to the planet's radius. The largest planetary vortex, the GRS of Jupiter, spans ${\sim }11^\circ$ (longitude) by $7.4^\circ$ (latitude), or $13\,000\ \text {km} \times 9200$ km (Simon et al. 2018; Wong et al. 2021), so ${\mathcal {L}} \simeq 13\,000$ km. Because ${\mathcal {L}}$ is small compared to the mean Jovian radius $r_0 \simeq 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, \theta _0, \phi _0)$ be the centre of the vortex of interest. Then the local Cartesian coordinates in figure 1 are

(2.1)$$\begin{gather} x=[r_0 \cos(\theta_0)] (\phi - \phi_0) [1 + {O}({\mathcal{L}}/r_0 )^2], \end{gather}$$

(2.2)$$\begin{gather}y=r_0 (\theta-\theta_0) [1 + {O}({\mathcal{L}}/r_0 )^2], \end{gather}$$

(2.3)$$\begin{gather}z=r -r_0 [1 + {O}({\mathcal{L}}/r_0 )^2]. \end{gather}$$

Note that $\hat {\boldsymbol {z}}$ is a radially outward unit vector and is along the vertical direction in which the flow is stratified; $\hat {\boldsymbol {x}}$ is eastwards; and $\hat {\boldsymbol {y}}=\hat {\boldsymbol {z}}\times \hat {\boldsymbol {x}}$ is northwards. When computing the GRS at $\theta _0 \simeq -23^{\circ }$, we can ignore Jupiter's small oblateness and neglect the curvature correction terms because $({\mathcal {L}}/r_0)^2 \simeq 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 $(\partial ^2 /\partial x^2 + \partial ^2 /\partial 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 ${\mathcal {L}}$. For example, the GRS has a relatively quiet interior surrounded by a high-velocity collar of thickness $\sim$2000 km, so although ${\mathcal {L}} \simeq 13\,000$ km, we have $L \simeq 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\sim 2000$ km to try to simulate qualitatively the velocity field of the GRS at cloud-top level. However, $L_R \equiv H_P \bar {N}/f$, where $\bar {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 $\bar {N}$ varies by three orders of magnitude over the height of the GRS, so does $L_R$. See § 3.2.)

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,\theta _0,\phi _0$). The angle shown is $-\theta _0$ because the GRS is in the southern hemisphere. For the Cartesian coordinates $(x,y,z)$ in red, $\hat {\boldsymbol {x}}$ is eastwards, $\hat {\boldsymbol {z}}$ is radially outwards, and $\hat {\boldsymbol {y}}=\hat {\boldsymbol {z}}\times \hat {\boldsymbol {x}}$ is along the local north direction. For the Cartesian coordinates $(x_c,y_c,z_c)$ in green, $\hat {\boldsymbol {x}}_c$ is eastwards, $\hat {\boldsymbol {z}}_c$ is along the direction of the planetary rotational axis, and $\hat {\boldsymbol {y}}_c=\hat {\boldsymbol {z}}_c\times \hat {\boldsymbol {x}}_c$. This illustration is chosen with $\theta _0 < 0$, which would be appropriate for the GRS of Jupiter, which lies in the southern hemisphere.

Figure 2. Comparison of the vertical profiles of the background thermodynamic functions based on the observed $\bar {T}(\bar {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 $\bar {P}=1$ bar, near the top of the visible clouds.

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 $\hat {\boldsymbol {z}}_{{c}}$ is in the direction of the planet's rotation vector, $x_c \equiv x$, and $\hat{\boldsymbol{y}}_c$ is orthogonal to $\hat{\boldsymbol{x}}_c$ and $\hat{\boldsymbol{z}}_c$ such that the right-hand rule is satisfied. To convert between the two local Cartesian coordinate systems:

(2.4)$$\begin{gather} x=x_c, \end{gather}$$

(2.5)$$\begin{gather}y=y_c \sin \theta_0+z_c \cos \theta_0, \end{gather}$$

(2.6)$$\begin{gather}z=-y_c \cos \theta_0+z_c \sin \theta_0, \end{gather}$$

(2.7)$$\begin{gather}x_c=x, \end{gather}$$

(2.8)$$\begin{gather}y_c=y \sin \theta_0-z \cos \theta_0, \end{gather}$$

(2.9)$$\begin{gather}z_c=y \cos \theta_0+z \sin \theta_0. \end{gather}$$

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

(2.10)$$\begin{gather} x=x', \end{gather}$$

(2.11)$$\begin{gather}y=y' +(z-z') \cot \theta_0. \end{gather}$$

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.

2.2. Anelastic approximation and the governing equations

2.2.1. Anelastic approximation

We denote the solutions of the steady hydrostatic (i.e. the velocity ${\boldsymbol {v}} = 0$) equations with overbars. The Euler equation and ideal gas equation become

(2.12)$$\begin{gather} \frac{\mathrm{d} \bar{P}(z)}{\mathrm{d} z}=-\bar{\rho}(z)\,g, \end{gather}$$

(2.13)$$\begin{gather}\bar{P}(z)=R\,\bar{\rho}(z)\,\bar{T}(z), \end{gather}$$

where $P$ is pressure, $\rho$ is mass density, $T$ is temperature, $g$ is the acceleration of gravity pointing in the $-\hat {\boldsymbol {z}}$ 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 \times 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 $\varTheta$, with

(2.14)$$\begin{gather} \varTheta\equiv T\left(\frac{P_{ref}}{P}\right)^{R/C_p}, \end{gather}$$

(2.15)$$\begin{gather}\bar{\varTheta}(z)=\bar{T}(z)\left(\frac{P_{ref}}{\bar{P}(z)}\right)^{R/C_p}, \end{gather}$$

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 $\bar {T}(z)$. Using the observed $\bar {T}(z)$, all of the other barred thermodynamic quantities can computed with (2.12)–(2.15).

We decompose the thermodynamic variables as

(2.16)$$\begin{gather} \tilde{\rho}(x,y,z,t)=\rho(x,y,z,t)-\bar{\rho}(z), \end{gather}$$

(2.17)$$\begin{gather}\tilde{P}(x,y,z,t)=P(x,y,z,t)-\bar{P}(z), \end{gather}$$

(2.18)$$\begin{gather}\tilde{T}(x,y,z,t)=T(x,y,z,t)-\bar{T}(z), \end{gather}$$

(2.19)$$\begin{gather}\tilde{\varTheta}(x,y,z,t)=\varTheta(x,y,z,t)-\bar{\varTheta}(z). \end{gather}$$

The anelastic approximation requires that the non-hydrostatic thermodynamic variables (denoted with the tildes) are much smaller than their hydrostatic counterparts (denoted with the overbars) at each height $z$. In addition, the Mach number must be small. These requirements are consistent with Jovian observations.

2.2.2. 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. 2015, 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

(2.20)\begin{gather} 0=\boldsymbol{\nabla}\boldsymbol{\cdot}{(\bar{\rho} \boldsymbol{v})},\end{gather}

(2.21) \begin{gather} \begin{aligned}&\frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{v}-\frac{\tilde{\varTheta}}{\bar{\varTheta}}\,g\hat{\boldsymbol{z}}+\boldsymbol{\nabla}{\biggl(\frac{\tilde{P}}{\bar{\rho}}\biggr)} = \boldsymbol{v}\times\boldsymbol{f} \nonumber\\ &\quad\simeq f(v_y \sin \theta_0 - v_z \cos \theta_0) \hat{\boldsymbol{x}} \nonumber\\ &\qquad + \beta y (v_y + v_z \tan \theta_0) \hat{\boldsymbol{x}} \nonumber\\ &\qquad -v_x (\,f \sin \theta_0 + \beta y) \hat{\boldsymbol{y}} \nonumber\\ &\qquad +v_x (\,f \cos \theta_0 - \beta y \tan \theta_0) \hat{\boldsymbol{z}}, \end{aligned}\nonumber\\[-1.5pc] \end{gather}

(2.22) \begin{gather} \begin{aligned}\frac{\partial }{\partial t}\Bigg(\frac{\tilde{\varTheta}}{\bar{\varTheta}}\Bigg)&=-\frac{(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}){\tilde{\varTheta}}}{\bar{\varTheta}}-\frac{\bar{N}^2}{g}\,v_z \nonumber\\ &=-\frac{(\boldsymbol{v}_{{\perp}} \boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}}) {\tilde{\varTheta}}}{\bar{\varTheta}}-\frac{{N}^2}{g}\,v_z, \end{aligned}\nonumber\\[-1.5pc] \end{gather}

(2.23)\begin{gather} \frac{\tilde{\rho}}{\bar{\rho}}= \frac{C_v}{C_p} \, \frac{\tilde{P}}{\bar{P}} - \frac{\tilde{\varTheta}}{\bar{\varTheta}}, \end{gather}

(2.24)\begin{gather} \frac{\tilde{T}}{\bar{T}}=\frac{\tilde{P}}{\bar{P}}-\frac{\tilde{\rho}}{\bar{\rho}}, \end{gather}

where $N$ is the Brunt–Väisälä frequency, with

(2.25a,b)\begin{equation} {N}^2(x, y, z)\equiv\frac{g}{\bar{\varTheta}}\,\frac{\partial {\varTheta}}{\partial z} \quad \text{and} \quad \bar{N}^2(z)\equiv\frac{g}{\bar{\varTheta}}\,\frac{\mathrm{d} \bar{\varTheta}}{\mathrm{d} z} = \frac{g}{\bar{T}(z)}\left(\frac{\mathrm{d} \bar{T}(z)}{\mathrm{d} z}+\frac{g}{C_p}\right)\!, \end{equation}

and $\boldsymbol{f}=f \hat {\boldsymbol {z}}_c$ is the full Coriolis vector, which is along the direction of $\hat {\boldsymbol {z}}_c$, or $\boldsymbol {f}=$ $f \cos \theta \,\hat {\boldsymbol {y}}+f \sin \theta \, \hat {\boldsymbol {z}} \simeq (\,f \cos \theta _0 - \beta \tan \theta _0\,y) \hat {\boldsymbol {y}}+(\,f \sin \theta _0 + \beta y) \hat {\boldsymbol {z}}$, with $\beta \equiv (\,f \cos \theta _0)/r_0$. In writing the gravitational acceleration as $- g \hat {\boldsymbol {z}}$, 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

(2.26)\begin{align} \frac{\partial (\omega_z + \beta y)}{\partial t}&=-(\boldsymbol{v}_\perp\boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) (\omega_z + \beta y) -v_z\,\frac{\partial \omega_z}{\partial z}+({\boldsymbol{\omega}}_\perp\boldsymbol{\cdot}\boldsymbol{\nabla}_\perp)v_z \nonumber\\ &\quad -\omega_z (\boldsymbol{\nabla}_\perp\boldsymbol{\cdot}\boldsymbol{v}_\perp) -(\,f \sin\theta_0 + \beta y) (\boldsymbol{\nabla}_\perp\boldsymbol{\cdot}\boldsymbol{v}_\perp) \nonumber\\ &\quad +(\,f \cos \theta_0 - \beta y \tan \theta_0) \frac{\partial v_z}{\partial y} - \beta v_z \tan \theta_0, \end{align}

(2.27)\begin{align} {\frac{\partial (\boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} {\boldsymbol{v}}_{{\perp}})}{\partial t}} &=-\boldsymbol{\nabla}_\perp\boldsymbol{\cdot} \left[(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{v}_{{\perp}} + \boldsymbol{v}\times\boldsymbol{f} \right] - \nabla^2_{{\perp}} \biggl( \frac{\tilde{P}}{\bar{\rho}} \biggr)\!, \end{align}

(2.28)\begin{align} {\frac{\partial v_z}{\partial t}} &=- ({\boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{\nabla}) v_z + g\, \frac{\tilde{\varTheta}}{\bar{\varTheta}} -\frac{\partial }{\partial z}\biggl(\frac{\tilde{P}}{\bar{\rho}}\biggr) +f v_x \cos \theta, \end{align}

where ${\boldsymbol {\omega }}$ is the vorticity vector, $\boldsymbol {v}_\perp$ is the horizontal velocity, and ${\boldsymbol {\omega }}_\perp$ is the horizontal vorticity. Also,

(2.29a,b)\begin{equation} \boldsymbol{\nabla}_\perp \equiv \hat{\boldsymbol{x}}\,\frac{\partial}{\partial x}+\hat{\boldsymbol{y}}\,\frac{\partial}{\partial y},\quad \boldsymbol{\nabla}_\perp^2=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2} \end{equation}

is the horizontal Laplacian operator, and

(2.30)\begin{equation} \boldsymbol{\nabla}_\perp\boldsymbol{\cdot}\boldsymbol{V} \equiv \frac{\partial }{\partial x}(\boldsymbol{V}\boldsymbol{\cdot}\hat{\boldsymbol{x}})+\frac{\partial }{\partial y}(\boldsymbol{V}\boldsymbol{\cdot}\hat{\boldsymbol{y}}) \end{equation}

is the horizontal divergence of any given vector $\boldsymbol {V}$.

3. Choice of hydrostatic fields

3.1. Our non-use of the hydrostatic energy equation

The reason why we use the observed $\bar {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

(3.1)$$\begin{gather} {\frac{\partial T}{\partial t}} + ({\boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{\nabla}) T = \frac{\boldsymbol{\nabla} \boldsymbol{\cdot}( k\,\boldsymbol{\nabla} {T})}{C_v \rho} - \frac{P (\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{v}})}{C_v \rho} \nonumber\\ {}+ \text{source}/\text{sink terms due to radiative transfer and viscous heating}, \end{gather}$$

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 $\bar {N} \simeq 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), and Yadav & 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 $\bar {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, $\bar {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 $\bar {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.)

3.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 $\bar {T}(\bar {P})$ at the GRS latitude $\theta _0=23\,^\circ$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 $\bar {P}(z)$, $\bar {P}(\bar {T})$ and $\bar {P}(\bar {\rho })$ as shown by the blue dots in figure 2. Although the temperature data appear to be smooth, when we differentiate to compute $\bar {N}(z)$ with (2.25a,b), we find that it is not smooth, especially at heights deeper than $\sim 200$ mbar (see the blue dots in figure 2d). A non-smooth $\bar {N}(z)$ has the potential of producing artificial convective instabilities and internal gravity waves in our calculations. In addition, the measured values of temperature extend to a depth of only $\sim$1 bar. Therefore we have taken the observed temperature measurements, and change them in three ways. (1) At heights above $800$ mbar, we replace the observed $\bar {N}^2(z)$ (blue dots in figure 2d) with the locally smoothed red curve. (2) We define the top of the convective zone to be at 10 bar, so at heights deeper than 10 bar, we set $\bar {N} = 0^+$. (3) We construct a smooth monotonic curve (red) to extrapolate $\bar {N}$ between $800$ mbar and 10 bar. We then use (2.12), (2.13), (2.15) and (2.25a,b) with the smooth $\bar {N}$ in figure 2(d) to compute the smoothed (red) functions $\bar {P}(z)$, $\bar {P}(\bar {T})$ and $\bar {P}(\bar {\rho })$ in figures 2(a–c). The smoothed functions are not very different from the original observational data.

4. Choice of zonal flow

The Jovian east–west zonal flow $v_{zonal}(y, z)\,\hat {\boldsymbol {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 $\theta$ 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 $^{\circ }$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 $\tilde {T}$ or density $\tilde {\rho }$ data as a function of $(y,z)$. The temperature and density measurements are sufficiently uncertain that the results are ambiguous (Fletcher et al. 2016, 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 $\sim$200 km of the Jovian vortices studied here.

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} = - \sigma y$, where $\sigma = -9.1\times 10^{-6}\ \mathrm {s}^{-1}$. The velocity is shown as a function of planetographic latitude $\theta$ and of the north–south length $y$, with $y=0$ chosen to be the approximate centre of the GRS.

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 $\sim$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 ($\sim$140 mbar) is similar to the cloud-top measurement at the latitude of our interest ($23\,^\circ$S).

Any zonal flow ${\boldsymbol {v}} = v_{zonal}(y, z)\,\hat {\boldsymbol {x}}$ is a solution to the steady anelastic equations subject to the requirements that

(4.1)$$\begin{gather} {\frac{\tilde{P}(x, y, z)}{\bar{\rho}(z)}} =- \int_{0}^y {{\rm d} y}'\,v_{zonal}(y', z)\,(\,f \sin \theta_0 + \beta y') + {\frac{\tilde{P}(x, 0, z)}{\bar{\rho}(z)}} , \end{gather}$$

(4.2)$$\begin{gather}{\frac{\tilde{\varTheta}}{\bar{\varTheta}}}= {\frac{1}{g}} \biggl[ {\frac{\partial}{\partial z}}\biggl({\frac{\tilde{P}}{\bar{P}}} \biggr) - v_{zonal} (\,f \cos \theta_0 - \beta y \tan \theta_0)\biggr]\!. \end{gather}$$

Therefore, the only constraint on our choice of zonal flow is that

(4.3)\begin{equation} v_{zonal}(y,z) = {\mathcal{F}}(y_c) \equiv {\mathcal{F}}(y \sin \theta_0 - z \cos \theta_0), \end{equation}

where ${\mathcal {F}}(y \sin \theta _0 - z \cos \theta _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 2014, 2019, 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.

5. Initial condition set-up

5.1. Strategy

The unforced and undissipated equations that govern the flow do not have a unique solution. For example, $\tilde {P} = \tilde {\rho } = \tilde {\varTheta } = \tilde {T} = {\boldsymbol {v}} =0$ is a solution with no zones and no vortices. Purely zonal flows $\boldsymbol {v}_{zonal}(y,z) \equiv v_{zonal}(y, z)\, \hat {\boldsymbol {x}}$ with no vortices also exist with the zonal vorticity equal to ${\boldsymbol {\omega }}_{zonal} \equiv \boldsymbol {\nabla }\times \boldsymbol {v}_{zonal}$, and with zonal pressure $\tilde {P}_{zonal}$, zonal temperature $\tilde {T}_{zonal}$ and zonal potential temperature $\tilde {\varTheta }_{zonal}$ determined by (4.1) and (4.2) with $\boldsymbol {v}=\boldsymbol {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 $\boldsymbol {v}_{vortex}\equiv \boldsymbol {v}-\boldsymbol {v}_{zonal}$, with similar decompositions for pressure, temperature and potential temperature. For the vertical vorticity, the contribution due to the vortex is $\omega _{z,vortex} \equiv ({\boldsymbol {\omega }} - {\boldsymbol {\omega }}_{zonal}) \boldsymbol {\cdot } \hat {\boldsymbol {z}}$, where ${\boldsymbol {\omega }}_{zonal} \boldsymbol {\cdot } \hat {\boldsymbol {z}} \equiv - \partial v_{zonal}/\partial y$. We also define $\omega _{z,vortex,0}(x,y) \equiv \omega _{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.)

5.2. 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 $\omega _{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 (Marcus, Kundu & Lee 2000; Shetty & Marcus 2010; Marcus & Shetty 2011) or shallow-water equation (Dowling & Ingersoll 1989; Cho & Polvani 1996; Stegner & Dritschel 2000; García-Melendo & Sánchez-Lavega 2017; Li et al. 2020) to look at vortex dynamics. To create near equilibria, we assume $v_z=0$ and that the flow is steady in time. With these assumptions, (2.20), (2.26)–(2.28) and (2.22) become

(5.1)$$\begin{gather} \boldsymbol{\nabla}_{{\perp}} \boldsymbol{\cdot} {\boldsymbol{v}}_{{\perp}} = 0, \end{gather}$$

(5.2)$$\begin{gather}(\boldsymbol{v}_\perp\boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) (\omega_z + \beta y) = 0, \end{gather}$$

(5.3)$$\begin{gather}\nabla^2_{{\perp}} \biggl(\frac{\tilde{P}}{\bar{\rho}} \biggr) = \boldsymbol{\nabla}_\perp\boldsymbol{\cdot} [-(\boldsymbol{v}_{{\perp}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{{\perp}})\boldsymbol{v}_{{\perp}}] + (\,f \sin \theta_0 + \beta y) \omega_z - \beta v_x , \end{gather}$$

(5.4)$$\begin{gather}g\,\frac{\tilde{\varTheta}}{\bar{\varTheta}} = \frac{\partial }{\partial z}\,\biggl(\frac{\tilde{P}}{\bar{\rho}}\biggr) -v_x\, (\,f \cos \theta_0 - \beta y \tan \theta_0), \end{gather}$$

(5.5)$$\begin{gather}(\boldsymbol{v}_\perp\boldsymbol{\cdot}\boldsymbol{\nabla}_\perp) \tilde{\varTheta} = 0. \end{gather}$$

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 ${\boldsymbol {v}}_{\perp }(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)\,\hat {\boldsymbol {x}}$ is given, then a multi-parameter family of steady or slowly drifting vortices, with $\omega _{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 ${\boldsymbol {v}}_{\perp }(x, y, z)$, we substitute it into the right-hand side of (5.3) and solve for $\tilde {P}(x, y, z)$. We then use (5.4) with $\tilde {P}(x, y, z)$ to find $\tilde {\varTheta }(x, y, z)$, and use $\tilde {\varTheta }(x, y, z)$ in (2.23)–(2.24) to find $\tilde {\rho }(x, y, z)$ and $\tilde {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 $\beta \equiv 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 $\tilde \omega _z=\omega _0$ inside the vortex. For $\boldsymbol {v}_{zonal}=-\sigma y \hat {\boldsymbol {x}}$, they found

(5.6)\begin{equation} \sigma/\omega_0 = \frac{\epsilon^2-\epsilon}{\epsilon+1}, \end{equation}

where $\epsilon$ is the horizontal aspect ratio (length of of the ellipse's major diameter in $x$ to its minor diameter in $y$). As $\sigma$ increases, the zonal flow shear stretches the vortex in the $x$ direction, and $\epsilon$ 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 $\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla }\omega _z=0$. In the limit of $\omega _0\rightarrow 0$, we have $\epsilon \rightarrow \infty$, 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 $\sigma$ and $\omega _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 $\beta =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 $\beta \equiv 0$, and even though their vertical vorticities $\omega _{z,vortex}$ are uniform and do not look like the GRS, which has a pronounced local minimum of $\omega _{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 $\omega _{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 $\sim$2000 km. To create our Moore–Saffman-like vortices at the cloud-top level, we define the following elliptical coordinates for convenience:

(5.7)$$\begin{gather} \eta=\sqrt{\left(\frac{x-x_0}{\epsilon}\right)^2+(y-y_0)^2}, \end{gather}$$

(5.8)$$\begin{gather}\chi=\tan^{-1}\left[\frac{\epsilon(y-y_0)}{x-x_0}\right]\!, \end{gather}$$

where $(x_0, y_0)$ is the centre of the vortex. We define the velocity field of the vortex as

(5.9)$$\begin{gather} v_{vortex}(x,y)= \begin{cases} \dfrac{\epsilon}{1+\epsilon}\,\omega_0 \eta & (\eta< R_v) ,\\ \dfrac{\epsilon}{1+\epsilon}\,\omega_0\,\dfrac{R_v^2} {\eta} \exp\left({\dfrac{-\eta+R_v}{L_r}}\right) & (\eta > R_v), \end{cases} \end{gather}$$

(5.10)$$\begin{gather}\boldsymbol{v}_{vortex}(x,y)=-v_{vortex} (\sin \chi)\,\hat{\boldsymbol{x}} + v_{vortex} (\cos \chi)\, \hat{\boldsymbol{y}}, \end{gather}$$

(5.11)$$\begin{gather}\omega_{z,vortex}(x,y)=\frac{\partial }{\partial x}v_{y,vortex}-\frac{\partial }{\partial y}v_{x,vortex}, \end{gather}$$

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\ \text {km}$ throughout this study. The value of $\omega _0$ is determined by $\epsilon$ via (5.6), and $\omega _{z,vortex}(x_0,y_0)=\omega _0$. In the limit $L_r \rightarrow \infty$, the vortex given by (5.9)–(5.11) is an exact Moore–Saffman vortex with vorticity $\omega _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.

5.3. 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.

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 $\theta _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 $\sim$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.

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 $\boldsymbol {v}_{vortex,0}(x, y)$ and vertical vorticity ${\omega }_{z,vortex,0}(x, y)$. In future studies, we will set ${\omega }_{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 $\omega _{z,vortex,0}(x,y)$ is given. We are particularly interested in whether the vortex's final quasi-equilibrium $\boldsymbol {v}_{vortex}(x, y, z=0)$ is approximately the same as its initial reference $\boldsymbol {v}_{vortex,0}(x, y)$. In using (5.9)–(5.11), the value of $\omega _0$ is slaved to the value of $\epsilon$ by (5.6) and the value of $\sigma$. In this paper, we use $\sigma = -9.1\times 10^{-6}\ \mathrm {s}^{-1}$ when calculating $\omega _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 \ne 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 $\omega _{z,vortex,0}(x,y) \equiv \omega _{z,vortex}(x,y, z=0)$ to other heights $z$. (The two-dimensional velocity at any $z$ is determined uniquely by $\omega _{z,vortex}(x,y, z)$ because $\boldsymbol {\nabla }_{\perp } \boldsymbol {\cdot } \boldsymbol {v}_{vortex} =0$ at all $z$.) For example, we could create a columnar stacked vortex aligned with the $z$ axis by defining $\omega _{z,vortex}(x,y,z)=\omega _{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 $\mathrm {CA}\unicode{x2013}\mathrm {Rad}$ family):

(5.12)\begin{equation} \mathrm{CA}\unicode{x2013}\mathrm{Rad}\ \mathrm{family} \quad \omega_{z,vortex}(x,y,z)=\omega_{z,vortex,0}(x,y)\,\frac{F(z)}{F(0)},\end{equation}

with

(5.13)\begin{equation} F(z)=\begin{cases} {\rm e}^{-[(z-z_0)/H_{top}]^2} & (z>{z_0}),\\ {\rm e}^{-[(z-z_0)/H_{bot}]^2} & (z<{z_0}), \end{cases} \end{equation}

which projects $\omega _{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 $\omega _{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 \theta _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 $\mathrm {CA}\unicode{x2013}\mathrm {Rad}$ because it projects $\omega _{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 $\omega _{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 $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ family, we set

(5.14)\begin{equation} \mathrm{CA}\unicode{x2013}\mathrm{Rot} \, \mathrm{family} \quad \omega_{z,vortex}(x,y,z)=\omega_{z,vortex,0}(x, y-z \cot \theta_0)\, \frac{F(z)}{F(0)}, \end{equation}

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 $\omega _{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,

(5.15)$$\begin{gather} \mathrm{CV}\unicode{x2013}\mathrm{Rad} \, \mathrm{Family} \quad \omega_{z,vortex}(x,y,z)=\omega_{z,vortex,0}(x', y'), \end{gather}$$

(5.16)$$\begin{gather}x'=\frac{F(0)}{F(z)}\,x, \end{gather}$$

(5.17)$$\begin{gather}y'=\frac{F(0)}{F(z)}\,y \end{gather}$$

produces a family of vortices, $\mathrm {CV}\unicode{x2013}\mathrm {Rad}$, oriented as in figure 4(a), where the magnitude of $\omega _{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, $\mathrm {CV}\textit {--}\mathrm {Rot}$, but oriented as in figure 4(b), can be created with

(5.18)$$\begin{gather} \mathrm{CV}\textit{--}\mathrm{Rot} \, \mathrm{family} \quad \omega_{z,vortex}(x,y,z)=\omega_{z,vortex,0}\ (x', y'), \end{gather}$$

(5.19)$$\begin{gather}x'=\frac{F(0)}{F(z)}\,x, \end{gather}$$

(5.20)$$\begin{gather}y'=\frac{F(0)}{F(z)}\,(y-z \cot \theta_0). \end{gather}$$

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 $\mathrm {CV}\textit {--}\mathrm {Rad}$ and $\mathrm {CV}\textit {--}\mathrm {Rot}$ families would be exact solutions to (5.1)–(5.3) with $\beta =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 $\mathrm {CA}\unicode{x2013}\mathrm {Rad}$ and $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ families are far from equilibrium and not approximate solutions to (5.1)–(5.2). The shear $\sigma$ of the zonal flow at the tops and bottoms of these vortices is the same as it is at $z=0$, but $\omega _{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 $\mathrm {CA}\unicode{x2013}\mathrm {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 $\omega _{z,vortex}(x, y=0, z)$ for vortices from the $\mathrm {CV}\textit {--}\mathrm {Rot}$ and $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ families are shown in figure 5.

Figure 5. Examples of $\omega _{z,vortex}(x,y=0,z)/(\,f\sin \theta _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 \theta _0$. Both vortices have $H_{top}=50$ km, $H_{bot} =120$ km, $\omega _0 =3.0\times 10^{-5}$ s$^{-1}$, $\epsilon = 1.5$, $R_v =4600$ km and $L_r =2000$ km. (a) A constant-vorticity, area-varying $\mathrm {CV}\textit {--}\mathrm {Rot}$ vortex. (b) A constant-area, vorticity-varying $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ vortex. Note that the planetary vortices are pancake vortices. For example, our vertical domain size is $\sim$300 km, while the horizontal domain size is $\sim$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.

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 $\tilde {P}_{vortex}/\bar {\rho }$ and $|\boldsymbol {v}_{vortex}|$. It is also the $x$–$y$ plane where the initial $\mathrm {CA}$ vortices have their maximum values of $|\tilde {\omega }_z|$, and the $x$–$y$ plane where the initial $\mathrm {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 $\tilde {P}_{vortex}/\bar {\rho }$, 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 $\mathrm {CA}$ family will have its biggest vertical vorticity, and where a vortex in the $\mathrm {CV}$ family will have its biggest horizontal area, so it is useful to track the flow fields in the mid-plane.

6. Numerical results with $\beta \equiv 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\ {\rm km}\le H_{bot} \le 180\ {\rm km}$, $50\ {\rm km}\le H_{top} \le 52.5\ {\rm km}$ and $0\le z_0 \le 15$ km (equivalently, between 500 mbar and 1 bar), and with parameter values for $\omega _{z,vortex,0}(x, y)$ (in (5.9)–(5.11)) in the range $-0.3 \le Ro \le -0.1$, or equivalently, $1.2 \le \epsilon \le 1.5$. We chose $R_v = 4600$ km and $L_r = 2000$ km. Here, the initial Rossby number is defined as $Ro \equiv \omega _0/(2f \sin \theta _0)$. Note that $\omega _0$ is slaved to the value of $\epsilon$ by (5.6) and the value of $\sigma$. In this paper, we use $\sigma = -9.1\times 10^{-6}\ \mathrm {s}^{-1}$ when calculating $\omega _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 $\theta _0 = 23\,^\circ$S and $f = 3.5\times 10^{-4}\ \mathrm {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.

6.1. Orientation: aligned with the $z$ or $z_c$ axis?

The planetary vortices of interest in this study have low Rossby number ($Ro\lesssim 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 $\mathrm {CV}\textit {--}\mathrm {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

(6.1)\begin{equation} y_{ori}(z)=\frac{\iint[\tilde{P}_{vortex}(x,y,z)/\bar{\rho}(z)] \,y\,{\rm d}\kern0.06em x\,{\rm d}y}{\iint[\tilde{P}_{vortex}(x,y,z)/\bar{\rho}(z)] \, {\rm d}\kern0.06em x\,{\rm d}y}. \end{equation}

We use $\tilde {P}_{vortex}/\bar {\rho }$ as the weighting factor because it is almost never negative and is a good proxy for the stream function of the horizontal velocity. If ${\rm d}z/{\rm d}y_{ori} = \pm \infty$, then the central axis is aligned with the $z$ axis, and (using (2.11)) if ${\rm d}z/{\rm d}y_{ori} = \tan \theta _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 ${\rm d}z/{\rm d}y_{ori} \simeq \tan \theta _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 (${\sim }0.01$), the change in the location of the central axis during the simulation is extremely small.

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 \theta _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} \equiv 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 $\mathrm {CV}\textit {--}\mathrm {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 $\omega _{z, vortex, 0}(x,y)$, we chose $Ro = -0.1$, or equivalently, $\epsilon = 1.5$.

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\ {\rm km} < z < 0\ {\rm 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 $\mathrm {CV}\textit {--}\mathrm {Rad}$ and $\mathrm {CA}\unicode{x2013}\mathrm {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 $\mathrm {CV}\textit {--}\mathrm {Rot}$ and $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ families.

6.2. $\mathrm {CV}\textit {--}\mathrm {Rot}$

Using (5.18), the maximum vertical vorticity $\omega _{z,vortex}$ of the initial vortices is the same at each height, but the area of the vortex changes. We examined several initial $\mathrm {CV}\textit {--}\mathrm {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 $\mathrm {CV}\textit {--}\mathrm {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, $\epsilon = 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 $\bar {P} \simeq 1$ bar). One reason why we chose this value $z_0$ was so that there are no initial, local convective instabilities (see § 8).

6.2.1. Vortex embedded within a zonal velocity with constant shear (case CV1)

Figure 7 shows $\omega _{z,vortex}/|\,f\sin \theta _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 $\omega _{z,vortex}/|\,f\sin \theta _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 $\omega _{z,vortex}$ at the central rotation axis of the vortex is visible at most values of $z$.

Figure 7. Plots of $\omega _{z,vortex}(x,y)/|\,f\sin \theta _0|$ at the mid-plane of the vortex (case CV1). See supplementary movie 1 available at https://doi.org/10.1017/jfm.2024.132 for more details. At all times, the vortex mid-plane is defined as the height $z$ where the anticyclone has its maximum $\tilde {P}_{vortex}/\bar {\rho }$. The time is in units of Earth days: (a) $t=0$, (b) $t\sim 40$, (c) $t\sim 80$, (d) $t\sim 120$. The colour bar shows the value of $\omega _{z,vortex}(x,y)/|\,f\sin \theta _0|$. The turn-around time of this vortex is $\sim$4 days; that of the GRS is $\sim$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 $\mathrm {CV}\textit {--}\mathrm {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 $\omega _{z,vortex}$ at its centre. The vortex remains shielded, but most of the negative $\omega _{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 $\mathrm {CV}\textit {--}\mathrm {Rot}$ vortices embedded in constant zonal shear. There is little qualitative change of the vortex after 120 days.

Figure 8. Plots of $\omega _{z,vortex}(x,y=0,z)/|\,f\sin \theta _0|$ in the $x$–$z$ plane of the vortex (case CV1). See supplementary movie 2 available at https://doi.org/10.1017/jfm.2024.132 for more details. The colour bar is as in figure 7. In this plane at day 40, the shed filaments have mostly negative values of $\omega _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 $\omega _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 \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

The truly surprising result is that the vortex becomes hollow. The initial $\omega _{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), $\omega _{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 $\tilde {\varTheta }_{vortex}/\bar {\varTheta }$ 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:

(6.2)\begin{equation} \frac{\partial }{\partial z}\biggl(\frac{\tilde{P}_{vortex}}{\bar{\rho}}\biggr)\sim \frac{\tilde{\varTheta}_{vortex}}{\bar{\varTheta}}\,g. \end{equation}

Figure 9. Plots of $\tilde {\varTheta }_{vortex}(x,y=0,z)/\bar {\varTheta }$ (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 $\tilde {\varTheta }_{vortex}(x,y,z) \simeq 0$. The high-$\tilde {\varTheta }_{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 $\tilde {\varTheta }_{vortex}/\bar {\varTheta }$ 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 \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

An anticyclone has a high pressure centre, so $\tilde {\varTheta }_{vortex}<0$ at heights above the mid-plane, and $\tilde {\varTheta }_{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. 2010, 2016).

Figure 10 shows $\tilde {\varTheta }_{vortex}(x,y)/\bar {\varTheta }$ 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 $\sim$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 $\tilde {\varTheta }_{vortex}/\bar {\varTheta }$ at the outer edge of the vortex, approximately coincident with the shield of negative $\omega _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 $\tilde {P}_{vortex}/\bar {P}$ at the location of this annulus (and because $\tilde {T}_{vortex}/\bar {T}=\tilde {\varTheta }_{vortex}/\bar {\varTheta }+(R/C_p) \tilde {P}_{vortex}/\bar {P}$), there is a high-temperature annulus of $\tilde {T}_{vortex}/\bar {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 $\tilde {\varTheta }_{vortex}$ was created by a downward vertical velocity that formed in the annulus. That is, the high-temperature annulus was created by adiabatic heating. de Pater et al. (2010) observed $5\ \mathrm {\mu }{\rm m}$ bright annular rings around the GRS and other Jovian anticyclones, and hypothesized: (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\ \mathrm {\mu }{\rm 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 10. Plots of $\tilde {\varTheta }_{vortex}(x,y)/\bar {\varTheta }$ 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-$\tilde {\varTheta }_{vortex}$ annular ring at the outer edge of the vortex and the ambient flow is approximately 1 K. The high-$\tilde {\varTheta }_{vortex}$ annulus in (b–d) is not present in (a) at $t=0$. The high-$\tilde {\varTheta }_{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 \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

Figure 11 shows the maximum values of $|\boldsymbol {v}_{vortex}|$, $\omega _{z,vortex}$, $\tilde {P}_{vortex}/\bar {\rho }$ 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 $(\bar {N}^2 - N^2)$ as a function of $z$ along the central axis of the vortex. Figures 11(a,c) show that the vertical profiles of $|\boldsymbol {v}_{vortex}|$ and $\tilde {P}_{vortex}/\bar {\rho }$ 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 $\mathrm {CV}\textit {--}\mathrm {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 $\bar {N}^2 - N^2 = - (g/\bar {\varTheta })(\partial \tilde {\varTheta }/\partial 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 \theta _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 $\tilde {N}^2_{vortex}/g < 0$, mixing has super-stratified the fluid, making it more stable to convection.

Figure 11. Vertical profiles of the vortex in case CV1. (a) Maximum $|\boldsymbol {v}_{vortex}|$ at each $x$–$y$ plane as a function of $z$. (b) Maximum $\omega _{z,vortex}$ at each $x$–$y$ plane as a function of $z$. (c) Maximum $\tilde {P}_{vortex}/\bar {\rho }$ 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$: $\|\boldsymbol {v}_{z,vortex}\|_{\infty }/\|\boldsymbol {v}_{\perp,vortex}\|_{\infty }\times 100$. (e) Plot of $(\bar {N}^2-N^2)/(\,f\sin \theta _0)^2$ along the $z$ axis of the vortex. Plots (a–c) are normalized by their maximum values; (d) is normalized by the maximum $|\boldsymbol {v}_{vortex}|$ of the vortex. The green dashed line in (e) indicates where $\bar {N}^2-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.

6.2.2. 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 $\omega _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.).

Figure 12. Plots of $\omega _{z,vortex}(x,y)/|\,f\sin \theta _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 \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

Figure 13. Plots of $\omega _{z,vortex}(x,y=0,z)/|\,f\sin \theta _0|$ of the CV2 vortex plotted as we did for the CV1 vortex in figure 8. Here, (a) $t = 0$ days, (b) $t \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

Figure 14. Plots of $\tilde {\varTheta }(x,y)/\bar {\varTheta }$ of the CV2 vortex plotted as we did for the CV1 vortex in figure 9. The black arrows indicates the locations of the westward/eastward zonal flow's maximum location. Here, (a) $t = 0$ days, (b) $t \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

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 $\omega _{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 $\omega _{z,vortex} >0$; the shields are the regions just outside interiors where $\omega _{z,vortex} <0$. The $\omega _{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 $\omega _{z,vortex}$ at $y=0$ surrounded in the north and south by large ‘shoulders’ of positive $\omega _{z,vortex}$. A measure of the hollowness of a vortex is the ratio of the largest value of $\omega _{z,vortex}$ in the shoulders divided by the local minimum value of $\omega _{z,vortex}$ at the vortex centre. Both the GRS and CV2 have similar ratios, $\sim$2. However, the GRS is larger than CV2, and its region of hollowness is also larger.

Figure 15. Comparison of the $\omega _{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 $\omega _{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 $\sim$2. The main difference between the GRS and CV2 is that the GRS is larger.

6.3. $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$

Using (5.14) to create an initial $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ vortex makes the horizontal area of the vertical vorticity $\omega _{z,vortex}$ of the initial vortex the same at each height, but the magnitude of $\omega _{z,vortex}$ decreases away from the mid-plane, vanishing at the tops and bottoms of the vortex. None of our simulations that began with $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ vortices produced a quasi-steady $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ 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 $\omega _{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 $\mathrm {CV}\textit {--}\mathrm {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, $\epsilon = 1.5$ (equivalently, $Ro =-0.1$) and $z_0 = 8.2$ km (equivalently, $\bar {P}=700$ mbar). Initially, the $\mathrm {CA}\unicode{x2013}\mathrm {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 $\mathrm {CV}\textit {--}\mathrm {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 $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ vortices are far from equilibrium and that they quickly evolve into $\mathrm {CV}\textit {--}\mathrm {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.

Figure 16. Evolution of $\omega _{z,vortex}(x,y)/|\,f\sin \theta _0|$ in the mid-plane of an initial $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ vortex embedded in a constant zonal shear (case CA1). The surviving vortices are all part of the $\mathrm {CV}\textit {--}\mathrm {Rot}$ family (see figure 17). Here, (a) $t = 0$ days, (b) $t \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

Figure 17. Evolution of $\omega _{z,vortex}(x,y=0,z)/|\,f\sin \theta _0|$ of the CA1 vortex in figure 16 shown in the $x$–$z$ plane passing through the central axis of the initial vortex. See supplementary movie 3 available at https://doi.org/10.1017/jfm.2024.132 for more details. The horizontal area of the initial vortex in (a) is approximately constant in $z$, but the magnitude of $\omega _{z,vortex}$ decreases with distance from the mid-plane, making it part of the $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ family. By day 40, the initial vortex and the large vortices spawned from it have changed from $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ vortices to $\mathrm {CV}\textit {--}\mathrm {Rot}$ vortices. Here, (a) $t = 0$ days, (b) $t \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

Figure 18. Plots of $\omega _{z,vortex}(x,y)/|\,f\sin \theta _0|$ in the mid-plane of the same initial vortex shown in figure 16, but embedded in the Jovian zonal flow. We call this the CA2 vortex. At 120 days, there are two large $\mathrm {CV}\textit {--}\mathrm {Rot}$ vortices and one smaller one. Based on simulations of two-dimensional quasi-geostrophic vortices, we believe that the difference in the number of vortices that survive for 120 days in this figure and in figure 16 is not physically significant. The initial filamentation of the out-of-equilibrium initial vortices is sensitive to small details of the initial flow, and when the filaments roll up into vortices, their latitudes are also sensitive to the details of the initial condition. If many of the rolled-up vortices have latitudes far from the latitude of the surviving vortex, then few of them will merge together. Here, (a) $t = 0$ days, (b) $t \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

Figure 19. Plots of $\omega _{z,vortex}(x,y=0,z)/|\,f\sin \theta _0|$ in the $x$–$z$ plane passing through the central axis of the initial vortex (case CA2). Although the initial vortex in (a) is in the $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ family, the surviving vortices in (d) are in the $\mathrm {CV}\textit {--}\mathrm {Rot}$ family. Here, (a) $t = 0$ days, (b) $t \sim 40$ days, (c) $t \sim 80$ days, (d) $t \sim 120$ days.

6.4. 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 \equiv 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 ${\boldsymbol {v}}$, $\tilde {\rho }$ and $\tilde {P}$ to change were long compared to the turn-around time of the vortex, but the initial characteristic time for $\tilde {\varTheta }$ to change was much faster. Of the two terms on the far right-hand side of (2.22), only $({\boldsymbol {v}}_{\perp } \boldsymbol {\cdot } \boldsymbol {\nabla }_{\perp }) \tilde {\varTheta }$ can initially change $\partial \tilde {\varTheta }/\partial t$ because $v_z$ is initially zero. When the vortex reaches its final steady state, the magnitude of the $({\boldsymbol {v}}_{\perp } \boldsymbol {\cdot } \boldsymbol {\nabla }_{\perp }) \tilde {\varTheta }$ term has decreased approximately by a factor of 3 from its initial magnitude, and is balanced primarily by the $\bar {N}^2 v_z \bar {\varTheta }/g$ term. Thus even though $v_z$ is small compared to $v_{\perp }$ (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 $\tilde {\varTheta }$ coincide with the two-dimensional horizontal streamlines.)

7. Scaling analyses and the shape of vortices

7.1. Scaling analysis of the vertical vorticity equation

In § 5.2, we showed that a steady vortex with $v_z \equiv 0$ has $(\boldsymbol {v}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla }_\perp )(\omega _z + \beta y) =0$. However, none of our computed vortices has $v_z \equiv 0$, and it is important to understand how small $(\boldsymbol {v}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla }_\perp )(\omega _z + \beta y)$ is, how it scales with the small dimensionless parameters of the flow, such as the Rossby number $Ro \equiv \|\omega _z\|_{\infty }/(2 f |\sin \theta _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 $\langle v_z \rangle / \langle v_{\perp } \rangle$ 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 \equiv \|v_{\perp }\|_2/\|\omega _z\|_2$ evaluated at the vortex mid-plane. We also define the characteristic vertical size of a vortex $H \equiv \|v_{\perp }\|_2/\|\omega _{\perp }\|_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 $\bar {\rho }$ and $\bar {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:

(7.1)$$\begin{gather} \langle v_z \rangle = Ro^2\,\frac{{H}}{L}\,\langle v_{{\perp}} \rangle, \end{gather}$$

(7.2)$$\begin{gather}\left\langle (\boldsymbol{\nabla}_\perp\boldsymbol{\cdot}\boldsymbol{v}_\perp) \right\rangle = \langle v_z \rangle/H = Ro^2 \, \langle v_{{\perp}} \rangle/L, \end{gather}$$

(7.3)$$\begin{gather}\left\langle (\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{v}) \right\rangle = \langle v_z \rangle/H = Ro^2\,\langle v_{{\perp}} \rangle/L, \end{gather}$$

(7.4)$$\begin{gather}\left\langle (\boldsymbol{v}_\perp\boldsymbol{\cdot} \boldsymbol{\nabla}_\perp)(\omega_z + \beta y)\right\rangle = Ro \left[{\frac{\langle v_\perp \rangle}{L}}\right]^2, \end{gather}$$

(7.5)$$\begin{gather}\left\langle\, f \, \sin\theta_0+ \beta y) (\boldsymbol{\nabla}_\perp\boldsymbol{\cdot}\boldsymbol{v}_\perp) \right\rangle = Ro \left[{\frac{\langle v_\perp \rangle}{L}}\right]^2, \end{gather}$$

(7.6)$$\begin{gather}\left\langle \omega_z (\boldsymbol{\nabla}_\perp\boldsymbol{\cdot}\boldsymbol{v}_\perp) \right\rangle = Ro^2 \left[{\frac{\langle v_\perp \rangle}{L}}\right]^2, \end{gather}$$

(7.7)$$\begin{gather}\left\langle v_z\,\frac{\partial \omega_z}{\partial z}\right\rangle = Ro^2 \left[{\frac{\langle v_\perp \rangle}{L}}\right]^2, \end{gather}$$

(7.8)$$\begin{gather}\left\langle ({\boldsymbol{\omega}}_\perp\boldsymbol{\cdot}\boldsymbol{\nabla}_\perp)v_z\right\rangle = Ro^2 \left[{\frac{\langle v_\perp \rangle}{L}}\right]^2, \end{gather}$$

(7.9)$$\begin{gather}\left\langle (\,f \cos_0 - \beta y \tan \theta_0)\,\frac{\partial v_z}{\partial y}\right\rangle = Ro\cot \theta_0\, {\frac{H}{L}} \left[{\frac{\langle v_\perp \rangle}{L}}\right]^2, \end{gather}$$

(7.10)$$\begin{gather}\left\langle \beta v_z \tan \theta_0 \right\rangle = Ro\,{\frac{H}{r_0}} \left[{\frac{\langle v_\perp \rangle}{L}}\right]^2. \end{gather}$$

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 $(\boldsymbol {v}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla }_\perp )(\omega _z + \beta y)$ is not identically zero as in (5.2), it is small compared to $[{{\langle v_\perp \rangle }/{L}}]^2$. Thus the two-dimensional streamlines and the isocontours of $(\omega _z + \beta 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 $(\boldsymbol {v}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla }_\perp )(\omega _z + \beta y)$ and $f \sin \theta _0\, (\boldsymbol {\nabla }_\perp \boldsymbol {\cdot }\boldsymbol {v}_\perp )$. This balance is illustrated in figure 20, which shows the values of $-(\boldsymbol {v}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla }_\perp )(\omega _z + \beta y)$ (figure 20a) and $-f \sin \theta \,(\boldsymbol {\nabla }_\perp \boldsymbol {\cdot }\boldsymbol {v}_\perp )$ (figure 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\sim 0.1$. For the simulations presented here, $H/L\sim 0.01$. The value of $H/L$ for Jovian vortices is not well constrained by observations. Also note that $\langle L \rangle /r_0 \simeq 0.03$ for the largest Jovian vortices.

Figure 20. The dominant two components of (2.26) at the vortex mid-plane (using the final state of case CV1). The plotted values are normalized by ${\|\boldsymbol {v}\|_\infty ^2}/{L^2}$: (a) $-(\boldsymbol {v}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla }_\perp ) \omega _z$, (b) $-f \sin \theta _0\,(\boldsymbol {\nabla }_\perp \boldsymbol {\cdot }\boldsymbol {v}_\perp )$.

7.2. Shape of the vortex

The shapes of $\omega _{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

(7.11)\begin{equation} N^2\equiv\frac{g}{\bar{\varTheta}}\,\frac{\partial \varTheta}{\partial z}=\bar{N}^2+\frac{g}{\bar{\varTheta}}\,\frac{\partial \tilde{\varTheta}}{\partial z}, \end{equation}

and using the hydrostatic approximation (6.2),

(7.12)\begin{align} N^2&=\bar{N}^2+\frac{g}{\bar{\varTheta}}\,\frac{\partial }{\partial z}\biggl[\frac{\bar{\varTheta}}{g}\,\frac{\partial }{\partial z}\biggl(\frac{\tilde{P}}{\bar{\rho}}\biggr)\biggr] \end{align}

(7.13)\begin{align} &=\bar{N}^2+\frac{\mathrm{d} \ln\bar{\varTheta}}{\mathrm{d} z}\,\frac{\partial }{\partial z}\biggl(\frac{\tilde{P}}{\bar{\rho}}\biggr)+\frac{\partial^2}{\partial z^2}\biggl(\frac{\tilde{P}}{\bar{\rho}}\biggr)\!. \end{align}

To eliminate $({\tilde {P}}/{\bar {\rho }})$ in (7.13), we use the small-Rossby number limit (quasi-geostrophic approximation) of (5.3) with $\beta =0$:

(7.14)\begin{equation} (\,f \sin \theta_0)\,\omega_z(x, y, z) = \nabla_\perp^2 \biggl(\frac{\tilde{P}(x, y, z)}{\bar{\rho}}\biggr)\!. \end{equation}

For each $z$ within the vortex, averaging (7.14) over the horizontal area within the vortex, we define ${\mathcal {A}}(z)$ with

(7.15) \begin{equation} (\,f \sin \theta_0)\,\omega_z(z)\biggl|_{averaged} = \biggl[\nabla_\perp^2 \biggl(\frac{\tilde{P}}{\bar{\rho}}\biggr)\biggr]\biggl|_{averaged} \equiv- \biggr[\biggl(\frac{\tilde{P}}{\bar{\rho}}\biggr)\biggl|_{averaged}\biggr] \biggl/ {\mathcal{A}}(z), \end{equation}

where ${\mathcal {A}}(z)$ has dimensions of length squared. Our numerical calculations show that $\sqrt {{\mathcal {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 $\mathrm {CV}\textit {--}\mathrm {Rot}$ vortices, it is reasonable to approximate $\omega _z(z)|_{averaged}$ as independent of $z$, and assign it a value $\omega _A$. Thus (7.13)–(7.15) become

(7.16)\begin{equation} \frac{{\rm d}^2{\mathcal{A}}}{{\rm d}z^2}+\frac{\mathrm{d} \ln\bar{\varTheta}}{\mathrm{d} z}\,\frac{\mathrm{d} {\mathcal{A}}}{\mathrm{d} z}=-\frac{N^2(z)|_{averaged}-\bar{N}^2(z)}{f \omega_A \sin \theta_0 }. \end{equation}

Given ${\rm d}\ln \bar {\varTheta }/{\rm d}z$, $\bar {N}^2(z)$, $N(z)|_{averaged}$ and $\omega _A$, we can solve (7.16) for $\sqrt {{\mathcal {A}}(z)}$, the mean radius of the vortex, using the boundary conditions that ${\mathcal {A}} =0$ at the top and bottom of the vortex. Figure 21 shows the results of solving (7.16) for $\sqrt {{\mathcal {A}}(z)}$ for the CV1 and CV2 vortices, and that the equation works well in describing the vortex shapes.

Figure 21. Comparing $2\sqrt {{\mathcal {A}}(z)}$ (dashed line) to the vertical vorticity $\omega _{z,vortex}(x,y=0,z)/|\,f\sin \theta _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.

8. Constraints due to local convective instabilities

We have found that if an initial vortex has any location in which $N(x, y, z)$ is imaginary, or equivalently, $\partial \varTheta (x, y, z) /\partial z \le 0$, then the flow is, as expected, locally unstable to convection. We can show that the convective instability limits how close the tops and bottoms of vortices can be to the vortex mid-planes.

For all $x$ and $y$ locations within a vortex, we define $z_{top}(x, y)$ and $z_{bottom}(x, y)$ to be the top and bottom of the vortex, where the tops and bottoms are the locations where $\tilde {P}(x, y, z)$ is approximately zero. We define $z_{cross}(x, y)$ to be the location near the mid-plane where $\tilde {\varTheta }(x, y, z) = 0$ (where the red and blue lines cross in figure 22a). The requirement that

(8.1)\begin{equation} N^2(x,y,z)=\frac{g}{\bar{\varTheta}}\,\frac{\partial \varTheta}{\partial z}\geq0 \end{equation}

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}(x, y)$ to the top of the vortex at $z_{top}(x, y)$:

(8.2)$$\begin{gather} \tilde{P}(x, y, z_{cross})/[g\,\bar{\rho}(z_{cross})] = \int_{z_{bottom}}^{z_{cross}} {\rm d}z\,\tilde{\varTheta}(x, y, z)/\bar{\varTheta}(z), \end{gather}$$

(8.3)$$\begin{gather}\tilde{P}(x, y, z_{cross})/[g\,\bar{\rho}(z_{cross})] =-\int_{z_{cross}}^{z_{top}} {\rm d}z\,\tilde{\varTheta}(x, y, z)/\bar{\varTheta}(z). \end{gather}$$

Figure 22. The potential temperature of an anticyclone (CV1). The horizontal black dashed line is the top of the convection zone at $z_{convection} =10$ bar. (a) The blue curve is potential temperature $\varTheta (z)$ along the central axis of the vortex. The dotted red curve is $\bar {\varTheta }(z)$, and $z_{cross}$ is the height near the mid-plane where the two curves cross. (b) The blue curve is potential temperature anomaly $\tilde {\varTheta }(z)$, and $z_{cross}$ is the height near the mid-plane where $\tilde {\varTheta }=0$ on the blue curve. The green dashed line is where $\tilde {\varTheta }=0$.

Figure 23(a) shows schematically the potential temperature $\varTheta$ of an anticyclone similar to the one in figure 22(a). Figure 23(c) shows the corresponding $\varTheta /\bar {\varTheta }$. 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 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 $\partial \varTheta /\partial z \ge 0$ to be stable, but to minimize $z_{top}$, it needs to have $\partial \varTheta /\partial z$ as small as possible. Thus the vortex with the minimum $z_{top}$ has $\partial \varTheta /\partial z =0$ at $z \ge 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 23. Schematic showing the minimum possible height of the top of the vortex $z_{top}$ without being convectively unstable; see text for details. The dashed horizontal lines show the top of the convective zone at $z_{convection}$; the dotted horizontal lines show $z_{cross}$. The values of $z_{cross}$ and $z_{convection}$ are the same in all four plots. (a,b) The red dotted curves are $\bar {\varTheta }(z)$, and the blue solid curves are $\varTheta (x, y, z)$ at some fixed $(x, y)$ within the vortex: (a) for a vortex with $z_{top}$ greater than the minimum value; (b) for the vortex with the same flow at $z < z_{cross}$ as in (a), but with a minimum $z_{top}$. Local convective stability requires that the slope is $\partial \varTheta /\partial z \ge 0$. (c,d) The vertical red dotted curves are $\varTheta (z)/\bar {\varTheta }(z) \equiv 1$, and the blue solid curves are $\varTheta (z)/\bar {\varTheta }(z)$ at $(x, y)$: (c) for the vortex in (a), and (d) for the vortex in (b). In (c,d), hydrostatic equilibrium requires that the area of the lobes at $z > z_c$ equals the area of the lobes at $z< z_c$.

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 $\varTheta < \bar {\varTheta }$. This means that $\partial \varTheta /\partial 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.

Figure 24. Schematic with curves for $\varTheta (z)$ and $\bar {\varTheta }(z)$, and lines for $z_{cross}$ and $z_{convection}$, plotted as in figure 23(a). If $z_{cross} < z_{convection}$, then $\varTheta (z_{convection}) < \varTheta (z_{cross})$, so $\partial \varTheta /\partial z < 0$ somewhere in the range $z_{cross} < z < z_{convection}$, and the vortex is locally unstable to convection.

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.

1. (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 $\bar {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.

2. (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 of accuracy or stability due to its novel time-stepping scheme that is an improvement upon the Krylov-like method of Barranco & Marcus (2006). (See the Appendix.)

3. (iii) There is a stable family (that we call $\mathrm {CV}\textit {--}\mathrm {Rot}$) of vortices in which the vertical vorticity is nearly uniform in the direction parallel to the planet's spin axis, but the horizontal cross-sectional area of the vortex varies with the distance in that direction. The tops and bottoms of the vortices are the locations where the cross-sectional area goes to zero, rather than locations where the vertical vorticity goes to zero.

4. (iv) An alternative proposed family ($\mathrm {CA}\unicode{x2013}\mathrm {Rot}$) of vortices (Galanti et al. 2019; Parisi et al. 2021) has its horizontal cross-sectional area nearly constant in the direction parallel to the planet's spin axis, but the magnitude of the vertical vorticity varies in that direction. The tops and bottoms of the vortices are the locations where the vertical vorticity goes to zero, rather than the locations where the horizontal cross-sectional area goes to zero. We have shown analytically that these proposed vortices are far from equilibrium. The initial-value code shows that the $\mathrm {CA}\unicode{x2013}\mathrm {Rot}$ vortices break apart and re-assemble as one or more $\mathrm {CV}\textit {--}\mathrm {Rot}$ vortices.

5. (v) We develop an algorithm for creating families of vortices that are close to equilibrium. Many properties of the vortices, such as the heights of their tops and bottoms, can be specified, as well as their horizontal velocities at the cloud-top levels. When the initial-value code is initialized with a vortex from the $\mathrm {CV}\textit {--}\mathrm {Rot}$ family, it often comes to a final quasi-steady state close to its initial conditions. Thus there is an easy way to produce long-lived, Jovian-like vortices.

6. (vi) The computed anticyclonic vortices have cool tops (approximately $5\ {\rm K}$), which is similar to the $\sim$8 K cool top observations (Hanel et al. 1979; Flasar et al. 1981; Orton et al. 1996; Sada, Beebe & Conrath 1996; Simon-Miller et al. 2002; Fletcher et al. 2010; Orton et al. 2022) and warm bottoms with respect to the surrounding atmosphere, similar to Jovian anticyclones.

7. (vii) An unexpected result is that initially non-hollow vortices evolve to quasi-steady hollow vortices, like the GRS (i.e. with a high-speed thin annulus at the outer edge of the vortex, and a quiet centre with a minimum of vorticity).

8. (viii) A second unexpected result is that at the outer edges of quasi-steady final anticyclonic vortices, there is a thin annular region that is warmer than the ambient fluid. The initial conditions do not have this annular temperature anomaly; it is created when a down-welling flow (i.e. $v_z < 0$) forms spontaneously within the annulus. These warm outer annuli might explain the rings of bright infrared emissions that surround the GRS and many of the other Jovian anticyclones observed by de Pater et al. (2010) and Fletcher et al. (2010).

9. (ix) We find new scaling laws for vortices in a vertically unbounded rotating, stratified atmosphere (rather than that confined vertically in a thin shell as in Chapter 6 of Pedlosky (1982), i.e. so that the vertical scale height of the horizontal flow is much larger than that of the vertical flow). We show that the vertical scale $H$ is the same for the horizontal and vertical components of the velocity. We determine how the ratio of the magnitudes of the vertical to horizontal velocities scales with Rossby number $Ro$ and $H/L$, where $L$ is the horizontal scale of the velocities. We determine how the divergence of the velocity and other properties of the flow scale with $Ro$ and $H/L$.

10. (x) Vortices whose central axes are initially aligned with the local gravity (i.e. aligned with the spherical radial direction of the planet) evolve so that their central axes are approximately aligned with the planet's axis of rotation.

11. (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.

12. (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.

13. (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 $\beta$ 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 $\sim$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 $\mathrm {CV}\textit {--}\mathrm {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 $\mathrm {CV}\textit {--}\mathrm {Rot}$ vortex that evolves into a long-lived vortex as large as the GRS, and with a horizontal velocity field at $\sim$1 bar that matches it quantitatively? Can we start with small-area $\mathrm {CV}\textit {--}\mathrm {Rot}$ vortices that evolve into long-lived vortices with sizes and horizontal velocity fields at $\sim$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 $\sim$4.5 year thermal dissipation time of the atmosphere near 1 bar to determine whether the numerically computed Jovian vortices are still long-lived.