Cloudy thermal bubble

This example sets up, runs, and visualizes simulations of "thermal bubbles" (just circular regions of warm air) rising through a neutral background. We run a dry simulation and two "cloudy" simulations, both with and without precipitation. In the cloudy cases, we simulate a pocket of warm air rising in a saturated, condensate-laden environment.

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

Dry thermal bubble

We first set up a dry thermal bubble simulation without moisture processes. This serves as a baseline for comparison with the moist case.

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

thermodynamic_constants = ThermodynamicConstants()
reference_state = ReferenceState(grid, thermodynamic_constants, surface_pressure=1e5, potential_temperature=300)
dynamics = AnelasticDynamics(reference_state)
advection = WENO(order=9)
model = AtmosphereModel(grid; dynamics, thermodynamic_constants, advection)

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

Potential temperature perturbation

We add a localized potential temperature perturbation for the dry bubble. In the dry case, this perturbation directly affects buoyancy without any moisture-related effects.

r₀ = 2e3
z₀ = 2e3
Δθ = 2 # K
θ₀ = model.dynamics.reference_state.potential_temperature
g = model.thermodynamic_constants.gravitational_acceleration

function θᵢ(x, z)
    r = sqrt((x / r₀)^2 + ((z - z₀) / r₀)^2)
    return θ₀ + Δθ * cos(π * min(1, r) / 2)^2
end

set!(model, θ=θᵢ)

Initial dry bubble visualization

Plot the initial potential temperature to visualize the dry thermal bubble.

θ = liquid_ice_potential_temperature(model)
E = total_energy(model)
∫E = Integral(E) |> Field

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

Simulation rising

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

function progress(sim)
    u, v, w = sim.model.velocities
    msg = @sprintf("Iter: % 4d, t: % 14s, Δt: % 14s, ⟨E⟩: %.8e J, extrema(θ): (%.2f, %.2f) K, max|w|: %.2f m/s",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt), mean(E), extrema(θ)..., maximum(abs, w))
    @info msg
    return nothing
end

add_callback!(simulation, progress, TimeInterval(100))

u, v, w = model.velocities
outputs = (; θ, w)

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

simulation.output_writers[:jld2] = writer

run!(simulation)

fig = Figure()
axθ = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)")
axw = Axis(fig[2, 1], aspect=2, xlabel="x (m)", ylabel="z (m)")

hmθ = heatmap!(axθ, θ)
hmw = heatmap!(axw, w)

Colorbar(fig[1, 2], hmθ, label = "θ (K) at t = $(prettytime(simulation.model.clock.time))")
Colorbar(fig[2, 2], hmw, label = "w (m/s) at t = $(prettytime(simulation.model.clock.time))")

fig

Just running to t=1000 is pretty boring, Let's run the simulation for a longer time, just for fun!

simulation.stop_time = 30minutes run!(simulation)

Visualization

Visualize the potential temperature and the vertical velocity through time and create an animation.

θt = FieldTimeSeries(filename, "θ")
wt = FieldTimeSeries(filename, "w")

times = θt.times
fig = Figure(size = (800, 800), fontsize = 12)
axθ = Axis(fig[1, 1], aspect=2, xlabel="x (m)", ylabel="z (m)")
axw = Axis(fig[2, 1], aspect=2, xlabel="x (m)", ylabel="z (m)")

n = Observable(length(θt))
θn = @lift θt[$n]
wn = @lift wt[$n]

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

θ_range = (minimum(θt), maximum(θt))
w_range = maximum(abs, wt)

hmθ = heatmap!(axθ, θn, colorrange = θ_range, colormap = :thermal)
hmw = heatmap!(axw, wn, colorrange = (-w_range, w_range), colormap = :balance)

Colorbar(fig[1, 2], hmθ, label = "θ (K)", vertical = true)
Colorbar(fig[2, 2], hmw, label = "w (m/s)", vertical = true)

CairoMakie.record(fig, "dry_thermal_bubble.mp4", 1:length(θt), framerate = 12) do nn
    n[] = nn
end

Moist thermal bubble with warm-phase saturation adjustment

Now we set up a moist thermal bubble simulation with warm-phase saturation adjustment, following the methodology described by Bryan and Fritsch (2002). This simulation includes moisture processes, where excess water vapor condenses to liquid water, releasing latent heat that enhances the buoyancy of the rising bubble.

For pedagogical purposes, we build a new model with warm-phase saturation adjustment microphysics. (We could have also used this model for the dry simulation):

microphysics = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())
moist_model = AtmosphereModel(grid; dynamics, thermodynamic_constants, advection, microphysics)

AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×128 RectilinearGrid{Float64, Bounded, Flat, Bounded} on CPU with 5×0×5 halo
├── dynamics: AnelasticDynamics(p₀=100000.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: WENO{5, Float64, Float32}(order=9)
│   ├── ρθ: WENO{5, Float64, Float32}(order=9)
│   └── ρqᵉ: WENO{5, Float64, Float32}(order=9)
├── forcing: @NamedTuple{ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵉ::Returns{Float64}, ρe::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: SaturationAdjustment

Moist thermal bubble initial conditions

For the moist bubble, we initialize both temperature and moisture perturbations. The bubble is warm and moist, leading to condensation and latent heat release as it rises and cools. First, we set the potential temperature to match the dry case, then we use the diagnostic saturation specific humidity field to set the moisture.

Set potential temperature to match the dry bubble initially

set!(moist_model, θ=θᵢ, qᵗ=0.025)

Compute saturation specific humidity using the diagnostic field, and adjust the buoyancy to match the dry bubble Note, this isn't quite right and needs to be fixed.

using Breeze.Thermodynamics: dry_air_gas_constant, vapor_gas_constant

qᵛ⁺ = SaturationSpecificHumidityField(moist_model, :equilibrium)
θᵈ = liquid_ice_potential_temperature(moist_model) # note, current state is dry
Rᵈ = dry_air_gas_constant(thermodynamic_constants)
Rᵛ = vapor_gas_constant(thermodynamic_constants)
Rᵐ = Rᵈ * (1 - qᵛ⁺) + Rᵛ * qᵛ⁺
θᵐ = θᵈ * Rᵈ / Rᵐ

set!(moist_model, θ=θᵐ)

Simulation

moist_simulation = Simulation(moist_model; Δt=2, stop_time=30minutes)
conjure_time_step_wizard!(moist_simulation, cfl=0.7)
Oceananigans.Diagnostics.erroring_NaNChecker!(moist_simulation)

E = total_energy(moist_model)
θ = liquid_ice_potential_temperature(moist_model)

function progress_moist(sim)
    ρqᵉ = sim.model.moisture_density
    u, v, w = sim.model.velocities

    msg = @sprintf("Iter: % 4d, t: % 14s, Δt: % 14s, ⟨E⟩: %.8e J, extrema(θ): (%.2f, %.2f) K \n",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt), mean(E), extrema(θ)...)

    msg *= @sprintf("   extrema(ρqᵉ): (%.2e, %.2e), max(qˡ): %.2e, max|w|: %.2f m/s, mean(ρqᵉ): %.2e",
                    extrema(ρqᵉ)..., maximum(qˡ), maximum(abs, w), mean(ρqᵉ))

    @info msg
    return nothing
end

add_callback!(moist_simulation, progress_moist, TimeInterval(3minutes))

θ = liquid_ice_potential_temperature(moist_model)
u, v, w = moist_model.velocities
qᵛ = specific_humidity(moist_model)
qˡ = moist_model.microphysical_fields.qˡ
qˡ′ = qˡ - Field(Average(qˡ, dims=1))
moist_outputs = (; θ, w, qˡ′)

moist_filename = "cloudy_thermal_bubble.jld2"
moist_writer = JLD2Writer(moist_model, moist_outputs; filename=moist_filename,
                          schedule = TimeInterval(10seconds),
                          overwrite_existing = true)

moist_simulation.output_writers[:jld2] = moist_writer

run!(moist_simulation)

[ Info: Initializing simulation...
┌ Info: Iter:    0, t:      0 seconds, Δt:  2.200 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.64, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 0.00 m/s, mean(ρqᵉ): 1.91e-02
[ Info:     ... simulation initialization complete (16.701 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.294 seconds).
┌ Info: Iter:   70, t:      3 minutes, Δt:  4.287 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.64, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 1.80 m/s, mean(ρqᵉ): 1.91e-02
┌ Info: Iter:  113, t:      6 minutes, Δt:  6.277 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.63, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 3.37 m/s, mean(ρqᵉ): 1.91e-02
┌ Info: Iter:  149, t:      9 minutes, Δt:  8.354 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.63, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 4.34 m/s, mean(ρqᵉ): 1.91e-02
┌ Info: Iter:  173, t:     12 minutes, Δt: 11.120 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.63, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 4.52 m/s, mean(ρqᵉ): 1.91e-02
┌ Info: Iter:  191, t:     15 minutes, Δt: 13.172 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.62, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 4.03 m/s, mean(ρqᵉ): 1.91e-02
┌ Info: Iter:  209, t:     18 minutes, Δt: 13.782 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.62, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 3.68 m/s, mean(ρqᵉ): 1.91e-02
┌ Info: Iter:  227, t:     21 minutes, Δt: 13.735 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.62, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 3.62 m/s, mean(ρqᵉ): 1.91e-02
┌ Info: Iter:  245, t:     24 minutes, Δt: 13.246 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.62, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 3.53 m/s, mean(ρqᵉ): 1.91e-02
┌ Info: Iter:  263, t:     27 minutes, Δt: 13.282 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.63, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 3.49 m/s, mean(ρqᵉ): 1.91e-02
[ Info: Simulation is stopping after running for 1.236 minutes.
[ Info: Simulation time 30 minutes equals or exceeds stop time 30 minutes.
┌ Info: Iter:  281, t:     30 minutes, Δt: 12.108 seconds, ⟨E⟩: 3.14202269e+05 J, extrema(θ): (295.63, 299.21) K 
└    extrema(ρqᵉ): (1.09e-02, 2.89e-02), max(qˡ): 2.08e-02, max|w|: 3.89 m/s, mean(ρqᵉ): 1.91e-02

Visualization of moist thermal bubble

θt = FieldTimeSeries(moist_filename, "θ")
wt = FieldTimeSeries(moist_filename, "w")
qˡ′t = FieldTimeSeries(moist_filename, "qˡ′")

times = θt.times
fig = Figure(size = (1800, 800), fontsize = 12)
axθ = Axis(fig[1, 2], aspect=2, xlabel="x (m)", ylabel="z (m)")
axw = Axis(fig[1, 3], aspect=2, xlabel="x (m)", ylabel="z (m)")
axl = Axis(fig[2, 2:3], aspect=2, xlabel="x (m)", ylabel="z (m)")

θ_range = (minimum(θt), maximum(θt))
w_range = maximum(abs, wt)
qˡ′_range = (minimum(qˡ′t), maximum(qˡ′t))

n = Observable(length(θt))
θn = @lift θt[$n]
wn = @lift wt[$n]
qˡ′n = @lift qˡ′t[$n]

hmθ = heatmap!(axθ, θn, colorrange = θ_range, colormap = :thermal)
hmw = heatmap!(axw, wn, colorrange = (-w_range, w_range), colormap = :balance)
hml = heatmap!(axl, qˡ′n, colorrange = qˡ′_range, colormap = :balance)

Colorbar(fig[1, 1], hmθ, label = "θ (K)", vertical = true)
Colorbar(fig[1, 4], hmw, label = "w (m/s)", vertical = true)
Colorbar(fig[2, 4], hml, label = "qˡ (kg/kg)", vertical = true)

CairoMakie.record(fig, "cloudy_thermal_bubble.mp4", 1:length(θt), framerate = 24) do nn
    n[] = nn
end

Moist thermal bubble with precipitating one-moment microphysics

Next, we extend the moist thermal bubble example to a precipitating case using OneMomentCloudMicrophysics, which adds prognostic rain via autoconversion (cloud droplets coalescing to form rain) and accretion (rain collecting cloud droplets). This follows the CM1 benchmark configuration (iinit=4, isnd=4).

Note: The one-moment microphysics requires the CloudMicrophysics.jl package to be loaded, which activates the BreezeCloudMicrophysicsExt extension.

using CloudMicrophysics
BreezeCloudMicrophysicsExt = Base.get_extension(Breeze, :BreezeCloudMicrophysicsExt)
using .BreezeCloudMicrophysicsExt: OneMomentCloudMicrophysics

Build a new model with one-moment microphysics. We use saturation adjustment for cloud formation, but now rain is a prognostic variable that evolves via microphysical processes. We also use the same initial conditions as the moist case, but with slightly lower total water (qᵗ = 0.020) following the CM1 benchmark.

precip_cloud_formation = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())
precip_microphysics = OneMomentCloudMicrophysics(; cloud_formation=precip_cloud_formation)
precip_model = AtmosphereModel(grid; dynamics, thermodynamic_constants, advection,
                               microphysics=precip_microphysics)

qᵗ_precip = 0.020  # CM1 qt_mb value for saturated neutrally-stable sounding
set!(precip_model, θ=θᵢ, qᵗ=qᵗ_precip)

Simulation

We run the simulation for 60 minutes to allow precipitation to develop. The one-moment scheme requires time for cloud liquid to accumulate and autoconversion to produce rain.

precip_simulation = Simulation(precip_model; Δt=2, stop_time=60minutes)
conjure_time_step_wizard!(precip_simulation, cfl=0.7)
Oceananigans.Diagnostics.erroring_NaNChecker!(precip_simulation)

θ_precip = liquid_ice_potential_temperature(precip_model)
u_p, v_p, w_precip = precip_model.velocities
qˡ_precip = precip_model.microphysical_fields.qˡ    # Total liquid (cloud + rain)
qᶜˡ_precip = precip_model.microphysical_fields.qᶜˡ  # Cloud liquid only
qʳ_precip = precip_model.microphysical_fields.qʳ    # Rain mixing ratio

function progress_precip(sim)
    qᶜˡmax = maximum(qᶜˡ_precip)
    qʳmax = maximum(qʳ_precip)
    wmax = maximum(abs, w_precip)

    msg = @sprintf("Iter: %4d, t: %14s, Δt: %14s, max|w|: %.2f m/s",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt), wmax)
    msg *= @sprintf(", max(qᶜˡ): %.2e, max(qʳ): %.2e", qᶜˡmax, qʳmax)

    @info msg
    return nothing
end

add_callback!(precip_simulation, progress_precip, TimeInterval(5minutes))

precip_outputs = (; θ=θ_precip, w=w_precip, qᶜˡ=qᶜˡ_precip, qʳ=qʳ_precip)

precip_filename = "precipitating_thermal_bubble.jld2"
precip_writer = JLD2Writer(precip_model, precip_outputs; filename=precip_filename,
                           schedule = TimeInterval(30seconds),
                           overwrite_existing = true)

precip_simulation.output_writers[:jld2] = precip_writer

run!(precip_simulation)

[ Info: Initializing simulation...
[ Info: Iter:    0, t:      0 seconds, Δt:  2.200 seconds, max|w|: 0.00 m/s, max(qᶜˡ): 1.77e-02, max(qʳ): 0.00e+00
[ Info:     ... simulation initialization complete (14.997 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.458 seconds).
[ Info: Iter:   94, t:      5 minutes, Δt:  5.187 seconds, max|w|: 3.36 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.53e-02
[ Info: Iter:  145, t:     10 minutes, Δt:  7.330 seconds, max|w|: 7.90 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.53e-02
[ Info: Iter:  208, t:     15 minutes, Δt:  3.905 seconds, max|w|: 14.72 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.54e-02
[ Info: Iter:  312, t:     20 minutes, Δt:  2.684 seconds, max|w|: 20.21 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.55e-02
[ Info: Iter:  432, t:     25 minutes, Δt:  2.607 seconds, max|w|: 20.79 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.58e-02
[ Info: Iter:  549, t:     30 minutes, Δt:  2.803 seconds, max|w|: 19.51 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.63e-02
[ Info: Iter:  659, t:     35 minutes, Δt:  3.003 seconds, max|w|: 17.74 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.64e-02
[ Info: Iter:  759, t:     40 minutes, Δt:  3.117 seconds, max|w|: 15.49 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.63e-02
[ Info: Iter:  872, t:     45 minutes, Δt:  2.111 seconds, max|w|: 25.75 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.63e-02
[ Info: Iter: 1029, t:     50 minutes, Δt:  1.720 seconds, max|w|: 32.19 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.63e-02
[ Info: Iter: 1235, t:     55 minutes, Δt:  1.501 seconds, max|w|: 35.80 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.62e-02
[ Info: Simulation is stopping after running for 5.811 minutes.
[ Info: Simulation time 1 hour equals or exceeds stop time 1 hour.
[ Info: Iter: 1400, t:         1 hour, Δt:  2.126 seconds, max|w|: 17.10 m/s, max(qᶜˡ): 1.76e-02, max(qʳ): 1.63e-02

Visualization of a precipitating thermal bubble

θts = FieldTimeSeries(precip_filename, "θ")
wts = FieldTimeSeries(precip_filename, "w")
qᶜˡts = FieldTimeSeries(precip_filename, "qᶜˡ")
qʳts = FieldTimeSeries(precip_filename, "qʳ")

times_precip = θts.times
Nt = length(times_precip)

θ_range_p = (minimum(θts), maximum(θts))
w_range_p = maximum(abs, wts)
qᶜˡ_range = (0, max(1e-6, maximum(qᶜˡts)))
qʳ_range = (0, max(1e-6, maximum(qʳts)))

fig = Figure(size=(1400, 700), fontsize=11)
axθ = Axis(fig[1, 2], aspect=2, xlabel="x (m)", ylabel="z (m)", title="θ (K)")
axw = Axis(fig[1, 3], aspect=2, xlabel="x (m)", ylabel="z (m)", title="w (m/s)")
axqᶜˡ = Axis(fig[2, 2], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Cloud liquid qᶜˡ (kg/kg)")
axqʳ = Axis(fig[2, 3], aspect=2, xlabel="x (m)", ylabel="z (m)", title="Rain qʳ (kg/kg)")

n = Observable(1)
θn = @lift θts[$n]
wn = @lift wts[$n]
qᶜˡn = @lift qᶜˡts[$n]
qʳn = @lift qʳts[$n]

hmθ = heatmap!(axθ, θn, colorrange=θ_range_p, colormap=:thermal)
hmw = heatmap!(axw, wn, colorrange=(-w_range_p, w_range_p), colormap=:balance)
hmqᶜˡ = heatmap!(axqᶜˡ, qᶜˡn, colorrange=qᶜˡ_range, colormap=:dense)
hmqʳ = heatmap!(axqʳ, qʳn, colorrange=qʳ_range, colormap=:amp)

Colorbar(fig[1, 1], hmθ, label="θ (K)", vertical=true, width=15)
Colorbar(fig[1, 4], hmw, label="w (m/s)", vertical=true, width=15)
Colorbar(fig[2, 1], hmqᶜˡ, label="qᶜˡ (kg/kg)", vertical=true, width=15)
Colorbar(fig[2, 4], hmqʳ, label="qʳ (kg/kg)", vertical=true, width=15)

colgap!(fig.layout, 10)
rowgap!(fig.layout, 10)

CairoMakie.record(fig, "precipitating_thermal_bubble.mp4", 1:Nt, framerate=12) 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.5
Commit 5fe89b8ddc1 (2026-02-09 16:05 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.5
  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.6 `.`
  [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.106.4
  [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.