Barotropic jet and Double Drake on a Yin-Yang grid

This example runs two simulations on a grid that covers the whole sphere, poles included. The first is the unstable midlatitude jet of , a standard test of shallow-water models on the sphere, in a single layer of fluid. The second is the "Double Drake" ocean of , in which two thin continents and a circumpolar channel shape wind-driven, baroclinically unstable fronts in both hemispheres.

The Yin-Yang grid

A latitude-longitude grid cannot include the poles, where its meridians converge. The Yin-Yang grid of avoids the poles with two identical latitude-longitude panels, Yin and Yang. Each panel spans 270° of longitude and 90° of latitude around its own equator, where the grid is nearly uniform. Yang is Yin rotated so that its equator passes through the geographic poles. The two panels overlap at their edges, and each one fills its halos by interpolating from the interior of the other.

A curve through the overlap divides the sphere between the panels, and integrals and averages count every point of the sphere once, on the side of the curve that owns it. The Yin panel ends at 45°N, so the jet of the first simulation runs along the curve where the panels meet.

Install dependencies

julia

using Pkg
pkg"add Oceananigans, CairoMakie, CUDA"

julia

using Oceananigans
using Oceananigans.Units
using Oceananigans.Architectures: on_architecture
using Oceananigans.Grids: λnode, φnode, znode
using Oceananigans.Operators: rotation_angle, ζ₃ᶠᶠᶜ, ℑxyᶜᶜᵃ
using Oceananigans.MultiRegion: owned_weights
using CUDA
using Printf
using CairoMakie

We run on the GPU. To run on the CPU, replace GPU() with CPU(): the jet then takes about 5 minutes, and Double Drake about 17 minutes, on a laptop with 12 threads.

julia

arch = GPU()

GPU{CUDABackend}(CUDABackend(false, true))

Barotropic jet

A zonal jet at 45°N starts in balance with the height of a single layer of fluid. A small bump in the height disturbs the jet, which meanders and rolls up into a chain of vortices within six days.

A single layer of fluid

A HydrostaticFreeSurfaceModel with one layer, no buoyancy and a MutableVerticalDiscretization solves the nonlinear shallow-water equations: the vertical coordinate stretches the layer with the free surface , so the fluid depth is . The layer is km deep at rest.

Each panel of the Earth-sized YinYangGrid has 128 cells in latitude and 384 in longitude, a spacing of 0.7°, or 78 km.

julia

H = 10kilometers

jet_grid = YinYangGrid(arch; size=(384, 128, 1), z=MutableVerticalDiscretization((-H, 0)))

jet_model = HydrostaticFreeSurfaceModel(jet_grid; buoyancy=nothing, tracers=(),
                                        coriolis = HydrostaticSphericalCoriolis(),
                                        closure = HorizontalScalarDiffusivity(ν=1e5),
                                        free_surface = ExplicitFreeSurface(),
                                        timestepper = :SplitRungeKutta3)

HydrostaticFreeSurfaceModel{CUDAGPU, MultiRegionGrid}(time = 0 seconds, iteration = 0)
├── grid: 390×134×1 YinYangGrid{Float32, Connected, Connected, Bounded} on CUDAGPU with 3×3×3 halo
├── timestepper: SplitRungeKuttaTimeStepper
├── tracers: ()
├── closure: HorizontalScalarDiffusivity{ExplicitTimeDiscretization}(ν=100000.0, κ=nothing)
├── buoyancy: Nothing
├── free surface: ExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
├── advection scheme: 
│   └── momentum: MultiRegionObject{Tuple{VectorInvariant{1, Float32, ExplicitTimeDiscretization, Oceananigans.Advection.EnstrophyConserving{Float32}, Oceananigans.Advection.VelocityStencil, Oceananigans.Advection.EnergyConserving{Float32}, Oceananigans.Advection.EnergyConserving{Float32}, Oceananigans.Advection.EnergyConserving{Float32}, Oceananigans.Advection.OnlySelfUpwinding{Oceananigans.Advection.EnergyConserving{Float32}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.divergence_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.divergence_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.u_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.v_smoothness)}}, false}, VectorInvariant{1, Float32, ExplicitTimeDiscretization, Oceananigans.Advection.EnstrophyConserving{Float32}, Oceananigans.Advection.VelocityStencil, Oceananigans.Advection.EnergyConserving{Float32}, Oceananigans.Advection.EnergyConserving{Float32}, Oceananigans.Advection.EnergyConserving{Float32}, Oceananigans.Advection.OnlySelfUpwinding{Oceananigans.Advection.EnergyConserving{Float32}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.divergence_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.divergence_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.u_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.v_smoothness)}}, false}}}
├── vertical_coordinate: ZStarCoordinate
└── coriolis: HydrostaticSphericalCoriolis{Oceananigans.Advection.EnstrophyConserving{Float32}, Float32, Oceananigans.Coriolis.HydrostaticFormulation}

Following , a horizontal viscosity of removes the enstrophy that cascades to the grid scale.

The jet

The zonal velocity is

₁ ₁

and zero elsewhere, with , , $₁$ , and $₁$ , so that the jet peaks at 45°N.

julia

uₘₐₓ = 80
φ₀ = π / 7
φ₁ = π / 2 - φ₀
eₙ = exp(-4 / (φ₁ - φ₀)^2)

jet(φ) = φ₀ < φ < φ₁ ? uₘₐₓ / eₙ * exp(1 / ((φ - φ₀) * (φ - φ₁))) : 0.0

jet (generic function with 1 method)

The height of the fluid is in gradient-wind balance with the jet. Since the jet is zonally symmetric, depends only on latitude and from geostrophic balance of the -momentum equation we have:

where is the radius of the Earth and . We integrate this balance across the jet with the midpoint rule, and remove the global mean of the result from the free surface, so that the mean depth is .

julia

g = Oceananigans.defaults.gravitational_acceleration
Ω = Oceananigans.defaults.planet_rotation_rate
R = Oceananigans.defaults.planet_radius

function integral(f, a, b; n=500)
    Δ = (b - a) / n
    return Δ * sum(f(a + (k - 1/2) * Δ) for k in 1:n)
end

dhdφ(φ) = - R * jet(φ) * (2Ω * sin(φ) + jet(φ) * tan(φ) / R) / g
height_departure(φ) = integral(dhdφ, φ₀, clamp(φ, φ₀, φ₁))
mean_departure = integral(φ -> height_departure(φ) * cos(φ), -π/2, π/2) / 2

-158.1764715256522

To destabilize the jet, we add the bump

with m, $₂$ , and , for .

julia

function bump(λ, φ; ĥ=120, φ₂=π/4, α=1/3, β=1/15)
    λ = mod(λ + π, 2π) - π
    return ĥ * cos(φ) * exp(-(λ / α)^2) * exp(-((φ₂ - φ) / β)^2)
end

bump (generic function with 1 method)

The functions above take radians, and set! passes longitude and latitude in degrees. set! also rotates the eastward velocity into the basis of each panel.

julia

uᵢ(λ, φ, z) = jet(deg2rad(φ))
ηᵢ(λ, φ, z) = height_departure(deg2rad(φ)) - mean_departure + bump(deg2rad(λ), deg2rad(φ))

set!(jet_model, u=uᵢ, η=ηᵢ)

The simulation

The free surface is explicit, so the time step must resolve gravity waves, which travel at . With a time step of 40 seconds, the fastest waves cross about a third of the smallest cells in one step.

julia

jet_simulation = Simulation(jet_model; Δt=40, stop_time=6days)

Simulation of HydrostaticFreeSurfaceModel{CUDAGPU, MultiRegionGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 40 seconds
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 6 days
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

The panels disagree slightly on the fluxes between them, so the volume of the fluid drifts. A conservation fixer removes this drift.

julia

add_conservation_fixer!(jet_simulation)

A progress message reports the wall time since the previous message.

julia

wall_clock = Ref(time_ns())

function wall_time()
    elapsed = 1e-9 * (time_ns() - wall_clock[])
    wall_clock[] = time_ns()
    return prettytime(elapsed)
end

function progress(sim)
    u, v, w = sim.model.velocities
    @info @sprintf("iteration: %d, time: %s, max|u|: %.2f m/s, max|w|: %.1e m/s, wall time: %s",
                   iteration(sim), prettytime(sim), maximum(abs, u), maximum(abs, w), wall_time())
    return nothing
end

add_callback!(jet_simulation, progress, IterationInterval(1000))

Every three hours, we save the vertical vorticity of each panel.

julia

@inline ζᶜᶜᶜ(i, j, k, grid, u, v) = ℑxyᶜᶜᵃ(i, j, k, grid, ζ₃ᶠᶠᶜ, u, v)

function vertical_vorticity(model)
    u, v, w = model.velocities
    return Field(KernelFunctionOperation{Center, Center, Center}(ζᶜᶜᶜ, model.grid, u, v))
end

jet_ζ = vertical_vorticity(jet_model)

jet_times = Float64[]
vorticity = (yin = Matrix{Float64}[], yang = Matrix{Float64}[])

function save_vorticity(sim)
    compute!(jet_ζ)
    push!(jet_times, time(sim))
    push!(vorticity.yin, Array(interior(jet_ζ[1], :, :, 1)))
    push!(vorticity.yang, Array(interior(jet_ζ[2], :, :, 1)))
    return nothing
end

add_callback!(jet_simulation, save_vorticity, TimeInterval(3hours))

We also track the volume of the fluid. The coordinate stretches the cells with the free surface, so the volume is the integral of one over the grid.

julia

fluid = CenterField(jet_grid)
set!(fluid, 1)
volume = Field(Integral(fluid))
total(integral) = first(on_architecture(CPU(), interior(compute!(integral))))
initial_volume = total(volume)

run!(jet_simulation)

@info @sprintf("The relative change in volume after %s is %.1e", prettytime(jet_simulation),
               total(volume) / initial_volume - 1)

[ Info: Initializing simulation...
[ Info: iteration: 0, time: 0 seconds, max|u|: 79.77 m/s, max|w|: 2.6e-06 m/s, wall time: 37.598 seconds
[ Info:     ... simulation initialization complete (14.671 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (41.387 seconds).
[ Info: iteration: 1000, time: 11.111 hours, max|u|: 78.69 m/s, max|w|: 2.4e-06 m/s, wall time: 46.071 seconds
[ Info: iteration: 2000, time: 22.222 hours, max|u|: 77.75 m/s, max|w|: 2.8e-06 m/s, wall time: 2.930 seconds
[ Info: iteration: 3000, time: 1.389 days, max|u|: 76.80 m/s, max|w|: 2.6e-06 m/s, wall time: 2.919 seconds
[ Info: iteration: 4000, time: 1.852 days, max|u|: 75.79 m/s, max|w|: 2.7e-06 m/s, wall time: 3.032 seconds
[ Info: iteration: 5000, time: 2.315 days, max|u|: 75.12 m/s, max|w|: 2.5e-06 m/s, wall time: 2.948 seconds
[ Info: iteration: 6000, time: 2.778 days, max|u|: 74.68 m/s, max|w|: 2.8e-06 m/s, wall time: 2.979 seconds
[ Info: iteration: 7000, time: 3.241 days, max|u|: 73.51 m/s, max|w|: 2.6e-06 m/s, wall time: 2.932 seconds
[ Info: iteration: 8000, time: 3.704 days, max|u|: 72.81 m/s, max|w|: 2.3e-06 m/s, wall time: 2.862 seconds
[ Info: iteration: 9000, time: 4.167 days, max|u|: 71.35 m/s, max|w|: 2.4e-06 m/s, wall time: 2.967 seconds
[ Info: iteration: 10000, time: 4.630 days, max|u|: 69.80 m/s, max|w|: 2.3e-06 m/s, wall time: 2.882 seconds
[ Info: iteration: 11000, time: 5.093 days, max|u|: 69.23 m/s, max|w|: 2.4e-06 m/s, wall time: 2.876 seconds
[ Info: iteration: 12000, time: 5.556 days, max|u|: 70.60 m/s, max|w|: 2.8e-06 m/s, wall time: 2.890 seconds
[ Info: Simulation is stopping after running for 1.589 minutes.
[ Info: Simulation time 6 days equals or exceeds stop time 6 days.
[ Info: The relative change in volume after 6 days is -1.2e-07

Visualization

We plot the vorticity on the sphere, seen from opposite sides, while the camera circles the globe once. Both views show Yin and Yang together. Each cell is drawn as a quad of one color between the corners of its panel, so the plot shows the solution as the panels hold it. Yin draws the cells it owns, and Yang, just beneath, the cells in which it owns any area, so that the overlap appears once.

julia

sphere_point(λ, φ, r) = Point3f(r * cosd(φ) * cosd(λ), r * cosd(φ) * sind(λ), r * sind(φ))

function cell_quads(panel, drawn; radius=1)
    Nz = size(panel, 3)
    corner(i, j) = sphere_point(λnode(i, j, Nz, panel, Face(), Face(), Center()),
                                φnode(i, j, Nz, panel, Face(), Face(), Center()), radius)
    cells = findall(drawn)
    vertices = [corner(c[1] + di, c[2] + dj) for c in cells for (di, dj) in ((0, 0), (1, 0), (1, 1), (0, 1))]
    return vertices, cells
end

cell_quads (generic function with 1 method)

CairoMakie orders the triangles of each plot by depth, but not the plots themselves. We therefore draw both panels as one mesh, with two triangles per quad, and give the four corners of each quad the value of its cell. Cells that land marks on a panel get NaN, which Double Drake uses below to leave its continents gray.

julia

function panel_mesh(grid; land=nothing)
    yin = on_architecture(CPU(), grid[1])
    yang = on_architecture(CPU(), grid[2])
    weights = owned_weights(grid)
    yin_vertices, yin_cells = cell_quads(yin, Array(interior(weights[1], :, :, 1)) .≥ 1/2)
    yang_vertices, yang_cells = cell_quads(yang, Array(interior(weights[2], :, :, 1)) .> 0; radius=0.995)

    vertices = vcat(yin_vertices, yang_vertices)
    quads = 1:length(vertices) ÷ 4
    triangles = [[4q-3 for q in quads] [4q-2 for q in quads] [4q-1 for q in quads];
                 [4q-3 for q in quads] [4q-1 for q in quads] [4q   for q in quads]]

    cell_values(yin, yang) = vcat(yin[yin_cells], yang[yang_cells])
    masked = isnothing(land) ? falses(length(yin_cells) + length(yang_cells)) : cell_values(land(yin), land(yang))
    on_quads(yin, yang) = repeat(ifelse.(masked, NaN, cell_values(yin, yang)), inner=4)

    return (; vertices, triangles, on_quads)
end

jet_mesh = panel_mesh(jet_grid)

The box of an Axis3 with viewmode=:fit keeps its size as the camera moves. We fit the box inside the sphere, so that the sphere fills the axis.

julia

function globe!(position, azimuth, elevation; width=400)
    ax = Axis3(position; width, height=width, aspect=:data, azimuth, elevation,
               limits=(-0.6, 0.6, -0.6, 0.6, -0.6, 0.6), viewmode=:fit, clip=false, protrusions=0)
    hidedecorations!(ax)
    hidespines!(ax)
    return ax
end

jet_frame = Observable(1)
ζₙ = @lift 1e5 * jet_mesh.on_quads(vorticity.yin[$jet_frame], vorticity.yang[$jet_frame])
azimuth = @lift 2π * ($jet_frame - 1) / length(jet_times)
title = @lift "Vertical vorticity, t = " * prettytime(jet_times[$jet_frame])

colormap = :balance
colorrange = (-15, 15)

jet_fig = Figure()
Label(jet_fig[0, 1:2], title, fontsize=20, tellwidth=false)

sides = ((1, "One side", 0),
         (2, "The opposite side", π))

for (column, name, offset) in sides
    Label(jet_fig[1, column], name, fontsize=16, tellwidth=false)
    ax = globe!(jet_fig[2, column], @lift($azimuth + offset), π/4; width=450)
    mesh!(ax, jet_mesh.vertices, jet_mesh.triangles; color=ζₙ, colormap, colorrange, shading=NoShading)
end

Colorbar(jet_fig[3, 1:2]; colormap, colorrange, label="Vertical vorticity (10⁻⁵ s⁻¹)",
         vertical=false, width=500)

resize_to_layout!(jet_fig)

CairoMakie.record(jet_fig, "yin_yang_barotropic_jet.mp4", eachindex(jet_times), framerate=8) do i
    jet_frame[] = i
end

By day 4 the jet meanders, and by day 6 it has rolled up into a chain of vortices.

Double Drake

The second simulation has wind-driven, baroclinically unstable fronts in an ocean that covers the whole sphere. Two thin continents run from the North Pole to 35°S along 0° and 90°E. They divide the ocean into a small and a large basin, which connect through a circumpolar channel in the south. This is the "Double Drake" geometry of , an idealized Atlantic and Pacific with a Drake Passage.

The grid

To resolve baroclinic eddies with a modest number of cells, we use a small planet with a radius of 1000 km that rotates at Earth's rate. Each panel has 64 cells in latitude and 192 in longitude, a spacing of 1.4°, or 25 km. Four layers span 4 km of ocean. We build the panels with YinYangGrid, and use a MutableVerticalDiscretization for a vertical coordinate.

julia

radius = 1000kilometers
z_faces = [-4000, -1500, -500, -100, 0]
Nz = length(z_faces) - 1

underlying_grid = YinYangGrid(arch; size=(192, 64, Nz), halo=(5, 5, 4), radius,
                              z=MutableVerticalDiscretization(z_faces))

200×72×4 YinYangGrid{Float32, Connected, Connected, Bounded} on CUDAGPU with 5×5×4 halo
├── core: 192×64 cells per panel spanning 270° × 90° with Δλ′ = Δφ′ = 1.40625°
├── overlap: 4 cells (5.625°) beyond the core on every side, 4×4 interpolation stencils
├── panels: λ′ ∈ [-140.625, 140.625], φ′ ∈ [-50.625, 50.625]
├── Yin (region 1):  panel frame = geographic frame
├── Yang (region 2): (x, y, z) ↦ (-x, z, y), panel (λ′, φ′) = (0, 0) at (λ, φ) = (180, 0)
└── z: Bounded  z ∈ [-4000.0, 0.0] variably and mutably spaced with min(Δr)=100.0, max(Δr)=2500.0

The continents are 16° wide strips along 0° and 90°E. They join a polar cap north of 82°N and end at 35°S. Their width is measured from the plane of each meridian, so they do not thin toward the pole.

julia

function continent(λ, φ)
    x, y = cosd(φ) * cosd(λ), cosd(φ) * sind(λ)
    along_0 = (abs(y) < sind(8)) & (x > 0)
    along_90E = (abs(x) < sind(8)) & (y > 0)
    return (φ > 82) | ((along_0 | along_90E) & (φ > -35))
end

bottom_height(λ, φ) = ifelse(continent(λ, φ), 0, -4000)

drake_grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height))

YinYangGrid{Float32, Connected, Connected, Bounded} partitioned on CUDAGPU: 
├── region_grids: 200×72×4 RotatedLatitudeLongitudeGrid{Float32, Connected, Connected, Bounded} on CUDAGPU with 5×5×4 halo 
├── partition: YinYangPartition with (1: Yin, 2: Yang) 
└── connectivity: YinYangConnectivity with 192×64 core, 4-cell overlap and 4×4 interpolation stencils

Wind stress

Westerlies blow over the mid-latitudes and trade winds over the tropics, with the strongest winds over the Southern Ocean:

with .

julia

@inline function zonal_wind_stress(φ)
    westerlies = 2exp(-((φ + 50) / 12)^2) + exp(-((φ - 45) / 12)^2)
    trades = 0.7exp(-((φ + 18) / 12)^2) + 0.6exp(-((φ - 18) / 12)^2)
    return 0.15 * cos(deg2rad(φ)) * (westerlies - trades)
end

zonal_wind_stress (generic function with 1 method)

Each panel stores velocities in its own basis, which is rotated by an angle from the geographic basis. We rotate the eastward stress into the panel basis. The momentum flux is minus the stress divided by the reference density.

julia

@inline function u_wind_flux(i, j, grid, clock, fields, ρ₀)
    φ = φnode(i, j, grid.Nz, grid, Face(), Center(), Center())
    θ = rotation_angle(i, j, grid, Face(), Center())
    return - zonal_wind_stress(φ) * cos(θ) / ρ₀
end

@inline function v_wind_flux(i, j, grid, clock, fields, ρ₀)
    φ = φnode(i, j, grid.Nz, grid, Center(), Face(), Center())
    θ = rotation_angle(i, j, grid, Center(), Face())
    return - zonal_wind_stress(φ) * sin(θ) / ρ₀
end

ρ₀ = 1020 # kg m⁻³
u_wind = FluxBoundaryCondition(u_wind_flux; discrete_form=true, parameters=ρ₀)
v_wind = FluxBoundaryCondition(v_wind_flux; discrete_form=true, parameters=ρ₀)

FluxBoundaryCondition: DiscreteBoundaryFunction v_wind_flux with parameters 1020

A quadratic BulkDrag acts on the sea floor.

julia

drag = BulkDrag(coefficient=2.5e-3)
u_bcs = FieldBoundaryConditions(top=u_wind, bottom=drag, immersed=ImmersedBoundaryCondition(bottom=drag))
v_bcs = FieldBoundaryConditions(top=v_wind, bottom=drag, immersed=ImmersedBoundaryCondition(bottom=drag))

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── bottom: FluxBoundaryCondition: BulkDragFunction(QuadraticFormulation(), Nothing, Cᴰ=0.0025)
├── top: FluxBoundaryCondition: DiscreteBoundaryFunction v_wind_flux with parameters 1020
└── immersed: ImmersedBoundaryCondition with west=Nothing, east=Nothing, south=Nothing, north=Nothing, bottom=Flux, top=Nothing

Fronts in both hemispheres

A warm thermocline fills the top kilometer between two fronts at , over a cold deep ocean and cold polar water:

with C, C, km and . The meanders , with amplitude and , seed the instability.

The surface layer is restored to the unperturbed fronts, with , on a time scale of 30 days.

julia

front = (; Td=2, ΔT=20, h=1000, φ₀=57, δ=0.1)

@inline function temperature_flux(i, j, grid, clock, fields, p)
    φ = φnode(i, j, grid.Nz, grid, Center(), Center(), Center())
    z = znode(i, j, grid.Nz, grid, Center(), Center(), Center())
    χ = (sin(deg2rad(φ))^2 - sin(deg2rad(p.φ₀))^2) / p.δ
    T★ = p.Td + p.ΔT / 2 * exp(z / p.h) * (1 - tanh(χ))
    return @inbounds p.rate * (fields.T[i, j, grid.Nz] - T★)
end

restoring_rate = 100 / 30days # surface layer thickness over the restoring time [m s⁻¹]
restoring_parameters = merge(front, (; rate=restoring_rate))
T_restoring = FluxBoundaryCondition(temperature_flux; discrete_form=true, parameters=restoring_parameters)
T_bcs = FieldBoundaryConditions(top=T_restoring)

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── top: FluxBoundaryCondition: DiscreteBoundaryFunction temperature_flux with parameters (Td=2, ΔT=20, h=1000, φ₀=57, δ=0.1, rate=3.85802e-5)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)

The model

Convective adjustment mixes temperature and momentum where the stratification is unstable, and a weak background diffusivity and viscosity act everywhere. Both are integrated implicitly in the vertical.

julia

vertical_mixing = ConvectiveAdjustmentVerticalDiffusivity(convective_κz=1, convective_νz=1,
                                                          background_κz=1e-5, background_νz=1e-3)

horizontal_mixing = HorizontalScalarDiffusivity(ν=100, κ=20)

α = 2e-4 # thermal expansion coefficient [K⁻¹]
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=α), constant_salinity=35)

SeawaterBuoyancy{Float32}:
├── gravitational_acceleration: 9.80665
├── constant_salinity: 35.0
└── equation_of_state: LinearEquationOfState(thermal_expansion=0.0002, haline_contraction=0.00078)

Where the fronts meet the continents, they drive coastal currents of about 6 m s⁻¹. With WENO alone, the sharp fronts these currents carry would undershoot . The initial state and the restoring lie between and , so we bound the advection to that range.

On a Yin-Yang grid, the split-explicit free surface substeps without extending its halos.

julia

drake_model = HydrostaticFreeSurfaceModel(drake_grid; buoyancy,
                                          coriolis = HydrostaticSphericalCoriolis(),
                                          closure = (vertical_mixing, horizontal_mixing),
                                          free_surface = SplitExplicitFreeSurface(substeps=12, extend_halos=false),
                                          momentum_advection = WENOVectorInvariant(order=5),
                                          tracer_advection = WENO(order=5, bounds=(front.Td, front.Td + front.ΔT)),
                                          tracers = :T,
                                          timestepper = :SplitRungeKutta3,
                                          vertical_coordinate = ZStarCoordinate(),
                                          boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs))

HydrostaticFreeSurfaceModel{CUDAGPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── grid: 200×72×4 YinYangGrid{Float32, Connected, Connected, Bounded} on CUDAGPU with 5×5×4 halo
├── timestepper: SplitRungeKuttaTimeStepper
├── tracers: T
├── closure: Tuple with 2 closures:
│   ├── ConvectiveAdjustmentVerticalDiffusivity{VerticallyImplicitTimeDiscretization}(background_κz=1.0e-5 convective_κz=1.0 background_νz=0.001 convective_νz=1.0)
│   └── HorizontalScalarDiffusivity{ExplicitTimeDiscretization}(ν=100.0, κ=(T=20.0,))
├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.0002, haline_contraction=0.00078) with ĝ = NegativeZDirection()
├── free surface: SplitExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── substepping: FixedSubstepNumber(8)
├── advection scheme: 
│   ├── momentum: MultiRegionObject{Tuple{WENOVectorInvariant{3, 3, Float32, ExplicitTimeDiscretization, WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, WENO{2, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, UpwindBiased{1, Float32, ExplicitTimeDiscretization, Nothing, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Oceananigans.Advection.VelocityStencil, WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, WENO{2, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, UpwindBiased{1, Float32, ExplicitTimeDiscretization, Nothing, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, WENO{2, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, UpwindBiased{1, Float32, ExplicitTimeDiscretization, Nothing, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, WENO{2, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, UpwindBiased{1, Float32, ExplicitTimeDiscretization, Nothing, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Oceananigans.Advection.OnlySelfUpwinding{Centered{2, Float32, ExplicitTimeDiscretization, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.divergence_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.divergence_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.u_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.v_smoothness)}}, false}, WENOVectorInvariant{3, 3, Float32, ExplicitTimeDiscretization, WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, WENO{2, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, UpwindBiased{1, Float32, ExplicitTimeDiscretization, Nothing, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Oceananigans.Advection.VelocityStencil, WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, WENO{2, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, UpwindBiased{1, Float32, ExplicitTimeDiscretization, Nothing, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, WENO{2, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, UpwindBiased{1, Float32, ExplicitTimeDiscretization, Nothing, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, WENO{2, Float32, Oceananigans.Utils.BackendOptimizedDivision, ExplicitTimeDiscretization, Nothing, UpwindBiased{1, Float32, ExplicitTimeDiscretization, Nothing, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Oceananigans.Advection.OnlySelfUpwinding{Centered{2, Float32, ExplicitTimeDiscretization, Centered{1, Float32, ExplicitTimeDiscretization, Nothing}}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.divergence_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.divergence_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.u_smoothness)}, Oceananigans.Advection.FunctionStencil{typeof(Oceananigans.Advection.v_smoothness)}}, false}}}
│   └── T: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5, bounds=(2.0, 22.0))
├── vertical_coordinate: ZStarCoordinate
└── coriolis: HydrostaticSphericalCoriolis{Oceananigans.Advection.EnstrophyConserving{Float32}, Float32, Oceananigans.Coriolis.HydrostaticFormulation}

A balanced initial state

The fronts start in balance with an initial current. The deepest layer is at rest: the free surface cancels the bottom pressure gradient of the warm water above it, so that the pressure anomaly at each depth is

where is the center of the deepest cell. We integrate with the same trapezoidal rule between cell centers as the model's hydrostatic pressure. The zonal velocity is in gradient-wind balance with , and the meridional velocity, which only the meanders induce, is geostrophic:

where, again, and is the radius of the sphere. With , the meridional pressure gradient is , where

The root of the quadratic for that reduces to geostrophic balance for weak flow, and the geostrophic , are

julia

A = g * α * front.ΔT / 2

zᶜ = on_architecture(CPU(), znodes(drake_grid, Center()))
eᶜ = exp.(zᶜ ./ front.h)
I(z) = sum((eᶜ[k] + eᶜ[k+1]) / 2 * clamp(z - zᶜ[k], 0, zᶜ[k+1] - zᶜ[k]) for k in 1:Nz-1)

I (generic function with 1 method)

The pressure at the top cell center includes buoyancy integrated over the distance from that center to its mirror image above the surface:

julia

Iˢ = I(0) - 2zᶜ[end] * eᶜ[end]

function initial_state(λ, φ, z; ϵ=0.035, φ₀=front.φ₀, δ=front.δ)
    s, c = sind(φ), cosd(φ)
    D = sind(φ₀)^2 * cosd(φ₀)^2
    λ = deg2rad(λ)
    P = cos(6λ) + cos(9λ) / 2 + sin(13λ) / 4
    ∂λP = - 6sin(6λ) - 9sin(9λ) / 2 + 13cos(13λ) / 4
    m = s^2 * c^2 / D
    dm = (c^2 - s^2) / D  ## dm/d(sin²φ)
    χ = (s^2 - sind(φ₀)^2 - ϵ * m * P) / δ

    a = A * I(z) * sech(χ)^2 / (Ω * radius * δ)
    U = a * (1 - ϵ * dm * P)
    u = 2c * U / (1 + sqrt(1 + 2U / (Ω * radius)))
    v = a * ϵ * s * c * ∂λP / (2D)
    T = front.Td + front.ΔT / 2 * exp(z / front.h) * (1 - tanh(χ))
    η = - A * Iˢ * tanh(χ) / g

    return (; u, v, T, η)
end

initial_state (generic function with 1 method)

set! takes the geographic components of the velocity and rotates them into each panel's basis.

julia

set!(drake_model, u = (λ, φ, z) -> initial_state(λ, φ, z).u,
                  v = (λ, φ, z) -> initial_state(λ, φ, z).v,
                  T = (λ, φ, z) -> initial_state(λ, φ, z).T,
                  η = (λ, φ, z) -> initial_state(λ, φ, z).η)

The simulation

We run for twenty days with a time step of 150 seconds.

julia

drake_simulation = Simulation(drake_model; Δt=150, stop_time=20days)

Simulation of HydrostaticFreeSurfaceModel{CUDAGPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 2.500 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 20 days
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

The panels disagree slightly on the fluxes between them, so a conservation fixer removes the volume and heat they add.

julia

add_conservation_fixer!(drake_simulation)

add_callback!(drake_simulation, progress, IterationInterval(1000))

For plotting, we save the surface vorticity and temperature of each panel every three hours.

julia

drake_ζ = vertical_vorticity(drake_model)

drake_times = Float64[]
surface_ζ = (yin = Matrix{Float64}[], yang = Matrix{Float64}[])
surface_T = (yin = Matrix{Float64}[], yang = Matrix{Float64}[])

function save_surface_fields(sim)
    compute!(drake_ζ)
    T = sim.model.tracers.T
    push!(drake_times, time(sim))
    push!(surface_ζ.yin, Array(interior(drake_ζ[1], :, :, Nz)))
    push!(surface_ζ.yang, Array(interior(drake_ζ[2], :, :, Nz)))
    push!(surface_T.yin, Array(interior(T[1], :, :, Nz)))
    push!(surface_T.yang, Array(interior(T[2], :, :, Nz)))
    return nothing
end

add_callback!(drake_simulation, save_surface_fields, TimeInterval(3hours))

We also track the heat content, , which changes only through surface restoring, since the conservation fixer removes the small errors of the interpolation between the panels.

julia

heat_content = Field(Integral(drake_model.tracers.T))
initial_heat_content = total(heat_content)

2.6694637f17

We restart the wall clock and run.

julia

wall_clock[] = time_ns()

run!(drake_simulation)

@info @sprintf("The heat content changed by %.3f%% in %s", 100 * (total(heat_content) / initial_heat_content - 1),
               prettytime(drake_simulation))

[ Info: Initializing simulation...
[ Info: iteration: 0, time: 0 seconds, max|u|: 1.80 m/s, max|w|: 5.5e-02 m/s, wall time: 24.788 seconds
[ Info:     ... simulation initialization complete (17.084 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (9.288 minutes).
[ Info: iteration: 1000, time: 1.736 days, max|u|: 6.17 m/s, max|w|: 7.2e-02 m/s, wall time: 9.644 minutes
[ Info: iteration: 2000, time: 3.472 days, max|u|: 5.54 m/s, max|w|: 9.9e-02 m/s, wall time: 25.706 seconds
[ Info: iteration: 3000, time: 5.208 days, max|u|: 5.07 m/s, max|w|: 3.2e-02 m/s, wall time: 27.357 seconds
[ Info: iteration: 4000, time: 6.944 days, max|u|: 4.79 m/s, max|w|: 2.5e-02 m/s, wall time: 20.300 seconds
[ Info: iteration: 5000, time: 8.681 days, max|u|: 4.56 m/s, max|w|: 2.6e-02 m/s, wall time: 22.719 seconds
[ Info: iteration: 6000, time: 10.417 days, max|u|: 4.86 m/s, max|w|: 2.6e-02 m/s, wall time: 17.485 seconds
[ Info: iteration: 7000, time: 12.153 days, max|u|: 4.98 m/s, max|w|: 2.7e-02 m/s, wall time: 17.318 seconds
[ Info: iteration: 8000, time: 13.889 days, max|u|: 4.98 m/s, max|w|: 3.5e-02 m/s, wall time: 20.415 seconds
[ Info: iteration: 9000, time: 15.625 days, max|u|: 4.67 m/s, max|w|: 2.3e-02 m/s, wall time: 17.456 seconds
[ Info: iteration: 10000, time: 17.361 days, max|u|: 4.66 m/s, max|w|: 2.0e-02 m/s, wall time: 17.298 seconds
[ Info: iteration: 11000, time: 19.097 days, max|u|: 4.08 m/s, max|w|: 1.5e-02 m/s, wall time: 19.074 seconds
[ Info: Simulation is stopping after running for 13.473 minutes.
[ Info: Simulation time 20 days equals or exceeds stop time 20 days.
[ Info: The heat content changed by -0.110% in 20 days

Visualization

We plot the surface fields on the sphere, seen from the north and from the south, while the camera circles the globe once, with the same mesh of quads as for the jet. Land is gray.

julia

drake_mesh = panel_mesh(drake_grid; land = panel -> interior(bottom_height_field(panel), :, :, 1) .≥ 0)

drake_frame = Observable(1)
ζₙ = @lift 1e5 * drake_mesh.on_quads(surface_ζ.yin[$drake_frame], surface_ζ.yang[$drake_frame])
Tₙ = @lift drake_mesh.on_quads(surface_T.yin[$drake_frame], surface_T.yang[$drake_frame])
azimuth = @lift π/4 + 2π * ($drake_frame - 1) / length(drake_times)
title = @lift "Surface vorticity and temperature, t = " * prettytime(drake_times[$drake_frame])

drake_fig = Figure()
Label(drake_fig[0, 1:2], title, fontsize=20, tellwidth=false)

views = ((1, "From the north", +π/6),
         (2, "From the south", -π/6))

for (column, name, elevation) in views
    Label(drake_fig[1, column], name, fontsize=16, tellwidth=false)
end

rows = ((ζₙ, :balance, (-10, 10), "Vertical vorticity (10⁻⁵ s⁻¹)"),
        (Tₙ, :thermal, (0, 24), "Temperature (ᵒC)"))

for (row, (values, colormap, colorrange, label)) in enumerate(rows)
    for (column, name, elevation) in views
        ax = globe!(drake_fig[row + 1, column], azimuth, elevation)
        mesh!(ax, drake_mesh.vertices, drake_mesh.triangles; color=values, colormap, colorrange, nan_color=:gray70, shading=NoShading)
    end
    Colorbar(drake_fig[row + 1, 3]; colormap, colorrange, label)
end

resize_to_layout!(drake_fig)

CairoMakie.record(drake_fig, "yin_yang_double_drake.mp4", eachindex(drake_times), framerate=12) do i
    drake_frame[] = 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.13.2
Commit 45b6d15e34b (2026-10-14 20:33 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.2+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.2+0.x64.linux.gnu/share/julia
  JULIA_VERSION_ENZYME = 1.13.1
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  LD_LIBRARY_PATH = 
  JULIA_MAX_NUM_PRECOMPILE_FILES = 24
  JULIA_VERSION = 1.13.2
  JULIA_CUDA_USE_COMPAT = false
  JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-35332/docs/
  JULIA_DEBUG = Literate

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

julia

import Pkg
Pkg.status()

Status `~/Oceananigans.jl-35332/docs/Project.toml`
  [79e6a3ab] Adapt v4.7.4
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.16
⌃ [8e50ac2c] ConservativeRegridding v0.2.17
⌅ [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.115.0 `..`
  [429524aa] Optim v2.3.2
  [f27b6e38] Polynomials v4.1.4
⌅ [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.

`NonhydrostaticModel`

`HydrostaticFreeSurfaceModel`

`ShallowWaterModel`

Barotropic jet and Double Drake on a Yin-Yang grid

The Yin-Yang grid

Install dependencies

Barotropic jet

A single layer of fluid

The jet

The simulation

Visualization

Double Drake

The grid

Wind stress

Fronts in both hemispheres

The model

A balanced initial state

The simulation

Visualization

Julia version and environment information

Barotropic jet and Double Drake on a Yin-Yang grid ​

The Yin-Yang grid ​

Install dependencies ​

Barotropic jet ​

A single layer of fluid ​

The jet ​

The simulation ​

Visualization ​

Double Drake ​

The grid ​

Wind stress ​

Fronts in both hemispheres ​

The model ​

A balanced initial state ​

The simulation ​

Visualization ​

Julia version and environment information ​

Barotropic jet and Double Drake on a Yin-Yang grid

The Yin-Yang grid

Install dependencies

Barotropic jet

A single layer of fluid

The jet

The simulation

Visualization

Double Drake

The grid

Wind stress

Fronts in both hemispheres

The model

A balanced initial state

The simulation

Visualization

Julia version and environment information