# Exploring Jupiter’s Polar Deformation Lengths with High Resolution Shallow
Water Modeling

Ali Hyder  Department of Astronomy, New Mexico State University, PO BOX 30001,
MSC 4500, Las Cruces, NM 88003-8001, USA.  Wladimir Lyra  Department of
Astronomy, New Mexico State University, PO BOX 30001, MSC 4500, Las Cruces, NM
88003-8001, USA.  Nancy Chanover  Department of Astronomy, New Mexico State
University, PO BOX 30001, MSC 4500, Las Cruces, NM 88003-8001, USA.  Raúl
Morales-Juberías  Department of Physics, New Mexico Institute of Mining and
Technology, Socorro, NM 87801, USA.  Jason Jackiewicz  Department of
Astronomy, New Mexico State University, PO BOX 30001, MSC 4500, Las Cruces, NM
88003-8001, USA.

###### Abstract

The polar regions of Jupiter host a myriad of dynamically interesting
phenomena including vortex configurations, folded-filamentary regions (FFRs),
and chaotic flows. Juno observations have provided unprecedented views of the
high latitudes, allowing for more constraints to be placed upon the
troposphere and the overall atmospheric energy cycle. Moist convective events
are believed to be the primary drivers of energetic storm behavior as observed
on the planet. Here, we introduce a novel single layer shallow water model to
investigate the effects of polar moist convective events at high resolution,
the presence of dynamical instabilities over long timescales, and the
emergence of FFRs at high latitudes. We use a flexible, highly parallelizable,
finite-difference hydrodynamic code to explore the parameter space set up by
previous models. We study the long term effects of deformation length
($L_{d}$), injected pulse size, and injected geopotential. We find that models
with $L_{d}$ beyond $1500$ km (planetary Burger number, Bu$=4.4\times
10^{-4}$) tend to homogenize their potential vorticity (PV) in the form of
dominant stable polar cyclones, while lower $L_{d}$ cases tend to show less
stability with regards to Arnol’d-type flows. We also find that large
turbulent forcing scales consistently lead to the formation of high latitude
FFRs. Our findings support the idea that moist convection, occurring at high
latitudes, may be sufficient to produce the dynamical variety seen at the
Jovian poles. Additionally, derived values of localized horizontal shear and
$L_{d}$ may constrain FFR formation and evolution.

Jupiter, vortex, polar dynamics, shallow water — numerical modeling

††software: The Pencil Code (Brandenburg & Dobler, 2002, 2010; Brandenburg &
Scannapieco, 2020), Numpy (van der Walt et al., 2011; Harris et al., 2020),
Scipy (Jones et al., 2001), Matplotlib (Hunter, 2007)

##  1 Introduction

The variety of cloud morphologies seen in Jupiter’s visible troposphere is
extensive, making it a rich laboratory for the study of a wide range of
dynamical processes. It contains planetary-scale wind fields that are
comprised of high speed jets between the belts and zones, along with long-
lived storm features such as the Great Red Spot (GRS) at the lower latitudes.
Although Jupiter’s jet streams extend up to the mid-latitudes, there is a
distinct difference in appearance beyond $\pm 60^{\circ}$ (Orton et al.,
2017), where there is a complete disintegration of the horizontal banded
structures of the lower to mid-latitudes, and the local turbulent flows self-
organize into more discretized vortex behavior.

An important discovery made by the Jovian InfraRed Auroral Mapper (JIRAM)
(Adriani et al., 2017) instrument on board the Juno spacecraft was that of
clustered polar cyclones at the high latitudes, also referred to as vortex
configurations (Adriani et al., 2018). Such behavior is in stark contrast to
Saturn’s polar hexagon (Godfrey, 1988) and tends to emerge in 2D Euler flows.
Jin & Dubin (1998) showed the emergence of such features in idealized cases
using regional maximum entropy theory, which states that vortex configurations
form due to strong vortices mixing the surrounding fluid, maximizing the local
entropy, thereby limiting the underlying chaotic motion and stabilizing the
system. However, the mechanisms through which such structures emerge in an
ever-changing atmosphere are unclear. It is important to understand what
mechanisms are involved in the formation, structure, and stability of cyclonic
configurations to better constrain the overall atmospheric energy cycle of
Jupiter.

Recently, 3D polar simulations were used to model the formation of such
configurations via turbulent convection, exploring the role of inertial
stability in a convectively unstable background (Cai et al., 2021). Although
this model successfully showed the emergence of 3D polar configurations, here
we show how a simplified 1.5-layer shallow water model can also capture
relevant dynamics pertinent to such instabilities. Such a model allows for an
exploration of a larger parameter space due to lower computational demands.

For the case of Jupiter, a host of simulations have been performed to
understand the global circulation of the system and characterize the energy
cycle. This includes, but is not limited to, full general circulation models
of the entire atmosphere (e.g. Williams, 1978, 1988; Kaspi et al., 2009), as
well as localized process modeling for predictions of storm behavior and
evolution to determine the quasi-steady state of such systems (e.g. García-
Melendo & Sánchez-Lavega, 2017). Shallow water models of Jupiter’s weather
layer have provided significant insights on the nature of the Jovian dynamics,
including the emergence of super rotating jets and long-lived storms (Dowling,
1993; Showman, 2007). However, models that focus on the high latitudes,
specifically the polar regions, have been limited due to geometric constrains
and low resolution.

O’Neill et al. (2015) performed a comprehensive analysis of gas giant polar
regions using a 2-layer shallow water system forced by ‘hetons’, which are
defined as 2-layer vortices with opposite rotation in the vertical direction
(Hogg & Stommel, 1985; Pedlosky, 1985; Yano & Flierl, 1992). They found that
the planetary Burger number, defined as

${\rm Bu}=(L_{d}/a)^{2},$ (1)

quantified the overall equilibrated dynamics and morphology of the atmosphere.
$L_{d}$ is defined as the deformation length and $a$ is the planet’s radius of
curvature at the pole. It is important to note that the planetary Burger
number is distinct from the local Burger number, Bu${}_{l}=(L_{d}/R_{\rm
storm})^{2}$, where $R_{\rm storm}$ is a storm radius. Here, we use the term
“morphology” as a qualitative measure of what is readily observed in the
models, specifically after equilibration of the storm forcing in the long-
term.

Brueshaber et al. (2019, hereafter BSD19) used the 1.5-layer shallow-water
capability of the Explicit Planetary Isentropic Coordinate General Circulation
Model (EPIC GCM; Dowling et al., 1998) to investigate mechanisms involved in
such morphological differences across the upper atmospheres of the giant
planets. They quantified the effects of cyclone to anticyclone ratio, storm
strength, and planetary Burger number, Bu, on the long-term evolution of their
models, and also found Bu to be the most important dimensionless parameter
that characterizes the primary mode of the dynamics. Brueshaber & Sayanagi
(2021) further explored Bu values relevant for Saturn and Saturn-like systems
to study the detailed behavior of emergent large-scale polar cyclones that
dominate most of the polar region.

In this study, we build upon previous work by developing a new shallow water
code that circumvents the polar singularity using the polar-$\beta$ plane
approximation, also called the $\gamma$-plane. The model injects storms
stochastically in the domain, resulting in 2D turbulence and long-term
dynamics. Using our model, we attempt to reproduce the range of observed
Jovian polar morphologies using previously employed parameterizations of moist
convection (O’Neill et al., 2015; Showman, 2007). In §2, we provide non-
dimensionalization of the main parameters and details of our model. In §3, we
present results from our model. In §4 we discuss our results regarding the
emergence of folded-filamentary behavior and low deformation length effects,
and in §5 we present our conclusions and plans for future work using this
model. We provide scalability tests and benchmarks that support the validity
of the model in Appendix A and B, respectively.

##  2 Model Description

Modeling efforts in the Jovian atmosphere community employ well-tested and
robust codes. These include (but are not limited to) the EPIC GCM, Isca
(Vallis et al., 2018), DYNAMICO (Boissinot et al., 2019; Bardet et al., 2019),
and the MITgcm (Marshall et al., 1997). These codes provide a multitude of
numerical tools necessary to model phenomenologically distinct scales of gas
giant atmospheres. For the case of Jupiter, EPIC has served as the workhorse
for those interested in the full primitive equations, while the MITgcm has
been used extensively for the anelastic approximation (Kaspi, 2008). Here, we
introduce an atmospheric module for the shallow water approximation using the
Pencil Code111The Pencil Code is available for download at http://pencil-
code.nordita.org/.

Pencil is a $6^{th}$-order spatial and $3^{rd}$-order temporal finite-
difference code, well-suited for modeling three dimensional weakly
compressible turbulence (Brandenburg & Dobler, 2002, 2010; Brandenburg &
Scannapieco, 2020). It uses a collocated Eulerian mesh instead of a
staggered-C grid that is widely used in other codes. Staggered grids offer
stability with second-order spatial accuracy. Pencil, on the other hand,
offers higher order spatial accuracy using a simple collocated grid with
explicit centered finite differences. Although that is not suitable for
second-order schemes, Pencil achieves high accuracy due to the high-order
nature of the stencil. It has been heavily employed in models of the solar
convective turbulence in the solar corona (e.g. Bourdin et al., 2013),
hydrodynamical instabilities in protoplanetary disks (e.g. Raettig et al.,
2013, 2015; Lyra et al., 2015, 2017, 2018), and is suitable for weakly
compressible hydrodynamical and hydromagnetic turbulence (Haugen et al.,
2004a). Here, we present the first application of Pencil to atmospheric
dynamics with a focus on the meteorological processes of Jupiter’s polar
region. Although Pencil has the capability to perform 3D numerical
simulations, we limit our work to the shallow water approximation to make use
of the high resolution that can be afforded to 2D models.

In order to accurately simulate the non-linear fluid dynamics, the first-
baroclinic Rossby radius of deformation must be resolved. For polar models,
typical values for the Rossby radius of deformation of around $2000$ km have
been used, even though deformation lengths of that scale correspond to
equatorial and mid-latitude regions (although BSD19 used $1000$ km along with
a significantly large planetary radius to simulate a low Bu environment). The
highest resolution achieved in their model was $512^{2}$ ($\Delta x\approx
150$ km), which is sufficient for stable evolution but is not enough to
resolve detailed behavior of the large cyclonic and anticyclonic regions with
increasingly low deformation length. Previous works have generally explored
the poles of Jupiter and Saturn with a resolution of $256^{2}-512^{2}$ (e.g.
Morales-Juberías et al., 2011, 2015; O’Neill et al., 2015; Brueshaber &
Sayanagi, 2021); for Jupiter size scales, such resolutions correspond to
$\Delta x=150-300{\rm\ km}$. Our model takes the base resolution further by
exploring Jupiter’s dynamics with a resolution of $1024^{2}\ (\Delta x=75{\rm\
km},4\times$ more resolution in area than previous shallow water models$)$,
which proves necessary to study the long term effects of the lowest
deformation lengths and resolved forcing scales that might be dynamically
relevant for atmospheric evolution. We test the model up to a resolution of
$2048^{2}$ ($\Delta x=38{\rm\ km},16\times$ more resolution than earlier
models). However, as we explore a vast parameter space, we limit model
resolution to $1024^{2}$.

###  2.1 Shallow Water Approximation

Although multiple layers can be used, we limit our analysis to a single-layer
forced turbulence approach, thereby restricting the parameter space to
barotropic instabilities. Perturbations representing moist convection are made
directly to this layer. The bottom layer, often referred to as the
“$1/2$-layer”, represents the deep atmosphere, which is assumed to be
quiescent. Such systems are conventionally called 1.5-layer models. In
multilayer systems, these perturbations can be generated deeper and their
effects understood in context of the reaction of the surrounding layers.
Single layer systems are limited in terms of the overall dynamics and energy
equilibration as they do not permit baroclinic modes. Although restricted,
they are sufficient for studying the long term evolution of non-stratified
systems. Since the system is limited to a single layer with a “quiescent” deep
layer, the momentum equations in the corotating frame reduce to the following:

$\displaystyle\frac{D(gh)}{Dt}$ $\displaystyle=$ $\displaystyle-
gh{\boldsymbol{\nabla}}\cdot{{\boldsymbol{u}}}+\sum^{n}_{i=1}S_{i,\rm
storm}+S_{\rm sub}+S_{\rm mass}$ (2)
$\displaystyle+D_{3}\nabla^{6}(gh)+{\boldsymbol{\nabla}}\cdot{{\boldsymbol{J}}},$
$\displaystyle\frac{D{\boldsymbol{u}}}{Dt}$ $\displaystyle=$
$\displaystyle-{\boldsymbol{\nabla}}{\left(gh\right)}-f\hat{{\boldsymbol{z}}}\times{\boldsymbol{u}}+\nu_{3}\nabla^{6}{\boldsymbol{u}}+{\boldsymbol{\nabla}}\cdot{{\boldsymbol{\xi}}}$
(3)

where $g$ is the reduced gravity, $h$ is the fluid height ($gh$ is the
“geopotential”), $f$ is the Coriolis parameter, and $\hat{{\boldsymbol{z}}}$
is the unit vector in the vertical (outward) direction.

$\frac{D{}}{Dt}\equiv\frac{\partial{}}{\partial{t}}+{\boldsymbol{u}}\cdot{\boldsymbol{\nabla}}$
(4)

is the Lagrangian derivative operator, and
${\boldsymbol{u}}=u\hat{{\boldsymbol{x}}}+v\hat{{\boldsymbol{y}}}$ is the
horizontal velocity. The coefficients $D_{3}$ and $\nu_{3}$ are $6^{th}$-order
hyperdiffusion and hyperviscosity coefficients, respectively, needed for
numerical stability such that dissipation occurs with a higher-order
differential rather than the standard Laplacian form, maximizing the inertial
range. The last terms in the mass and momentum equations represent divergence
damping, also called divergence limiting (Skamarock & Klemp, 1992), and take
the form of mass diffusion and a simplified bulk viscosity

$\displaystyle{\boldsymbol{J}}\equiv D_{\rm bulk}{\boldsymbol{\nabla}}{(gh)},$
(5) $\displaystyle\xi_{ij}\equiv\nu_{\rm bulk}\,\partial_{j}u_{i}.$ (6)

Divergence damping is employed to filter out high frequency acoustic waves
that are not relevant for long timescale atmospheric dynamics (Klemp et al.,
2018). The coefficient $\nu_{\rm bulk}$ is an artificial bulk viscosity. We
also define $D_{\rm bulk}$ as an artificial “bulk diffusion”. They are equal
in value, defined numerically as

$D_{\rm bulk}=\nu_{\rm bulk}\equiv c_{\rm
bulk}\left<\max_{5}[(-{\boldsymbol{\nabla}}\cdot{\boldsymbol{u}})^{+}]\right>{\left[\min(\Delta
x)\right]}^{2}{.}$ (7)

where the superscript plus sign indicates the positive part of the quantity.
This formulation requires that the divergence limiters are proportional to the
maximum of positive flow convergence, as evaluated over five grid cells in
each direction for a total of 25 zones in 2D (the given cell plus its
immediate neighbors). The angled brackets represent a quadratic smoothing
function that smooths the divergence over seven zones in each direction, with
weights (1, 6, 15, 20, 15, 6, 1)/64. The result is then scaled by the square
of the smallest grid spacing. The quantity $c_{\rm bulk}$ is a constant
defining the strength of the bulk diffusivities, set to unity in our
simulations. The divergence limiters in the equations of motion have the
effect of smoothing a numerical discontinuity until it is resolved by the
stencil. Notably, the formulation in terms of divergences of diffusion fluxes
implies that these extra terms, although artificial, conserve mass and
momentum. This type of divergence limiter, taking the mathematical form of a
bulk viscosity, is routinely used in compressible hydrodynamical simulations
to treat shocks (e.g. Haugen et al., 2004b; Richert et al., 2015; Lyra et al.,
2016; Hord et al., 2017).

$S_{i,\rm storm}$ represents the forcing model and drives the geopotential
gradients via stochastic perturbations (iterated over $i$) to which the
velocity fields must adjust geostrophically. The sum is conducted over each
perturbation. The total number of active injected storms remains fixed
throughout this work and is maintained at $n=38$. This is sufficiently high
for large scale dynamics and interactions to occur over long timescales but is
low enough that emergent features are not immediately disrupted by concurrent
moist convective events. The effect of varying n is not studied in this work.
The choice of n is motivated by scaling the $100$ storm injections used by
BSD19 for a planetary radius of $200,000$ km. As we use Jupiter’s radius
explicitly, we lowered the number of active storms from 100 to 38. $S_{\rm
sub}$ adds or removes the overall mass outside of an injection site to ensure
that the domain-averaged geopotential and deformation length does not vary
(O’Neill et al., 2015, 2016). Finally, $S_{\rm mass}$ represents the mass
relaxation that removes mass by steadily bringing the instantaneous average
thickness to the steady-state over the relaxation timescale. The forms of
equations 2 and 3 closely follow the formulations provided in BSD19 and
Brueshaber & Sayanagi (2021) as the fundamental equations of motion for both
models are the same, apart from the $6^{th}$-order hyperviscosity term that is
implemented in our case.

We provide perturbations to the system as Gaussians in space and time
following BSD19:

$S_{\rm storm}\equiv s_{\rm max}\exp{\left(-\frac{r^{2}}{R_{\rm
storm}^{2}}-\frac{(t-t_{\rm p})^{2}}{\tau^{2}}\right)},$ (8)

where $s_{\rm max}$ is the maximal injection rate, $r$ is the mass injection
site offset from the pole, $R_{\rm storm}$ is the peak injection pulse size
(or perturbation radius), and $\tau$ is the characteristic mass injection
time, which is set to $10^{5}\,$s. $t_{\rm p}$ is the peak mass injection
time, which is set to half of the storm duration. The interval between storm
injection is set to $2\tau$. The mass relaxation term, $S_{mass}$, follows the
same moist convection formulation as in Showman (2007). The subsidence term,
$S_{sub}$, adds or removes the mass injected from a given vortex over the
whole domain, ensuring that geopotential perturbations do not cause a
variation in the steady-state geopotential (BSD19).

Typically, the shallow water continuity equation is written solely in terms of
the fluid height, $h$. However, since the reduced gravity, $g$, remains
constant for the 1.5-layer model, we analyze our simulations in terms of the
geopotential, $gh$ (Dowling & Ingersoll, 1989). We employ a forced turbulence
model in our system, where the perturbations only occur on the geopotential
(O’Neill et al., 2015). This allows for horizontal velocities to respond to
perturbations in the flow and stabilize accordingly.

The mass injections are parameterizations of moist convective events, the
details of which cannot be incorporated in a shallow water approximation. Our
model captures the effective resultant behavior due to such storms being
churned up, perturbing the dynamics, and the consequent winds adjusting
geostrophically. Such parameterization captures the bulk properties of mass
and energy perturbations due to moist convection, even though complexities of
the process are ignored. The perturbation imparts surface gravity waves that
dissipate some of the energy away but most of it remains at the injection
site, driving winds whose maximum velocities depend upon the storm intensity
applied. We thus use estimates of the wind speeds from observations (Orton et
al., 2017; Tabataba-Vakili et al., 2020) to sufficiently constrain the
injection strength to ensure we are not overforcing. The injection strengths
that produce realistic wind speeds are used as upper limits of the forcing
model.

###  2.2 Non-dimensionalization of Physical Parameters

Non-dimensionalization can reduce the overall complexity of the problem and
provide insight into the key length and time scales of a system. In order to
properly reproduce the polar dynamics in an appropriately scaled manner, we
must ensure that we use the correct non-dimensional values of the parameters
we are interested in. We begin our non-dimensionalization with time. A
characteristic timescale for the dynamics of Jupiter can be obtained via its
rotation frequency, $\Omega_{J}=1.76\times 10^{-4}\,$s-1. However, since we
are using the Rossby deformation radius as a given known parameter, we may
also determine the Coriolis frequency, $f_{0}=2\Omega$, using the deformation
length,

$L_{d}\equiv\sqrt{\langle gh\rangle}/f_{0},$ (9)

where $\langle gh\rangle$ is the steady-state geopotential. The unit of time,
$[t]$, is defined as,

$[t]\equiv\frac{1}{\Omega}=\frac{2}{f_{0}}=\frac{2L_{d}}{\sqrt{\langle
gh\rangle}}$ (10)

The choice of steady-state geopotential, $\langle gh\rangle$, and deformation
length sets the rotation rate for the planet. Here, the angled brackets around
the geopotential denote a domain-average. This results in two distinct methods
for choosing the rotation: using the observed rotation rate of Jupiter,
thereby setting $\langle gh\rangle$ uniquely for a given model, or by varying
the rotation rate of the planet explicitly. The two approaches are equivalent
in their non-dimensional form but allow for a distinction between the relevant
parameters being tested. The bottom layer of the fluid can have spatial
variations across the domain, often referred to as ‘dynamical topography’
(Dowling, 2020). For the purposes of this study, we do not implement any such
spatial variations as that would significantly increase the parameter space we
explore in this work.

The most appropriate unit of length for these polar simulations is the Jovian
osculating radius, $a=(R_{e}^{2}/R_{p})=76,452$ km, where $R_{e}$ and $R_{p}$
are Jupiter’s equatorial and polar radii, respectively222$R_{e}=71,492$ km;
$R_{p}=66,854$ km (Bagenal et al., 2007).. Thus, varying the deformation
length directly changes the planetary Burger number, which is a function of
planetary radius and the Rossby deformation length. For each corresponding
planetary Bu value, we produce models with mass injection sizes that satisfy
$R_{\rm storm}$$/L_{d}<1$ (local Burger number, Bu${}_{l}>1$), effectively
exploring a storm-to-deformation length ratio. The deformation radius
represents the length scale where the pressure gradient of a storm, acting
outwards (inwards), is balanced by the Coriolis force acting inward (outward).
At this length scale, rotational effects become as important as local
geostrophic behavior. The deformation radius further modulates the overall
distance two vortices can remain from each other without non-linear
interactions causing disruption or coalescence. The planetary Burger number,
therefore acts as a proxy for rotational dominance in an effectively non-
stratified environment, suitable for the shallow water approximation.

Following BSD19, we set a steady-state geopotential, $\langle
gh\rangle=2\times 10^{5}$ m2 s-2. This parameterization represents a model
suite with varying differential rotation, ranging from $0.63\,\Omega_{J}$ to
$1.69\,\Omega_{J}$. Since we maintain a fixed domain-averaged geopotential
through a given suite of simulations, the timescale solely depends on the
value of the planetary Burger number via equation 1, with smaller Bu implying
faster rotation, as ${\rm Bu}\propto L_{d}^{2}\propto t^{2}$. The
corresponding non-dimensionalizations relevant to our model suite are
presented in Table 1. For the timescale, we present both non-dimensional
scalings in the code depending on the choice of fixed geopotential or fixed
rotation. Finally, as we extend the application of previous work, we define
the nondimensional eddy potential vorticity as in BSD19 to facilitate an
easier comparison:

$Q_{e}^{*}=\left[\frac{\zeta+f}{h}-\frac{f}{\langle
h\rangle}\right]\cdot\frac{\langle h\rangle}{f_{0}},$ (11)

where $\zeta$ is the local relative vorticity, $\langle h\rangle$ is the
steady-state fluid thickness, and $f_{0}=2\Omega$.

Table 1: Nondimensionalization of relevant parameters. These are the scales used to convert from the code units used herein to physically meaningful values. Physical parameter | Model conversion factor
---|---
Length | $a$
Time |  $2L_{d}/\sqrt{\langle gh\rangle}$ or $1/\Omega_{J}$
Energy | $a^{6}\Omega^{4}$
Mass injection | $a^{2}\Omega^{3}$

###  2.3 Geometry and the Gamma-Plane Approximation

Figure 1: Illustration of the $\gamma$-plane. The approximation is equivalent
to a $\beta$-plane at the lower latitudes but applied to the pole. The polar
coordinates are converted to a Cartesian grid for the numerical computation to
avoid the polar singularity.

We use a Cartesian approximation for the model to circumvent pole-centered
singularities that emerge due to a polar geometry. Although other methods are
applicable at the poles, the Cartesian approximation sufficiently captures the
main morphological features of the atmosphere in our model.

The $\gamma$-plane approximation, as used in BSD19 and also employed in the
form of the polar-$\beta$ plane by O’Neill et al. (2015), allows an
exploration of the shear effects dominating the high latitudes. It quantifies
the planetary vorticity gradient and offers parameterization of the length and
time scales directly. In short, $\gamma$ carries the effect of variations in
the planetary Burger number. With a fixed planetary radius, it traces the
behavior of the deformation length. The functionality is implemented using the
following correction to the Coriolis parameter at the pole:

$f=2\Omega-\gamma r^{2},$ (12)

where $\gamma=\Omega/a^{2}$ (Bridger & Stevens, 1980, BSD19). Our model uses
the Jovian osculating radius, $a$, as the domain size and use a doubly
periodic domain. We model latitudes poleward of $60^{\circ}$ North. Figure 1
illustrates the setup of the $\gamma$-plane with the natural polar
coordinates, $(r,\theta)$ where $\theta$ is the azimuthal angle, converted to
Cartesian coordinates $(x,y)$ for our model. The domain of interest extends
from the center of the box to a circle of $0.5a$ in radius, after which a
damping parameter is applied such that a traversing gravity wave is quickly
damped outside the domain. Strongest damping occurs farthest away from the
outer circle to minimize shear generation due to the boundary. We apply this
damping to a variable $X$ in the form of a relaxation term over a timescale
$\tau_{\rm dis}$, so that

$\frac{\partial X}{\partial t}=-\frac{X-X_{0}}{\tau_{\rm dis}}\sc{S(r)},$ (13)

where $X_{0}$ is the desired value. $\sc{S(r)}$ is a $5^{th}$-order step
function that smoothly goes from zero to unity from $0.5a$ to $0.55a$,
respectively. The timescale $\tau_{\rm dis}$ is chosen such that a gravity
wave dissipates within this buffer zone. Beyond $0.55a$, a freezing boundary
condition is used where all derivatives are set to zero (Lyra et al., 2008).
Thus, similar to BSD19, no vortex feature or propagating gravity wave can go
across the damping region and reenter the domain.

##  3 Results

The Jovian polar cyclones have an approximate range of sizes from $4000-7000$
km and exhibit velocities between $\sim 55-115$ m/s as measured from JunoCam
(Tabataba-Vakili et al., 2020) and JIRAM (Grassi et al., 2018) images. We
restrict our model to produce wind speeds comparable to observed velocities,
providing a way to constrain the turbulent forcing strength necessary to
sufficiently perturb the steady-state geopotential. Given $\langle
gh\rangle=2\times 10^{5}$ m2 s-2, we used mass injections that would produce
wind speeds on the order of 100 m/s. Previous models (e.g., O’Neill et al.
(2015); BSD19; Brueshaber & Sayanagi (2021)) used a range of perturbations but
exact values of injection strengths are not well-constrained. Our forced
turbulence model has two forms for two distinct applications: varying Burger
numbers (planetary and local) and varying storm strengths, described in §3.1.1
and 3.1.2 as Case A and Case B, respectively. For the case that utilizes
previously published values for a steady-state geopotential (Case A), we use
$s_{max}=0.25$ m2 s-3, which proves sufficient for our model as the resulting
wind speeds are in agreement with derived velocities. For the case of
Jupiter’s rotation (Case B), we vary the storm strength such that total
injected geopotential is reasonably small compared to the fluid thickness for
that simulation. Otherwise, the same perturbation strength for a lower
geopotential will result in overforcing. This method defines a unique
geopotential injection parameter that allows us to explore the effect of
forcing strength in an atmosphere with Jupiter’s rotation. The details of the
injection strength are provided in §3.1.2.

###  3.1 Parameter Space

As there are two methods for choosing the timescale for the models, we explore
both in order to study the different parameters. Previous models of the polar
region focused more on determining the general characteristics of gas giant
polar dynamics, and thus used timescales derived from a choice of deformation
length and steady-state geopotential (as opposed to using Jupiter’s
differential rotation), as well as significantly larger planetary radii. In
contrast, we used Jupiter’s osculating radius, $a$, explicitly as the length
scale for all of our models and explore the parameter space by studying both
time scales independently. This results in two cases. In Case A, we vary the
planetary and local Burger number by varying deformation lengths and storm
size. The specified steady-state (domain-averaged) geopotential and choice of
deformation length determine the timescale, $\Omega^{-1}$, for these
simulations (Table 2). This implies a unique rotation for each choice of
deformation length. In Case B, we vary the storm strength and use the observed
System III rotation rate of Jupiter, thereby varying the underlying steady-
state geopotential for a choice of deformation length. Our simulations produce
a variety of morphologically distinct and interesting structures. These
behaviors are a result of the higher resolution and computational accuracy of
our model. Our model shows similarity with previous work, in particular
O’Neill et al. (2015), BSD19, and Brueshaber & Sayanagi (2021), and
additionally shows the emergence of distinct novel behavior as described in
detail in the following subsections.

####  3.1.1 Case A - Varying Planetary & Local Burger Numbers

We use the ratio $R_{\rm storm}$$/L_{d}$ as one of the parameters in our
models. Although this ratio is the inverse square root of the local Burger
number, Bul, we vary the ratio itself to have a more direct relation to the
storm size. Brueshaber & Sayanagi (2021) employed this parameter to explore
the dynamical regimes relevant to Saturn and Saturn-like systems, i.e. higher
planetary Burger numbers. Here, we provide a completion of the parameter space
by studying the effects on Jupiter with lower planetary Burger numbers by
using smaller deformation lengths with higher resolution. The values employed
in the model suite are physically motivated due to the overall dynamical
effect of the Rossby deformation length. If Bu${}_{l}<1$, the forced
turbulence model would represent nonphysical storm injections as the imparted
pulse would be larger than the deformation length. As our mass injections
represent physically relevant moist convective events, a deep moist convective
plume would be undergoing deformation due to the Coriolis force as it emerges
in the upper atmosphere. Consequently, it would not be larger than the overall
deformation length. With Bu${}_{l}>1$, all emergent storms are such that they
grow by coalescence or aggregation of surrounding fluid until they attain a
phase of geostrophic adjustment. This does not preclude the existence of
storms with Bu${}_{l}<1$ as such features may form via mergers, but only
limits the type of moist convective plume length scales relevant in the
forcing model.

For a given deformation length, the $R_{\rm storm}$values are tested from
$522-1750$ km in increments of $250$ km (except for the first increment). The
lowest $R_{\rm storm}$value is kept at $522$ km as we maintain a
$6^{th}$-order spatial derivative across the radius, fully resolving the
injected pulse. This is a strict constraint that is maintained throughout the
study. The deformation lengths themselves are varied from $750-2000$ km (Bu
$\sim 9.6\times 10^{-5}-6.8\times 10^{-4}$). This defines a single simulation
for the $L_{d}=750$ km case, two for $L_{d}=1000$ km, three for $L_{d}=1250$
km, four for $L_{d}=1500$ km, five for $L_{d}=1750$, and six for $L_{d}=2000$
km. The results for this model suite are presented in §3.2.

####  3.1.2 Case B - Varying Storm Strength

We ran the simulations in Case A using a fixed domain-averaged geopotential
that was maintained for the whole parameter space, similar to what was done in
BSD19 and Brueshaber & Sayanagi (2021). However, the consequence of such a
system is that the suite represents simulations with varying rotations. To
ensure the applicability of the model to Jupiter’s case, we relax this
constraint. Using Jupiter’s rotation, $\Omega_{J}=1.76\times 10^{-4}$ s-1, we
must vary the storm strengths to prevent over-forcing of the fluid layer. The
time-dependent component of $S_{\rm storm}$ (equation 8) may be integrated as
we do not vary the temporal components of the forcing model. Integrating over
injection duration, $t_{d}=2.2\tau$ (BSD19), the injected storm geopotential,
$gh_{\rm inj}$, is given by

$\begin{split}gh_{\rm inj}&=s_{\rm
max}\int_{0}^{t_{d}}\exp{\left[-\frac{(t-t_{p})^{2}}{\tau^{2}}\right]}dt,\\\
&=s_{\rm max}\frac{\sqrt{\pi}}{2}\tau\left[{\rm
erf}\left(\frac{t_{p}}{\tau}\right)-{\rm
erf}\left(\frac{t_{p}-t_{d}}{\tau}\right)\right],\end{split}$ (14)

where we have evaluated the spatial component at the injection site ($r=0$).
The values of $s_{\rm max}$ are chosen such that the mass injection never
exceeds $\sim 20\%$ of the total fluid thickness, which ensures that the
produced winds will remain within observed values. We define a threshold
parameter for the forcing using the ratio of the geopotentials as

$\Phi_{r}\equiv\frac{gh_{\rm inj}}{\langle gh\rangle}\times 100\%,$ (15)

where $\langle gh\rangle$ is the domain-averaged steady-state geopotential
given for the choice of $L_{d}$. We express it as a percentage of the total
geopotential and vary it incrementally to observe the overall effect of storm
strength over fixed rotation for the cases of $L_{d}=(750,1250,2000)$ km,
producing a suite of $9$ models. The details of the models are presented in
§3.4.

###  3.2 Dependence on Deformation Length

Figure 2: Final snapshots of Case A simulations. Non-dimensional eddy
potential vorticity, $Q_{e}^{*}$, as a function of deformation length and
injected pulse size at day 20,000. We use the same color scheme as O’Neill et
al. (2015, 2016) where red indicates cyclonic motion and blue indicates
anticyclonic motion. Top to bottom: Increasing deformation length, $L_{d}$,
from $750$ km to $2000$ km with $250$ km increment. Left to right: Increasing
the injected pulse size, $R_{storm}$, from $\approx 500$ km to $1750$ km with
$250$ km increment. Each colorbar pertains to the column of simulations above
it. Resolution: $1024^{2}$ ($\Delta x=75$ km). Movies for these Case A
simulations are available on Zenodo: https://doi.org/10.5281/zenodo.6642986
(catalog doi:10.5281/zenodo.6642986) (Hyder et al., 2022). The movies show
that primary distinction between the deformation lengths manifest early on
even though the overall dynamics continue to evolve.

We produced a suite of 21 simulations that parameterize the entire domain of
the observed deformation lengths on Jupiter, and their relation to the size of
the injected pulse. The model suite is outlined in Table 2. Although moderate
resolutions are sufficient for lower latitudes where the Rossby deformation
length has values upwards of $1500$ km, high resolution is required to model
polar dynamics where the value tends to decrease (Read et al., 2006). In this
suite, the lowest deformation length achieved was $750$ km. The smallest
possible perturbation that could be made on such a domain was $522$ km,
ensuring that the perturbation radius was fully spatially resolved on the
stencil.

Table 2: Modeled parameter space for Case A runs shown in Figure 2. Here, the geopotential is maintained at $\langle gh\rangle=200000$ m2 s-2, and all runs have $s_{\rm max}$= 0.25 m2 s-3.  Model | $L_{d}$ | $R_{\rm storm}$ | $\Omega$ | Bul | Bu |  |  |  | Model | $L_{d}$ | $R_{\rm storm}$ | $\Omega$ | Bul | Bu
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
# | (km) | (km) | $\times\,\Omega_{J}$ | - | $\times 10^{-4}$ |  |  |  | # | (km) | (km) | $\times\,\Omega_{J}$ | - | $\times 10^{-4}$
A1a | 750 | 522 | 1.69 | 2.06 | 0.96 |  |  |  | A5b | 1750 | 750 | 0.73 | 5.43 | 5.24
A2a | 1000 | 522 | 1.27 | 3.67 | 1.71 |  |  |  | A5c | 1750 | 1000 | 0.73 | 3.07 | 5.24
A2b | 1000 | 750 | 1.27 | 1.78 | 1.71 |  |  |  | A5d | 1750 | 1250 | 0.73 | 1.96 | 5.24
A3a | 1250 | 522 | 1.02 | 5.72 | 2.67 |  |  |  | A5e | 1750 | 1500 | 0.73 | 1.36 | 5.24
A3b | 1250 | 750 | 1.02 | 2.78 | 2.67 |  |  |  | A6a | 2000 | 522 | 0.64 | 14.7 | 6.84
A3c | 1250 | 1000 | 1.02 | 1.56 | 2.67 |  |  |  | A6b | 2000 | 750 | 0.64 | 7.11 | 6.84
A4a | 1500 | 522 | 0.85 | 8.26 | 3.85 |  |  |  | A6c | 2000 | 1000 | 0.64 | 4.00 | 6.84
A4b | 1500 | 750 | 0.85 | 4.00 | 3.85 |  |  |  | A6d | 2000 | 1250 | 0.64 | 2.56 | 6.84
A4c | 1500 | 1000 | 0.85 | 2.25 | 3.85 |  |  |  | A6e | 2000 | 1500 | 0.64 | 1.78 | 6.84
A4d | 1500 | 1250 | 0.85 | 1.44 | 3.85 |  |  |  | A6f | 2000 | 1750 | 0.64 | 1.31 | 6.84
A5a | 1750 | 522 | 0.73 | 11.3 | 5.24 |  |  |  |  |  |  |  |  |

The suite shows the overall effect on nondimensional potential vorticity,
$Q_{e}^{*}$ (equation 11), of varying deformation length while keeping the
injected pulse size fixed for each column (Figure 2). Deformation length
$>1500$ km (and $1250$ km to a lesser extent) with the smallest perturbation
exhibits the formation of dominant polar cyclones, as expected since the
deformation length is directly related to the planetary Burger number. Large
values of $L_{d}$ produce Saturn-like behavior, where large scale polar
cyclones dominate most of the flow, supporting the findings of BSD19.

In simulation A1a, (where $L_{d}=750$ km), turbulent self-organization of the
emergent features results in a far more incoherent dynamical appearance as
compared to the dominant polar cyclones at high $L_{d}$. The deformation
length modulates the proximity at which such storms can remain stable in each
other’s presence. With perturbations set to the smallest of scales, the flow
field can evolve the mass pulses and allow them to interact, coalesce and
stabilize over time. A small deformation length prevents them from
aggregating, so most of the small storms that emerge remain relatively small
and coherent. Those that grow beyond the deformation length are affected by
the Coriolis force. The potential vorticity homogenizes over the domain and a
cyclonic pattern similar to an $m=3$ wave emerges in simulation A1a. Although
this feature evolves in time due to the forcing model, the overall behavior
remains the same. We note that there is a dependence on the changing zonal
wavenumber on deformation lengths. This is seen very clearly in the movies
where the wave number behavior emerges early on in the simulations (Hyder et
al., 2022). For large planetary Bu, the low wave number begins to dominate
much earlier in the evolution than low planetary Bu cases. In §6 below, we
explain how the main zonal wavenumber emerging in our different simulations is
controlled by the size of the largest balanced vortex that emerges from our
forcing in each case.

Figure 2 shows that simulations with $L_{d}=1500$ km and above (A4a, A4b, A5a,
A5b, etc.) exhibit more intense cyclonic features in the domain due to
increased deformation length and increased perturbation size. The maximum
dimensional potential vorticity from the suite reaches $\sim 4\times 10^{-3}$
s m-2, which is much higher than other models. However, potential vorticity is
expected to increase at the pole (Showman, 2007), and these simulations
represent an extreme parameter space (high perturbation radii). The peak
amplitude of the mass injection is the same throughout the suite, but varying
$R_{\rm storm}$ varies the FWHM of the Gaussian injection. Larger values for
$R_{\rm storm}$ input significantly more energy into the system. This drives
stronger winds and causes a higher degree of forcing on the overall domain.
Although the winds are realistic, the emergent morphological features are
different from what is observed on Jupiter’s poles. Although such values of
$R_{\rm storm}$$/L_{d}$ do not yield Jupiter-like dynamics, they are similar
in appearance to cyclonic folded-filamentary regions observed on Jupiter
across various latitudes, which are discussed in more detail in §3.3.

It is worth noting that solely employing small $R_{\rm storm}$$/L_{d}$ values
is insufficient to produce a Jupiter-like pole, where the vortices are stable
and do not merge. Indeed, simulation A6a (where $L_{d}=2000$ and Bu${}_{l}\sim
15$) produces a singular polar cyclone that dominates most of the polar region
upwards of $80^{\circ}$. This behavior breaks down with decreasing $L_{d}$
values. Therefore, to approach a realistic Jovian pole, one must model
increasingly small deformation lengths with varying injected pulse size. This
mechanism is not exhaustive as other parameters certainly play a role in the
emergence of vortex configurations on Jupiter and their overall stability
against the perpetual forcing from deeper atmospheric convection as shown by
Cai et al. (2021).

###  3.3 Folded-Filamentary Regions

Jupiter’s midlatitudes are dominated by anticyclones. Most of the cyclonic
features are expected to travel polewards via the Beta-Gyre Effect (or Beta
Drift) (Scott, 2011). However, cyclonicity can survive at low latitudes in the
form of folded-filamentary regions (FFRs) and brown barges. FFRs are elongated
cyclonic structures that exist throughout the Jovian atmosphere and were first
modelled numerically in Marcus (2004). In that work, Marcus (2004) notes that
clouds of elongated forms are consistent with long-lived cyclonic behavior,
and are dissimilar to clouds that form in anticyclonic regions.

There is observational evidence for such cyclonic filamentary structures at
the high latitudes in JunoCam images returned by the Juno spacecraft (Orton et
al., 2017). Our shallow water model successfully shows the emergence of FFRs
at the high latitudes, given certain shearing conditions above $R_{\rm storm}$
$\approx 1200$ km. To our knowledge, this is the first such model to produce
high latitude cyclonic FFRs and quantify a length scale beyond which such
structures tend to emerge. Formation of FFR-like beta-skirts, which are
regions of elevated potential vorticity surrounding a primary vortex
(Montgomery et al., 2006), are observed in the models by O’Neill et al. (2015)
and BSD19 at lower Bu. They note that such features tend to emerge in the
beta-skirts surrounding the primary cyclones. In our models, although elevated
beta-skirts that exhibit filamentary behavior are seen around primary cyclones
in the low $L_{d}$ cases, large scale FFR formation with significantly high
mean horizontal shear is seen solely for the high deformation lengths with a
perturbation scale beyond $1200$ km. Furthermore, the high $R_{\rm storm}$
FFR-like features tend to be the main cyclonic features of those simulations
rather than secondary features emerging in the vicinity of a cyclone.

Figure 2 shows the normalized form of non-dimensional potential vorticity. It
captures the skew-symmetric component of the velocity gradient tensor as it
contains the local vorticity. However, to understand the dynamical behavior of
FFRs that emerge in the high $L_{d}$ models, we use the rate of strain tensor,
which captures the symmetric component of the velocity gradient tensor and
takes the form

$S_{i,j}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial
u_{j}}{\partial x_{i}}\right),$ (16)

where $i,j$ represent the Cartesian coordinates. Physically, the rate of
strain tensor represents the elongation of fluid elements and exhibits the
total effect of local divergence and shear. We compute the magnitude of the
tensor using the Frobenius norm (Golub & Van Loan, 1996):

$\mid{\bf{S}}\mid=\sqrt{\sum_{i=1}^{m}\sum_{j=1}^{n}\mid S_{i,j}\mid^{2}}.$
(17)

where $m=n=2$, as our models are limited to two dimensions. The norm of the
rate of strain tensor encapsulates the combined effect of divergence that
occurs where the FFRs fold over, and shear that occurs where fluid is deformed
due to elongation.

Figure 3: Rate of strain tensor magnitudes for models A4d, A5d, A5e, A6d, A6e, and A6f. The strain rates are significantly higher at these regions relative to their environment. Figure 4: Mean shear magnitude at equilibrium for Case A simulations. The points within the red rectangle pertain to the six simulations from Figure 3 that exhibit FFR-type behavior. Table 3: Dynamical quantities of interest for models that exhibit FFR-type behavior as shown in Figure 3. The angle brackets denote a spatial average over the domain. Model | $\textit{v}_{\rm{max}}$ | $\langle v\rangle$ | ${(\nabla\cdot{\bf{v}})}_{\rm{max}}$ | $\langle\mid\nabla\cdot{\bf{v}}\mid\rangle$ | $\mid{\bf{S}}\mid_{\rm{max}}$ | $\langle\mid{\bf{S}}\mid\rangle$
---|---|---|---|---|---|---
# | ms-1 | ms-1 |  $\times 10^{-2}$ (s-1) |  $\times 10^{-4}$ (s-1) | (s-1) |  $\times 10^{-3}$ (s-1)
A4d | 121 | 23 | 6.953 | 1.680 | 0.432 | 7.147
A5d | 137 | 27 | 3.329 | 1.890 | 0.289 | 7.237
A5e | 156 | 29 | 6.640 | 2.996 | 0.467 | 7.899
A6d | 143 | 28 | 9.691 | 1.877 | 0.360 | 6.720
A6e | 152 | 32 | 2.038 | 2.278 | 0.220 | 7.464
A6f | 189 | 35 | 4.344 | 3.004 | 0.352 | 9.039

Figure 3 shows a number of interesting features. The plot shows the rate of
strain tensor magnitudes using equation 17 of a subset of models from Figure
2, which reflect the behavior of increased $R_{\rm storm}$. The strain-rates
are able to capture the high shear behavior in folded regions well. Such FFR-
like behavior is morphologically distinct from FFRs in beta-skirts around
cyclones, as the FFRs in our simulations represent the primary mode of
cyclonicity in the flow. The corresponding divergence, shear, and maximum wind
speed values are presented in Table 3.

Although cyclonic filamentary regions have been observed at the high latitudes
by Juno, their underlying wind speeds remain unconstrained. Thus, we limit the
emergent mean wind speed of models with such elongated structures to the
observed zonal speeds of the ambient winds and storms, ensuring that we remain
under the overforcing regime. The strength of the forcing, and more precisely,
the injected geopotential relative to the total geopotential, may have a
significant effect on the emergence of such morphology. Our model shows the
formation of FFRs, but is strongly dependent on the mass injection size
(perturbation radii). Figure 4 supports the idea that emergent FFRs are
regions of dominant shear. The presence of such features at Jupiter’s polar
regions as seen in Juno observations suggests that FFRs form either due to
intense vortex interactions with large moist convective events enabled due to
the planetary vorticity gradient, or they emerge due to other mechanisms not
captured by the shallow water approximation.

###  3.4 Forcing Strength

BSD19 and Brueshaber & Sayanagi (2021) found that storm energetics did not
play an important role in the equilibrated dynamical state of the system. In
their models, they used a fixed geopotential (our Case A). Here, we use
Jupiter’s observed rotation rate (Case B), and vary the forcing strengths that
may be relevant to the Jovian pole in particular. In order to vary the
injected geopotential relative to the total geopotential, we use $\Phi_{r}$,
as defined in equation 15. Deformation lengths of $750$, $1250$, and $2000$ km
were used for this model suite. For each choice of $L_{d}$ $(Bu)$, we varied
the injection strength such that $\Phi_{r}\approx 7\%,14\%$, or $20\%$. The
simulations are presented in Figure 5 and the parameters are provided in Table
4.

Table 4: Modeled parameter space for Case B runs as shown in Figure 5. Here, the storm strength is varied while Bul is kept constant for the choice of Bu.  Bu$=9.6\times 10^{-5}$ |  |  |  | Bu$=2.67\times 10^{-4}$ |  |  |  | Bu$=6.84\times 10^{-4}$
---|---|---|---|---|---|---|---|---
Bu${}_{l}=2.06$ |  |  |  | Bu${}_{l}=5.72$ |  |  |  | Bu${}_{l}=14.7$
Model | $s_{\rm max}$ | $\Phi_{r}$ |  |  |  | Model | $s_{\rm max}$ | $\Phi_{r}$ |  |  |  | Model | $s_{\rm max}$ | $\Phi_{r}$
# | (m2 s-3) | % |  |  |  | # | (m2 s-3) | % |  |  |  | # | (m2 s-3) | %
B1a | 0.028 | 7 |  |  |  | B2a | 0.079 | 7 |  |  |  | B3a | 0.202 | 7
B1b | 0.057 | 14 |  |  |  | B2b | 0.158 | 14 |  |  |  | B3b | 0.404 | 14
B1c | 0.085 | 20 |  |  |  | B2c | 0.237 | 20 |  |  |  | B3c | 0.606 | 20

We note that the main differences in the observed variations in dynamical
behavior occur due to the deformation length. However, the storm strengths
have an important impact on the low deformation length models, particularly
with regard to how potential vorticity homogenizes over the domain. As
expected, simulation B1a produces low values for $Q_{e}^{*}$ and results in
low overall wind speeds. Higher wind speeds are produced for higher $\Phi_{r}$
cases. Although all simulations show the formation of a polar cyclone to some
degree, the cyclonic beta-skirts of the low $L_{d}$ simulations, in general,
homogenize distinctly from the high $L_{d}$ cases. This holds true for all
values of $\Phi_{r}$. The cyclonic potential vorticity aggregates at the poles
as expected, but it tends to form an almost triangular pattern, similar to the
$m=3$ pattern seen in low $L_{d}$ models in Case A. Although, the final
snapshot of simulation B1a shows an almost bimodal structure in $Q_{e}^{*}$,
it also exhibits the $m=3$ wave behavior throughout its evolution. If an
anticyclonic region is surrounded by cyclonic fluid, as is seen in the final
snapshot of simulation B2a, the anticyclonic patch is ejected over time due to
the beta-gyre effect, leaving a net accumulation of cyclonicity at the pole.

Figure 5: $Q_{e}^{*}$-dependence on forcing strength and deformation length.
Top to bottom: Increasing deformation length, $L_{d}=\\{750,1250,2000\\}$ km.
Left to right: Increasing geopotential ratio, $\Phi_{r}=\\{7\%,14\%,20\%\\}$.
Movies for these Case B simulations are available on Zenodo:
https://doi.org/10.5281/zenodo.6642986 (catalog doi:10.5281/zenodo.6642986)
(Hyder et al., 2022). Similar to Case A, the movies show a stronger dependence
on deformation length compared to the forcing strength. However, for smaller
deformation lengths the effect of forcing strength is much larger.

Although the dynamics evolve in time, the energy stabilizes relatively early
in the simulations. Shallow water available potential energy (APE) and kinetic
energy (KE) are defined as

$\displaystyle APE$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int(gh)^{2}-{\langle gh\rangle}^{2}\ dA,$ (18)
$\displaystyle KE$ $\displaystyle=$ $\displaystyle\frac{1}{2}\int
gh(u^{2}+v^{2})\ dA,$ (19)

where $dA$ is an areal element of the domain, and the integrals are taken over
the entire box. The APE tends to be an order of magnitude higher than the
kinetic energy even as the simulations equilibrate. Figure 6 shows the
energies relevant to Case B. The energies have been normalized by the final
values of the $\Phi_{r}=7\%$ simulations, and averaged over $L_{d}$. The
shaded regions show the $1\sigma$ variation that is due to changing
deformation lengths. In terms of the energy, we note that the effect of
varying forcing strength is much larger than varying deformation length.
$\Phi_{r}$ directly affects the KE to APE conversion point. Simulations with
stronger forcing tend to convert their KE into APE earlier in their evolution,
whereas simulations with weak forcing exhibit a less efficient conversion
between KE and APE. This can be seen in the Case B movies where large
planetary Bu simulations manifest the low wave number behavior early in their
evolution compared to the low planetary Bu cases (Hyder et al., 2022). The
models stabilize by $10^{4}$ Earth days in KE while APE tends to equilibrate
more slowly, in agreement with Showman (2007). APE remains about an order of
magnitude higher than KE. Our APE remains above KE by a larger factor than the
low latitude shallow water study performed by Showman (2007), but agrees with
polar modeling results from BSD19.

Figure 6: Available potential energy (APE) and kinetic energy (KE) for
simulations in Case B, where $L_{d}$ has been averaged over. The energies have
been normalized by the final values of the $\Phi_{r}=7\%$ simulations. The
shaded regions show the $1\sigma$ variation across the various deformation
lengths. Changing $\Phi_{r}$ has a significant effect on when the APE becomes
dominant over KE.

###  3.5 Evolution of the Characteristic Length

The ratio of APE to KE is directly related to the length scale of emergent
vortices via

$\frac{APE}{KE}=C\left(\frac{L}{L_{d}}\right)^{2},$ (20)

where $L$ is the size of the largest balanced vortex/eddy, and $C$ is a
constant of proportionality (Pedlosky, 1979; Marshall, 2008). We determine $C$
by using simulation B3a. The characteristic size of the largest balanced
vortex emerging in this simulation is $\sim 0.55a$ (where $a$ is the domain
size). Using the final $APE/KE$ value for simulation B3a and $L\sim 0.55a$,
$C$ is found to be about $0.025$. We show $L/a$ as a function of time for the
simulations in Case B in Figure 7. Given $L$, the maximum number of eddies,
$N$, in a domain with size $a$ would be proportional to $a^{2}/L^{2}$.

Figure 7: Characteristic length of the dominant balanced vortex, $L/a$,
plotted as a function of time for all Case B simulations. The profiles plateau
as the simulations stabilize their APE to KE ratios.

The wavenumber behavior observed in our simulations comes as result of eddies
balancing out to occupy a circular domain. It is dependent on how the $N$
vortices self-organize to minimize the total energy. The $m=1$ behavior of
simulations B3a, B3b, and B3c corresponds to how the dominant polar vortex
centralizes in the domain, given $L/a\sim 0.55$. This leaves no room for any
other major storm to share the domain. For simulations B1c and B2a, an $m=3$
pattern emerges that makes these two simulations look almost identical in the
last simulated day. As seen in Figure 7, these two simulations reach
approximately the same value of $L/a\approx 0.3$ by the end of the simulation.
This value is similar to the one that would be obtained from circle-packing
arguments (L/a=0.333) (see e.g., Pirl, 1969) for a total number of vortices
equal to 7, which is the number of large balanced vortices seen filling the
domain in Figure 5 (one in the center of the domain and six of alternating
signs surrounding it). The same argument can be used to explain the other
patterns that emerge from our simulations. For example, for simulations B1a
and B1b the number of balanced vortices that emerge by the end of the
simulations is $N\sim\\{13,9\\}$, which from circle packing arguments would
correspond to values of L/a of 0.236 and 0.277 respectively (Pirl, 1969;
Fodor, 2003), while the values we obtain from our analysis are 0.225 and 0.25.
While our L/a values are in good agreement with those derived from circle
packing arguments, they are always slightly smaller, since the largest
balanced vortices that emerge from our simulations cannot be tightly packed.

These same arguments can also be used to explain the patterns that emerge in
the simulations in case A. The energy analysis for simulations A4a, A5a, and
A6a reveals that the ratio of L/a$\sim 0.5$, which would explain why most
panels in Figure 2 look like there are only two large balanced vortices
sharing the domain. The main difference with case B is that in this case none
of these eddies occupy the center of the domain, therefore allowing for an
additional dominant vortex, with a characteristic length scale of $\sim 0.5$,
to remain stable.

###  3.6 Dynamical Instability

Dynamical instabilities are a consequence of a system’s inability to withstand
perturbations to the mean flow. As small perturbations grow in an unstable
system, they tend to have long term effects on the overall behavior of the
flow (Chandrasekhar, 1961). In forced-dissipative simulations, such
perturbations continuously perturb the system, resulting in variations to the
underlying stability. As our quasi-2D model uses the shallow water
approximation, there are no solenoidal contributions to the system that may
excite baroclinic modes, thereby limiting the domain to barotropic
instabilities exclusively. We examine the long term behavior of our
simulations to see if they are Arnol’d- stable.

Figure 8: “Ma”-1 as a function of distance from the pole for all simulations
used in Case B. The dashed black line demarcates “Ma”${}^{-1}=1$.

The Arnol’d stability theorems utilize the pseudoenergy, $H$, which is a
combination of the kinetic energy and enstrophy (squared vorticity) and is
defined as

$\begin{split}H&\equiv\frac{1}{2}\int\left[(\nabla\psi^{\prime})^{2}+\frac{d\Psi}{dQ}(q^{\prime})^{2}\right]\
dA,\end{split}$ (21)

where $\Psi$ and $Q$ represent the geophysical steady-state stream function
(Mcintyre & Shepherd, 1987) and absolute vorticity, respectively.
Correspondingly, $\psi^{\prime}$ and $q^{\prime}$ are their perturbed values
(Vallis, 2006). Stability of the flow requires that the pseudoenergy is
positive-definite or negative-definite, thus the emergence of instabilities
rests on the behavior of $d\Psi/dQ$. Arnol’d’s first stability condition
requires that $d\Psi/dQ>0$, while the second requires that it be sufficiently
negative such that $H$ remains negative-definite across the domain. Dowling
(1993) showed that if $d\Psi/dQ<-L_{d}^{2}$, Arnol’d’s second criterion is
sufficiently satisfied. In our case, we may write

$\begin{split}\frac{d\Psi}{dQ}&=\frac{d\Psi/dr}{dQ/dr}=-v_{\theta}\left(\frac{d\zeta}{dr}+\frac{df}{dr}\right)^{-1},\end{split}$
(22)

where $v_{\theta}$ is the zonal wind. Dowling (2020) showed that the two
Arnol’d-stable branches can be concatenated in to a singular condition using
the inverse of the “vortical” Mach number, “Ma”-1, where “Ma”${}^{-1}<1$
implies an Arnol’d-stable flow. The vortical Mach number is defined as.

$\begin{split}``\rm{Ma}"^{-1}=-L_{d}^{2}\ \frac{dQ}{d\Psi}.\end{split}$ (23)

Andrews (1984) found that Arnol’d-type stability can only occur for zonally
symmetric cases. However, Carnevale & Shepherd (1990) and Mu & Wu (2001)
provided further extensions to the theorem, although the restrictions are non-
trivial (Read et al., 2020). Here, we only include this analysis as a means to
study the stable or unstable behavior of the longterm trends in our
simulations. Limiting the applicability to satisfy the results of Andrews
(1984) means that the analysis can only be applied to high deformation length
cases (simulations B3a, B3b, and B3c). However, for completion, we provide the
results for all Case B simulations to highlight the distinct profiles of
“Ma”-1.

We use the steady-state velocity and vorticity fields over the last 5000 days.
This results in azimuthally symmetric fields for most of the Case B
simulations. We then take an azimuthal average to produce a radial profile for
“Ma”-1. Since we are interested in the stability of the emergent structures
that are affected by the choice of deformation length and not the
perturbations at the grid scale, we convolve the profile with a Gaussian
kernel with a full width half maximum (FWHM) set to the corresponding $L_{d}$
for each simulation.

Figure 8 shows the profile of “Ma”-1 as a function of distance from the pole,
$r$. By day 20,000, all dynamics have equilibrated (see Figure 6). Simulations
B3a, B3b, and B3c show that the polar cyclones are Arnol’d-stable out to at
least $r\sim 0.3$, beyond which the “Ma”-1 profiles cross the “Ma”${}^{-1}=1$
transition. The profiles show stability where the polar cyclones dominate.
There is a noticeable dependence on the forcing strength, which is clear in
simulations B3a, B3b, and B3c. As the forcing strength is increased, the
crossing point is pushed equatorward from $r\sim 0.3$ (B3a) to $r\sim 0.4$
(B3c). This is because as the forcing strength increases, the polar cyclone
gets stronger (Figure 5) resulting in a larger domain of stability. The
interpretations of other simulations in Case B are difficult with regards to
Arnol’d-stability. However, the distinction between the high and low $L_{d}$
cases is clear. Lower $L_{d}$ simulations show less homogeneous behavior in
their PV fields, which is reflected in their “Ma”-1 profiles.

##  4 Discussion

Our model provides insight into the emergence of FFRs, application of moist
convective theory, and an exploration of the low deformation length regimes in
gas giant atmospheres. Here we discuss the primary results from our
simulations in context of current modeling efforts.

###  4.1 Emergence of FFRs

Our simulations show the emergence of FFRs at high $R_{\rm storm}$ values,
indicating that a larger turbulent forcing scale, which may inject more
localized available potential energy, is responsible for higher values of
horizontal shear (see Figure 4). Even though not all of the injected energy is
able to convert to kinetic energy, a large turbulent forcing scale injects
more geopotential in a region at the same forcing strength, $s_{\rm max}$.
Marcus (2004) showed that twisted filamentary structures are consistent with
long-lived cyclonic behavior. We observe that the high horizontal shear in our
simulations is exclusively localized to the cyclonic regions of the domain, in
agreement with their result.

Figure 9: Left: JunoCam RGB processed image of the South Polar Region at PJ20
(Rogers et al., 2022, image credit: NASA/SwRI/MISSS/Gerald Eichstädt/John
Rogers). The FFRs (demarcated in yellow) are offset from the pole and
constitute a belt of large scale FFRs beyond about $85^{\circ}$S. Right:
Simulation A5e from Figure 3. The FFR spans a large area in the domain and is
also offset from the pole, similar to what is seen by JunoCam.

The FFR behavior tends to occur consistently above a specific forcing scale in
our simulations. Above $R_{\rm storm}$$=1200$ km, horizontal shear begins to
dominate over the divergence by about an order of magnitude (see Table 3).
Thus, the bulk of the available potential energy is converted to shear as the
fluid element elongates and folds over. This ultimately results in extremely
high localized velocities, e.g. simulation A6f has $v_{\rm max}=189$ m s-1 and
is well beyond any observed wind speed at the polar regions, however, the
domain-averaged velocity magnitudes remain below $\sim 35$ m s-1, in agreement
with observations (Orton et al., 2017). Furthermore, some high latitude FFRs
have been identified as sources of lightning, suggesting a moist convective
origin (Borucki & Magalhaes, 1992; Fletcher et al., 2017). The formation of
FFRs above a certain forcing scale suggests that large scale moist convective
events might be the dynamical source of amorphous, filamentary cyclonicity at
the poles as they deposit more APE to the upper atmosphere that can be
converted to horizontal shear. In Figure 9, we show a comparison of a JunoCam
image of the South polar region during Perijove 20 (PJ20) with one of the
simulations that contain a large scale FFR-type feature. Qualitatively, the
modeled FFR is larger than any of the observed FFRs in the JunoCam image.
However, the simulation exhibits a similar structure in its morphology, and
position relative to the pole. The belt of FFRs $(70^{\circ}-85^{\circ}$ S)
hosts an assembly of FFRs that are similar in morphology to the large scale
FFR behavior we see in our simulations.

###  4.2 Barotropic Modes and Moist Convection

Our model shows that exclusively employing barotropic modes in the simulations
is sufficient to produce the dynamical variety of the polar atmosphere,
consistent with previous works (O’Neill et al., 2015, 2016; Brueshaber et al.,
2019; Brueshaber & Sayanagi, 2021). The consistency of these results provides
strong support to the idea that moist convective behavior of the deep
atmosphere functions as an important source of observed barotropic vortices of
the higher latitudes. Our findings lend support to the hypothesis put forth by
O’Neill et al. (2015) that the planetary Burger number functions as the
primary non-dimensional parameter that modulates the long term evolutionary
behavior of polar cyclones. It lends further credence to the findings of
BSD19, i.e. that exploration of the parameter space by varying planetary Bu,
injection strengths, and cyclone-to-anticyclone ratio provides dynamical
variation that may be relevant to the atmospheres of gas and ice giants.
Indeed, this distinction is further explored in Brueshaber & Sayanagi (2021),
which focused exclusively on the overall wind speeds and evolutionary
structure of polar cyclones that may elucidate the primary contrast between
Saturn-like (“S”) and Ice Giant-like (“I”) regimes in giant planets.

In Figure 2, the low deformation length case (simulation A1a) exhibits
behavior similar to the Jupiter-like (“J”) systems and begins to exhibit more
transitional (“T”) morphology as $L_{d}$ is increased. The intermediate
$L_{d}$ simulations in Figure 2 show more of a homogenized distribution of
cyclonic PV, however, no dominant polar cyclone emerges until an $L_{d}$ of
$1500$ km is reached, where the rotational frequency becomes lower than
$\Omega_{J}$. However, this is also dependent upon the turbulent forcing scale
of the model. For simulations with intermediate $L_{d}$, increasing the
injection scale results in the formation of a much more dominant polar
cyclone, which is reminiscent of high $L_{d}$ cases. Thus, there is a non-
linear relation between the interplay of $R_{\rm storm}$ and $L_{d}$, which
directly supports the findings of O’Neill et al. (2015). Both parameters
modulate the flow via conversion of the injected APE into localized winds, and
therefore KE, as well as imparting strong gravity waves. The choice of the
local Burger number, Bul, therefore, modulates how much of the injected
potential energy undergoes geostrophic adjustment and how the resulting KE
interacts with the ambient flow. The model, however, is limited to barotropic
modes as perturbations are made solely made to single layer geopotential
field. This sufficiently captures the effect of mass and energy changes due to
vertical mass transport via moist convective events (Sánchez-Lavega et al.,
2011; García-Melendo et al., 2013; García-Melendo & Sánchez-Lavega, 2017), but
is not sufficient to reveal deep atmospheric behavior of the moist convective
plumes themselves. A baroclinic anomaly, as used by O’Neill et al. (2015) to
simulate moist convective events, would result in vertical shearing of the
hetons, leaving behind a strongly barotropic vertical structure (Skinner &
Cho, 2021). However, for a comprehensive approach to moist convection of
stratified baroclinic anomalies, 3D models are necessary to constrain the deep
atmospheric behavior and its effect on the upper troposphere (see Garcia et
al., 2020, for a full 3D treatment of deep atmospheric flows).

###  4.3 Completion of Parameter Space - Exploring low planetary Burger number

This work introduces an accurate and high resolution model that is necessary
to simulate the high latitude dynamics on gas giants, particularly Jupiter.
Providing an extension of the parameter space setup in BSD19 and Brueshaber &
Sayanagi (2021), our model is able to resolve the lower end of the deformation
length parameter that functions as the primary variant that modulates the
planetary Burger number. Although BSD19 explores similar Bu regimes, their
model is restricted to extremely large planetary scales. The benefit of
addressing the case for Jupiter explicitly supports the applicability of the
moist convective theory to the Jovian high latitudes, and due to computational
capabilities afforded by the Pencil framework, allows for a much richer
picture of the shallow dynamics that may be driven by moist convective events
occurring in the deeper atmosphere.

We limited our work to a resolution of $1024^{2}$ in order to produce the main
aspects of the parameter space that is explored by BSD19. Furthermore,
Brueshaber & Sayanagi (2021) expand upon the details of the “S” and “I”
dynamical regimes to furnish the details of polar cyclone behavior, which
correspond to Bu $\sim 1.6\times 10^{-3}$ and Bu $\sim 10^{-2}$, respectively.
Here, we provided a necessary investigation of the Jovian “J” dynamical regime
$({\rm Bu}\sim 10^{-5})$. Our findings support the idea that barotropic modes
are sufficient to capture the dynamical variety observed on Jupiter’s polar
regions, even though the simulations are not able to produce the vortex
configurations (Adriani et al., 2018). The vortex configurations require a
balance between the meridional vorticity gradient generated by the polar
cyclone and the planetary vortcity gradient. This balance results in a region
of stability that is able to host multiple circumpolar cyclones, as long as
they remain at a distance of a few $L_{d}$ from each other in order to
minimize disruption (Gavriel & Kaspi, 2021). Our simulations successfully show
the emergence polar cyclones, even in the low $L_{d}$ cases, but are unable to
produce similar sized circumpolar structures, suggesting that the balance
between the polar cyclone vorticity gradient and the planetary vorticity
gradient is easily disrupted by continuous forcing. Using a freely decaying
model, Li et al. (2020) showed that once the cyclones are stabilized against
the planetary vorticity gradient, initialized vortex configurations evolve in
equilibrium. In a forced turbulence model such as ours, meridional vorticity
gradient balance is continuously disrupted and no configuration is able to
form.

###  4.4 Shallow Water Applications of the Pencil Code

This work applies the computational framework of Pencil to the atmosphere of
gas giant planets for the first time. The code is open source and allows for
users to model a wide variety of astrophysically relevant dynamics and
thermochemistry in with modern computational techniques. This allows for
increased time savings due to significantly faster model integration.
Specifically, the code allows for full $20,000$ days of evolution for the
lowest perturbation size and lowest deformation length simulation, A1a (most
computationally expensive simulation due to relevant size scales), within a
wall-clock time of $45$ hours. Our model also requires the turbulent forcing
scale to be resolved fully by the stencil, a condition that required far more
wall-clock time in previous works (e.g. Brueshaber & Sayanagi (2021)).

Our model produces results similar to those of previous works and further adds
to the increasing body of knowledge that motivates 3D modeling of the
atmosphere in the future (Garcia et al., 2020). We are able to show distinct
and realistic morphologies that emerge due to a consequence of our
parameterizations using the Pencil Code. Although we maintain the resolution
at $1024^{2}$, we test the code for convergence up to $2048^{2}$. Thus, the
parameter space can be expanded further with the inclusion of higher
resolution simulations and deformation lengths as low as $500$ km.

##  5 Conclusions

In this work, we showcase a novel application of a well-tested hydrodynamical
code to the atmosphere of gas giant planets capable of addressing wide
parameter space searches by utilizing the framework provided by PencilCode.
Applying this atmospheric module to the case of Jupiter’s polar regions using
high resolution, we find high latitude behavior reminiscent of FFRs. The
simulations show that these FFRs are regions of intense horizontal shear and
are generated by storms above a length scale of $\sim 1200$ km. We further
explore an immense parameter space defined by O’Neill et al. (2015), BSD19,
and Brueshaber & Sayanagi (2021), and produce results that support the idea
that moist convection functions as a dominant source of the dynamical variety
that is seen at the Jovian poles.

We find that the deformation length, $L_{d}$, functions as a primary indicator
of long term polar behavior, and specifically that increasing the turbulent
forcing scale for large $L_{d}$ simulations result in increased horizontal
shear. Our findings suggest that large scale moist convection may be the
dominant cause of the formation of high latitude FFRs, which is in support of
previous observations (e.g. Borucki & Magalhaes, 1992). Furthermore, we find
that high $L_{d}$ simulations show Arnol’d-stable behavior using the vortical
Mach number condition, “Ma”${}^{-1}=1$ (Dowling, 2020). Higher deformation
lengths are able to stabilize the domain via homogenization of local potential
vorticity in the form of polar cyclone formation. Low $L_{d}$ simulations are
difficult to interpret with regards to Arnol’d-stability, but show a clear
distinction from the high $L_{d}$ simulations.

We also note that the emergence of patterns with different dominant
wavenumbers in our simulations depends primarily on the size (relative to that
of our domain) of the largest balanced vortex that can emerge in each
simulation. Increasing the deformation length (decreasing planetary Bu)
reduces the number of balanced eddies that stabilize at the polar regions. The
wavenumber that emerges is dependent on number of the balanced eddies, as well
as how they self-organize within the domain. Morphologically, this behavior is
reminiscent of polar vortex Rossby waves (Chen & Yau, 2001; Houghton, 2002)
and should be studied further.

Our work can be expanded upon with the inclusion of stratification. The model
can be made multilayered to include baroclinicity and thermodynamical effects.
The parameter space may be reduced to study low deformation length cases at
significantly higher resolution while forcing the domain using baroclinic
modes. As multilayered shallow water systems crudely approximate
stratification, allowing baroclinic modes in a multilayered system will allow
a better quantification of the effect of anisotropic turbulence on the upper
atmosphere and how it modulates the long term behavior of these flows. Further
modifications may include variations in the underlying dynamical topography.
For simplicity, we use no variation in the underlying geopotential for a given
simulation. Added complexities may lead to further dynamically interesting
behavior as it may impede meridional eddy momentum flux transfer over long
timescales.

We thank two anonymous referees for their useful comments that helped improve
our manuscript. We also thank Dr. John Rogers and Gerald Eichstädt for the
contribution of the JunoCam image. A.H. also thanks Dr. Shawn Brueshaber and
Dr. Morgan O’Neill for important discussions regarding gas giant polar
modeling and providing important insight into the parameter space. This work
was supported by NASA’s New Frontiers Data Analysis Program (grant#
80NSSC20K0561) and by the NASA Fellowship Activity (grant# 80NSSC20K1456). The
authors acknowledge the Texas Advanced Computing Center (TACC) at The
University of Texas at Austin for providing HPC resources that have
contributed to the research results reported within this paper
(http://www.tacc.utexas.edu).

##  Appendix A Scaling Study

As we are presenting a new shallow water atmospheric module in Pencil, it is
important to assess the efficiency of our code using a scaling analysis that
characterizes the effectiveness of parallelization. We employed the
XSEDE/Stampede2 supercomputer at Texas Advanced Computing Center (TACC) for
our simulations. We perform a strong scaling analysis to determine how the
code handles large scale communication between active nodes on the cluster. We
fix the resolution at $1024^{2}$ while keeping all of the parameters for a
full simulation active, ensuring that the timestep average would be
representative of simulations presented in this paper. In Figure 10, we see
that a processor load of $32^{2}$ grid points minimizes total computation and
communication most effectively. This load signifies the turnover point where
communication time is already causing computational inefficiency, however, we
elect to use it as the increase in time from processor loads of $64^{2}$ or
$32\times 64$ is still important for simulations with long term evolution.

Figure 10: Strong scaling study using Stampede2 at a base resolution of
$1024^{2}$. With our fixed resolution, $32^{2}$ optimizes the balance between
computation and communication time.

The efficiency of the parallelization in Pencil using weak scaling is shown in
Figure 11. As expected, maintaining the optimal processor load, $32^{2}$, and
increasing the resolution via increasing total computing nodes used yields a
highly linear relation, indicative of efficient parallelization. Tests were
conducted up to resolutions of $2048^{2}$; at a resolution of $1024^{2}$, the
$32^{2}$ processor load optimizes the simulations over $64^{2}$. Thus, our
simulations can be used for studying the polar dynamics with $\times 16$ more
resolution than previous simulations if a resolution of $2048^{2}$ is
employed.

Figure 11: Weak scaling study using Stampede2 with varying processing loads.
The red dash-dotted line is for a processor load of $32^{2}$ while the cyan is
for $64^{2}$.

##  Appendix B Benchmark

Since our model uses approximations and methodologies explored in previous
works, we opt to test the validity of our model using a phenomenological
benchmark. We use the vortex solution developed by DeMaria & Chan (1984). This
solution has been utilized in finite-volume methods to address the stability
of the Jovian polar vortex configurations and merger events (Li et al., 2020).
Here, we use it as a measure of the stability of our code. The radial velocity
profile, $v(r)$, for the vortex can be written as:

$\begin{split}v(r)=v_{m}\left(\frac{r}{r_{m}}\right)\exp{\left\\{\frac{1}{b}\left[1-\left(\frac{r}{r_{m}}\right)^{b}\right]\right\\}},\end{split}$
(B1)

where $v_{m}$ is the maximum tangential velocity, $r$ is the distance from the
center of the vortex, $r_{m}$ is the radius of maximum wind speed, and $b$ is
a shielding parameter, which determines the strength of the surrounding
(opposite-signed) vorticity. Equation B1 is a numerical approximation of the
behavior of a tropical cyclone interacting with another cyclone in its
vicinity. A full conceptual argument of this interaction is provided in
(Chang, 1983) and (DeMaria & Chan, 1984). Likewise, the corresponding
geopotential profile is calculated numerically. This results in a minor phase
of adjustment during which the geopotential and the velocity profile must
relax into their equilibrium values. The adjustment phase only lasts for a
brief period and imparts gravity waves until the velocity field and
geopotential appropriately converge.

Figure 12: Behavior of diagnostic quantities of interest given the initial
vortex solution of Li et al. (2020); DeMaria & Chan (1984). The values
stabilize by about $200$ Earth days and show sufficient convergence from
resolutions $1024^{2}$ (dashed blue) to $2048^{2}$ (solid red).

In Figure 12, we show the early evolution of the vortex solution from (DeMaria
& Chan, 1984). We offset the vortex from the pole by a small amount to
highlight any azimuthal asymmetries that may exist in the simulation. The
vortex stabilizes at about $200$ Earth days and varies very slowly and no
azimuthal asymmetries are seen. To test convergence, the model is run at
resolutions $1024^{2}$ and $2048^{2}$. All diagnostics show code convergence
as the differences are only on the order of $10^{-6}$. The maximum relative
error is exhibited in the mean-square divergence. However, as we employ
incompressibility of the flow in question, the overall effect is on the order
of $10^{-10}$. Both KE (in the form of $\langle u^{2}\rangle$) and APE tend to
oscillate and converge in the same manner for both resolutions, and are
sufficiently stable by day $200$. The adjustment phase is dependent upon the
numerical integration scheme used for the initial geopotential as larger
errors will impart stronger gravity waves as the flow adjusts to the
initialized velocity profile given by equation B1. However, low-order
numerical integration only effects low resolutions and is not a concern for
the resolutions used in this work. Thus, employing $1024^{2}$ for our
simulations is sufficient to capture relevant dynamical behavior of a
$1\frac{1}{2}$-layer shallow water system.

## References

  * Adriani et al. (2017) Adriani, A., Filacchione, G., Di Iorio, T., et al. 2017, Space Science Reviews, 213, 393, doi: 10.1007/s11214-014-0094-y
  * Adriani et al. (2018) Adriani, A., Mura, A., Orton, G., et al. 2018, Nature, 555, 216, doi: 10.1038/nature25491
  * Andrews (1984) Andrews, D. G. 1984, Geophysical and Astrophysical Fluid Dynamics, 28, 243, doi: 10.1080/03091928408230366
  * Bagenal et al. (2007) Bagenal, F., Dowling, T. E., & McKinnon, W. B. 2007, Jupiter
  * Bardet et al. (2019) Bardet, D., Spiga, A., Guerlet, S., et al. 2019, in EGU General Assembly Conference Abstracts, EGU General Assembly Conference Abstracts, 14516
  * Boissinot et al. (2019) Boissinot, A., Spiga, A., Guerlet, S., Cabanes, S., & Hourdin, F. 2019, in EGU General Assembly Conference Abstracts, EGU General Assembly Conference Abstracts, 15821
  * Borucki & Magalhaes (1992) Borucki, W. J., & Magalhaes, J. A. 1992, Icarus, 96, 1, doi: 10.1016/0019-1035(92)90002-O
  * Bourdin et al. (2013) Bourdin, P. A., Bingert, S., & Peter, H. 2013, A&A, 555, A123, doi: 10.1051/0004-6361/201321185
  * Brandenburg & Dobler (2002) Brandenburg, A., & Dobler, W. 2002, Computer Physics Communications, 147, 471, doi: 10.1016/S0010-4655(02)00334-X
  * Brandenburg & Dobler (2010) —. 2010, Pencil: Finite-difference Code for Compressible Hydrodynamic Flows, Astrophysics Source Code Library.  http://ascl.net/1010.060
  * Brandenburg & Scannapieco (2020) Brandenburg, A., & Scannapieco, E. 2020, ApJ, 889, 55, doi: 10.3847/1538-4357/ab5e7f
  * Bridger & Stevens (1980) Bridger, A. F. C., & Stevens, D. E. 1980, Journal of Atmospheric Sciences, 37, 534 , doi: 10.1175/1520-0469(1980)037<0534:LAWATP>2.0.CO;2
  * Brueshaber & Sayanagi (2021) Brueshaber, S. R., & Sayanagi, K. M. 2021, Icarus, 361, 114386, doi: 10.1016/j.icarus.2021.114386
  * Brueshaber et al. (2019) Brueshaber, S. R., Sayanagi, K. M., & Dowling, T. E. 2019, Icarus, 323, 46, doi: 10.1016/j.icarus.2019.02.001
  * Cai et al. (2021) Cai, T., Chan, K. L., & Mayr, H. G. 2021, The Planetary Science Journal, 2
  * Carnevale & Shepherd (1990) Carnevale, G. F., & Shepherd, T. G. 1990, Geophysical and Astrophysical Fluid Dynamics, 51, 1, doi: 10.1080/03091929008219847
  * Chandrasekhar (1961) Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability
  * Chang (1983) Chang, S. W.-J. 1983, Monthly Weather Review, 111, 1806 , doi: 10.1175/1520-0493(1983)111<1806:ANSOTI>2.0.CO;2
  * Chen & Yau (2001) Chen, Y., & Yau, M. K. 2001, Journal of the Atmospheric Sciences, 58, 2128 , doi: 10.1175/1520-0469(2001)058<2128:SBIASH>2.0.CO;2
  * DeMaria & Chan (1984) DeMaria, M., & Chan, J. C. L. 1984, Monthly Weather Review, 112, 1643 , doi: 10.1175/1520-0493(1984)112<1643:CONSOT>2.0.CO;2
  * Dowling (1993) Dowling, T. E. 1993, Journal of Atmospheric Sciences, 50, 14 , doi: 10.1175/1520-0469(1993)050<0014:ARBPVA>2.0.CO;2
  * Dowling (2020) —. 2020, The Planetary Science Journal, 1, 6, doi: 10.3847/PSJ/ab789d
  * Dowling et al. (1998) Dowling, T. E., Fischer, A. S., Gierasch, P. J., et al. 1998, Icarus, 132, 221, doi: 10.1006/icar.1998.5917
  * Dowling & Ingersoll (1989) Dowling, T. E., & Ingersoll, A. P. 1989, Journal of Atmospheric Sciences, 46, 3256 , doi: 10.1175/1520-0469(1989)046<3256:JGRSAA>2.0.CO;2
  * Fletcher et al. (2017) Fletcher, L. N., Orton, G. S., Rogers, J. H., et al. 2017, Icarus, 286, 94, doi: 10.1016/j.icarus.2017.01.001
  * Fodor (2003) Fodor, F. 2003, Beiträge zur Algebra und Geometrie, 44
  * Garcia et al. (2020) Garcia, F., Chambers, F. R. N., & Watts, A. L. 2020, MNRAS, 499, 4698, doi: 10.1093/mnras/staa2962
  * García-Melendo et al. (2013) García-Melendo, E., Hueso, R., Sánchez-Lavega, A., et al. 2013, Nature Geoscience, 6, 525, doi: 10.1038/ngeo1860
  * García-Melendo & Sánchez-Lavega (2017) García-Melendo, E., & Sánchez-Lavega, A. 2017, Icarus, 286, 241, doi: https://doi.org/10.1016/j.icarus.2016.10.006
  * Gavriel & Kaspi (2021) Gavriel, N., & Kaspi, Y. 2021, Nature Geoscience, 14, 559, doi: 10.1038/s41561-021-00781-6
  * Godfrey (1988) Godfrey, D. A. 1988, Icarus, 76, 335, doi: 10.1016/0019-1035(88)90075-9
  * Golub & Van Loan (1996) Golub, G. H., & Van Loan, C. F. 1996, Matrix Computations (3rd Ed.) (USA: Johns Hopkins University Press)
  * Grassi et al. (2018) Grassi, D., Adriani, A., Moriconi, M. L., et al. 2018, Journal of Geophysical Research (Planets), 123, 1511, doi: 10.1029/2018JE005555
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Haugen et al. (2004a) Haugen, N. E., Brandenburg, A., & Dobler, W. 2004a, Phys. Rev. E, 70, 016308, doi: 10.1103/PhysRevE.70.016308
  * Haugen et al. (2004b) Haugen, N. E. L., Brandenburg, A., & Mee, A. J. 2004b, MNRAS, 353, 947, doi: 10.1111/j.1365-2966.2004.08127.x
  * Hogg & Stommel (1985) Hogg, N. G., & Stommel, H. M. 1985, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 397, 1, doi: 10.1098/rspa.1985.0001
  * Hord et al. (2017) Hord, B., Lyra, W., Flock, M., Turner, N. J., & Mac Low, M.-M. 2017, ApJ, 849, 164, doi: 10.3847/1538-4357/aa8fcf
  * Houghton (2002) Houghton, J. 2002, The Physics of Atmospheres
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hyder et al. (2022) Hyder, A., Lyra, W., Chanover, N., Morales-Juberías, R., & Jackiewicz, J. 2022, Exploring Jupiter’s Polar Deformation Lengths with High Resolution Shallow Water Modeling, Zenodo, doi: 10.5281/zenodo.6642986
  * Jin & Dubin (1998) Jin, D. Z., & Dubin, D. H. E. 1998, Phys. Rev. Lett., 80, 4434, doi: 10.1103/PhysRevLett.80.4434
  * Jones et al. (2001) Jones, E., Oliphant, T., & Peterson, P. 2001
  * Kaspi (2008) Kaspi, Y. 2008, PhD thesis, Massachusetts Institute of Technology
  * Kaspi et al. (2009) Kaspi, Y., Flierl, G. R., & Showman, A. P. 2009, Icarus, 202, 525, doi: 10.1016/j.icarus.2009.03.026
  * Klemp et al. (2018) Klemp, J. B., Skamarock, W. C., & Ha, S. 2018, Monthly Weather Review, 146, 1911 , doi: 10.1175/MWR-D-17-0384.1
  * Li et al. (2020) Li, C., Ingersoll, Andrew, P., Klipfel, Alexandra, P., & Brettle, H. 2020, Proceedings of the National Academy of Science, 117, 24082, doi: 10.1073/pnas.2008440117
  * Lyra et al. (2008) Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008, A&A, 479, 883, doi: 10.1051/0004-6361:20077948
  * Lyra et al. (2017) Lyra, W., McNally, C. P., Heinemann, T., & Masset, F. 2017, AJ, 154, 146, doi: 10.3847/1538-3881/aa8811
  * Lyra et al. (2018) Lyra, W., Raettig, N., & Klahr, H. 2018, Research Notes of the American Astronomical Society, 2, 195, doi: 10.3847/2515-5172/aaeac9
  * Lyra et al. (2016) Lyra, W., Richert, A. J. W., Boley, A., et al. 2016, ApJ, 817, 102, doi: 10.3847/0004-637X/817/2/102
  * Lyra et al. (2015) Lyra, W., Turner, N. J., & McNally, C. P. 2015, A&A, 574, A10, doi: 10.1051/0004-6361/201424919
  * Marcus (2004) Marcus, P. S. 2004, Nature, 428, 828, doi: 10.1038/nature02470
  * Marshall (2008) Marshall, J. 2008, Atmosphere, ocean, and climate dynamics an introductory text, International geophysics series ; v. 93 (Amsterdam ; [Burlington, MA]: Elsevier Academic Press)
  * Marshall et al. (1997) Marshall, J., Adcroft, A., Hill, C., Perelman, L., & Heisey, C. 1997, Journal of Geophysical Research: Oceans, 102, 5753, doi: https://doi.org/10.1029/96JC02775
  * Mcintyre & Shepherd (1987) Mcintyre, M., & Shepherd, T. 1987, Journal of Fluid Mechanics, 181, 527 , doi: 10.1017/S0022112087002209
  * Montgomery et al. (2006) Montgomery, M. T., Nicholls, M. E., Cram, T. A., & Saunders, A. B. 2006, Journal of Atmospheric Sciences, 63, 355, doi: 10.1175/JAS3604.1
  * Morales-Juberías et al. (2011) Morales-Juberías, R., Sayanagi, K. M., Dowling, T. E., & Ingersoll, A. P. 2011, Icarus, 211, 1284, doi: 10.1016/j.icarus.2010.11.006
  * Morales-Juberías et al. (2015) Morales-Juberías, R., Sayanagi, K. M., Simon, A. A., Fletcher, L. N., & Cosentino, R. G. 2015, ApJ, 806, L18, doi: 10.1088/2041-8205/806/1/L18
  * Mu & Wu (2001) Mu, M., & Wu, Y.-H. 2001, Surveys in Geophysics, 22, 383, doi: 10.1023/A:1014229917728
  * O’Neill et al. (2015) O’Neill, M. E., Emanuel, K. A., & Flierl, G. R. 2015, Nature Geoscience, 8, 523, doi: 10.1038/ngeo2459
  * Orton et al. (2017) Orton, G. S., Hansen, C., Caplinger, M., et al. 2017, Geophys. Res. Lett., 44, 4599, doi: 10.1002/2016GL072443
  * O’Neill et al. (2016) O’Neill, M. E., Emanuel, K. A., & Flierl, G. R. 2016, Journal of the Atmospheric Sciences, 73, 1841 , doi: 10.1175/JAS-D-15-0314.1
  * Pedlosky (1979) Pedlosky, J. 1979, Geophysical fluid dynamics / Joseph Pedlosky (Springer Verlag New York), xii, 624 p. :
  * Pedlosky (1985) Pedlosky, J. 1985, Journal of Atmospheric Sciences, 42, 1477 , doi: 10.1175/1520-0469(1985)042<1477:TIOCHC>2.0.CO;2
  * Pirl (1969) Pirl, U. 1969, Mathematische Nachrichten, 40, 111
  * Raettig et al. (2015) Raettig, N., Klahr, H., & Lyra, W. 2015, ApJ, 804, 35, doi: 10.1088/0004-637X/804/1/35
  * Raettig et al. (2013) Raettig, N., Lyra, W., & Klahr, H. 2013, ApJ, 765, 115, doi: 10.1088/0004-637X/765/2/115
  * Read et al. (2020) Read, P., Kennedy, D., Lewis, N., et al. 2020, Nonlinear Processes in Geophysics, 27, 147, doi: 10.5194/npg-27-147-2020
  * Read et al. (2006) Read, P. L., Gierasch, P. J., Conrath, B. J., et al. 2006, Quarterly Journal of the Royal Meteorological Society, 132, 1577, doi: 10.1256/qj.05.34
  * Richert et al. (2015) Richert, A. J. W., Lyra, W., Boley, A., Mac Low, M.-M., & Turner, N. 2015, ApJ, 804, 95, doi: 10.1088/0004-637X/804/2/95
  * Rogers et al. (2022) Rogers, J., Eichstädt, G., Hansen, C., et al. 2022, Icarus, 372, 114742, doi: https://doi.org/10.1016/j.icarus.2021.114742
  * Sánchez-Lavega et al. (2011) Sánchez-Lavega, A., del Río-Gaztelurrutia, T., Hueso, R., et al. 2011, Nature, 475, 71, doi: 10.1038/nature10203
  * Scott (2011) Scott, R. 2011, Geophysical & Astrophysical Fluid Dynamics, 105, 409, doi: 10.1080/03091929.2010.509927
  * Showman (2007) Showman, A. P. 2007, Journal of Atmospheric Sciences, 64, 3132, doi: 10.1175/JAS4007.1
  * Skamarock & Klemp (1992) Skamarock, W. C., & Klemp, J. B. 1992, Monthly Weather Review, 120, 2109 , doi: 10.1175/1520-0493(1992)120<2109:TSOTSN>2.0.CO;2
  * Skinner & Cho (2021) Skinner, J. W., & Cho, J. Y. K. 2021, MNRAS, doi: 10.1093/mnras/stab2809
  * Tabataba-Vakili et al. (2020) Tabataba-Vakili, F., Rogers, J. H., Eichstädt, G., et al. 2020, Icarus, 335, 113405, doi: 10.1016/j.icarus.2019.113405
  * Vallis (2006) Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics, doi: 10.2277/0521849691
  * Vallis et al. (2018) Vallis, G. K., Colyer, G., Geen, R., et al. 2018, Geoscientific Model Development, 11, 843, doi: 10.5194/gmd-11-843-2018
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Williams (1978) Williams, G. P. 1978, Journal of Atmospheric Sciences, 35, 1399, doi: 10.1175/1520-0469(1978)035<1399:PCBROJ>2.0.CO;2
  * Williams (1988) —. 1988, Climate Dynamics, 3, 45, doi: 10.1007/BF01080901
  * Yano & Flierl (1992) Yano, J.-I., & Flierl, G. R. 1992, Dynamics of Atmospheres and Oceans, 16, 439, doi: https://doi.org/10.1016/0377-0265(92)90001-A