Splitting supercell — anelastic vs compressible
This example simulates the development of a splitting supercell thunderstorm with both Breeze dynamical cores side-by-side, following the idealized test case described by Klemp et al. (2015) and the DCMIP2016 supercell intercomparison by Zarzycki et al. (2019). This benchmark evaluates the model's ability to capture deep moist convection with warm-rain microphysics and strong updrafts.
For microphysics we use the Kessler scheme, which includes prognostic cloud water and rain water with autoconversion, accretion, rain evaporation, and sedimentation processes. This is the same scheme used in the DCMIP2016 supercell intercomparison (Zarzycki et al., 2019).
We run two simulations from identical initial conditions: one with the anelastic solver and one with the fully compressible solver (split-explicit acoustic substepping). At the end we compare horizontal (xy at $z \approx 5 \, {\rm km}$) and vertical (xz at $y = L_y/2$) slices side-by-side, and plot maximum vertical-velocity time series for both runs on the same axes.
Physical setup
The simulation initializes a conditionally unstable atmosphere with a warm bubble perturbation that triggers deep convection. The environment includes:
- A realistic tropospheric potential temperature profile with a tropopause at 12 km
- Relative humidity that decreases with height, with the resulting water vapor mixing ratio capped at 0.014 kg/kg "to approximate a well-mixed boundary layer in the lowest kilometer" (Klemp et al. (2015)).
- Wind shear in the lower 5 km to promote storm rotation and supercell development
Potential temperature profile
The background potential temperature follows a piecewise profile (Equation 14 in Klemp et al. (2015)):
\[θ(z) = \begin{cases} θ_0 + (θ_{\rm tr} - θ_0) \left(\dfrac{z}{z_{\rm tr}}\right)^{5/4} & z \leq z_{\rm tr} \\ θ_{\rm tr} \exp\left[\dfrac{g}{c_p^d T_{\rm tr}} (z - z_{\rm tr})\right] & z > z_{\rm tr} \end{cases}\]
where $θ_0 = 300 \, {\rm K}$ is the surface potential temperature, $θ_{\rm tr} = 343 \, {\rm K}$ is the tropopause potential temperature, $z_{\rm tr} = 12 \, {\rm km}$ is the tropopause height, and $T_{\rm tr} = 213 \, {\rm K}$ is the tropopause temperature.
Warm bubble perturbation
A localized warm bubble triggers convection (Equations 17–18 in Klemp et al. (2015)):
\[θ'(x, y, z) = \begin{cases} Δθ \cos^2\left(π R / 2 \right) & R < 1 \\ 0 & R \geq 1 \end{cases}\]
where $R = \sqrt{(r/r_h)^2 + [(z-z_c)/r_z]^2}$ is the normalized radius, $r = \sqrt{(x-x_c)^2 + (y-y_c)^2}$ is the horizontal distance from the bubble center, $Δθ = 3 \, {\rm K}$ is the perturbation amplitude, $r_h = 10 \, {\rm km}$ is the horizontal radius, and $r_z = 1.5 \, {\rm km}$ is the vertical radius.
Wind shear profile
The zonal wind increases linearly with height up to the shear layer $z_s = 5 \, {\rm km}$, with a smooth transition zone, providing the environmental shear necessary for supercell development and mesocyclone formation (Equations 15-16 in Klemp et al. (2015)).
using Breezeusing Breeze: DCMIP2016KesslerMicrophysics, TetensFormulausing Breeze.Thermodynamics: hydrostatic_density, hydrostatic_temperature, pressure_balanced_densityusing Oceananigans: Oceananigansusing Oceananigans.Unitsusing Oceananigans.Grids: znodesusing CairoMakieusing CUDAusing PrintfDomain and grid
The domain is 168 km × 168 km × 20 km with 168 × 168 × 40 grid points, giving 1 km horizontal resolution and 500 m vertical resolution. The grid uses periodic lateral boundary conditions and bounded top/bottom boundaries.
Oceananigans.defaults.FloatType = Float32Nx, Ny, Nz = 168, 168, 40Lx, Ly, Lz = 168kilometers, 168kilometers, 20kilometersgrid = RectilinearGrid(GPU(), size = (Nx, Ny, Nz), x = (0, Lx), y = (0, Ly), z = (0, Lz), halo = (5, 5, 5), topology = (Periodic, Periodic, Bounded))168×168×40 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── Periodic x ∈ [0.0, 168000.0) regularly spaced with Δx=1000.0
├── Periodic y ∈ [0.0, 168000.0) regularly spaced with Δy=1000.0
└── Bounded z ∈ [0.0, 20000.0] regularly spaced with Δz=500.0Background profiles
Thermodynamic constants and the surface/standard pressures shared by both runs:
constants = ThermodynamicConstants(saturation_vapor_pressure = TetensFormula())p₀ = 100000pˢᵗ = 100000100000Stratification parameters define the troposphere–stratosphere transition:
θ₀ = 300 # K - surface potential temperatureθᵖ = 343 # K - tropopause potential temperaturezᵖ = 12000 # m - tropopause heightTᵖ = 213 # K - tropopause temperatureqᵛ_max = 0.014 # kg/kg - cap on water vapor mixing ratio from Klemp et al. (2015)Wind shear parameters:
zˢ = 5kilometers # m - shear layer heightuˢ = 30 # m/s - maximum shear wind speeduᶜ = 15 # m/s - storm motion (Galilean translation speed)Thermodynamic constants used inside the profile functions:
g = constants.gravitational_accelerationcᵖᵈ = constants.dry_air.heat_capacityBackground potential temperature profile (Equation 14 in Klemp et al. (2015)):
function θ_background(z) θᵗ = θ₀ + (θᵖ - θ₀) * (z / zᵖ)^(5/4) θˢ = θᵖ * exp(g / (cᵖᵈ * Tᵖ) * (z - zᵖ)) return (z ≤ zᵖ) * θᵗ + (z > zᵖ) * θˢendθ_background (generic function with 1 method)Relative humidity profile (Equations 11–12 by Klemp et al. (2015)) combined with the water vapor cap $qᵛ_{max}$. The local temperature and density are obtained by numerically integrating the hydrostatic balance with the actual $θ(z)$ profile:
function qᵛ_bg(z) ℋ = (1 - 3/4 * (z / zᵖ)^(5/4)) * (z ≤ zᵖ) + 1/4 * (z > zᵖ) T = hydrostatic_temperature(z, p₀, θ_background, pˢᵗ, constants) ρ = hydrostatic_density(z, p₀, θ_background, pˢᵗ, constants) qᵛ⁺ = saturation_specific_humidity(T, ρ, constants, PlanarLiquidSurface()) return min(ℋ * qᵛ⁺, qᵛ_max)endqᵛ_bg (generic function with 1 method)Zonal wind profile with linear shear below $zˢ$ and smooth transition (Equations 15-16):
function u_background(z) uˡ = uˢ * (z / zˢ) - uᶜ uᵗ = (-4/5 + 3 * (z / zˢ) - 5/4 * (z / zˢ)^2) * uˢ - uᶜ uᵘ = uˢ - uᶜ return (z < (zˢ - 1000)) * uˡ + (abs(z - zˢ) ≤ 1000) * uᵗ + (z > (zˢ + 1000)) * uᵘendu_background (generic function with 1 method)Warm bubble perturbation
The warm bubble parameters following Equations 17–18 in Klemp et al. (2015):
Δθ = 3 # K - perturbation amplituderᵇʰ = 10kilometers # m - bubble horizontal radiusrᵇᵛ = 1500 # m - bubble vertical radiuszᵇ = 1500 # m - bubble center heightxᵇ = Lx / 2 # m - bubble center x-coordinateyᵇ = Ly / 2 # m - bubble center y-coordinateThe total initial potential temperature combines the background profile with the cosine-squared warm bubble perturbation:
function θᵢ(x, y, z) θ̄ = θ_background(z) r = sqrt((x - xᵇ)^2 + (y - yᵇ)^2) R = sqrt((r / rᵇʰ)^2 + ((z - zᵇ) / rᵇᵛ)^2) θ′ = ifelse(R < 1, Δθ * cos(π * R / 2)^2, 0.0) return θ̄ + θ′enduᵢ(x, y, z) = u_background(z)uᵢ (generic function with 1 method)Initial-condition fields
We evaluate the background $qᵛ$ and the bubble-augmented $θ$ once on column / 3D fields and reuse them for both the anelastic and compressible runs. The hydrostatic integration in qᵛ_bg then runs only once per vertical level rather than once per horizontal grid point.
qᵛ_column = Field{Nothing, Nothing, Center}(grid)set!(qᵛ_column, qᵛ_bg)θ_background_column = Field{Nothing, Nothing, Center}(grid)set!(θ_background_column, θ_background)θ_initial_field = CenterField(grid)set!(θ_initial_field, θᵢ)168×168×40 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── grid: 168×168×40 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 178×178×50 OffsetArray(::CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}, -4:173, -4:173, -4:45) with eltype Float32 with indices -4:173×-4:173×-4:45
└── max=489.259, min=300.34, mean=357.189Visualization of initial conditions and warm bubble perturbation
We visualize the background potential temperature, water vapor mixing ratio, and wind shear profiles that define the environmental stratification:
θ_profile = set!(Field{Nothing, Nothing, Center}(grid), z -> θ_background(z))qᵛ_profile = set!(Field{Nothing, Nothing, Center}(grid), 1000 * qᵛ_column) # convert kg/kg -> g/kgu_profile = set!(Field{Nothing, Nothing, Center}(grid), z -> u_background(z))fig = Figure(size=(1000, 400), fontsize=14)axθ = Axis(fig[1, 1], xlabel="θ (K)", ylabel="z (km)", title="Potential temperature")lines!(axθ, θ_profile, linewidth=2, color=:magenta)hlines!(axθ, [zᵖ / 1000], color=:gray, linestyle=:dash)axqᵛ = Axis(fig[1, 2], xlabel="qᵛ (g/kg)", ylabel="z (km)", title="Water vapor mixing ratio")lines!(axqᵛ, qᵛ_profile, linewidth=2, color=:dodgerblue)hlines!(axqᵛ, [zᵖ / 1000], color=:gray, linestyle=:dash)axu = Axis(fig[1, 3], xlabel="u (m/s)", ylabel="z (km)", title="Wind profile")lines!(axu, u_profile, linewidth=2, color=:orangered)hlines!(axu, [zˢ / 1000], color=:gray, linestyle=:dash)vlines!(axu, [0], color=:black, linestyle=:dot)fig
Visualize the warm bubble perturbation on a vertical slice through the domain center:
θ′_slice = set!(Field{Center, Nothing, Center}(grid), (x, z) -> θᵢ(x, yᵇ, z) - θ_background(z))fig = Figure(size=(700, 400), fontsize=14)ax = Axis(fig[1, 1], xlabel="x (km)", ylabel="z (km)", title="Warm bubble perturbation θ′")hm = heatmap!(ax, θ′_slice, colormap=:thermal, colorrange=(0, Δθ))Colorbar(fig[1, 2], hm, label="θ′ (K)")fig
Reference state — anelastic
Breeze dynamics subtract a hydrostatically-balanced reference column from the prognostic state so that the time-tendency variables carry only the deviation from rest. For the anelastic core, density is a fixed background $\bar ρ(z)$ rather than a prognostic variable, and we construct the reference from a single surface potential temperature θ₀. Anelastic dynamics are insensitive to the exact reference stratification as long as buoyancy perturbations remain small compared with it.
reference_state_anelastic = ReferenceState(grid, constants; surface_pressure = p₀, potential_temperature = θ₀)ReferenceState{Float32}(p₀=100000.0, θ₀=300.0, pˢᵗ=100000.0)Dynamics — anelastic
AnelasticDynamics integrates an incompressible-with-stratification equation set that filters acoustic waves by construction through a pressure-Poisson solve at every RK substep:
dynamics_anelastic = AnelasticDynamics(reference_state_anelastic)AnelasticDynamics(p₀=100000.0, θ₀=300.0)
└── pressure_anomaly: not materializedMicrophysics
Kessler warm-rain microphysics carries prognostic cloud water qᶜˡ and rain water qʳ alongside vapor qᵛ, with autoconversion, accretion, rain evaporation, and sedimentation:
microphysics = DCMIP2016KesslerMicrophysics()Breeze.Microphysics.DCMIP2016KesslerMicrophysics{Float32}(237.3f0, 36.34f0, 0.001f0, 0.1364f0, 0.001f0, 0.001f0, 2.2f0, 0.875f0, 1.6f0, 124.9f0, 0.2046f0, 0.525f0, 2.55f8, 540000.0f0, 0.8f0)Advection
We use WENO advection at order 9. Klemp et al. (2015) note that supercell intensity and structure are highly sensitive to numerical diffusion; high-order WENO keeps it low without adding an explicit diffusion operator.
advection = WENO(order=9)WENO{5, Float32, Nothing}(order=9)
├── buffer_scheme: WENO{4, Float32, Nothing}(order=7)
│ └── buffer_scheme: WENO{3, Float32, Nothing}(order=5)
│ └── buffer_scheme: WENO{2, Float32, Nothing}(order=3)
│ └── buffer_scheme: Centered(order=2)
└── advecting_velocity_scheme: Centered(order=8)Building the anelastic model
AtmosphereModel ties together the grid, dynamics, microphysics, advection, and thermodynamic constants. The same constructor signature works for both dynamical cores — only the dynamics argument changes.
model_anelastic = AtmosphereModel(grid; dynamics = dynamics_anelastic, microphysics, advection, thermodynamic_constants = constants)AtmosphereModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 168×168×40 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── dynamics: AnelasticDynamics(p₀=100000.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float32}
├── timestepper: SSPRungeKutta3
├── advection scheme:
│ ├── momentum: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│ ├── ρθ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│ ├── ρqᵛ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│ ├── ρqᶜˡ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│ └── ρqʳ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
├── forcing: @NamedTuple{ρu::Returns{Float32}, ρv::Returns{Float32}, ρw::Returns{Float32}, ρθ::Returns{Float32}, ρqᵛ::Returns{Float32}, ρqᶜˡ::Returns{Float32}, ρqʳ::Returns{Float32}, ρe::Returns{Float32}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: DCMIP2016KesslerMicrophysicsInitializing the anelastic model
set! accepts pointwise functions (uᵢ) and pre-built fields (θ_initial_field, qᵛ_column). It calls update_state! internally so a single invocation refreshes all auxiliary diagnostics.
set!(model_anelastic, θ=θ_initial_field, qᵛ=qᵛ_column, u=uᵢ)Slice indices for output
We save horizontal slices at $z \approx 5 \, {\rm km}$ (mid-troposphere, where the rotating updraft is well-developed) and vertical slices through the bubble center at $y = L_y/2$. Both runs use the same indices.
z_centers = znodes(grid, Center())const k_5km = searchsortedfirst(z_centers, 5000)const j_center = Ny ÷ 2 + 1@info "Saving xy at z = $(z_centers[k_5km]) m (k = $k_5km); xz at y = $(Ly/2) m (j = $j_center)"[ Info: Saving xy at z = 5250.0 m (k = 11); xz at y = 84000.0 m (j = 85)
Simulation driver
run_simulation accepts an already-built AtmosphereModel, runs it for 2 hours with a CFL-controlled time-step wizard, periodically writes horizontal (xy at $z \approx 5 \, {\rm km}$) and vertical (xz at $y = L_y/2$) slices, and collects the maximum vertical-velocity time series. We use this same driver for both the anelastic and compressible runs.
function run_simulation(model, label) @info "=== Running case: $label ===" θ = liquid_ice_potential_temperature(model) θ′ = Field(θ - deepcopy(θ)) ## lazy field, recomputed each output qᶜˡ = model.microphysical_fields.qᶜˡ qʳ = model.microphysical_fields.qʳ qᵛ = model.microphysical_fields.qᵛ u, v, w = model.velocities simulation = Simulation(model; Δt=2, stop_time=2hours) conjure_time_step_wizard!(simulation, cfl=0.7) Oceananigans.Diagnostics.erroring_NaNChecker!(simulation) wall_clock = Ref(time_ns()) function progress(sim) elapsed = 1e-9 * (time_ns() - wall_clock[]) compute!(θ′) msg = @sprintf("[%s] Iter: %d, t: %s, Δt: %s, wall: %s, max|u|: %.2f, max w: %.2f, min w: %.2f, extrema(θ'): (%.2f, %.2f)", label, iteration(sim), prettytime(sim), prettytime(sim.Δt), prettytime(elapsed), maximum(abs, u), maximum(w), minimum(w), minimum(θ′), maximum(θ′)) msg *= @sprintf(", max(qᵛ): %.2e, max(qᶜˡ): %.2e, max(qʳ): %.2e", maximum(qᵛ), maximum(qᶜˡ), maximum(qʳ)) @info msg return nothing end add_callback!(simulation, progress, IterationInterval(100)) max_w_ts = Float64[] max_w_times = Float64[] function collect_max_w(sim) push!(max_w_times, time(sim)) push!(max_w_ts, maximum(w)) return nothing end add_callback!(simulation, collect_max_w, TimeInterval(1minutes)) slice_outputs = ( wxy = view(w, :, :, k_5km), qᶜˡxy = view(qᶜˡ, :, :, k_5km), qʳxy = view(qʳ, :, :, k_5km), wxz = view(w, :, j_center, :), θ′xz = view(θ′, :, j_center, :), qᶜˡxz = view(qᶜˡ, :, j_center, :), qʳxz = view(qʳ, :, j_center, :), ) slices_filename = "splitting_supercell_$(label)_slices.jld2" simulation.output_writers[:slices] = JLD2Writer(model, slice_outputs; filename=slices_filename, schedule=TimeInterval(2minutes), overwrite_existing=true) run!(simulation) wall_seconds = simulation.run_wall_time @info @sprintf("[%s] DONE. wall time = %.1f s (%.2f min) over %d iterations", label, wall_seconds, wall_seconds / 60, iteration(simulation)) return (; label, wall_seconds, slices_filename, max_w_ts, max_w_times, iterations = iteration(simulation))endrun_simulation (generic function with 1 method)Run the anelastic simulation:
results = Dict{String, Any}()results["anelastic"] = run_simulation(model_anelastic, "anelastic")(label = "anelastic", wall_seconds = 121.06437841800022, slices_filename = "splitting_supercell_anelastic_slices.jld2", max_w_ts = [0.0, 0.6458489894866943, 1.2690489292144775, 1.8292303085327148, 2.313955545425415, 2.7805943489074707, 3.455655574798584, 4.44831657409668, 5.74930477142334, 7.691911697387695, 10.316685676574707, 13.835740089416504, 17.5748233795166, 21.25345802307129, 26.14861488342285, 31.712913513183594, 35.89370346069336, 42.29914855957031, 45.563697814941406, 51.075439453125, 51.842620849609375, 54.47051239013672, 60.406219482421875, 66.54891204833984, 62.63206100463867, 64.96622467041016, 70.39079284667969, 67.93971252441406, 69.59407043457031, 70.38851165771484, 68.14762115478516, 66.47178649902344, 82.73602294921875, 76.72563171386719, 72.49556732177734, 73.287841796875, 72.44225311279297, 75.18572235107422, 75.24552917480469, 74.74150848388672, 75.3452377319336, 76.48583984375, 76.05244445800781, 69.89737701416016, 72.27714538574219, 76.59481811523438, 76.74520111083984, 77.40391540527344, 75.12325286865234, 78.39433288574219, 76.85580444335938, 78.89058685302734, 78.40741729736328, 73.13424682617188, 69.67193603515625, 69.94641876220703, 73.36933898925781, 73.35005950927734, 71.22443389892578, 68.33670806884766, 70.29177856445312, 67.93128204345703, 64.79580688476562, 64.27740478515625, 59.82387924194336, 56.97315979003906, 55.37782287597656, 56.638206481933594, 56.43620300292969, 59.55274963378906, 58.50358581542969, 60.149330139160156, 61.785438537597656, 61.86437225341797, 64.61840057373047, 66.19480895996094, 62.586421966552734, 62.490325927734375, 64.22573852539062, 56.826171875, 52.588539123535156, 54.99455261230469, 54.586524963378906, 51.8397216796875, 59.9561767578125, 62.20354080200195, 60.22412109375, 56.85688400268555, 57.22355270385742, 53.79008865356445, 55.91199493408203, 53.540794372558594, 53.08961486816406, 55.5501708984375, 50.87851333618164, 52.61809158325195, 57.645843505859375, 60.609073638916016, 58.54133605957031, 58.129844665527344, 57.7296028137207, 60.91950225830078, 60.6228141784668, 62.536373138427734, 59.93550109863281, 61.18586730957031, 59.82748794555664, 59.64795684814453, 55.816795349121094, 61.52462387084961, 68.08128356933594, 62.77248001098633, 59.22852325439453, 58.783973693847656, 54.26276397705078, 54.10453414916992, 55.61895751953125, 55.59037399291992, 53.59627914428711, 52.93680191040039, 54.456565856933594], max_w_times = [0.0, 60.0, 120.0, 180.0, 240.0, 300.0, 360.0, 420.0, 480.0, 540.0, 600.0, 660.0, 720.0, 780.0, 840.0, 900.0, 960.0, 1020.0, 1080.0, 1140.0, 1200.0, 1260.0, 1320.0, 1380.0, 1440.0, 1500.0, 1560.0, 1620.0, 1680.0, 1740.0, 1800.0, 1860.0, 1920.0, 1980.0, 2040.0, 2100.0, 2160.0, 2220.0, 2280.0, 2340.0, 2400.0, 2460.0, 2520.0, 2580.0, 2640.0, 2700.0, 2760.0, 2820.0, 2880.0, 2940.0, 3000.0, 3060.0, 3120.0, 3180.0, 3240.0, 3300.0, 3360.0, 3420.0, 3480.0, 3540.0, 3600.0, 3660.0, 3720.0, 3780.0, 3840.0, 3900.0, 3960.0, 4020.0, 4080.0, 4140.0, 4200.0, 4260.0, 4320.0, 4380.0, 4440.0, 4500.0, 4560.0, 4620.0, 4680.0, 4740.0, 4800.0, 4860.0, 4920.0, 4980.0, 5040.0, 5100.0, 5160.0, 5220.0, 5280.0, 5340.0, 5400.0, 5460.0, 5520.0, 5580.0, 5640.0, 5700.0, 5760.0, 5820.0, 5880.0, 5940.0, 6000.0, 6060.0, 6120.0, 6180.0, 6240.0, 6300.0, 6360.0, 6420.0, 6480.0, 6540.0, 6600.0, 6660.0, 6720.0, 6780.0, 6840.0, 6900.0, 6960.0, 7020.0, 7080.0, 7140.0, 7200.0], iterations = 1652)Now the compressible core
CompressibleDynamics integrates the fully compressible Euler equations and resolves acoustic waves explicitly via split-explicit substepping: a small Δt advances sound and buoyancy oscillations, while the standard CFL Δt advances advection and physics. Density is now a prognostic variable, so the reference state must closely match the actual atmosphere. We pass the same θ_background(z) and qᵛ_bg(z) profiles used in the initial condition; CompressibleDynamics builds an ExnerReferenceState internally that satisfies discrete hydrostatic balance to machine precision, so the slow vertical-momentum tendency vanishes on a rest atmosphere.
dynamics_compressible = CompressibleDynamics(SplitExplicitTimeDiscretization(); surface_pressure = p₀, standard_pressure = pˢᵗ, reference_potential_temperature = θ_background, reference_vapor_mass_fraction = qᵛ_bg)CompressibleDynamics{SplitExplicitTimeDiscretization}
├── dry_density: not materialized
├── pressure: not materialized
├── terrain_metrics: Breeze.TerrainFollowingDiscretization.SlopeOutsideInterpolation
├── time_discretization: SplitExplicitTimeDiscretization
└── reference_state: @NamedTuple{reference_potential_temperature::typeof(Main.var"##277".θ_background), reference_vapor_mass_fraction::typeof(Main.var"##277".qᵛ_bg)}Build the model — same constructor, different dynamics:
model_compressible = AtmosphereModel(grid; dynamics = dynamics_compressible, microphysics, advection, thermodynamic_constants = constants)AtmosphereModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 168×168×40 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── dynamics: CompressibleDynamics{SplitExplicitTimeDiscretization}
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float32}
├── timestepper: AcousticRungeKutta3
├── advection scheme:
│ ├── momentum: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│ ├── ρθ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│ ├── ρqᵛ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│ ├── ρqᶜˡ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│ └── ρqʳ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
├── forcing: @NamedTuple{ρᵈ::Returns{Float32}, ρu::Returns{Float32}, ρv::Returns{Float32}, ρw::Returns{Float32}, ρθ::Returns{Float32}, ρqᵛ::Returns{Float32}, ρqᶜˡ::Returns{Float32}, ρqʳ::Returns{Float32}, ρe::Returns{Float32}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: DCMIP2016KesslerMicrophysicsInitial density. Naively setting ρ to the reference density while perturbing θ would change ρθ and seed a spurious acoustic pulse from the warm bubble. Instead we rescale the reference density via pressure_balanced_density so that ρθ (and therefore the equation-of-state pressure) is unchanged at t = 0:
ρ_initial = CenterField(grid)set!(ρ_initial, pressure_balanced_density(model_compressible.dynamics.reference_state.density, θ_background_column, θ_initial_field))168×168×40 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── grid: 168×168×40 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 178×178×50 OffsetArray(::CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}, -4:173, -4:173, -4:45) with eltype Float32 with indices -4:173×-4:173×-4:45
└── max=1.12725, min=0.095289, mean=0.480205Set the same θ, qᵛ, u initial state plus the balanced ρ:
set!(model_compressible, θ=θ_initial_field, qᵛ=qᵛ_column, u=uᵢ, ρ=ρ_initial)Run the compressible simulation:
results["compressible"] = run_simulation(model_compressible, "compressible")(label = "compressible", wall_seconds = 164.81950471499974, slices_filename = "splitting_supercell_compressible_slices.jld2", max_w_ts = [0.0, 0.6768927574157715, 1.3134430646896362, 1.8922226428985596, 2.3776144981384277, 2.727368116378784, 3.1537787914276123, 3.895456552505493, 4.7892985343933105, 6.252162456512451, 8.02805233001709, 10.770435333251953, 13.934611320495605, 17.90035629272461, 23.103870391845703, 26.289026260375977, 29.372283935546875, 30.003835678100586, 32.524009704589844, 36.969810485839844, 41.579856872558594, 43.058258056640625, 43.94227981567383, 41.9036979675293, 49.632347106933594, 56.447837829589844, 60.34217071533203, 62.92866897583008, 60.1523323059082, 61.782989501953125, 62.98213577270508, 64.43629455566406, 62.93006134033203, 64.89041900634766, 68.26017761230469, 61.69541931152344, 66.83203125, 65.58805847167969, 65.78307342529297, 64.5566177368164, 63.67531967163086, 73.617431640625, 72.09612274169922, 74.4419937133789, 71.78388977050781, 71.78034973144531, 66.44620513916016, 63.81211853027344, 69.87799835205078, 72.47650146484375, 69.48582458496094, 69.9297103881836, 66.65293884277344, 67.37059783935547, 67.79581451416016, 70.44039916992188, 68.7064208984375, 66.49777221679688, 66.09059143066406, 66.69266510009766, 66.01874542236328, 65.32340240478516, 64.69502258300781, 61.28749084472656, 61.62250518798828, 54.353878021240234, 51.94252395629883, 50.76690673828125, 50.54977035522461, 50.45380783081055, 51.32086181640625, 53.81473159790039, 53.00284194946289, 56.002418518066406, 57.167930603027344, 54.87980270385742, 55.23634338378906, 51.76028060913086, 51.35083770751953, 50.48495864868164, 44.977909088134766, 43.393470764160156, 43.07719421386719, 44.28581619262695, 51.44919204711914, 59.06858825683594, 59.029903411865234, 57.5195198059082, 56.53203201293945, 52.082401275634766, 52.32869338989258, 53.1827278137207, 53.8145866394043, 54.36850357055664, 49.02791213989258, 50.468101501464844, 55.65583419799805, 62.84329605102539, 64.68345642089844, 62.0468864440918, 58.135643005371094, 48.72262954711914, 49.50233840942383, 46.428035736083984, 46.1258430480957, 47.248897552490234, 47.88462448120117, 46.70795440673828, 47.66108322143555, 55.55215072631836, 54.43608093261719, 53.03281021118164, 52.42317581176758, 53.48455047607422, 54.39635467529297, 54.61353302001953, 56.53194046020508, 55.06249237060547, 52.81375503540039, 52.43895721435547, 55.35820007324219], max_w_times = [0.0, 60.0, 120.0, 180.0, 240.0, 300.0, 360.0, 420.0, 480.0, 540.0, 600.0, 660.0, 720.0, 780.0, 840.0, 900.0, 960.0, 1020.0, 1080.0, 1140.0, 1200.0, 1260.0, 1320.0, 1380.0, 1440.0, 1500.0, 1560.0, 1620.0, 1680.0, 1740.0, 1800.0, 1860.0, 1920.0, 1980.0, 2040.0, 2100.0, 2160.0, 2220.0, 2280.0, 2340.0, 2400.0, 2460.0, 2520.0, 2580.0, 2640.0, 2700.0, 2760.0, 2820.0, 2880.0, 2940.0, 3000.0, 3060.0, 3120.0, 3180.0, 3240.0, 3300.0, 3360.0, 3420.0, 3480.0, 3540.0, 3600.0, 3660.0, 3720.0, 3780.0, 3840.0, 3900.0, 3960.0, 4020.0, 4080.0, 4140.0, 4200.0, 4260.0, 4320.0, 4380.0, 4440.0, 4500.0, 4560.0, 4620.0, 4680.0, 4740.0, 4800.0, 4860.0, 4920.0, 4980.0, 5040.0, 5100.0, 5160.0, 5220.0, 5280.0, 5340.0, 5400.0, 5460.0, 5520.0, 5580.0, 5640.0, 5700.0, 5760.0, 5820.0, 5880.0, 5940.0, 6000.0, 6060.0, 6120.0, 6180.0, 6240.0, 6300.0, 6360.0, 6420.0, 6480.0, 6540.0, 6600.0, 6660.0, 6720.0, 6780.0, 6840.0, 6900.0, 6960.0, 7020.0, 7080.0, 7140.0, 7200.0], iterations = 1488)Wall-time summary
println("\n========== Wall-time summary ==========")@printf("%-14s %12s %12s %12s\n", "case", "wall (s)", "wall (min)", "iterations")for label in ("anelastic", "compressible") r = results[label] @printf("%-14s %12.1f %12.2f %12d\n", label, r.wall_seconds, r.wall_seconds / 60, r.iterations)end
========== Wall-time summary ==========
case wall (s) wall (min) iterations
anelastic 121.1 2.02 1652
compressible 164.8 2.75 1488
Horizontal slice comparison (z ≈ 5 km)
Top row: anelastic, bottom row: compressible. Columns show vertical velocity $w$, cloud water $qᶜˡ$, and rain water $qʳ$. The simulated supercell exhibits splitting behavior, with the initial storm dividing into right- and left-moving cells, consistent with the DCMIP2016 intercomparison results (Zarzycki et al., 2019).
xy_ts = Dict(label => ( wxy = FieldTimeSeries(results[label].slices_filename, "wxy"), qᶜˡxy = FieldTimeSeries(results[label].slices_filename, "qᶜˡxy"), qʳxy = FieldTimeSeries(results[label].slices_filename, "qʳxy"), ) for label in ("anelastic", "compressible"))wlim = maximum(maximum(abs, xy_ts[l].wxy) for l in keys(xy_ts)) / 2qᶜˡlim = maximum(maximum(xy_ts[l].qᶜˡxy) for l in keys(xy_ts)) / 4qʳlim = maximum(maximum(xy_ts[l].qʳxy) for l in keys(xy_ts)) / 4times_xy = xy_ts["anelastic"].wxy.timesNt_xy = min(length(times_xy), length(xy_ts["compressible"].wxy.times))fig = Figure(size=(1200, 750), fontsize=12)fig[1, 1] = Label(fig, "anelastic", rotation=π/2, fontsize=14, tellheight=false)fig[2, 1] = Label(fig, "compressible", rotation=π/2, fontsize=14, tellheight=false)axw_anel = Axis(fig[1, 2], aspect=1, xlabel="x (m)", ylabel="y (m)", title="w (m/s)")axqᶜˡ_anel = Axis(fig[1, 4], aspect=1, xlabel="x (m)", ylabel="y (m)", title="qᶜˡ (kg/kg)")axqʳ_anel = Axis(fig[1, 6], aspect=1, xlabel="x (m)", ylabel="y (m)", title="qʳ (kg/kg)")axw_comp = Axis(fig[2, 2], aspect=1, xlabel="x (m)", ylabel="y (m)")axqᶜˡ_comp = Axis(fig[2, 4], aspect=1, xlabel="x (m)", ylabel="y (m)")axqʳ_comp = Axis(fig[2, 6], aspect=1, xlabel="x (m)", ylabel="y (m)")n_xy = Observable(1)wxy_anel_n = @lift xy_ts["anelastic"].wxy[$n_xy]qᶜˡxy_anel_n = @lift xy_ts["anelastic"].qᶜˡxy[$n_xy]qʳxy_anel_n = @lift xy_ts["anelastic"].qʳxy[$n_xy]wxy_comp_n = @lift xy_ts["compressible"].wxy[$n_xy]qᶜˡxy_comp_n = @lift xy_ts["compressible"].qᶜˡxy[$n_xy]qʳxy_comp_n = @lift xy_ts["compressible"].qʳxy[$n_xy]title_xy = @lift "Splitting supercell, xy at z ≈ 5 km, t = " * prettytime(times_xy[$n_xy])hmw_anel = heatmap!(axw_anel, wxy_anel_n, colormap=:balance, colorrange=(-wlim, wlim))hmqᶜˡ_anel = heatmap!(axqᶜˡ_anel, qᶜˡxy_anel_n, colormap=:dense, colorrange=(0, qᶜˡlim))hmqʳ_anel = heatmap!(axqʳ_anel, qʳxy_anel_n, colormap=:amp, colorrange=(0, qʳlim)) heatmap!(axw_comp, wxy_comp_n, colormap=:balance, colorrange=(-wlim, wlim)) heatmap!(axqᶜˡ_comp, qᶜˡxy_comp_n, colormap=:dense, colorrange=(0, qᶜˡlim)) heatmap!(axqʳ_comp, qʳxy_comp_n, colormap=:amp, colorrange=(0, qʳlim))Colorbar(fig[1:2, 3], hmw_anel)Colorbar(fig[1:2, 5], hmqᶜˡ_anel)Colorbar(fig[1:2, 7], hmqʳ_anel)fig[0, :] = Label(fig, title_xy, fontsize=14, tellwidth=false)CairoMakie.record(fig, "splitting_supercell_xy_comparison.mp4", 1:Nt_xy; framerate = 10, compression = 23) do nn n_xy[] = nnendVertical slice comparison (y = Ly/2)
Vertical (xz) slice cut through the bubble center. Columns show $w$ and the potential-temperature perturbation $θ' = θ - θ_{\rm initial}$; rows show anelastic (top) vs compressible (bottom).
xz_ts = Dict(label => ( wxz = FieldTimeSeries(results[label].slices_filename, "wxz"), θ′xz = FieldTimeSeries(results[label].slices_filename, "θ′xz"), ) for label in ("anelastic", "compressible"))wlim_xz = maximum(maximum(abs, xz_ts[l].wxz) for l in keys(xz_ts)) / 2θ′lim_xz = maximum(maximum(abs, xz_ts[l].θ′xz) for l in keys(xz_ts)) / 2times_xz = xz_ts["anelastic"].wxz.timesNt_xz = min(length(times_xz), length(xz_ts["compressible"].wxz.times))fig = Figure(size=(1100, 700), fontsize=12)fig[1, 1] = Label(fig, "anelastic", rotation=π/2, fontsize=14, tellheight=false)fig[2, 1] = Label(fig, "compressible", rotation=π/2, fontsize=14, tellheight=false)axw_anel_xz = Axis(fig[1, 2], xlabel="x (m)", ylabel="z (m)", title="w (m/s)")axθ_anel_xz = Axis(fig[1, 4], xlabel="x (m)", ylabel="z (m)", title="θ' (K)")axw_comp_xz = Axis(fig[2, 2], xlabel="x (m)", ylabel="z (m)")axθ_comp_xz = Axis(fig[2, 4], xlabel="x (m)", ylabel="z (m)")n_xz = Observable(1)wxz_anel_n = @lift xz_ts["anelastic"].wxz[$n_xz]θ′xz_anel_n = @lift xz_ts["anelastic"].θ′xz[$n_xz]wxz_comp_n = @lift xz_ts["compressible"].wxz[$n_xz]θ′xz_comp_n = @lift xz_ts["compressible"].θ′xz[$n_xz]title_xz = @lift "Splitting supercell, xz at y = Ly/2, t = " * prettytime(times_xz[$n_xz])hmw_anel_xz = heatmap!(axw_anel_xz, wxz_anel_n, colormap=:balance, colorrange=(-wlim_xz, wlim_xz))hmθ_anel_xz = heatmap!(axθ_anel_xz, θ′xz_anel_n, colormap=:balance, colorrange=(-θ′lim_xz, θ′lim_xz)) heatmap!(axw_comp_xz, wxz_comp_n, colormap=:balance, colorrange=(-wlim_xz, wlim_xz)) heatmap!(axθ_comp_xz, θ′xz_comp_n, colormap=:balance, colorrange=(-θ′lim_xz, θ′lim_xz))Colorbar(fig[1:2, 3], hmw_anel_xz)Colorbar(fig[1:2, 5], hmθ_anel_xz)fig[0, :] = Label(fig, title_xz, fontsize=14, tellwidth=false)CairoMakie.record(fig, "splitting_supercell_xz_comparison.mp4", 1:Nt_xz; framerate = 10, compression = 23) do nn n_xz[] = nnendMaximum vertical velocity time series
The maximum updraft velocity is a key diagnostic for supercell intensity. Strong supercells typically develop updrafts exceeding 30–50 m/s. As noted by Klemp et al. (2015), the simulated storm intensity and structure are highly sensitive to numerical diffusion; no explicit numerical diffusion is applied here. Plotting both runs together highlights how closely the two dynamical cores agree (or where they diverge) under identical microphysics and initial conditions.
fig = Figure(size=(700, 400), fontsize=14)ax = Axis(fig[1, 1], xlabel="Time (s)", ylabel="Maximum w (m/s)", title="Maximum vertical velocity", xticks=0:1800:7200)lines!(ax, results["anelastic"].max_w_times, results["anelastic"].max_w_ts, linewidth=2, color=:dodgerblue, label="anelastic")lines!(ax, results["compressible"].max_w_times, results["compressible"].max_w_ts, linewidth=2, color=:orangered, label="compressible")axislegend(ax, position=:lt)fig
Julia version and environment information
This example was executed with the following version of Julia:
using InteractiveUtils: versioninfoversioninfo()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 PkgPkg.status()Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
⌃ [86bc3604] AtmosphericProfilesLibrary v0.1.8
[660aa2fb] Breeze v0.8.0 `.`
[052768ef] CUDA v6.2.1
[13f3f980] CairoMakie v0.15.13
⌅ [6a9e3e04] CloudMicrophysics v0.37.1
[e30172f5] Documenter v1.17.0
[daee34ce] DocumenterCitations v1.4.1
[b6400b83] DocumenterCodeBlocks v1.1.0
[7da242da] Enzyme v0.13.198
[46192b85] GPUArraysCore v0.2.0
[63c18a36] KernelAbstractions v0.9.42
[98b081ad] Literate v2.21.0
[85f8d34a] NCDatasets v0.14.15
[9e8cae18] Oceananigans v0.110.14
⌅ [a01a1ee8] RRTMGP v0.21.9
[3c362404] Reactant v0.2.278
[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.