Contents lists available at
ScienceDirect
Icarus
journal homepage:
www.elsevier.com/locate/icarus
Research Paper
Vortex crystals at Jupiter’s poles: Emergence controlled by initial small-scale
turbulence
Sihe Chen
a
,
∗
,
Andrew P. Ingersoll
a
,
Cheng Li
b
a
Planetary Science Department, California Institute of Technology, 1200 E California Blvd, Pasadena, 91125, CA, USA
b
Department of Climate and Space Sciences and Engineering, University of Michigan, 2455 Hayward St, Ann Arbor, 48109, MI, USA
A R T I C L E I N F O
Dataset link:
https://github.com/chengcli/cano
e
Keywords:
Jupiter
Vortex crystal
Shallow-water system
A B S T R A C T
At the poles of Jupiter, cyclonic vortices are clustered together in patterns made up of equilateral triangles
called vortex crystals. Such patterns are seen in laboratory flows but never before in a planetary atmosphere,
where the planet’s rotation and gravity add new physics. Siegelman (2022b) used a one-layer quasi-geostrophic
(QG) model with an infinite radius of deformation to study the emergence of vortex crystals from small-scale
turbulence, and Li (2020) showed that shielding of the vortices is important for the stability of the vortex
crystals. Here we use the shallow water (SW) equations at the pole of a rotating planet to study the emergence
and evolution of vortices starting from an initial random pattern of small-scale turbulence. The flow is in a
single layer with a free surface whose slope produces the horizontal pressure gradient force. With the planet’s
radius and rotation used to define the units, only three input parameters are needed to define the system: the
mean kinetic energy of the initial turbulence, the horizontal scale of the initial turbulence, and the radius of
deformation of the undisturbed fluid layer. We identified a non-dimensional number,
훥ℎ
∕
ℎ
, which is related
to the relative layer thickness variation of the initial turbulence and determines whether the vortex crystal
or chaotic patterns emerge: Small
훥ℎ
∕
ℎ
values lead to vortex crystals, and large
훥ℎ
∕
ℎ
values lead to chaotic
patterns. The value
훥ℎ
∕
ℎ
is related to the radius of deformation as
퐿
−2
푑
. This means that a large polar radius
of deformation is positively correlated to the emergence of vortex crystals, and this implies either a polar
atmosphere enriched with water or deeper roots for the vortices than previously estimated.
1.
Introduction
At
each
pole
of
Jupiter,
a
central
cyclonic
vortex
is
surrounded
by
a
ring
of
cyclonic
vortices
at
a
latitude
of
83
◦
,
corresponding
to
a
circle
8700
km
in
radius
(
Adriani
et
al.
,
2018
).
Since
they
were
discovered
in
2017,
the
polar
vortex
patterns
at
the
poles
of
Jupiter
have
been
remarkably
stable
(
Adriani
et
al.
,
2020
;
Mura
et
al.
,
2021
).
Individual
cyclones
have
radii,
as
defined
by
their
azimuthal
velocity
maxima,
of
1000
km.
In
the
north,
the
ring
has
8
circumpolar
cy-
clones;
in
the
south,
the
ring
has
5
circumpolar
cyclones.
The
south
circumpolar
vortices
drift
westward
with
1
.
5
±
0
.
2
◦
during
each
of
the
spacecraft’s
53-day
orbits,
while
the
north
vortices
do
not
present
significant
drifts
(
Tabataba-Vakili
et
al.
,
2020
).
In
2019,
a
cyclone
joined
the
pattern
in
the
south
but
subsequently
disappeared
without
affecting
the
stability
of
the
pattern
(
Mura
et
al.
,
2021
).
The
central
cyclone
and
the
circumpolar
cyclones
resemble
approximately
a
two-
dimensional
structure
comprised
of
equilateral
triangles.
Such
patterns
have
been
studied
in
the
laboratory
and
are
appropriately
called
‘‘vor-
tex
crystals’’
(
Fine
et
al.
,
1995
;
Schecter
et
al.
,
1999
;
Siegelman
et
al.
,
2022b
).
∗
Corresponding
author.
E-mail
address:
sihechen@caltech.edu
(S.
Chen)
.
Saturn
has
an
eastward-moving
jet
stream
at
approximately
78
◦
N
that
flows
to
form
a
giant
hexagon
centered
on
the
pole.
The
hexagon
was
discovered
by
Godfrey
(
1988
)
in
Voyager
data
taken
in
1981
and
has
endured
to
the
present.
The
hexagon’s
corners
are
changes
in
the
jet’s
direction
but
are
not
closed-streamline
regions
like
the
polar
cyclones
of
Jupiter.
There
were
several
successful
attempts
to
model
the
hexagon
(
Aguiar
et
al.
,
2010
;
Morales-Juberías
et
al.
,
2015
),
which
led
to
a
deeper
understanding
of
the
differences
between
the
dynamics
of
the
polar
atmosphere
and
the
beta
plane
of
mid-latitude.
O’Neill
et
al.
(
2015
,
2016
),
Brueshaber
et
al.
(
2019
),
and
Hyder
et
al.
(
2022
)
used
forced-dissipative
models
to
study
the
Jupiter
vor-
tices.
O’Neill
et
al.
(
2015
)
suggested
that
the
size
ratio
between
the
geostrophic
turbulence
and
planetary
radius
is
important
in
the
genesis
of
a
polar
cyclone.
O’Neill
et
al.
(
2016
)
found
that
the
energy
density
is
the
major
predictor
of
the
polar
vortex
behavior.
Brueshaber
et
al.
(
2019
)
suggested
that
a
small
radius
of
deformation
is
related
to
the
discrete
vortex
structures,
containing
both
cyclones
and
anticyclones,
but
not
the
observed
vortex
crystal
structures.
Meanwhile,
Hyder
et
al.
https://doi.org/10.1016/j.icarus.2024.116438
Received
20
September
2024;
Received
in
revised
form
7
December
2024;
Accepted
11
December
2024
Icarus
429
(2025)
116438
Available
online
18
December
2024
0019-1035/©
2024
The
Authors.
Published
by
Elsevier
Inc.
This
is
an
open
access
article
under
the
CC
BY
license
(
http://creativecommons.org/licenses/by/4.0/
).
S.
Chen
et
al.
(
2022
)
found
that
the
radius
of
deformation
and
the
turbulent
forcing
scale
affect
the
emergence
of
the
folded
filamentary
region
(FFR),
and
a
larger
radius
of
deformation
tends
to
homogenize
the
field
of
PV
at
the
pole.
While
none
of
these
efforts
produced
vortex
crystal
patterns
like
those
seen
on
Jupiter
(
Adriani
et
al.
,
2018
;
Li
et
al.
,
2020b
)
offered
a
single-layer
shallow
water
(SW)
model
to
explain
them.
They
assumed
a
finite
radius
of
deformation,
corresponding
to
a
layer
of
limited
thickness
and
limited
density
contrast
with
the
deep
atmosphere
below.
They
found
that
rings
of
anticyclonic
vorticity
–
shielding
–
around
each
cyclone
were
necessary
to
maintain
the
shape
of
the
polygons
and
keep
the
vortices
from
merging.
They
showed
that
the
pattern
was
stable
to
attach
by
an
intruder
drifting
from
lower
latitudes,
resembling
the
sixth
cyclone
observed
at
the
south
pole
in
2019
(
Mura
et
al.
,
2021
).
Li
et
al.
(
2020b
)
is
an
extension
of
DeMaria
and
Chan
(
1984
),
where
the
effects
of
shielding
on
the
interaction
between
two
vortices
are
discussed.
Subsequently,
Gavriel
and
Kaspi
(
2022
)
used
a
model
in
which
the
vortex
shielding
is
treated
as
the
source
of
a
repulsive
force
to
explain
the
observed
patterns
of
the
polar
vortex
motions.
Siegelman
et
al.
(
2022b
)
assumed
an
infinite
radius
of
deformation
in
a
quasi-
geostrophic
(QG)
model
and
was
the
only
work
that
simulated
the
emergence
of
stable
vortex
polygons.
They
studied
an
initial
condition
problem
of
decaying
turbulence
and
demonstrated
the
spontaneous
formation
of
vortex
crystals
with
minimal
shielding.
However,
the
question
is
still
open:
What
is
the
essential
difference
between
the
results
of
Li
et
al.
(
2020b
)
and
Siegelman
et
al.
(
2022b
)?
Is
it
the
shielding
of
Li
et
al.
(
2020b
),
the
infinite
radius
of
defor-
mation
of
Siegelman
et
al.
(
2022b
),
the
initial
small-scale
turbulence
of
Siegelman
et
al.
(
2022b
),
or
the
approximations
inherent
in
the
QG
model
compared
to
the
SW
model?
Here,
we
use
the
SW
equations
at
the
pole
of
a
rotating
planet
to
study
the
emergence
and
evolution
of
vortices
starting
from
an
initial
random
pattern
of
small-scale
turbulence.
The
flow
is
in
a
single
layer
of
variable
thickness,
ℎ
,
with
a
free
surface
capable
of
propagating
gravity
waves.
We
use
a
reduced-gravity
model
(
Chassignet
and
Cushman-
Roisin
,
1991
),
where
the
active
layer
is
floating
hydrostatically
on
a
much
deeper
fluid
layer
at
rest.
We
use
reduced
gravity
푔
푟
=
푔
훥휌
∕
휌
,
a
constant,
where
푔
is
the
gravitational
acceleration
at
the
planet’s
surface,
and
훥휌
∕
휌
is
the
fractional
density
difference
between
the
two
layers.
The
geopotential,
훷
=
푔
푟
ℎ
,
is
one
of
the
dependent
parameters
of
the
SW
system,
the
other
being
the
velocity
퐮
.
The
upper
boundary
of
the
active
layer
is
a
free
surface
with
zero
pressure,
so
−∇
훷
is
the
horizontal
pressure
gradient
force
per
unit
mass
within
the
layer.
Since
훷
is
proportional
to
the
thickness,
it
obeys
the
horizontal
con-
tinuity
equation.
This,
along
with
the
horizontal
momentum
equation,
constitutes
the
full
set
of
equations
of
the
SW
system:
휕
훷
휕
푡
+
(
퐮
⋅
∇)
훷
+
훷
∇
⋅
퐮
=
0
(1)
휕
퐮
휕
푡
+
(
퐮
⋅
∇)
퐮
+
푓
̂
퐤
×
퐮
=
−∇
훷
(2)
This
system
materially
conserves
the
SW
version
of
the
potential
vor-
ticity
(PV),
푞
=
휁
+
푓
훷
(3)
Here
휁
=
̂
퐤
⋅
(∇
×
퐮
)
is
the
relative
vorticity
of
the
flow.
We
are
modeling
a
thin
horizontal
layer
on
a
rotating
sphere,
so
only
the
vertical
component
of
the
Coriolis
parameter
matters.
Thus,
̂
퐤
is
the
vertical
unit
vector,
and
푓
=
2
훺
sin
휙
is
the
Coriolis
parameter,
where
훺
is
the
planetary
rotation
rate,
and
휙
is
latitude.
Although
we
did
not
explicitly
include
viscosity,
the
numerical
scheme
inherently
introduces
artificial
viscosity,
which
stabilizes
the
flow
but
can
dissipate
energy
at
smaller
scales.
In
this
system,
except
for
the
planet
radius
푎
and
rotation
rate
훺
,
there
are
only
3
independent
parameters
—
the
root
mean
square
velocity
푈
of
the
initial
turbulence,
the
horizontal
wavelength
퐿
푖
of
the
initial
turbulence,
and
the
radius
of
deformation
퐿
푑
=
√
̄
훷
∕2
훺
evalu-
ated
at
the
pole,
where
̄
훷
is
the
horizontal
mean
of
훷
,
and
√
̄
훷
is
the
mean
speed
of
horizontally
propagating
gravity
waves.
Since
훷
is
the
product
of
the
thickness
of
the
upper
layer
and
the
density
difference
between
the
lower
and
upper
layers,
it
is
where
the
vertical
structure
enters
the
SW
equations.
The
planetary
radius
푎
is
71,492
km,
and
the
rotation
period
is
9
ℎ
55
푚
29
.
71
푠
.
These
two
numbers
do
not
change
and
serve
to
define
the
units
of
length
and
time,
respectively.
Our
results
hold
for
other
planets
when
the
same
values
of
dimensionless
numbers
푈
∕(2
훺
푎
)
,
퐿
푖
∕
푎
,
and
퐿
푑
∕
푎
are
used.
Other
approaches
vary
the
planet
size
(
O’Neill
et
al.
,
2015
;
Brueshaber
et
al.
,
2019
),
which
is
another
way
to
change
the
three
basic
dimensionless
numbers.
We
acknowledge
that
this
single-layer
SW
model
is
not
realistic
enough
to
capture
all
the
important
parameters
of
Jupiter’s
atmo-
sphere:
the
nature
of
small-scale
turbulence,
e.g.
whether
it
arises
from
convection,
baroclinic
instability,
breaking
gravity
waves,
or
breaking
planetary
waves;
the
thickness
and
vertical
structure
of
the
layer
where
the
cyclones
exist;
the
nature
of
dissipation
and
whether
there
are
upscale
and
downscale
cascades
of
energy,
etc.
However,
we
cannot
achieve
perfect
realism.
There
are
few
observational
data
for
the
poles,
especially
about
the
vertical
structure.
Current
knowledge
is
related
to
small-scale
turbulence
and
inverse
cascades,
at
least
near
cloud-top
altitudes
(
Ingersoll
et
al.
,
2022
;
Siegelman
et
al.
,
2022a
).
In
the
future,
with
more
information
from
the
Juno
mission
about
the
roots
of
the
vortices
and
stable
layers
at
depth,
a
better
constraint
for
the
model
assumptions
can
be
obtained.
This
paper
aims
to
explore
the
range
of
parameters
that
produce
vor-
tex
crystals,
providing
insights
for
such
measurements
and
constraints
to
the
future
modeling
work
of
the
polar
region.
We
vary
the
three
parameters
and
show
that
some
regions
of
the
parameter
space
lead
to
regular
patterns,
namely
vortex
crystals,
and
others
lead
to
chaotic
patterns.
2.
Parameters
and
length
scales
We
explore
the
initial
turbulence
scales
퐿
푖
up
to
the
scale
similar
to
the
circumpolar
cyclones,
210,
350,
and
700
km,
giving
퐿
푖
∕
푎
values
of
3
×
10
−3
,
5
×
10
−3
,
and
1
×
10
−2
.
Jovian
vortices
were
measured
to
have
velocity
profiles
that
reach
a
maximum
of
80
to
110
m/s
at
approximately
1000
km
from
their
centers
(
Grassi
et
al.
,
2018
).
Hence,
we
chose
푈
values
of
24,
70,
120,
and
170
m/s,
bracketing
and
exploring
beyond
the
observed
values.
The
Rossby
number,
푈
∕(
푓
퐿
푖
)
,
then
takes
values
of
0.10
to
2.31.
The
value
of
퐿
푑
in
the
polar
region
is
not
measured,
but
the
mid-latitude
value
is
estimated
to
be
1000
km
(
Wong
et
al.
,
2011
).
For
the
central
vortex
in
the
north,
the
requirement
that
the
thickness
cannot
be
negative
leads
to
the
smallest
value
of
퐿
푑
being
at
least
750
km
(
Ingersoll
et
al.
,
2022
).
Li
et
al.
(
2020b
)
used
a
fixed
value
퐿
=
1000
km
for
the
observed
radius
of
the
cyclones
and
varied
the
burger
number
퐵
푢
=
퐿
2
푑
∕
퐿
2
from
1
to
1000,
corresponding
to
퐿
푑
from
1000
km
to
3
.
1
×
10
4
km.
Here,
we
choose
퐿
푑
values
to
be
larger
than
the
750-km
lower
bound,
ranging
from
1750
to
8400
km.
In
particular,
we
focused
on
regions
of
transition
between
an
unstructured
pattern
and
vortex
crystals,
which
leads
to
different
퐿
푑
choices
with
different
퐿
푖
values.
The
input
parameters
and
the
associated
dimensionless
counterparts
of
all
60
simulations
are
summarized
in
the
supplementary
material,
Table
A.2
.
In
this
paper,
we
define
a
run
of
the
model
with
a
particular
choice
of
the
three
parameters
as
a
simulation.
A
snapshot
of
a
simulation
at
a
particular
time
is
an
image
that
has
an
overall
pattern
that
may
be
regular,
i.e.,
a
vortex
crystal,
or
irregular,
i.e.,
individual
vortex
patches
of
various
sizes
and
uneven
separation
distances.
The
three
basic
nondimensional
parameters
푈
∕(2
훺
푎
)
,
퐿
푖
∕
푎
,
and
퐿
푑
∕
푎
can
be
combined
in
a
variety
of
ways,
the
goal
being
a
combi-
nation
that
controls
the
final
output
of
the
model.
One
combination
is
the
Rossby
number
푈
∕(2
훺
퐿
푖
)
based
on
the
wavelength
of
the
initial
turbulence,
which
gives
the
importance
of
fluid
accelerations
relative
Icarus
429
(2025)
116438
2
S.
Chen
et
al.
to
the
Coriolis
acceleration,
where
the
former
arises
from
motions
on
the
scale
of
the
initial
disturbance.
With
the
choices
of
parameters
in
Table
A.2
,
we
have
푈
∕(2
훺
퐿
푖
)
values
ranging
from
0.1
to
2.31.
The
individual
cyclones
at
Jupiter’s
poles
have
peak
azimuthal
velocities
of
80
m/s
and
peak
relative
vorticities
휁
=
̂
퐤
⋅
(∇
×
퐮
)
≈
3
.
0
×
10
−4
s
−
1
,
the
latter
being
almost
equal
to
the
polar
planetary
vorticity
2
훺
of
3
.
5
×
10
−4
s
−
1
(
Ingersoll
et
al.
,
2022
).
The
second
important
combination
of
the
basic
initial
parameters
is
훥ℎ
∕
ℎ
.
It
is
the
fractional
departure
of
the
mean
surface
height
from
equilibrium,
and
when
the
value
becomes
of
order
1,
the
shallow
water
thickness
can
become
negative.
Then
the
sound
speed
becomes
undefined,
and
the
simulation
crashes.
To
estimate
the
value
of
훥ℎ
∕
ℎ
in
the
SW
system,
we
use
a
scaling
based
on
the
gradient
wind
balance:
푢
2
푟
+
2
훺
푢
=
−
푔
푟
푑
ℎ
푑
푟
,
(4)
Here
we
lack
a
characteristic
length
scale.
This
length
scale
is
highly
dependent
on
the
results
of
the
simulation,
but
we
want
it
to
be
related
to
the
initial
parameters
only.
Hence,
in
finding
a
parameter
that
controls
the
formation
of
vortex
crystals,
we
use
퐿
푖
as
the
length
scale.
We
consider
cyclones
here,
where
the
azimuthal
velocity
is
positive
and
the
Coriolis
term
2
훺
푢
has
the
same
sign
as
the
centrifugal
term
푢
2
∕
푟
.
Taking
푑
ℎ
∕
푑
푟
∼
훥ℎ
∕
푟
,
푢
∼
푈
,
푟
∼
퐿
푖
,
and
퐿
2
푑
∼
푔
푟
ℎ
∕(2
훺
)
2
,
we
obtain
훥ℎ
∕
ℎ
∼
푈
2
+
2
훺
푈
퐿
푖
(2
훺
퐿
푑
)
2
=
(
푈
2
훺
퐿
푑
)
2
[
1
+
(
푈
2
훺
퐿
푖
)
−1
]
.
(5)
This
covers
both
small
and
large
values
of
the
Rossby
number
푈
∕(2
훺
퐿
푖
)
.
When
this
Rossby
number
is
small,
the
initial
scale
of
the
turbu-
lence
matters,
which
gives
훥ℎ
∕
ℎ
∼
푈
퐿
푖
2
훺
퐿
2
푑
;
(6)
When
this
Rossby
number
based
on
퐿
푖
becomes
large,
the
velocity
dominates,
and
the
initial
scale
of
turbulence
has
minimal
influence
on
the
thickness
perturbation:
훥ℎ
∕
ℎ
∼
푈
2
(2
훺
퐿
푑
)
2
.
(7)
In
this
paper,
we
use
the
expression
presented
in
Eq. (
5
) to
calculate
the
values
of
훥ℎ
∕
ℎ
.
The
third
important
combination
of
initial
parameters
is
the
polar
beta
effect,
which
is
different
from
the
mid-latitude
value.
Here,
훽
=
푑
푓
∕
푑
푦
=
2
훺
cos
휙
∕
푎
,
where
푦
is
the
poleward
coordinate,
and
푎
is
the
radius
of
the
planet.
At
the
mid-latitudes,
the
quantity
√
푈
∕
훽
is
a
length
scale
below
which
eddies
dominate
and
above
which
zonal
jets
dominate
(
Rhines
,
1975
).
With
the
turbulence
cascading
from
small
scales
to
large
scales,
the
Rhines
scale
acts
as
a
cross-over
scale:
when
it
reaches
this
scale,
the
훽
-term
dominates,
and
the
cascade
stops.
However,
since
훽
goes
to
zero
at
the
poles,
the
cascade
can
never
reach
that
scale
at
the
poles,
and
eddies
should
always
dominate.
The
question
then
arises:
at
what
latitude
does
the
transition,
with
zonal
jets
below
that
latitude
and
eddies
above
it,
take
place?
If
we
expand
훽
at
radial
distance
푟
away
from
the
pole
using
its
derivative
at
the
pole:
훽
≈
2
훺
푟
∕
푎
2
,
and
assume
that
the
eddy
scale
inside
푟
is
푟
itself,
as
on
Jupiter.
Then
√
푈
훽
≈
√
푈
푎
2
2
훺
푟
∼
푟.
(8)
Solving
for
푟
gives
푟
∼
퐿
훾
=
(
푈
푎
2
2
훺
)
1
3
.
(9)
The
answer
has
been
given
before
(
Ingersoll
et
al.
,
2022
;
Siegelman
et
al.
,
2022b
),
and
this
gives
a
polar
lengthscale
퐿
훾
.
Here,
we
apply
an
assumption
about
the
eddies
to
provide
the
physical
intuition
of
the
lengthscale.
With
푈
=
80
m/s,
the
expression
gives
10,500
km,
which
is
a
fairly
good
match
to
the
8700
km
radius
of
the
circumpolar
ring.
It
corresponds
to
another
dimensionless
number,
the
critical
colatitude,
휃
푐
,
휃
푐
=
퐿
훾
∕
푎
=
(
푈
2
푎훺
)
1
3
.
(10)
3.
Model
description
We
solve
shallow-water
equations
using
an
Athena++-based
shallow-water
model
(
Li
and
Chen
,
2019
;
Stone
et
al.
,
2020
)
to
study
the
genesis
of
the
vortices
and
their
subsequent
evolution.
The
flow
takes
place
in
a
plane
tangent
to
the
surface
of
the
planet
at
the
pole.
We
choose
a
2048
×
2048
grid,
with
푥
푚푎푥
=
푦
푚푎푥
=
6
×
10
4
km
and
the
size
of
domain
퐿
=
2
푥
푚푎푥
.
This
domain
covers
from
the
pole
to
about
40
degrees
of
latitude.
The
setup
of
the
initial
condition
is
similar
to
Siegelman
et
al.
(
2022b
).
We
start
with
a
narrow
ring
in
the
wavenumber
space,
̃
훷
(
푘
푥
,
푘
푦
)
,
centered
at
the
origin
with
a
radius
of
2
휋
∕
퐿
푖
,
randomly
seeded
and
uniformly
distributed
in
direction.
Then
the
initial
geopo-
tential
is
obtained
from
the
inverse
Fourier
transform
of
̃
훷
(
푘
푥
,
푘
푦
)
.
The
initial
velocity
field
is
geostrophically
balanced
by
taking
horizontal
derivatives
of
the
initial
훷
field
푓
̂
퐤
×
퐮
=
−∇
훷
.
(11)
This
region
extends
to
50
◦
latitude,
making
it
much
larger
than
the
vortices.
After
the
initial
time
step,
the
flow
evolves
according
to
Eqs. (
1
) and (
2
).
We
also
impose
a
sponge
layer
that
starts
at
50
◦
latitude,
where
the
velocities
are
dissipated
by
Rayleigh
friction.
The
dissipation
rate
increases
linearly
from
0
at
50
◦
latitude
to
0.5
day
−
1
at
40
◦
latitude.
The
amplitude
of
the
initial
turbulence
is
chosen
so
that
the
root
mean
square
velocity
in
the
field
between
60
◦
latitude
and
the
pole
is
푈
.
In
other
words,
the
mean
specific
kinetic
energy
between
60
◦
latitude
and
the
pole
is
푈
2
∕2
.
The
choice
of
60
◦
is
made
as
most
of
the
energy
in
the
final
patterns
is
always
confined
within
this
region.
To
quantitatively
determine
how
irregular
a
pattern
is,
we
locate
all
the
vortices
in
the
image
—
patches
of
fluid
whose
absolute
value
of
vorticity
is
larger
than
0
.
2
훺
,
and
we
calculate
the
area
of
each
patch.
Then
we
calculate
the
mean
area
퐴
and
standard
deviation
푑
퐴
of
all
the
patches
in
this
image.
The
irregularity
measure
of
the
image
is
then
푑
퐴
∕
퐴
.
For
aid
in
graphical
depiction
in
Figs.
3
–
6
,
we
use
a
logarithmic
function
to
scale
the
푑
퐴
∕
퐴
values
to
irregularity
levels,
ranging
from
0
to
9.
The
irregularity
levels
are
calculated
as
5
.
8
×
(
log
10
(
푑
퐴
∕
퐴
)
+
0
.
9
)
.
(12)
The
irregularity
levels
are
then
rounded
to
the
nearest
integer
for
a
single-digit
representation.
Under
this
scaling,
irregularity
levels
smaller
than
5
generally
correspond
to
vortex
crystals,
and
irregular
patterns
correspond
to
numbers
6
to
9.
Each
simulation’s
dA/A
value
and
irregularity
level
at
day
600
are
calculated
and
shown
in
the
supplementary
material
Table
A.2
.
4.
Results
Fig.
1
shows
the
evolution
of
the
flow
undergoing
an
inverse
energy
cascade
under
different
parameters.
In
all
simulations,
rapid
evolution
occurs
in
the
first
10
days,
where
vortices
develop
from
the
initial
turbulence.
Simulation
8
evolves
into
a
regular
pattern
of
similarly
sized
vortices.
In
comparison,
simulations
58,
16,
and
6
end
up
with
larger
vortices,
meaning
more
extensive
merging
occurs.
By
comparing
the
chaotic
simulations
58,
16,
and
6,
each
of
which
varies
one
of
the
parameters
relative
to
the
vortex
crystal
simulation
8,
it
is
clear
that
merging
happens
when
we
increase
퐿
푖
or
푈
or
when
we
decrease
퐿
푑
.
These
variations
all
lead
to
a
larger
훥ℎ
∕
ℎ
.
We
have
provided
movies
in
the
online
materials
for
these
four
simulations
in
the
supplementary
materials,
each
with
a
short
10-day
Icarus
429
(2025)
116438
3
S.
Chen
et
al.
Fig.
1.
The
vortex
development
as
a
function
of
time,
for
four
simulations:
an
example
where
the
vortices
organize
to
a
vortex
crystal
(simulation
8),
and
examples
where
the
vortices
do
not
form
regular
patterns
(simulations
58,
16,
and
6),
where
each
has
one
of
the
parameters
changed
to
obtain
chaotic
patterns.
The
parameter
changed
relative
to
simulation
8
is
shown
in
bold
fonts.
The
patterns
at
days
0,
10,
50,
100,
and
600
are
presented.
The
circles
are
at
80,
70,
60,
and
50
degrees
latitudes.
The
푑
퐴
∕
퐴
value,
as
defined
in
the
text
in
conjunction
with
Eq. (
12
),
is
shown
for
each
simulation.
Fig.
2.
The
conversion
of
kinetic
energy
(noted
as
the
red
shade)
to
potential
energy
(blue
shade)
and
the
dissipation
of
kinetic
energy
(the
gray
shade).
The
푥
-axis
(time
axis)
is
scaled
to
exaggerate
the
first
10
days’
variations.
It
is
linear
on
days
0–10,
and
10–600,
respectively.
The
same
four
example
simulations
as
in
Fig.
1
are
shown,
with
kinetic
energy
normalized
by
the
initial
kinetic
energy,
and
the
potential
energy
is
chosen
to
be
zero
at
the
initial
state.
movie
with
each
frame
captured
at
one-hour
intervals
and
a
long
600-
day
movie
with
each
frame
captured
at
one-day
intervals.
In
the
first
10
days,
the
flow
quickly
organizes
into
vortices.
Compared
to
simulation
8,
simulations
58
and
16
show
significant
differences
in
the
first
10
days,
one
due
to
a
different
initial
wavelength
of
the
turbulence
and
another
due
to
a
stronger
stirring.
Despite
the
difference
from
the
starting
turbulence,
the
further
evolution
of
simulations
58,
16,
and
6
have
merging
events
after
day
50:
vortices
approach
each
other,
form
a
fast-rotating
pair,
and
eventually
merge.
Some
of
the
merged
vortices
show
an
irregular
shape.
In
contrast,
the
vortices
in
simulation
8
mostly
vibrate
in
a
vortex
crystal
structure
but
do
not
merge
into
each
other.
The
simulations
can
be
viewed
from
the
perspective
of
energy
conversion.
The
same
four
examples
shown
in
Fig.
1
are
shown
in
Fig.
2
.
Initially,
most
of
the
energy
is
stored
in
the
smallest
scales
as
kinetic
energy.
The
kinetic
energy
can
be
converted
to
potential
energy,
radiated
away
as
gravity
waves
and
dissipated
in
the
sponge
layer,
or
by
artificial
viscosity.
In
the
first
several
days
of
simulation,
the
kinetic
energy
is
quickly
lost
by
converting
to
potential
energy
and
radiating
gravity
waves
into
the
sponge
layer.
The
strength
of
this
dissipation
increases
with
a
larger
initial
velocity
or
a
smaller
initial
turbulence
scale.
After
the
first
10
days,
the
kinetic
energy
loss
for
the
vortex
crystal
pattern
becomes
minimal
(simulation
8).
Other
simulations
continue
to
have
vortices
merging,
and
the
kinetic
energy
continues
to
be
converted
to
potential
energy
or
dissipated
in
the
sponge
layer.
In
Figs.
3
–
5
,
we
organize
60
images
at
day
600
and
label
them
by
훥ℎ
∕
ℎ
.
Notice
that
simulations
with
small
훥ℎ
∕
ℎ
show
a
more
regular
pattern,
closer
to
a
vortex
crystal,
as
compared
to
simulations
with
large
훥ℎ
∕
ℎ
.
The
levels
of
irregularity
defined
as
the
single-digit
numbers
in
Figs.
3
–
5
agree
well
with
our
subjective
impression.
Numbers
less
than
5
look
like
regular
vortex
crystals.
Numbers
5
to
6
look
somewhat
irregular,
being
a
transitioning
regime
to
more
chaotic
patterns,
marked
by
digits
7
to
9.
It
is
not
a
perfect
metric,
and
the
connection
to
퐿
푑
and
other
input
parameters
is
not
perfect,
perhaps
because
there
are
random
elements
in
the
merging
process.
Blank
spaces
in
the
left
panel
are
where
훥ℎ
∕
ℎ
was
so
large
that
the
program
crashed
before
reaching
600
days.
Crashing
was
always
due
to
the
geopotential
훷
becoming
negative.
If
one
scans
across
any
of
the
4
images
with
constant
퐿
푑
and
퐿
푖
in
Figs.
3
–
5
,
it
is
clear
that
irregularity
increases
with
푈
.
Similarly,
if
one
scans
down
any
of
the
4
images
with
constant
푈
and
퐿
푖
,
it
is
clear
that
irregularity
decreases
as
퐿
푑
increases.
Both
of
these
observations
are
consistent
with
irregularity
increasing
with
훥ℎ
∕
ℎ
.
As
shown
in
the
movies
and
Fig.
1
,
either
the
vortex
crystal
structure
appears
during
the
first
10
days
and
the
vortices
stay
clear
of
each
other,
or
it
never
appears
—
the
flow
keeps
evolving,
and
the
vortices
keep
merging.
Fig.
6
is
a
regime
diagram
in
the
3D
space
defined
by
the
input
parameters
푈
,
퐿
푑
,
and
퐿
푖
.
In
this
figure,
we
project
the
irregularity
levels
defined
in
Eq. (
12
) onto
the
3
axes,
and
based
on
the
value
ranges
Icarus
429
(2025)
116438
4