Cloudy Kelvin-Helmholtz instability

This example sets up a two-dimensional ($x$–$z$) Kelvin–Helmholtz instability in a moist, stably stratified atmosphere.

The configuration is intentionally simple but reasonably "meteorological":

• We impose horizontal wind $U(z)$ with a shear layer.
• We impose a stably stratified potential temperature profile $θ(z)$ with a specified dry Brunt–Väisälä frequency $N$.
• We embed a Gaussian moisture layer $q(z)$ centered on the shear layer.

As the shear layer rolls up, the moist layer is advected and deformed, producing billow-like patterns reminiscent of observed "wave clouds". Breeze encapsulates much of this thermodynamics for us via the AtmosphereModel and saturation adjustment.

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

Random.seed!(301)

Random.TaskLocalRNG()

Domain and grid

We use a 2D $x$–$z$ slice with periodic boundaries in $x$ and rigid, impermeable boundaries at the top and bottom.

Grid resolution is modest but enough to clearly resolve the Kelvin-Helmholtz billows and rolled-up moisture filament.

Nx, Nz = 384, 128   # resolution
Lx, Lz = 10e3, 3e3  # domain extent

grid = RectilinearGrid(CPU(); size = (Nx, Nz), x = (0, Lx), z = (0, Lz),
                       topology = (Periodic, Flat, Bounded))

384×1×128 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [0.0, 10000.0) regularly spaced with Δx=26.0417
├── Flat y                      
└── Bounded  z ∈ [0.0, 3000.0]  regularly spaced with Δz=23.4375

Model and microphysics

We construct the AtmosphereModel model with saturation adjustment microphysics.

microphysics = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())
model = AtmosphereModel(grid; advection=WENO(order=5), microphysics)

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

Background thermodynamic state

We set a reference potential temperature $θ₀$ and a linear $θ$ gradient that corresponds to a desired dry Brunt–Väisälä frequency $N$. For a dry atmosphere,

\[N² = \frac{g}{θ₀} \frac{∂θ}{∂z} ,\]

We initialize with a potential temperature that gives constant Brunt–Väisälä frequency, representative of mid-tropospheric stability. The (dry) Brunt–Väisälä frequency is

\[N² = \frac{g}{θ} \frac{∂θ}{∂z}\]

and thus, for constant $N²$, the above implies that $θ = θ₀ \exp{(N² z / g)}$.

thermo = ThermodynamicConstants()
g = thermo.gravitational_acceleration
θ₀ = model.dynamics.reference_state.potential_temperature
N = 0.01                  # target dry Brunt–Väisälä frequency (s⁻¹)
θᵇ(z) = θ₀ * exp(N^2 * z / g)

θᵇ (generic function with 1 method)

Shear and moisture profiles

We want:

• A shear layer centered at height $z₀$ with the zonal flow transitioning from a lower speed $U_{\rm bot}$ to an upper speed $U_{\rm top}$.
• A moist layer centered at the same height with a Gaussian relative humidity profile.

The above mimics a moist, stably stratified layer embedded in stronger flow above and weaker flow below.

First, we set up the shear layer using a $\tanh$ profile:

z₀  = 1e3  # center of shear & moist layer (m)
Δzu = 150  # shear layer half-thickness (m)
U₀  =  5   # base wind speed (m/s)
ΔU  = 20   # upper-layer wind (m/s)
uᵇ(z) = U₀ + ΔU * (1 + tanh((z - z₀) / Δzu)) / 2

uᵇ (generic function with 1 method)

For the moisture layer, we specify a Gaussian relative humidity profile centered at $z₀$. The peak relative humidity is supersaturated ($ℋ₀ > 1$), which triggers immediate cloud formation via saturation adjustment.

ℋ₀  = 1.6    # peak relative humidity (supersaturated)
Δzℋ = 200    # moist layer half-width (m)
ℋᵇ(x, z) = ℋ₀ * exp(-(z - z₀)^2 / 2Δzℋ^2)

ℋᵇ (generic function with 1 method)

We initialize the model via Oceananigans set!, adding also a bit of random noise. Note that we use the ℋ keyword argument to set moisture via relative humidity.

δθ = 0.01
δu = 1e-3
δℋ = 0.05

ϵ() = rand() - 1/2
θᵢ(x, z) = θᵇ(z) + δθ * ϵ()
uᵢ(x, z) = uᵇ(z) + δu * ϵ()
ℋᵢ(x, z) = ℋᵇ(x, z) + δℋ * ϵ()

set!(model; u=uᵢ, θ=θᵢ, ℋ=ℋᵢ)

The Kelvin-Helmholtz instability

The Miles–Howard criterion tells us that Kelvin–Helmholtz instability occurs where the Richardson number,

\[Ri = \frac{N²}{(∂uᵇ/∂z)²}\]

is less than 1/4 (Miles, 1961; Howard, 1961). With the parameters chosen above this is the case.

Let's plot the initial state as well as the Richardson number and relative humidity.

U = Field(Average(model.velocities.u, dims=(1, 2)))
Ri = N^2 / ∂z(U)^2

Qᵛ = Field(Average(specific_humidity(model), dims=1))
θ = Field(Average(liquid_ice_potential_temperature(model), dims=1))
ℋ = Field(Average(RelativeHumidity(model), dims=1))

fig = Figure(size=(1000, 500))

axu = Axis(fig[1, 1], xlabel = "uᵇ (m/s)", ylabel = "z (m)", title = "Zonal velocity")
axq = Axis(fig[1, 2], xlabel = "qᵛ (kg/kg)", title="Specific humidity")
axℋ = Axis(fig[1, 3], xlabel = "ℋ", title="Relative humidity")
axθ = Axis(fig[1, 4], xlabel = "θ (K)", title="Potential temperature")
axR = Axis(fig[1, 5], xlabel = "Ri", title="Richardson number")

lines!(axu, U)
lines!(axq, Qᵛ)
lines!(axℋ, ℋ)
lines!(axθ, θ)
lines!(axR, Ri)
lines!(axR, [1/4, 1/4], [0, Lz], linestyle = :dash, color = :black)

xlims!(axR, 0, 0.8)
axR.xticks = 0:0.25:1

for ax in (axq, axℋ, axθ, axR)
    ax.yticksvisible = false
    ax.yticklabelsvisible = false
    ax.ylabelvisible = false
end

fig

Set up and run the simulation

We construct a simulation and use the time-step wizard to keep the CFL number under control.

stop_time = 12minutes   # total simulation time
simulation = Simulation(model; Δt=1, stop_time)
conjure_time_step_wizard!(simulation; cfl = 0.7)
Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)

We also add a progress callback:

function progress(sim)
    u, v, w = model.velocities
    max_w = maximum(abs, w)
    @info @sprintf("iteration: %d, time: %s, Δt: %s, max|w|: %.2e m/s",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt), max_w)
    return nothing
end

add_callback!(simulation, progress, IterationInterval(200))

Output

We save the model velocities, the cross-stream component of vorticity, $ξ = ∂_z u - ∂_x w$, the potential temperatures and the specific humidities (vapour, liquid, ice).

u, v, w = model.velocities
ξ = ∂z(u) - ∂x(w)
θ = liquid_ice_potential_temperature(model)
outputs = merge(model.velocities, model.microphysical_fields, (; ξ, θ))

filename = "wave_clouds.jld2"

output_writer = JLD2Writer(model, outputs; filename,
                           schedule = TimeInterval(4),
                           overwrite_existing = true)

simulation.output_writers[:fields] = output_writer

JLD2Writer scheduled on TimeInterval(4 seconds):
├── filepath: wave_clouds.jld2
├── 8 outputs: (u, v, w, qᵛ, qˡ, qᵉ, ξ, θ)
├── array_type: Array{Float32}
├── including: [:thermodynamic_constants]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

Run!

Now we are ready to run the simulation.

run!(simulation)

[ Info: Initializing simulation...
[ Info: iteration: 0, time: 0 seconds, Δt: 729.144 ms, max|w|: 4.19e-04 m/s
[ Info:     ... simulation initialization complete (33.468 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.543 seconds).
[ Info: iteration: 200, time: 2.224 minutes, Δt: 728.140 ms, max|w|: 4.18e-01 m/s
[ Info: iteration: 400, time: 4.447 minutes, Δt: 702.729 ms, max|w|: 1.26e+00 m/s
[ Info: iteration: 600, time: 6.499 minutes, Δt: 642.906 ms, max|w|: 2.73e+00 m/s
[ Info: iteration: 800, time: 8.312 minutes, Δt: 545.514 ms, max|w|: 8.74e+00 m/s
[ Info: iteration: 1000, time: 9.875 minutes, Δt: 470.163 ms, max|w|: 1.21e+01 m/s
[ Info: iteration: 1200, time: 11.356 minutes, Δt: 441.857 ms, max|w|: 1.21e+01 m/s
[ Info: Simulation is stopping after running for 3.889 minutes.
[ Info: Simulation time 12 minutes equals or exceeds stop time 12 minutes.

Read output and visualize

We load the saved output as Oceananigans' FieldTimeSeries

ξt = FieldTimeSeries(filename, "ξ")
θt = FieldTimeSeries(filename, "θ")
qˡt = FieldTimeSeries(filename, "qˡ")

times = ξt.times
Nt = length(ξt)

181

and then use CairoMakie to plot and animate the output.

n = Observable(Nt)

ξn = @lift ξt[$n]
θn = @lift θt[$n]
qˡn = @lift qˡt[$n]

fig = Figure(size=(800, 800), fontsize=14)

axξ = Axis(fig[1, 1], ylabel="z (m)", title = "Vorticity", titlesize = 20)
axl = Axis(fig[2, 1], ylabel="z (m)", title = "Liquid mass fraction", titlesize = 20)
axθ = Axis(fig[3, 1], xlabel="x (m)", ylabel="z (m)", title = "Potential temperature", titlesize = 20)

hmξ = heatmap!(axξ, ξn, colormap = :balance, colorrange = (-0.25, 0.25))
hml = heatmap!(axl, qˡn, colormap = Reverse(:Blues_4), colorrange = (0, 0.003))
hmθ = heatmap!(axθ, θn, colormap = :thermal, colorrange = extrema(θt))

Colorbar(fig[1, 2], hmξ, label = "s⁻¹", vertical = true)
Colorbar(fig[2, 2], hml, label = "kg/kg", vertical = true)
Colorbar(fig[3, 2], hmθ, label = "Κ", vertical = true)

fig

We can also make a movie:

CairoMakie.record(fig, "wave_clouds.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.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 Pkg
Pkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
  [86bc3604] AtmosphericProfilesLibrary v0.1.7
  [660aa2fb] Breeze v0.4.7 `.`
⌃ [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.107.3
  [a01a1ee8] RRTMGP v0.21.7
  [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.