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.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
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: WENO{3, Float64, Float32}(order=5)
│   ├── ρθ: WENO{3, Float64, Float32}(order=5)
│   └── ρqᵗ: WENO{3, Float64, Float32}(order=5)
├── 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(model.specific_moisture, 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="Total moisture")
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)

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
├── 7 outputs: (u, v, w, qᵛ, qˡ, ξ, θ)
├── array_type: Array{Float32}
├── including: [:grid]
├── file_splitting: NoFileSplitting
└── file size: 56.5 KiB

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 (22.763 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.545 seconds).
[ Info: iteration: 200, time: 2.224 minutes, Δt: 728.146 ms, max|w|: 4.15e-01 m/s
[ Info: iteration: 400, time: 4.447 minutes, Δt: 703.651 ms, max|w|: 1.25e+00 m/s
[ Info: iteration: 600, time: 6.510 minutes, Δt: 645.425 ms, max|w|: 2.74e+00 m/s
[ Info: iteration: 800, time: 8.322 minutes, Δt: 546.826 ms, max|w|: 8.72e+00 m/s
[ Info: iteration: 1000, time: 9.890 minutes, Δt: 469.980 ms, max|w|: 1.21e+01 m/s
[ Info: iteration: 1200, time: 11.370 minutes, Δt: 442.280 ms, max|w|: 1.21e+01 m/s
[ Info: Simulation is stopping after running for 3.985 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.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.6
  [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
  [85f8d34a] NCDatasets v0.14.10
  [9e8cae18] Oceananigans v0.103.1
  [a01a1ee8] RRTMGP v0.21.6
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.0
  [9a3f8284] Random v1.11.0

This page was generated using Literate.jl.