Inertia-gravity waves

This example simulates the propagation of inertia-gravity waves in a stably stratified atmosphere, following the classical benchmark test case described by Skamarock and Klemp (1994). This test evaluates the accuracy of numerical pressure solvers by introducing a small-amplitude temperature perturbation into a stratified environment with constant Brunt-Väisälä frequency, triggering propagating inertia-gravity waves.

The test case is particularly useful for validating anelastic and compressible solvers, as discussed at the CM1 inertia-gravity wave test page.

Physical setup

The background state is a stably stratified atmosphere with constant Brunt-Väisälä frequency $N$, which gives a potential temperature profile

\[θ^{\rm bg}(z) = θ_0 \exp\left( \frac{N^2 z}{g} \right)\]

where $θ_0 = 300 \, {\rm K}$ is the surface potential temperature and $g$ is the gravitational acceleration.

The initial perturbation is a localized temperature anomaly centered at $x = x_0$:

\[θ'(x, z) = Δθ \frac{\sin(π z / L_z)}{1 + (x - x_0)^2 / a^2}\]

with amplitude $Δθ = 0.01 \, {\rm K}$, half-width parameter $a = 5000 \, {\rm m}$, and perturbation center $x_0 = L_x / 3$. A uniform mean wind $U = 20 \, {\rm m \, s^{-1}}$ advects the waves.

using Breeze
using Oceananigans.Units
using Statistics
using Printf
using CairoMakie

Problem parameters

We define the thermodynamic base state and mean wind following Skamarock and Klemp (1994):

p₀ = 100000  # Pa - surface pressure
θ₀ = 300     # K - reference potential temperature
U  = 20      # m s⁻¹ - mean wind
N  = 0.01    # s⁻¹ - Brunt-Väisälä frequency
N² = N^2

0.0001

Grid configuration

The domain is 300 km × 10 km with 300 × 10 grid points, matching the nonhydrostatic case configuration in the paper by Skamarock and Klemp (1994).

Nx, Nz = 300, 10
Lx, Lz = 300kilometers, 10kilometers

grid = RectilinearGrid(CPU(), size = (Nx, Nz), halo = (5, 5),
                       x = (0, Lx), z = (0, Lz),
                       topology = (Periodic, Flat, Bounded))

300×1×10 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 5×0×5 halo
├── Periodic x ∈ [0.0, 300000.0) regularly spaced with Δx=1000.0
├── Flat y                       
└── Bounded  z ∈ [0.0, 10000.0]  regularly spaced with Δz=1000.0

Atmosphere model setup

We use the anelastic formulation with liquid-ice potential temperature thermodynamics:

constants = ThermodynamicConstants()
reference_state = ReferenceState(grid, constants; surface_pressure=p₀, potential_temperature=θ₀)
dynamics = AnelasticDynamics(reference_state)
advection = WENO(minimum_buffer_upwind_order=3)
model = AtmosphereModel(grid; dynamics, advection)

AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 300×1×10 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 5×0×5 halo
├── dynamics: AnelasticDynamics(p₀=100000.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: 
│   ├── momentum: WENO{3, Float64, Float32}(order=5)
│   ├── ρθ: WENO{3, Float64, Float32}(order=5)
│   └── ρqᵗ: WENO{3, Float64, Float32}(order=5)
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

Initial conditions

The perturbation parameters from the paper by Skamarock and Klemp (1994):

Δθ = 0.01               # K - perturbation amplitude
a  = 5000               # m - perturbation half-width parameter
x₀ = Lx / 3             # m - perturbation center in x

100000.0

The background potential temperature profile with a constant Brunt-Väisälä frequency:

g = model.thermodynamic_constants.gravitational_acceleration
θᵇᵍ(z) = θ₀ * exp(N² * z / g)

θᵇᵍ (generic function with 1 method)

The initial condition combines the background profile with the localized perturbation:

θᵢ(x, z) = θᵇᵍ(z) + Δθ * sin(π * z / Lz) / (1 + (x - x₀)^2 / a^2)

set!(model, θ=θᵢ, u=U)

Simulation

We run for 3000 seconds with a fixed time step equal 24 seconds, matching the simulation time in (Skamarock and Klemp, 1994):

simulation = Simulation(model; Δt=24, stop_time=3000)

Simulation of AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 24 seconds
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 50 minutes
├── 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

Progress callback:

θ = PotentialTemperature(model)
θᵇᵍf = CenterField(grid)
set!(θᵇᵍf, (x, z) -> θᵇᵍ(z))
θ′ = θ - θᵇᵍf

function progress(sim)
    u, v, w = sim.model.velocities
    msg = @sprintf("Iter: % 4d, t: % 14s, max(θ′): %.4e, max|w|: %.4f",
                   iteration(sim), prettytime(sim), maximum(θ′), maximum(abs, w))
    @info msg
    return nothing
end

add_callback!(simulation, progress, IterationInterval(20))

Output

We save the potential temperature for visualization, including an animation of the wave propagation:

outputs = merge(model.velocities, (; θ′))

filename = "inertia_gravity_wave.jld2"
simulation.output_writers[:jld2] = JLD2Writer(model, outputs; filename,
                                              schedule = TimeInterval(100),
                                              overwrite_existing = true)

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter:    0, t:      0 seconds, max(θ′): 9.7791e-03, max|w|: 0.0000
[ Info:     ... simulation initialization complete (33.350 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.811 seconds).
[ Info: Iter:   20, t:  6.667 minutes, max(θ′): 4.5377e-03, max|w|: 0.0060
[ Info: Iter:   40, t: 13.333 minutes, max(θ′): 3.8578e-03, max|w|: 0.0049
[ Info: Iter:   60, t:     20 minutes, max(θ′): 3.5220e-03, max|w|: 0.0042
[ Info: Iter:   80, t: 26.667 minutes, max(θ′): 3.2983e-03, max|w|: 0.0036
[ Info: Iter:  100, t: 33.333 minutes, max(θ′): 3.1264e-03, max|w|: 0.0032
[ Info: Iter:  120, t:     40 minutes, max(θ′): 2.9939e-03, max|w|: 0.0029
[ Info: Iter:  140, t: 46.667 minutes, max(θ′): 2.8867e-03, max|w|: 0.0027
[ Info: Simulation is stopping after running for 39.006 seconds.
[ Info: Simulation time 50 minutes equals or exceeds stop time 50 minutes.

Results: potential temperature perturbation

Following Skamarock and Klemp (1994), we visualize the potential temperature perturbation $θ' = θ - θ^{\rm bg}$. The final state at $t = 3000 \, {\rm s}$ can be compared directly to Figure 3b in the paper by Skamarock and Klemp (1994), which shows the analytic solution for incompressible flow.

The CM1 model test page provides additional comparisons between compressible, anelastic, and incompressible solvers.

θ′t = FieldTimeSeries(filename, "θ′")
times = θ′t.times
Nt = length(times)

31

Plot the final potential temperature perturbation (compare to Figure 3b by Skamarock and Klemp (1994)):

θ′N = θ′t[Nt]

fig = Figure(size=(800, 300))
ax = Axis(fig[1, 1], xlabel = "x (m)", ylabel = "z (m)",
          title = "Potential temperature perturbation θ′ at t = $(prettytime(times[Nt]))")

levels = range(-Δθ/2, stop=Δθ/2, length=20)
contourf!(ax, θ′N, colormap=:balance; levels)

fig

Animation of wave propagation

The animation shows the evolution of the potential temperature perturbation as the inertia-gravity waves propagate away from the initial disturbance:

n = Observable(1)
θ′n = @lift θ′t[$n]

fig = Figure(size=(800, 300))
ax = Axis(fig[1, 1], xlabel = "x (m)", ylabel = "z (m)")

title = @lift "Potential temperature perturbation θ′ at t = $(prettytime(times[$n]))"
fig[0, :] = Label(fig, title, fontsize=16, tellwidth=false)

heatmap!(ax, θ′n, colormap = :balance, colorrange = (-Δθ/2, Δθ/2))

record(fig, "inertia_gravity_wave.mp4", 1:Nt, framerate=8) do nn
    n[] = nn
end

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD EPYC 7R13 Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_GPG = 3673DF529D9049477F76B37566E3C7DC03D6E495
  JULIA_LOAD_PATH = :@breeze
  JULIA_VERSION = 1.12.2
  JULIA_DEPOT_PATH = /usr/local/share/julia:
  JULIA_PATH = /usr/local/julia
  JULIA_PROJECT = @breeze

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

import Pkg
Pkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
  [86bc3604] AtmosphericProfilesLibrary v0.1.7
  [660aa2fb] Breeze v0.2.1 `.`
  [052768ef] CUDA v5.9.5
  [13f3f980] CairoMakie v0.15.8
  [6a9e3e04] CloudMicrophysics v0.29.1
  [e30172f5] Documenter v1.16.1
  [daee34ce] DocumenterCitations v1.4.1
  [98b081ad] Literate v2.21.0
⌅ [9e8cae18] Oceananigans v0.102.5
  [a01a1ee8] RRTMGP v0.21.6
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.0
  [9a3f8284] Random v1.11.0
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.