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 Breezeusing Oceananigans: Oceananigansusing Oceananigans.Unitsusing Statisticsusing Printfusing 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)
│   ├── ρs: 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}, ρs::Returns{Float64}, ρqᵛ::Returns{Float64}, ρE::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 # KN² = 1e-6θ₀ = model.dynamics.reference_state.potential_temperatureg = model.thermodynamic_constants.gravitational_accelerationfunction θᵢ(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 θ̄ + θ′endset!(model, θ = θᵢ)ρs = static_energy_density(model)ρS = Field(Average(ρs, dims=1))ρs′ = Field(ρs - ρS)

128×1×128 Field{Center, Center, Center} on 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 StaticStridesOffsetArray(::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 ρs′ (J / kg)")hm = heatmap!(ax, ρs′)Colorbar(fig[1, 2], hm, label = "ρs′ (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)    ρs = static_energy_density(sim.model)    u, v, w = sim.model.velocities    msg = @sprintf("Iter: %d, t: %s, Δt: %s, extrema(ρs): (%.2f, %.2f) J/kg, max|u|: %.2f m/s, max|w|: %.2f m/s",                   iteration(sim), prettytime(sim), prettytime(sim.Δt),                   minimum(ρs), maximum(ρs),                   maximum(abs, u), maximum(abs, w))    @info msg    return nothingendadd_callback!(simulation, progress, TimeInterval(1minute))u, v, w = model.velocitiesT = model.temperatureoutputs = merge(model.velocities, model.tracers, (; ρs′, ρs, T))filename = "thermal_bubble.jld2"writer = JLD2Writer(model, outputs; filename,                    schedule = TimeInterval(10seconds),                    overwrite_files = true)simulation.output_writers[:jld2] = writerrun!(simulation)

[ Info: Initializing simulation...
[ Info: Iter: 0, t: 0 seconds, Δt: 2.200 seconds, extrema(ρs): (127306.19, 353642.72) J/kg, max|u|: 0.00 m/s, max|w|: 0.00 m/s
[ Info:     ... simulation initialization complete (19.321 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (954.829 ms).
[ Info: Iter: 28, t: 1 minute, Δt: 2.662 seconds, extrema(ρs): (127306.19, 353642.74) J/kg, max|u|: 3.60 m/s, max|w|: 8.89 m/s
[ Info: Iter: 52, t: 2 minutes, Δt: 3.297 seconds, extrema(ρs): (127306.18, 353642.79) J/kg, max|u|: 7.25 m/s, max|w|: 16.65 m/s
[ Info: Iter: 76, t: 3 minutes, Δt: 2.681 seconds, extrema(ρs): (127306.17, 353642.85) J/kg, max|u|: 10.87 m/s, max|w|: 20.96 m/s
[ Info: Iter: 103, t: 4 minutes, Δt: 2.326 seconds, extrema(ρs): (127306.14, 353642.90) J/kg, max|u|: 14.02 m/s, max|w|: 23.54 m/s
[ Info: Iter: 133, t: 5 minutes, Δt: 2.000 seconds, extrema(ρs): (127306.08, 353642.95) J/kg, max|u|: 16.38 m/s, max|w|: 27.48 m/s
[ Info: Iter: 168, t: 6 minutes, Δt: 1.808 seconds, extrema(ρs): (127305.95, 353642.98) J/kg, max|u|: 20.03 m/s, max|w|: 30.15 m/s
[ Info: Iter: 204, t: 7 minutes, Δt: 1.714 seconds, extrema(ρs): (127305.70, 353643.01) J/kg, max|u|: 22.19 m/s, max|w|: 30.08 m/s
[ Info: Iter: 243, t: 8 minutes, Δt: 1.632 seconds, extrema(ρs): (127305.21, 353643.04) J/kg, max|u|: 21.33 m/s, max|w|: 31.33 m/s
[ Info: Iter: 285, t: 9 minutes, Δt: 1.564 seconds, extrema(ρs): (127304.43, 353643.08) J/kg, max|u|: 25.19 m/s, max|w|: 29.49 m/s
[ Info: Iter: 327, t: 10 minutes, Δt: 1.658 seconds, extrema(ρs): (127303.29, 353643.12) J/kg, max|u|: 29.97 m/s, max|w|: 26.38 m/s
[ Info: Iter: 363, t: 11 minutes, Δt: 1.666 seconds, extrema(ρs): (127302.86, 353643.16) J/kg, max|u|: 35.54 m/s, max|w|: 23.78 m/s
[ Info: Iter: 400, t: 12 minutes, Δt: 1.745 seconds, extrema(ρs): (127294.96, 353643.22) J/kg, max|u|: 35.26 m/s, max|w|: 21.73 m/s
[ Info: Iter: 436, t: 13 minutes, Δt: 1.969 seconds, extrema(ρs): (127302.53, 353643.29) J/kg, max|u|: 35.78 m/s, max|w|: 20.73 m/s
[ Info: Iter: 472, t: 14 minutes, Δt: 1.629 seconds, extrema(ρs): (127277.62, 353643.40) J/kg, max|u|: 28.96 m/s, max|w|: 25.33 m/s
[ Info: Iter: 513, t: 15 minutes, Δt: 1.604 seconds, extrema(ρs): (127310.83, 353643.52) J/kg, max|u|: 26.87 m/s, max|w|: 30.40 m/s
[ Info: Iter: 556, t: 16 minutes, Δt: 1.368 seconds, extrema(ρs): (127501.41, 353643.66) J/kg, max|u|: 27.81 m/s, max|w|: 38.09 m/s
[ Info: Iter: 606, t: 17 minutes, Δt: 1.131 seconds, extrema(ρs): (127458.91, 353643.81) J/kg, max|u|: 27.13 m/s, max|w|: 48.44 m/s
[ Info: Iter: 666, t: 18 minutes, Δt: 1.090 seconds, extrema(ρs): (127268.68, 353643.98) J/kg, max|u|: 24.33 m/s, max|w|: 48.20 m/s
[ Info: Iter: 721, t: 19 minutes, Δt: 1.191 seconds, extrema(ρs): (127379.68, 353644.17) J/kg, max|u|: 23.43 m/s, max|w|: 44.33 m/s
[ Info: Iter: 775, t: 20 minutes, Δt: 1.235 seconds, extrema(ρs): (127640.48, 353644.41) J/kg, max|u|: 24.18 m/s, max|w|: 41.89 m/s
[ Info: Iter: 824, t: 21 minutes, Δt: 1.364 seconds, extrema(ρs): (127690.86, 353644.77) J/kg, max|u|: 24.93 m/s, max|w|: 37.91 m/s
[ Info: Iter: 868, t: 22 minutes, Δt: 1.521 seconds, extrema(ρs): (127933.56, 353645.82) J/kg, max|u|: 22.54 m/s, max|w|: 33.50 m/s
[ Info: Iter: 907, t: 23 minutes, Δt: 1.808 seconds, extrema(ρs): (127938.81, 356179.89) J/kg, max|u|: 20.14 m/s, max|w|: 28.17 m/s
[ Info: Iter: 943, t: 24 minutes, Δt: 1.944 seconds, extrema(ρs): (127951.75, 356280.07) J/kg, max|u|: 29.41 m/s, max|w|: 26.90 m/s
[ Info: Simulation is stopping after running for 40.301 seconds.
[ Info: Simulation time 25 minutes equals or exceeds stop time 25 minutes.
[ Info: Iter: 976, t: 25 minutes, Δt: 2.020 seconds, extrema(ρs): (127969.81, 356137.19) J/kg, max|u|: 28.94 m/s, max|w|: 26.17 m/s

Visualization

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

@info "Creating visualization..."ρs′t = FieldTimeSeries(filename, "ρs′")wt = FieldTimeSeries(filename, "w")times = ρs′t.timesNt = length(ρs′t)fig = Figure(size = (800, 800), fontsize = 12)axρ = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Energy perturbation ρs′ (J / kg)")axw = Axis(fig[2, 1], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Vertical velocity w (m / s)")n = Observable(Nt)ρs′n = @lift ρs′t[$n]wn = @lift wt[$n]title = @lift "Thermal bubble evolution — t = $(prettytime(times[$n]))"fig[0, :] = Label(fig, title, fontsize = 16, tellwidth = false)ρs′_range = (minimum(ρs′t), maximum(ρs′t))w_range = maximum(abs, wt)hmρ = heatmap!(axρ, ρs′n, colorrange = ρs′_range, colormap = :balance)hmw = heatmap!(axw, wn, colorrange = (-w_range, w_range), colormap = :balance)Colorbar(fig[1, 2], hmρ, label = "ρs′ (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, compression = 23) do nn    n[] = nnend

[ Info: Creating visualization...

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfoversioninfo()

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: 8 × AMD EPYC 7R13 Processor
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 5 GC (on 8 virtual cores)
Environment:
  JULIA_NUM_GC_THREADS = 4,1
  JULIA_GPG = 64B779A570972FFF7BFC2B54EAD471E1A1F2C10A
  JULIA_LOAD_PATH = :@breeze
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_VERSION = 1.13.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 PkgPkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
⌃ [86bc3604] AtmosphericProfilesLibrary v0.1.8
  [660aa2fb] Breeze v0.15.1 `.`
⌃ [052768ef] CUDA v6.4.1
  [13f3f980] CairoMakie v0.15.17
  [6a9e3e04] CloudMicrophysics v0.49.1
  [e30172f5] Documenter v1.20.0
  [daee34ce] DocumenterCitations v1.5.0
  [b6400b83] DocumenterCodeBlocks v1.5.1
  [1c52b33b] DocumenterLandingPage v0.2.2
  [7da242da] Enzyme v0.13.219
⌅ [46192b85] GPUArraysCore v0.2.1
⌅ [63c18a36] KernelAbstractions v0.9.47
  [98b081ad] Literate v2.21.1
  [85f8d34a] NCDatasets v0.14.15
  [cd8119b0] NumericalRadiation v0.2.1
  [9e8cae18] Oceananigans v0.118.1
  [a01a1ee8] RRTMGP v1.1.1
  [3c362404] Reactant v0.2.308
  [276daf66] SpecialFunctions v2.9.2
  [10745b16] Statistics v1.11.5
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.13.0
  [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.