Thermal bubble

This example sets up, runs, and visualizes a "thermal bubble" (just a circular region of warm air) rising through a stably-stratified background.

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

A simple model on a RectilinearGrid

grid = RectilinearGrid(CPU(); size = (128, 128), halo = (5, 5),
                       x = (-10e3, 10e3), z = (0, 10e3),
                       topology = (Periodic, Flat, Bounded))

128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── Periodic x ∈ [-10000.0, 10000.0) regularly spaced with Δx=156.25
├── Flat y                           
└── Bounded  z ∈ [0.0, 10000.0]      regularly spaced with Δz=78.125

This example uses StaticEnergy thermodynamics, which is an alternative to the default LiquidIcePotentialTemperature thermodynamics. StaticEnergy is useful for dry simulations that don't require potential temperature diagnostics.

reference_state = ReferenceState(grid, ThermodynamicConstants(eltype(grid)))
dynamics = AnelasticDynamics(reference_state)
advection = WENO(order=9)
model = AtmosphereModel(grid; dynamics, formulation=:StaticEnergy, advection)

AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=288.0)
├── formulation: StaticEnergyFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: WENO{5, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│   ├── ρe: WENO{5, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│   └── ρqᵛ: WENO{5, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
├── forcing: @NamedTuple{ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρe::Returns{Float64}, ρqᵛ::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

Moist static energy perturbation

We add a localized potential temperature perturbation that translates into a moist static energy anomaly.

r₀ = 2e3
Δθ = 10 # K
N² = 1e-6
θ₀ = model.dynamics.reference_state.potential_temperature
g = model.thermodynamic_constants.gravitational_acceleration

function θᵢ(x, z;
            x₀ = mean(xnodes(grid, Center())),
            z₀ = 0.3*grid.Lz,
            N² = N²)

    θ̄ = θ₀ * exp(N² * z / g)
    r = sqrt((x - x₀)^2 + (z - z₀)^2)
    θ′ = Δθ * max(0, 1 - r / r₀)

    return θ̄ + θ′
end

set!(model, θ = θᵢ)

ρe = static_energy_density(model)
ρE = Field(Average(ρe, dims=1))
ρe′ = Field(ρe - ρE)

128×1×128 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
├── grid: 128×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: BinaryOperation at (Center, Center, Center)
├── status: time=0.0
└── data: 138×1×138 OffsetArray(::Array{Float64, 3}, -4:133, 1:1, -4:133) with eltype Float64 with indices -4:133×1:1×-4:133
    └── max=7288.25, min=-850.782, mean=-3.41416e-12

Initial energy perturbation visualization

Plot the initial moist static energy perturbation to ensure the bubble looks as expected.

fig = Figure()
ax = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Initial energy perturbation ρe′ (J / kg)")
hm = heatmap!(ax, ρe′)
Colorbar(fig[1, 2], hm, label = "ρe′ (J/kg)")
fig

Simulation rising

simulation = Simulation(model; Δt=2, stop_time=25minutes)
conjure_time_step_wizard!(simulation, cfl=0.7)
Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)

function progress(sim)
    ρe = static_energy_density(sim.model)
    u, v, w = sim.model.velocities

    msg = @sprintf("Iter: %d, t: %s, Δt: %s, extrema(ρe): (%.2f, %.2f) J/kg, max|u|: %.2f m/s, max|w|: %.2f m/s",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt),
                   minimum(ρe), maximum(ρe),
                   maximum(abs, u), maximum(abs, w))

    @info msg
    return nothing
end

add_callback!(simulation, progress, TimeInterval(1minute))

u, v, w = model.velocities
T = model.temperature

outputs = merge(model.velocities, model.tracers, (; ρe′, ρe, T))

filename = "thermal_bubble.jld2"
writer = JLD2Writer(model, outputs; filename,
                    schedule = TimeInterval(10seconds),
                    overwrite_existing = true)

simulation.output_writers[:jld2] = writer

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter: 0, t: 0 seconds, Δt: 2.200 seconds, extrema(ρe): (127306.19, 353642.72) J/kg, max|u|: 0.00 m/s, max|w|: 0.00 m/s
[ Info:     ... simulation initialization complete (35.044 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.500 seconds).
[ Info: Iter: 28, t: 1 minute, Δt: 2.662 seconds, extrema(ρe): (127306.19, 353642.74) J/kg, max|u|: 3.61 m/s, max|w|: 8.96 m/s
[ Info: Iter: 52, t: 2 minutes, Δt: 3.210 seconds, extrema(ρe): (127306.19, 353642.79) J/kg, max|u|: 7.37 m/s, max|w|: 17.13 m/s
[ Info: Iter: 76, t: 3 minutes, Δt: 2.554 seconds, extrema(ρe): (127306.19, 353642.85) J/kg, max|u|: 11.34 m/s, max|w|: 22.15 m/s
[ Info: Iter: 105, t: 4 minutes, Δt: 2.185 seconds, extrema(ρe): (127306.19, 353642.90) J/kg, max|u|: 15.11 m/s, max|w|: 25.34 m/s
[ Info: Iter: 138, t: 5 minutes, Δt: 1.875 seconds, extrema(ρe): (127306.19, 353642.95) J/kg, max|u|: 18.55 m/s, max|w|: 29.97 m/s
[ Info: Iter: 176, t: 6 minutes, Δt: 1.624 seconds, extrema(ρe): (127306.19, 353642.98) J/kg, max|u|: 22.04 m/s, max|w|: 33.57 m/s
[ Info: Iter: 218, t: 7 minutes, Δt: 1.529 seconds, extrema(ρe): (127306.18, 353643.01) J/kg, max|u|: 25.04 m/s, max|w|: 33.91 m/s
[ Info: Iter: 260, t: 8 minutes, Δt: 1.405 seconds, extrema(ρe): (127306.16, 353643.04) J/kg, max|u|: 29.81 m/s, max|w|: 35.00 m/s
[ Info: Iter: 308, t: 9 minutes, Δt: 1.392 seconds, extrema(ρe): (127306.11, 353643.08) J/kg, max|u|: 38.04 m/s, max|w|: 32.43 m/s
[ Info: Iter: 356, t: 10 minutes, Δt: 1.375 seconds, extrema(ρe): (127306.18, 353643.12) J/kg, max|u|: 46.27 m/s, max|w|: 28.97 m/s
[ Info: Iter: 401, t: 11 minutes, Δt: 1.541 seconds, extrema(ρe): (127301.81, 353643.18) J/kg, max|u|: 45.49 m/s, max|w|: 25.90 m/s
[ Info: Iter: 445, t: 12 minutes, Δt: 1.342 seconds, extrema(ρe): (127270.29, 353643.26) J/kg, max|u|: 36.67 m/s, max|w|: 30.07 m/s
[ Info: Iter: 493, t: 13 minutes, Δt: 1.347 seconds, extrema(ρe): (127570.64, 353643.37) J/kg, max|u|: 33.46 m/s, max|w|: 38.09 m/s
[ Info: Iter: 546, t: 14 minutes, Δt: 1.093 seconds, extrema(ρe): (127615.88, 353643.50) J/kg, max|u|: 31.40 m/s, max|w|: 45.52 m/s
[ Info: Iter: 608, t: 15 minutes, Δt: 955.613 ms, extrema(ρe): (127183.61, 353643.64) J/kg, max|u|: 33.48 m/s, max|w|: 55.37 m/s
[ Info: Iter: 670, t: 16 minutes, Δt: 1.079 seconds, extrema(ρe): (128298.35, 353643.81) J/kg, max|u|: 29.62 m/s, max|w|: 47.47 m/s
[ Info: Iter: 733, t: 17 minutes, Δt: 919.310 ms, extrema(ρe): (128208.66, 353644.01) J/kg, max|u|: 27.13 m/s, max|w|: 57.16 m/s
[ Info: Iter: 799, t: 18 minutes, Δt: 988.444 ms, extrema(ρe): (128348.94, 353644.29) J/kg, max|u|: 26.78 m/s, max|w|: 53.38 m/s
[ Info: Iter: 865, t: 19 minutes, Δt: 990.684 ms, extrema(ρe): (128738.46, 353644.87) J/kg, max|u|: 26.92 m/s, max|w|: 52.99 m/s
[ Info: Iter: 928, t: 20 minutes, Δt: 1.038 seconds, extrema(ρe): (128778.59, 353679.96) J/kg, max|u|: 23.43 m/s, max|w|: 49.72 m/s
[ Info: Iter: 982, t: 21 minutes, Δt: 1.524 seconds, extrema(ρe): (128780.68, 355400.96) J/kg, max|u|: 52.77 m/s, max|w|: 26.37 m/s
[ Info: Iter: 1019, t: 22 minutes, Δt: 2.015 seconds, extrema(ρe): (128736.55, 355179.36) J/kg, max|u|: 53.91 m/s, max|w|: 22.74 m/s
[ Info: Iter: 1055, t: 23 minutes, Δt: 1.873 seconds, extrema(ρe): (128308.10, 355044.98) J/kg, max|u|: 53.47 m/s, max|w|: 24.02 m/s
[ Info: Iter: 1095, t: 24 minutes, Δt: 1.524 seconds, extrema(ρe): (127908.83, 354851.93) J/kg, max|u|: 52.88 m/s, max|w|: 26.05 m/s
[ Info: Simulation is stopping after running for 2.320 minutes.
[ Info: Simulation time 25 minutes equals or exceeds stop time 25 minutes.
[ Info: Iter: 1137, t: 25 minutes, Δt: 1.582 seconds, extrema(ρe): (127764.02, 354742.11) J/kg, max|u|: 52.31 m/s, max|w|: 30.10 m/s

Visualization

Visualize the moist static energy perturbation and the vertical velocity through time and create an animation.

@info "Creating visualization..."

ρe′t = FieldTimeSeries(filename, "ρe′")
wt = FieldTimeSeries(filename, "w")

times = ρe′t.times
Nt = length(ρe′t)

fig = Figure(size = (800, 800), fontsize = 12)
axρ = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Energy perturbation ρe′ (J / kg)")
axw = Axis(fig[2, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Vertical velocity w (m / s)")

n = Observable(Nt)

ρe′n = @lift ρe′t[$n]
wn = @lift wt[$n]

title = @lift "Thermal bubble evolution — t = $(prettytime(times[$n]))"
fig[0, :] = Label(fig, title, fontsize = 16, tellwidth = false)

ρe′_range = (minimum(ρe′t), maximum(ρe′t))
w_range = maximum(abs, wt)

hmρ = heatmap!(axρ, ρe′n, colorrange = ρe′_range, colormap = :balance)
hmw = heatmap!(axw, wn, colorrange = (-w_range, w_range), colormap = :balance)

Colorbar(fig[1, 2], hmρ, label = "ρe′ (J/kg)", vertical = true)
Colorbar(fig[2, 2], hmw, label = "w (m/s)", vertical = true)

CairoMakie.record(fig, "thermal_bubble.mp4", 1:Nt, framerate = 12) do nn
    n[] = nn
end

[ Info: Creating visualization...

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.6
Commit 15346901f00 (2026-04-09 19:20 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_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_VERSION = 1.12.6
  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.4.7 `.`
⌃ [052768ef] CUDA v5.11.0
  [13f3f980] CairoMakie v0.15.9
⌅ [6a9e3e04] CloudMicrophysics v0.32.0
  [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.4.1
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.107.3
  [a01a1ee8] RRTMGP v0.21.7
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.1
  [9a3f8284] Random v1.11.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.