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 Siegelman, Young & Ingersoll (2022) 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.

julia

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 GridFittedBoundary) at 1500 km from the pole forming the polar disk.

julia

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

grid = LambertConformalConicGrid(Float64;
                                 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_earth = grid.radius
R_bowl  = 1500kilometers
bowl_mask(λ, φ, z) = R_earth * (π/2 - deg2rad(φ)) > R_bowl
ibg = ImmersedBoundaryGrid(grid, GridFittedBoundary(bowl_mask))

128×128×1 ImmersedBoundaryGrid{Float64, Bounded, Bounded, Bounded} on CPU with 7×7×7 halo:
├── immersed_boundary: GridFittedBoundary{Field{Center, Center, Center, Nothing, Oceananigans.Grids.ZRegOrthogonalSphericalShellGrid{Float64, Bounded, Bounded, Bounded, Oceananigans.Grids.StaticVerticalDiscretization{OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64, Float64}, LambertConformalConic{Float64}, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, CPU, Float64}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Bool, 3, Array{Bool, 3}}, Bool, FieldBoundaryConditions{Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Nothing, @NamedTuple{bottom_and_top::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(16, 16)}, Oceananigans.Utils.OffsetStaticSize{(1:128, 1:128)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, south_and_north::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(128, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:128, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_south_and_north_halo!)}, west_and_east::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(128, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:128, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_west_and_east_halo!)}}, @NamedTuple{bottom_and_top::Tuple{Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition}, south_and_north::Tuple{Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition}, west_and_east::Tuple{Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition}}}, Nothing, Nothing}}
├── underlying_grid: 128×128×1 OrthogonalSphericalShellGrid{Float64, Bounded, Bounded, Bounded} on CPU with 7×7×7 halo
├── centered at: North Pole, (λ, φ) = (0.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.

julia

Δ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 radius r_ring = 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 projected (xp, yp) coordinates of the polar stereographic limit lets us write η and the geostrophic velocities (u, v) as plain Gaussians of distance from each vortex centre. set! is told these are in the grid's intrinsic frame.

julia

N_ring   = 6
r_ring   = 900kilometers
σ_vortex = 200kilometers
η₀       = -13

ring_λ      = [360k/N_ring for k in 0:N_ring-1]
ring_φ      = fill(90 - rad2deg(r_ring/R_earth), N_ring)
vortex_λ    = vcat(ring_λ, [0])
vortex_φ    = vcat(ring_φ, [90])
vortex_xpyp = [lcc_forward(grid.conformal_mapping, λv, φv)
               for (λv, φv) in zip(vortex_λ, vortex_φ)]

g_const  = Oceananigans.defaults.gravitational_acceleration
f_pole   = 2 * Oceananigans.defaults.planet_rotation_rate
geo_coef = g_const * η₀ / (f_pole * σ_vortex^2)

gaussian(xp, yp, xv, yv) = exp(-((xp - xv)^2 + (yp - yv)^2) / (2σ_vortex^2))

function η_init(λ, φ, z)
    xp, yp = lcc_forward(grid.conformal_mapping, λ, φ)
    return sum(η₀ * gaussian(xp, yp, xv, yv) for (xv, yv) in vortex_xpyp)
end

function u_init(λ, φ, z)
    xp, yp = lcc_forward(grid.conformal_mapping, λ, φ)
    return sum( geo_coef * (yp - yv) * gaussian(xp, yp, xv, yv) for (xv, yv) in vortex_xpyp)
end

function v_init(λ, φ, z)
    xp, yp = lcc_forward(grid.conformal_mapping, λ, φ)
    return sum(-geo_coef * (xp - xv) * gaussian(xp, yp, xv, yv) for (xv, yv) in vortex_xpyp)
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.

julia

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(time(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_existing = 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 (1.108 minutes)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (41.316 seconds).
iter  1000, t = 13.889 days, max|u| = 2.738 m/s, max|η| = 13.73 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.742 m/s, max|η| = 13.67 m, CFL = 0.196
iter  4000, t = 55.556 days, max|u| = 2.817 m/s, max|η| = 13.78 m, CFL = 0.200
iter  5000, t = 69.444 days, max|u| = 3.126 m/s, max|η| = 14.14 m, CFL = 0.199
iter  6000, t = 83.333 days, max|u| = 2.894 m/s, max|η| = 14.14 m, CFL = 0.215
iter  7000, t = 97.222 days, max|u| = 2.834 m/s, max|η| = 14.30 m, CFL = 0.213
iter  8000, t = 111.111 days, max|u| = 2.860 m/s, max|η| = 14.51 m, CFL = 0.223
[ Info: Simulation is stopping after running for 5.061 minutes.
[ Info: Simulation time 120 days equals or exceeds stop time 120 days.

Visualization

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

julia

using CairoMakie

CairoMakie.activate!(type = "png")

ηts = FieldTimeSeries(filename, "η")
ζts = FieldTimeSeries(filename, "ζ")
sts = FieldTimeSeries(filename, "s")
times = ηts.times
Nt = length(times)

η_lim = maximum(maximum(abs, interior(ηts[n])) for n in 1:Nt)
ζ_lim = maximum(maximum(abs, interior(ζts[n])) for n in 1:Nt) * 0.5
s_lim = maximum(maximum(abs, interior(sts[n])) for n in 1:Nt)

n = Observable(1)
title = @lift @sprintf("Polar vortex crystal — t = %.2f d", times[$n] / 86400)
ηₙ = @lift ηts[$n]
ζₙ = @lift ζts[$n]
sₙ = @lift sts[$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:

julia

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.4
Commit 01a2eadb047 (2026-01-06 16:56 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-18.1.7 (ORCJIT, znver4)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
  LD_LIBRARY_PATH = 
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.12.4+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.12.4+0.x64.linux.gnu/share/julia
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-32087/docs/
  JULIA_VERSION = 1.12.4
  JULIA_CUDA_USE_COMPAT = false
  JULIA_LOAD_PATH = @:@v#.#:@stdlib
  JULIA_VERSION_ENZYME = 1.10.10
  JULIA_DEBUG = Literate

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

julia

import Pkg
Pkg.status()

Status `~/Oceananigans.jl-32087/docs/Project.toml`
  [79e6a3ab] Adapt v4.6.1
⌅ [052768ef] CUDA v6.1.0
  [13f3f980] CairoMakie v0.15.12
  [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.4.1
  [4710194d] DocumenterVitepress v0.3.4
  [033835bb] JLD2 v0.6.4
  [63c18a36] KernelAbstractions v0.9.41
  [98b081ad] Literate v2.21.0
  [da04e1cc] MPI v0.20.26
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.110.4 `..`
  [f27b6e38] Polynomials v4.1.1
  [6038ab10] Rotations v1.7.1
  [d496a93d] SeawaterPolynomials v0.3.10
  [09ab397b] StructArrays v0.7.3
  [bdfc003b] TimesDates v0.3.3
  [0a941bbe] Zarr v0.10.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.12.0
  [44cfe95a] Pkg v1.12.1
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.

`NonhydrostaticModel`

`HydrostaticFreeSurfaceModel`

`ShallowWaterModel`

Polar vortex crystal

Grid

Model

Initial condition

Simulation

Visualization

Julia version and environment information

Polar vortex crystal ​

Grid ​

Model ​

Initial condition ​

Simulation ​

Visualization ​

Julia version and environment information ​

Polar vortex crystal

Grid

Model

Initial condition

Simulation

Visualization

Julia version and environment information