Polar vortex crystal

Self-organisation of cyclones into a stable ring around a central cyclone in rotating shallow water on a polar disk — a regime studied by as a model for Jupiter's polar vortex clusters observed by Juno's JIRAM instrument.

The grid is a LambertConformalConicGrid centred exactly on the North Pole with standard_parallel = 90 — the polar stereographic limit, where the cone is tangent to the sphere at the pole and the projection has no antemeridian wedge.

using Oceananigans
using Oceananigans.OrthogonalSphericalShellGrids
using Oceananigans.Units
using Printf

Grid

A 128² horizontal × 1 vertical-cell grid covering a 3200 km box on the pole, with a circular wall (an ImmersedBoundaryGrid with GridFittedBottom whose bottom height reaches the top of the domain) at 1500 km from the pole forming the polar disk.

Nx = Ny = 128
Δ  = 25kilometers
H  = 1000meters

grid = LambertConformalConicGrid(; size = (Nx, Ny, 1),
                                   center = (0, 90),
                                   spacing = Δ,
                                   standard_parallel = 90,
                                   latitude_of_origin = 90,
                                   central_longitude = 0,
                                   z = (-H, 0),
                                   halo = (7, 7, 7))

R = grid.radius
disk_radius = 1500kilometers
bowl_bottom(λ, φ) = R * deg2rad(90 - φ) > disk_radius ? 0 : -H
ibg = ImmersedBoundaryGrid(grid, GridFittedBottom(bowl_bottom))

128×128×1 ImmersedBoundaryGrid{Float64, Bounded, Bounded, Bounded} on CPU with 7×7×7 halo:
├── immersed_boundary: GridFittedBottom(mean(z)=-697.51, min(z)=-1000.0, max(z)=0.0)
├── underlying_grid: 128×128×1 LambertConformalConicGrid{Float64, Bounded, Bounded, Bounded} on CPU with 7×7×7 halo
├── centered at: North Pole, (λ, φ) = (0, 90.0)
├── longitude: Bounded  extent 28.6285 degrees variably spaced with min(Δλ)=0.217957, max(Δλ)=0.22483
├── latitude:  Bounded  extent 28.6285 degrees variably spaced with min(Δφ)=0.217957, max(Δφ)=0.22483
└── z:         Bounded  z ∈ [-1000.0, 0.0]     regularly spaced with Δz=1000.0

Model

Single-layer barotropic shallow-water-style dynamics via a HydrostaticFreeSurfaceModel with one vertical level. We use HydrostaticSphericalCoriolis (so f ≈ 2Ω near the pole), a SplitExplicitFreeSurface whose number of substeps is set from the barotropic-wave CFL, and WENOVectorInvariant momentum advection.

Δt = 20minutes

model = HydrostaticFreeSurfaceModel(ibg;
                                    coriolis           = HydrostaticSphericalCoriolis(),
                                    free_surface       = SplitExplicitFreeSurface(ibg; cfl = 0.7, fixed_Δt = Δt),
                                    momentum_advection = WENOVectorInvariant())

HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×1 ImmersedBoundaryGrid{Float64, Bounded, Bounded, Bounded} on CPU with 7×7×7 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: ()
├── closure: Nothing
├── buoyancy: Nothing
├── free surface: SplitExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── substepping: FixedSubstepNumber(14)
├── advection scheme: 
│   └── momentum: WENOVectorInvariant{5, Float64}(vorticity_order=9, vertical_order=5)
├── vertical_coordinate: ZCoordinate
└── coriolis: HydrostaticSphericalCoriolis{Oceananigans.Advection.EnstrophyConserving{Float64}, Float64, Oceananigans.Coriolis.HydrostaticFormulation}

Initial condition

Six Gaussian cyclones at a distance ring_radius = 900 km from the pole, evenly spaced in longitude, plus one central cyclone at the pole. The depression amplitude η₀ = -13 m and width σ = 200 km give an initial Rossby number Ro ≈ -2gη₀/(f²σ²) ≈ 0.3 at each vortex centre.

Working in the projected coordinates (x, y) of the polar stereographic limit lets us write η as a sum of Gaussians centred on (xᵥ, yᵥ), and the geostrophic velocities u = -(g/f) ∂η/∂y and v = (g/f) ∂η/∂x in closed form. set! is told these velocities are in the grid's intrinsic frame.

ring_size   = 6
ring_radius = 900kilometers
σ  = 200kilometers
η₀ = -13

ring_latitude = 90 - rad2deg(ring_radius / R)
ring_longitudes = [360k / ring_size for k in 0:ring_size-1]
vortex_coordinates = [[(λ, ring_latitude) for λ in ring_longitudes]; (0, 90)]
vortex_centres = [lcc_forward(grid.conformal_mapping, λ, φ) for (λ, φ) in vortex_coordinates]

g = Oceananigans.defaults.gravitational_acceleration
f = 2 * Oceananigans.defaults.planet_rotation_rate

ηᵥ(x, y, xᵥ, yᵥ) = η₀ * exp(-((x - xᵥ)^2 + (y - yᵥ)^2) / 2σ^2)

function η_init(λ, φ, z)
    x, y = lcc_forward(grid.conformal_mapping, λ, φ)
    return sum(ηᵥ(x, y, xᵥ, yᵥ) for (xᵥ, yᵥ) in vortex_centres)
end

function u_init(λ, φ, z)
    x, y = lcc_forward(grid.conformal_mapping, λ, φ)
    return sum(g / f * (y - yᵥ) / σ^2 * ηᵥ(x, y, xᵥ, yᵥ) for (xᵥ, yᵥ) in vortex_centres)
end

function v_init(λ, φ, z)
    x, y = lcc_forward(grid.conformal_mapping, λ, φ)
    return sum(- g / f * (x - xᵥ) / σ^2 * ηᵥ(x, y, xᵥ, yᵥ) for (xᵥ, yᵥ) in vortex_centres)
end

set!(model, η = η_init, u = u_init, v = v_init; intrinsic_velocities = true)

Simulation

120-day run at Δt = 20 minutes (outer timestep limited by the advective CFL; the split-explicit substeps handle the much faster barotropic gravity waves internally). The progress callback reports the advective CFL alongside the maximum velocity and free-surface displacement. Snapshots are saved every 12 hours.

simulation = Simulation(model; Δt, stop_time = 120days)

advective_cfl = AdvectiveCFL(simulation.Δt)

progress(sim) = @printf("iter %5d, t = %s, max|u| = %.3f m/s, max|η| = %.2f m, CFL = %.3f\n",
                        iteration(sim), prettytime(sim),
                        maximum(abs, sim.model.velocities.u),
                        maximum(abs, sim.model.free_surface.displacement),
                        advective_cfl(sim.model))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(1000))

u, v, w = model.velocities
η = model.free_surface.displacement
ζ = ∂x(v) - ∂y(u)
s = sqrt(u^2 + v^2)

filename = "polar_vortex_crystal.jld2"
simulation.output_writers[:fields] = JLD2Writer(model, (; η, ζ, s);
                                                filename,
                                                schedule = TimeInterval(12hours),
                                                overwrite_files = true)

run!(simulation)

[ Info: Initializing simulation...
iter     0, t = 0 seconds, max|u| = 2.643 m/s, max|η| = 12.97 m, CFL = 0.189
[ Info:     ... simulation initialization complete (36.837 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (9.171 seconds).
iter  1000, t = 13.889 days, max|u| = 2.741 m/s, max|η| = 13.72 m, CFL = 0.195
iter  2000, t = 27.778 days, max|u| = 2.784 m/s, max|η| = 13.61 m, CFL = 0.194
iter  3000, t = 41.667 days, max|u| = 2.755 m/s, max|η| = 13.67 m, CFL = 0.197
iter  4000, t = 55.556 days, max|u| = 2.816 m/s, max|η| = 13.76 m, CFL = 0.200
iter  5000, t = 69.444 days, max|u| = 3.145 m/s, max|η| = 14.20 m, CFL = 0.201
iter  6000, t = 83.333 days, max|u| = 2.932 m/s, max|η| = 14.14 m, CFL = 0.212
iter  7000, t = 97.222 days, max|u| = 2.834 m/s, max|η| = 14.38 m, CFL = 0.213
iter  8000, t = 111.111 days, max|u| = 2.882 m/s, max|η| = 14.47 m, CFL = 0.218
[ Info: Simulation is stopping after running for 3.584 minutes.
[ Info: Simulation time 120 days equals or exceeds stop time 120 days.

Visualization

Animate η, ζ, and |u| over the 120-day evolution.

using CairoMakie

CairoMakie.activate!(type = "png")

η_timeseries = FieldTimeSeries(filename, "η")
ζ_timeseries = FieldTimeSeries(filename, "ζ")
s_timeseries = FieldTimeSeries(filename, "s")
times = η_timeseries.times
Nt = length(times)

η_lim = maximum(abs, η_timeseries)
ζ_lim = maximum(abs, ζ_timeseries) / 2
s_lim = maximum(s_timeseries)

n = Observable(1)
title = @lift "Polar vortex crystal — t = " * prettytime(times[$n])
ηₙ = @lift η_timeseries[$n]
ζₙ = @lift ζ_timeseries[$n]
sₙ = @lift s_timeseries[$n]

fig = Figure(size = (1500, 540))
Label(fig[0, 1:6], title; fontsize = 18, tellwidth = false)

ax_η = Axis(fig[1, 1], aspect = 1, title = "η (m)")
hm_η = heatmap!(ax_η, ηₙ; colormap = :balance, colorrange = (-η_lim, η_lim))
Colorbar(fig[1, 2], hm_η)

ax_ζ = Axis(fig[1, 3], aspect = 1, title = "ζ (1/s)")
hm_ζ = heatmap!(ax_ζ, ζₙ; colormap = :balance, colorrange = (-ζ_lim, ζ_lim))
Colorbar(fig[1, 4], hm_ζ)

ax_s = Axis(fig[1, 5], aspect = 1, title = "|u| (m/s)")
hm_s = heatmap!(ax_s, sₙ; colormap = :speed, colorrange = (0, s_lim))
Colorbar(fig[1, 6], hm_s)

for ax in (ax_η, ax_ζ, ax_s)
    hidedecorations!(ax)
end

record(fig, "polar_vortex_crystal.mp4", 1:Nt; framerate = 12) do i
    n[] = i
end

---

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.13.0
Commit d1c37793dd2 (2026-09-09 19:00 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 9374F 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/share/julia
  JULIA_VERSION_ENZYME = 1.13.0
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  LD_LIBRARY_PATH = 
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_VERSION = 1.13.0
  JULIA_CUDA_USE_COMPAT = false
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-34931/docs/
  JULIA_DEBUG = Literate

These were the top-level packages installed in the environment:

import Pkg
Pkg.status()

Status `~/Oceananigans.jl-34931/docs/Project.toml`
  [79e6a3ab] Adapt v4.7.4
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.16
  [8e50ac2c] ConservativeRegridding v0.2.16
⌅ [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.5.0
  [4710194d] DocumenterVitepress v0.3.8
  [7da242da] Enzyme v0.13.216
  [033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.45
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.29
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.114.0 `..`
  [429524aa] Optim v2.3.2
  [f27b6e38] Polynomials v4.1.3
⌅ [3c362404] Reactant v0.2.292
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.11
  [09ab397b] StructArrays v0.7.4
  [bdfc003b] TimesDates v0.3.3
  [0a941bbe] Zarr v0.11.2
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.13.0
  [44cfe95a] Pkg v1.13.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

---

This page was generated using Literate.jl.