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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 5×0×5 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=288.0)
├── formulation: StaticEnergyFormulation
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: 
│   ├── momentum: WENO{5, Float64, Float32}(order=9)
│   ├── ρe: WENO{5, Float64, Float32}(order=9)
│   └── ρqᵗ: WENO{5, Float64, Float32}(order=9)
├── 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 Oceananigans.Architectures.CPU
├── grid: 128×1×128 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.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=7304.8, min=-852.751, mean=-1.1191e-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): (127015.20, 354968.39) J/kg, max|u|: 0.00 m/s, max|w|: 0.00 m/s
[ Info:     ... simulation initialization complete (28.805 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.175 seconds).
[ Info: Iter: 28, t: 1 minute, Δt: 2.662 seconds, extrema(ρe): (127015.20, 354968.49) J/kg, max|u|: 3.60 m/s, max|w|: 8.94 m/s
[ Info: Iter: 52, t: 2 minutes, Δt: 3.208 seconds, extrema(ρe): (127015.19, 354968.74) J/kg, max|u|: 7.37 m/s, max|w|: 17.15 m/s
[ Info: Iter: 76, t: 3 minutes, Δt: 2.540 seconds, extrema(ρe): (127015.17, 354969.01) J/kg, max|u|: 11.37 m/s, max|w|: 22.30 m/s
[ Info: Iter: 105, t: 4 minutes, Δt: 2.155 seconds, extrema(ρe): (127015.14, 354969.23) J/kg, max|u|: 15.23 m/s, max|w|: 25.74 m/s
[ Info: Iter: 138, t: 5 minutes, Δt: 1.825 seconds, extrema(ρe): (127015.08, 354969.39) J/kg, max|u|: 18.81 m/s, max|w|: 30.88 m/s
[ Info: Iter: 178, t: 6 minutes, Δt: 1.562 seconds, extrema(ρe): (127014.98, 354969.51) J/kg, max|u|: 22.51 m/s, max|w|: 35.07 m/s
[ Info: Iter: 220, t: 7 minutes, Δt: 1.457 seconds, extrema(ρe): (127014.82, 354969.60) J/kg, max|u|: 25.62 m/s, max|w|: 35.03 m/s
[ Info: Iter: 264, t: 8 minutes, Δt: 1.381 seconds, extrema(ρe): (127014.61, 354969.68) J/kg, max|u|: 31.99 m/s, max|w|: 36.02 m/s
[ Info: Iter: 312, t: 9 minutes, Δt: 1.387 seconds, extrema(ρe): (127014.37, 354969.74) J/kg, max|u|: 40.72 m/s, max|w|: 33.08 m/s
[ Info: Iter: 360, t: 10 minutes, Δt: 1.407 seconds, extrema(ρe): (127014.16, 354969.80) J/kg, max|u|: 46.43 m/s, max|w|: 29.12 m/s
[ Info: Iter: 404, t: 11 minutes, Δt: 1.496 seconds, extrema(ρe): (127014.31, 354969.85) J/kg, max|u|: 46.43 m/s, max|w|: 26.20 m/s
[ Info: Iter: 449, t: 12 minutes, Δt: 1.279 seconds, extrema(ρe): (126908.23, 354969.90) J/kg, max|u|: 38.48 m/s, max|w|: 34.33 m/s
[ Info: Iter: 499, t: 13 minutes, Δt: 1.157 seconds, extrema(ρe): (127174.51, 354969.94) J/kg, max|u|: 34.66 m/s, max|w|: 46.21 m/s
[ Info: Iter: 560, t: 14 minutes, Δt: 927.148 ms, extrema(ρe): (127238.33, 354969.97) J/kg, max|u|: 33.29 m/s, max|w|: 57.53 m/s
[ Info: Iter: 633, t: 15 minutes, Δt: 895.781 ms, extrema(ρe): (127433.76, 354970.00) J/kg, max|u|: 31.97 m/s, max|w|: 59.23 m/s
[ Info: Iter: 705, t: 16 minutes, Δt: 840.672 ms, extrema(ρe): (127424.28, 354970.03) J/kg, max|u|: 28.31 m/s, max|w|: 62.66 m/s
[ Info: Iter: 777, t: 17 minutes, Δt: 869.780 ms, extrema(ρe): (127393.96, 354970.05) J/kg, max|u|: 28.24 m/s, max|w|: 60.80 m/s
[ Info: Iter: 849, t: 18 minutes, Δt: 874.392 ms, extrema(ρe): (127363.13, 354970.07) J/kg, max|u|: 25.45 m/s, max|w|: 59.99 m/s
[ Info: Iter: 919, t: 19 minutes, Δt: 997.584 ms, extrema(ρe): (127339.65, 356153.25) J/kg, max|u|: 51.94 m/s, max|w|: 46.39 m/s
[ Info: Iter: 969, t: 20 minutes, Δt: 1.466 seconds, extrema(ρe): (127324.35, 356238.32) J/kg, max|u|: 63.92 m/s, max|w|: 31.04 m/s
[ Info: Iter: 1010, t: 21 minutes, Δt: 1.695 seconds, extrema(ρe): (127315.46, 356159.89) J/kg, max|u|: 63.33 m/s, max|w|: 31.74 m/s
[ Info: Iter: 1046, t: 22 minutes, Δt: 1.730 seconds, extrema(ρe): (127310.53, 356059.86) J/kg, max|u|: 63.40 m/s, max|w|: 24.58 m/s
[ Info: Iter: 1082, t: 23 minutes, Δt: 1.746 seconds, extrema(ρe): (127307.74, 355962.42) J/kg, max|u|: 63.31 m/s, max|w|: 26.03 m/s
[ Info: Iter: 1118, t: 24 minutes, Δt: 1.680 seconds, extrema(ρe): (127306.26, 355905.56) J/kg, max|u|: 62.63 m/s, max|w|: 28.89 m/s
[ Info: Simulation is stopping after running for 2.308 minutes.
[ Info: Simulation time 25 minutes equals or exceeds stop time 25 minutes.
[ Info: Iter: 1155, t: 25 minutes, Δt: 1.648 seconds, extrema(ρe): (127305.69, 355852.34) J/kg, max|u|: 61.62 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.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.