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.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, Float32}(order=9)
│   ├── ρe: WENO{5, Float64, Float32}(order=9)
│   └── ρqᵗ: WENO{5, Float64, Float32}(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=7277.39, min=-849.551, mean=1.71951e-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)

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): (126538.65, 353636.56) J/kg, max|u|: 0.00 m/s, max|w|: 0.00 m/s
[ Info:     ... simulation initialization complete (36.331 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.601 seconds).
[ Info: Iter: 28, t: 1 minute, Δt: 2.662 seconds, extrema(ρe): (126538.65, 353636.58) J/kg, max|u|: 3.62 m/s, max|w|: 8.96 m/s
[ Info: Iter: 52, t: 2 minutes, Δt: 3.204 seconds, extrema(ρe): (126538.64, 353636.62) J/kg, max|u|: 7.39 m/s, max|w|: 17.17 m/s
[ Info: Iter: 76, t: 3 minutes, Δt: 2.543 seconds, extrema(ρe): (126538.63, 353636.67) J/kg, max|u|: 11.38 m/s, max|w|: 22.26 m/s
[ Info: Iter: 105, t: 4 minutes, Δt: 2.166 seconds, extrema(ρe): (126538.62, 353636.70) J/kg, max|u|: 15.20 m/s, max|w|: 25.58 m/s
[ Info: Iter: 138, t: 5 minutes, Δt: 1.847 seconds, extrema(ρe): (126538.59, 353636.73) J/kg, max|u|: 18.73 m/s, max|w|: 30.47 m/s
[ Info: Iter: 177, t: 6 minutes, Δt: 1.592 seconds, extrema(ρe): (126538.54, 353636.75) J/kg, max|u|: 22.34 m/s, max|w|: 34.34 m/s
[ Info: Iter: 219, t: 7 minutes, Δt: 1.508 seconds, extrema(ρe): (126538.46, 353636.76) J/kg, max|u|: 25.16 m/s, max|w|: 34.49 m/s
[ Info: Iter: 263, t: 8 minutes, Δt: 1.392 seconds, extrema(ρe): (126538.36, 353636.78) J/kg, max|u|: 31.00 m/s, max|w|: 35.60 m/s
[ Info: Iter: 311, t: 9 minutes, Δt: 1.401 seconds, extrema(ρe): (126538.24, 353636.79) J/kg, max|u|: 39.59 m/s, max|w|: 32.88 m/s
[ Info: Iter: 359, t: 10 minutes, Δt: 1.384 seconds, extrema(ρe): (126538.13, 353636.80) J/kg, max|u|: 45.72 m/s, max|w|: 28.56 m/s
[ Info: Iter: 404, t: 11 minutes, Δt: 1.507 seconds, extrema(ρe): (126538.20, 353636.80) J/kg, max|u|: 45.83 m/s, max|w|: 25.66 m/s
[ Info: Iter: 448, t: 12 minutes, Δt: 1.293 seconds, extrema(ρe): (126561.27, 353636.81) J/kg, max|u|: 36.59 m/s, max|w|: 32.96 m/s
[ Info: Iter: 496, t: 13 minutes, Δt: 1.264 seconds, extrema(ρe): (126539.23, 353636.82) J/kg, max|u|: 34.08 m/s, max|w|: 43.95 m/s
[ Info: Iter: 553, t: 14 minutes, Δt: 1.006 seconds, extrema(ρe): (126894.44, 353636.82) J/kg, max|u|: 32.50 m/s, max|w|: 52.17 m/s
[ Info: Iter: 622, t: 15 minutes, Δt: 935.791 ms, extrema(ρe): (127122.06, 353636.83) J/kg, max|u|: 31.42 m/s, max|w|: 56.62 m/s
[ Info: Iter: 688, t: 16 minutes, Δt: 960.347 ms, extrema(ρe): (127086.20, 353636.83) J/kg, max|u|: 28.23 m/s, max|w|: 56.68 m/s
[ Info: Iter: 758, t: 17 minutes, Δt: 919.735 ms, extrema(ρe): (127097.55, 353636.83) J/kg, max|u|: 26.87 m/s, max|w|: 57.04 m/s
[ Info: Iter: 824, t: 18 minutes, Δt: 927.693 ms, extrema(ρe): (127091.78, 353636.83) J/kg, max|u|: 25.77 m/s, max|w|: 56.69 m/s
[ Info: Iter: 890, t: 19 minutes, Δt: 964.816 ms, extrema(ρe): (127136.82, 353636.84) J/kg, max|u|: 27.71 m/s, max|w|: 54.23 m/s
[ Info: Iter: 949, t: 20 minutes, Δt: 1.285 seconds, extrema(ρe): (127108.77, 355264.98) J/kg, max|u|: 56.67 m/s, max|w|: 30.29 m/s
[ Info: Iter: 994, t: 21 minutes, Δt: 1.736 seconds, extrema(ρe): (127081.10, 355058.52) J/kg, max|u|: 58.09 m/s, max|w|: 22.55 m/s
[ Info: Iter: 1030, t: 22 minutes, Δt: 1.906 seconds, extrema(ρe): (127050.27, 354907.25) J/kg, max|u|: 57.37 m/s, max|w|: 19.79 m/s
[ Info: Iter: 1066, t: 23 minutes, Δt: 1.804 seconds, extrema(ρe): (127003.85, 354789.91) J/kg, max|u|: 56.75 m/s, max|w|: 19.41 m/s
[ Info: Iter: 1106, t: 24 minutes, Δt: 1.588 seconds, extrema(ρe): (126961.75, 354700.92) J/kg, max|u|: 55.87 m/s, max|w|: 25.22 m/s
[ Info: Simulation is stopping after running for 2.461 minutes.
[ Info: Simulation time 25 minutes equals or exceeds stop time 25 minutes.
[ Info: Iter: 1148, t: 25 minutes, Δt: 1.486 seconds, extrema(ρe): (126928.42, 354636.20) J/kg, max|u|: 54.96 m/s, max|w|: 31.19 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.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: 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.4
  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.3.3 `.`
  [052768ef] CUDA v5.9.6
  [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
  [85f8d34a] NCDatasets v0.14.11
  [9e8cae18] Oceananigans v0.104.5
  [a01a1ee8] RRTMGP v0.21.7
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.1
  [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.