12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

QG SIMULATION OF JUPITER’S GREAT RED SPOT WITH DISCRETE
EXTERIOR CALCULUS

Pankaj Jagad
Mechanical Engineering, PSE Division
King Abdullah University of Science and Technology
Thuwal, Saudi Arabia
pankaj.jagad@kaust.edu.sa

Ravi Samtaney
Mechanical Engineering, PSE Division
King Abdullah University of Science and Technology
Thuwal, Saudi Arabia
ravi.samtaney@kaust.edu.sa

ABSTRACT

ing GRS using a DEC scheme extending our previous work(Jagad et al., 2021).

Simulating long lived coherent structures such as the great
red spot (GRS) on Jupiter is a challenge from a modeling and anumerical standpoint. Presently, we employ a reduced-gravity
quasi geostrophic (QG) formulation to simulate the GRS. Furthermore, we employ a discrete exterior calculus simulation
technique. DEC is known to be a robust numerical methodwith strong structure preserving properties. To setup the physical conditions we decompose the potential vorticity into threeparts: one governs the potential vorticity in the zonal flow,
while a second part comprises of the potential vorticity in
the GRS vortex, and finally an intermediate potential vorticity field which accounts for the deflection of the flow aroundthe GRS. An iterative technique is employed that constraints
the potential vorticity to observations of the GRS and zonal
flow velocity fields. Preliminary simulations of the GRS on afull spherical surface are presented.

PHYSICAL SETUP AND NUMERICAL METHOD
We employ 1.5 layer reduced-gravity quasi geostrophic(QG) flow assumption (Ingersoll & Cuong, 1981), wherein ashallow upper layer (containing visible clouds and vortices)of a constant density overlies a much deeper lower layer of aconstant density containing the shearing zonal flow. The twolayers are dynamically equivalent to a single layer with a rigid
bottom topography (parameterizing the flow in the lower layer)
and an effective gravity. The computational domain consistsof a spherical surface rotating about its north-south axis. The
governing equation for the system conserves potential vorticity
written as follows.
∂q+ u · ∇q = 0,∂t

INTRODUCTION
The persistence of the giant vortex, Jupiter’s great red spot(GRS), for over three hundred years has been a research topicof interest for not only astrophysicists but also fluid dynamists.
Modeling and simulating the GRS can shed some light on its
origin as well as longevity. Previous attempts at simulating theGRS have been limited to the β -plane and to a domain size
limited to the GRS plus a small neighborhood surrounding theGRS (Shetty et al., 2007). Here, we consider spherical geometry and a domain comprising of the entire spherical surface.
Discrete exterior calculus (DEC) is a robust numericalmethod with desirable conservation and mimetic properties(Hirani, 2003). Several DEC discretizations of the incompressible Navier-Stokes/Euler equations are known to exactly
conserve mass, vorticity and kinetic energy (Mullen et al.,2009; Mohamed et al., 2016). These structure preserving
attributes make DEC an appropriate choice for investigatingflows dominated by long-lived coherent structures. Moreover, the coordinate independence property of DEC makes isconvenient for investigating flows over arbitrary curved surfaces/manifolds. Presently, we develop a method for simulat-(1)where q is the potential vorticity and u is the velocity vector.
The potential vorticity is expressed ash (θ ) g
1q (θ , φ ,t) = ∇2 − 2 ψ (θ , φ ,t) + f (θ ) + b 2 , (2)LRfo LR
where ψ is the stream function, LR is the Rossby deformation
radius, hb is a function of latitude/colatitude representing thebottom topography, g is the reduced/effective gravity, f is theCoriolis parameter, and f0 is the Coriolis parameter√
at a reference location. The Rossby deformation radius LR = gHo / fo ,(H0 is the mean height of the layer) is a length scale at which
the Coriolis force becomes comparable to pressure forces associated with the hydrostatic equilibrium. The Coriolis parameter f (θ ) ≡ 2Ω cos θ , where Ω is the rate of rotation of
the domain and θ is the colatitude. The GRS drifts westward
with nearly a constant velocity, and in a frame of referencetranslating with GRS, the QG equations (1) - (2) have a steady1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022solution q (θ , φ ). Therefore, in a similar approach as Shettyet al. (2007), we decompose the total potential vorticity into
three components as follows:q (θ , φ ) = q∞ (θ ) + qGRS (θ , φ ) + qINT (θ , φ ) ,Here, we compute the potential energy term (gh) as

gh = f0 ψ − ghb ,(3)where q∞ is the potential vorticity in the zonal flow, qGRS is thepotential vorticity in the GRS vortex, and qINT is the potential
vorticity which accounts for the deflection of the flow around
GRS, and vanishes asymptotically into and away from GRS.
These components can further be expressed as:1
qGRS = ∇2 − 2 ψ GRS ,LR(4)h (θ ) g
1
q∞ = ∇2 − 2 ψ ∞ + f + b 2 ,LRf o LR(5)

1qINT = ∇2 − 2 ψ INT .
LR(6)where the bottom topography term ghb is determined as already discussed before. We use the aforementioned initial condition (see figure 4), and let the flow field evolve to a steadystate.

RESULTS AND DISCUSSIONFigure 1 shows the distribution of qGRS . It reveals thatqGRS is the superposition of four nested patches of nearly uniform potential vorticity, which is in contrast to the assumption of two nested patches for a model of qGRS in the previous
work of Shetty et al. (2007). Our computation shows two additional narrow patches, the inner most and the outer most,whereas the two intermediate patches are similar to that inShetty et al. (2007). Figure 2 shows the distribution of q∞ , andfigure 3 shows its profile as a function of the colatitude. The
computed zonal flow potential vorticity comprises of nearlyhomogenized bands, with each band separated from the adjacent ones by a steep meridional gradient of potential vorticity.
The profile of zonal flow potential vorticity has the appearance of a staircase, and q∞ varies from the north to the southpole non-monotonically. Figure 4 shows the initial potential
vorticity distribution computed from the superposition of theGRS on the zonal flow. Figure 5 shows the distribution of total potential vorticity after 5.5 hours (equivalent to 0.26 eddy
turnover time), which is still evolving. We will examine theflow physics of GRS for much longer simulations times in thefuture.

We employ the observations of the GRS and zonal flow
velocity fields (Mitchell et al., 1981; Limaye, 1986), and with
the GRS superimposed on the zonal flow, we compute the corresponding potential vorticity distribution as an initial guess.
To this end, we compute qGRS directly from equation (4), andcompute q∞ as follows. For the zonal flow, we compute theBernoulli function (B) distribution as∇2 ψ ∞ + f eˆr × u∞ = −∇B,

(10)

(7)

REFERENCES
Dowling, Timothy E & Ingersoll, Andrew P 1989 Jupiter’s
great red spot as a shallow water system. Journal of Atmospheric Sciences 46 (21), 3256–3278.
Hirani, Anil N 2003 Discrete exterior calculus. PhD thesis,California Institute of Technology.
Ingersoll, Andrew P & Cuong, PG 1981 Numerical model oflong-lived jovian vortices. Journal of Atmospheric Sciences
38 (10), 2067–2076.
Jagad, Pankaj, Abukhwejah, Abdullah, Mohamed, Mamdouh& Samtaney, Ravi 2021 A primitive variable discrete exterior calculus discretization of incompressible navier–stokes
equations over surface simplicial meshes. Physics of Fluids
33 (1), 017114.
Kwon, O. K. & Pletcher, R. H. 1981 Prediction of the incompressible flow over a rearward-facing step. Technical Report
HTL-26, CFD-4. Iowa State Univ., Ames, IA.
Lee, Y., Korpela, S. A. & Horne, R. N. 1982 Structure of multicellular natural convection in a tall vertical annulus. In Proceedings of 7th International Heat Transfer Conference (ed.
U. Grigul et al.), , vol. 2, pp. 221–226. Washington, D.C.:
Hemisphere Publishing Corp.
Limaye, Sanjay S 1986 Jupiter: New estimates of the mean
zonal flow at the cloud level. Icarus 65 (2-3), 335–352.
Mitchell, Jim L, Beebe, Reta F, Ingersoll, Andrew P & Garneau, Glenn W 1981 Flow fields within jupiter’s great red
spot and white oval bc. Journal of Geophysical Research:
Space Physics 86 (A10), 8751–8757.

where u∞ is the zonal flow velocity field. We subtract the kinetic energy from it to arrive at the potential energy distribuL2 f 2 (ω+ f )tion gh. Next, we compute q∞ = R 0 gh
, where ω is therelative vorticity. Then, we compute the product of the bottom topography function and the reduced gravity distribution
via equation (5). We assume the GRS center as the reference
location for the computation of f0 . The potential vorticity
distribution computed from the superposition of the GRS on
the zonal flow does not conserve advectively, because it does
not account for the deflection of flow around the GRS. Therefore, it is updated using an iterative procedure as follows. Ina frame translating with the GRS, for a current velocity field,
we compute the updated potential vorticity that is advectivelyconserved. Finally, we update the stream function by inverting the Helmholtz equation (2), and repeat the procedure until
the flow field converges to a steady solution. An alternative
approach is to solve the governing momentum and continuity
equations in primitive variables as follows.

∂u+ u · ∇u = −g∇h + f (θ ) u × eˆr ,∂t(8)∇ · u = 0.

(9)212th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 1. GRS potential vorticity distribution as computed from equation (4) and the observed GRS velocity field, and zoomed inview of GRS showing the details.

Figure 2. Zonal flow potential vorticity distribution as computed from the aforementioned procedure and the observed zonal flow
velocity field
0.00.5
1.0
1.5
2.0
2.5
3.00.0007 0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0.0000

qFigure 3. Zonal flow potential vorticity as a function of colatitude. There are nearly homogenized bands of potential vorticity organized into staircase.

Figure 4.

Potential vorticity distribution as computed from the superposition of GRS on the zonal flow.

Mohamed, Mamdouh S, Hirani, Anil N & Samtaney, Ravi2016 Discrete exterior calculus discretization of incompressible navier–stokes equations over surface simplicialmeshes. Journal of Computational Physics 312, 175–191.
Mullen, Patrick, Crane, Keenan, Pavlov, Dmitry, Tong, Yiying & Desbrun, Mathieu 2009 Energy-preserving integrators for fluid animation. ACM Transactions on Graphics(TOG) 28 (3), 1–8.
Shetty, Sushil, Asay-Davis, Xylar S & Marcus, Philip S 2007
On the interaction of jupiters great red spot and zonal
jet streams. Journal of the Atmospheric Sciences 64 (12),4432–4444.
Sparrow, E. M. 1980a Fluid-to-fluid conjugate heat transfer fora vertical pipe - internal forced convection and external natural convection. ASME Journal of Heat Transfer 102, 402–407.
Sparrow, E. M. 1980b Forced-convection heat transfer in aduct having spanwise-periodic rectangular protuberances.
Numerical Heat Transfer 3, 149–167.
Tung, C. Y. 1982 Evaporative heat transfer in the contact lineof a mixture. PhD thesis, Rensselaer Polytechnic Institute,Troy, NY.

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 5.

Potential vorticity distribution after 5.5 hours (equivalent to 0.26 eddy turnover time).

4