Internal wave example

In this example, we initialize an internal wave packet in two-dimensions and watch it propagate. This example illustrates how to set up a two-dimensional model, set initial conditions, and how to use BackgroundFields.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, CairoMakie"

The physical domain

First, we pick a resolution and domain size. We use a two-dimensional domain that's periodic in and is Flat in :

using Oceananigans

grid = RectilinearGrid(size=(128, 128), x=(-π, π), z=(-π, π), topology=(Periodic, Flat, Periodic))

128×1×128 RectilinearGrid{Float64, Periodic, Flat, Periodic} on CPU with 3×0×3 halo
├── Periodic x ∈ [-3.14159, 3.14159) regularly spaced with Δx=0.0490874
├── Flat y                           
└── Periodic z ∈ [-3.14159, 3.14159) regularly spaced with Δz=0.0490874

Internal wave parameters

Inertia-gravity waves propagate in fluids that are both (i) rotating, and (ii) density-stratified. We use Oceananigans' Coriolis abstraction to implement a background rotation rate:

coriolis = FPlane(f=0.2)

FPlane{Float64}(f=0.2)

On an FPlane, the domain is idealized as rotating at a constant rate with rotation period . coriolis is passed to NonhydrostaticModel below. Our units are arbitrary.

We use Oceananigans' background_fields abstraction to define a background buoyancy field  , where is the vertical coordinate and is the "buoyancy frequency". This means that the modeled buoyancy field perturbs the basic state .

Background fields are functions of the coordinates (here x, z since y is Flat), time t, and optional parameters. Here we have one parameter, the buoyancy frequency.

N = 1
stratification(x, z, t, N) = N^2 * z
B = BackgroundField(stratification, parameters=N)

BackgroundField{typeof(Main.var"##8399".stratification), Int64}
├── func: stratification (generic function with 1 method)
└── parameters: 1

We are now ready to instantiate our model. We pass grid, coriolis, and B to the NonhydrostaticModel constructor. We add a small amount of ScalarDiffusivity to keep the model stable during time-stepping, and specify that we're using a single tracer called b that we identify as buoyancy by setting buoyancy=BuoyancyTracer().

model = NonhydrostaticModel(grid; coriolis,
                            advection = Centered(order=4),
                            closure = ScalarDiffusivity(ν=1e-6, κ=1e-6),
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            background_fields = (; b=B))

NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Periodic} on CPU with 3×0×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme:
│   ├── momentum: Centered(order=4)
│   └── b: Centered(order=4)
├── tracers: b
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1.0e-6, κ=(b=1.0e-6,))
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.2)

A Gaussian wavepacket

Next, we set up an initial condition that excites an internal wave that propagates through our rotating, stratified fluid. This internal wave has the pressure field

where is the vertical wavenumber, is the horizontal wavenumber, is the wave frequency, and is a Gaussian envelope. The internal wave dispersion relation links the wave numbers and , the Coriolis parameter , and the buoyancy frequency :

m = 16
k = 8
f = coriolis.f

ω = sqrt((N^2 * k^2 + f^2 * m^2) / (k^2 + m^2))

We define a Gaussian envelope with amplitude and width , centered at  , so that we can observe wave propagation.

A = 1e-9
δ = grid.Lx / 15

a(x, z) = A * exp(-(x^2 + z^2) / 2δ^2)

An inertia-gravity wave is a linear solution to the Boussinesq equations. In order that our initial condition excites an inertia-gravity wave, we initialize the velocity and buoyancy perturbation fields to be consistent with the pressure field      at  . These relations are sometimes called the "polarization relations". At  , the polarization relations yield

u₀(x, z) = a(x, z) * k * ω   / (ω^2 - f^2) * cos(k * x + m * z)
v₀(x, z) = a(x, z) * k * f   / (ω^2 - f^2) * sin(k * x + m * z)
w₀(x, z) = a(x, z) * m * ω   / (ω^2 - N^2) * cos(k * x + m * z)
b₀(x, z) = a(x, z) * m * N^2 / (ω^2 - N^2) * sin(k * x + m * z)

set!(model, u=u₀, v=v₀, w=w₀, b=b₀)

Recall that the buoyancy b is a perturbation, so that the total buoyancy field is  .

A wave packet on the loose

We're ready to release the packet. We build a simulation with a constant time-step,

simulation = Simulation(model, Δt = 0.1 * 2π/ω, stop_iteration = 20)

Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 1.304 seconds
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: Inf days
├── stop_iteration: 20.0
├── 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

and add an output writer that saves the velocity fields every iteration:

filename = "internal_wave.jld2"
simulation.output_writers[:velocities] = JLD2Writer(model, model.velocities; filename,
                                                    schedule = IterationInterval(1),
                                                    overwrite_files = true)

JLD2Writer scheduled on IterationInterval(1):
├── filepath: internal_wave.jld2
├── 3 outputs: (u, v, w)
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

With initial conditions set and an output writer at the ready, we run the simulation

run!(simulation)

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (1.399 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (935.275 ms).
[ Info: Simulation is stopping after running for 2.534 seconds.
[ Info: Model iteration 20 equals or exceeds stop iteration 20.

Animating a propagating packet

To animate the propagating wavepacket we just simulated, we load CairoMakie and make a Figure and an Axis for the animation,

using CairoMakie
set_theme!(Theme(fontsize = 20))

fig = Figure(size = (600, 600))

ax = Axis(fig[2, 1]; xlabel = "x", ylabel = "z",
          limits = ((-π, π), (-π, π)), aspect = AxisAspect(1))

Next, we load the vertical velocity w as a FieldTimeSeries and make contour plots of it. We use Makie's Observable to animate the data. To dive into how Observables work, refer to Makie.jl's Documentation.

n = Observable(1)

w_timeseries = FieldTimeSeries(filename, "w")
w = @lift w_timeseries[$n]
wmax = 1e-8

contourf!(ax, w;
          levels = range(-wmax, wmax, length=10),
          colormap = :balance,
          extendlow = :auto,
          extendhigh = :auto)

title = @lift "ωt = $(round(ω * w_timeseries.times[$n], digits=2))"
fig[1, 1] = Label(fig, title, fontsize=24, tellwidth=false)

Label()

And, finally, we record a movie.

frames = 1:length(w_timeseries.times)

@info "Animating a propagating internal wave..."

record(fig, "internal_wave.mp4", frames, framerate=8) do i
    n[] = i
end

[ Info: Animating a propagating internal wave...

---

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.