Tropical Cyclone World (Cronin & Chavas, 2019)

This example implements the rotating radiative-convective equilibrium (RCE) experiment from Cronin and Chavas (2019). The experiment demonstrates that tropical cyclones can form and persist even in completely dry atmospheres, challenging the conventional wisdom that moisture is essential for TC dynamics.

The key innovation is the surface wetness parameter β, which controls the transition from completely dry (β = 0, no evaporation) to fully moist (β = 1) conditions. Cronin and Chavas (2019) found that TCs form in both limits, with a "no-storms-land" at intermediate β where spontaneous genesis does not occur. This script defaults to β = 1 (moist), which produces robust spontaneous TC genesis at moderate resolution. The simulation approximates the paper's 100-day nonrotating RCE spinup with an equilibrated initial temperature profile (dry adiabat in the troposphere, isothermal stratosphere) and uses warm-phase saturation adjustment microphysics for the moist case.

using Breeze
using Breeze.Thermodynamics: compute_reference_state!
using Oceananigans: Oceananigans
using Oceananigans.Units

using CairoMakie
using CUDA
using Printf
using Random

Random.seed!(2019)
Oceananigans.defaults.FloatType = Float32

Float32

Domain and grid

Cronin and Chavas (2019) used a 1152 km × 1152 km domain with 2 km horizontal resolution. To reduce computational costs for the purpose of this example, we use a 288 km × 288 km domain – 4x smaller in both horizontal directions – with a 2x coarser 4 km horizontal resolution. We keep the 28 km model top, but with 40 m spacing in the lowest kilometers rather than ~16 m, and 1000 m spacing above 3.5 km rather than 500 m (and a smooth transition in between).

arch = GPU()
paper_Lx = 1152e3
paper_Nx = 576
Lx = Ly = paper_Lx / 4
Nx = Ny = paper_Nx / 8 |> Int
H = 28e3

Δz_fine = 40 # m
Δz_coarse = 1000 # m

z = PiecewiseStretchedDiscretization(
    z  = [0, 1000, 3500, H],
    Δz = [Δz_fine, Δz_fine, Δz_coarse, Δz_coarse])

Nz = length(z) - 1

grid = RectilinearGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 5),
                       x = (0, Lx), y = (0, Ly), z,
                       topology = (Periodic, Periodic, Bounded))

72×72×60 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── Periodic x ∈ [0.0, 288000.0) regularly spaced with Δx=4000.0
├── Periodic y ∈ [0.0, 288000.0) regularly spaced with Δy=4000.0
└── Bounded  z ∈ [0.0, 28000.0]  variably spaced with min(Δz)=40.0, max(Δz)=1000.0

Reference state and dynamics

We use the anelastic formulation with a reference state initialized from the surface potential temperature T₀ = 300 K and standard surface pressure. The reference state is then adjusted to match the initial temperature and moisture profiles. This adjustment is critical for tall domains: without it, the constant-θ adiabat reference state diverges from the actual atmosphere in the stratosphere (Tref ≈ 26 K vs Tactual = 210 K at 28 km), producing catastrophic buoyancy forces.

T₀ = 300
constants = ThermodynamicConstants()

reference_state = ReferenceState(grid, constants;
                                 surface_pressure = 101325,
                                 potential_temperature = T₀,
                                 vapor_mass_fraction = 0)

ReferenceState{Float32}(p₀=101325.0, θ₀=300.0, pˢᵗ=100000.0)

Define equilibrium temperature and moisture profiles for adjustment and initialization

Tᵗˢ = 210
cᵖᵈ = constants.dry_air.heat_capacity
g = constants.gravitational_acceleration
Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(constants)
κ = Rᵈ / cᵖᵈ
pˢᵗ = reference_state.standard_pressure
Π₀ = (101325 / pˢᵗ)^κ

1.0037661f0

Analytical Exner function for a hydrostatic constant-θ atmosphere

Π(z) = Π₀ - g * z / (cᵖᵈ * T₀)

β = 1
q₀ = 15e-3 # surface specific humidity (kg/kg)
Hq = 3000   # moisture scale height (m)

Tᵇᵍ(z) = max(Tᵗˢ, T₀ * Π(z))
qᵇᵍ(z) = max(0, β * q₀ * exp(-z / Hq))

qᵇᵍ (generic function with 1 method)

Adjust reference state to match actual profiles

compute_reference_state!(reference_state, Tᵇᵍ, qᵇᵍ, constants)

dynamics = AnelasticDynamics(reference_state)
coriolis = FPlane(f = 3e-4)

FPlane{Float32}(f=0.0003)

Surface fluxes

Following the paper's bulk formulas (Eqs. 2-4), with drag coefficient Cᴰ = 1.5 × 10⁻³ and gustiness v★ = 1 m/s. The surface wetness parameter β scales the moisture flux coefficient.

Cᴰ = Cᵀ = 1.5e-3
Uᵍ = 1

ρu_bcs = FieldBoundaryConditions(bottom = BulkDrag(coefficient = Cᴰ, gustiness = Uᵍ))
ρv_bcs = FieldBoundaryConditions(bottom = BulkDrag(coefficient = Cᴰ, gustiness = Uᵍ))

ρe_bcs = FieldBoundaryConditions(bottom = BulkSensibleHeatFlux(coefficient = Cᵀ,
                                                               gustiness = Uᵍ,
                                                               surface_temperature = T₀))

ρqᵗ_bcs = FieldBoundaryConditions(bottom = BulkVaporFlux(coefficient = β*Cᵀ,
                                                         gustiness = Uᵍ,
                                                         surface_temperature = T₀))

boundary_conditions = (; ρu=ρu_bcs, ρv=ρv_bcs, ρe=ρe_bcs, ρqᵗ=ρqᵗ_bcs)

Radiative forcing

The paper (Eq. 1) prescribes a piecewise radiative tendency: constant cooling at Ṫ = 1 K/day for T > Tᵗˢ (troposphere), and Newtonian relaxation toward Tᵗˢ with timescale τᵣ = 20 days for T ≤ Tᵗˢ (stratosphere). We apply this as an energy forcing on ρe, so that Breeze handles the conversion to ρθ tendency.

Ṫ  = 1 / day
τᵣ = 20days
ρᵣ = reference_state.density
parameters = (; Tᵗˢ, Ṫ, τᵣ, ρᵣ, cᵖᵈ)

@inline function ρe_forcing_func(i, j, k, grid, clock, model_fields, p)
    @inbounds T = model_fields.T[i, j, k]
    @inbounds ρ = p.ρᵣ[i, j, k]
    ∂T∂t = ifelse(T > p.Tᵗˢ, -p.Ṫ, (p.Tᵗˢ - T) / p.τᵣ)
    return ρ * p.cᵖᵈ * ∂T∂t
end

ρe_forcing = Forcing(ρe_forcing_func; discrete_form=true, parameters)

DiscreteForcing{@NamedTuple{Tᵗˢ::Int64, Ṫ::Float64, τᵣ::Float64, ρᵣ::Oceananigans.Fields.Field{Nothing, Nothing, Oceananigans.Grids.Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{OffsetArrays.OffsetVector{Float32, CUDA.CuArray{Float32, 1, CUDA.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDA.CuArray{Float32, 1, CUDA.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDA.CuArray{Float32, 1, CUDA.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDA.CuArray{Float32, 1, CUDA.DeviceMemory}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Architectures.GPU{CUDA.CUDAKernels.CUDABackend}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 3, CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Nothing, Nothing, Nothing, Nothing, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Value, Float32}, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Nothing, @NamedTuple{bottom_and_top::KernelAbstractions.Kernel{CUDA.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, typeof(Oceananigans.BoundaryConditions.gpu__fill_bottom_and_top_halo!)}, south_and_north::Nothing, west_and_east::Nothing}, @NamedTuple{bottom_and_top::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Value, Float32}, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition}, south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Nothing, Nothing}}}, Nothing, Nothing}, cᵖᵈ::Float32}}
├── func: ρe_forcing_func (generic function with 1 method)
└── parameters: (Tᵗˢ = 210, Ṫ = 1.1574074074074073e-5, τᵣ = 1.728e6, ρᵣ = 1×1×60 Field{Nothing, Nothing, Oceananigans.Grids.Center} reduced over dims = (1, 2) on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── grid: 72×72×60 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Value, top: ZeroFlux, immersed: Nothing
└── data: 1×1×70 OffsetArray(::CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}, 1:1, 1:1, -4:65) with eltype Float32 with indices 1:1×1:1×-4:65
    └── max=1.16002, min=0.0238631, mean=0.745775, cᵖᵈ = 1005.0f0)

Sponge layer

Rayleigh damping with a Gaussian profile centered at 26 km (width 2 km) prevents spurious wave reflections from the rigid lid.

sponge_mask = GaussianMask{:z}(center=26000, width=2000)
ρw_sponge = Relaxation(rate=1/30, mask=sponge_mask)

forcing = (; ρe=ρe_forcing, ρw=ρw_sponge)

Model

We use 9th-order WENO advection and warm-phase saturation adjustment microphysics.

momentum_advection = WENO(order=9)
scalar_advection = (ρθ = WENO(order=5),
                    ρqᵗ = WENO(order=5, bounds=(0, 1)))

microphysics = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())

model = AtmosphereModel(grid; dynamics, coriolis, momentum_advection, scalar_advection,
                        microphysics, forcing, boundary_conditions)

AtmosphereModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 72×72×60 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float32}
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: WENO{5, Float32, Float32}(order=9)
│   ├── ρθ: WENO{3, Float32, Float32}(order=5)
│   └── ρqᵗ: WENO{3, Float32, Float32}(order=5, bounds=(0.0f0, 1.0f0))
├── forcing: ρw=>ContinuousForcing, ρe=>DiscreteForcing
├── tracers: ()
├── coriolis: FPlane{Float32}(f=0.0003)
└── microphysics: SaturationAdjustment

Initial conditions

We initialize with an equilibrated temperature profile: a dry adiabat in the troposphere transitioning to an isothermal stratosphere at Tᵗˢ = 210 K. This approximates the paper's 100-day nonrotating RCE spinup. Small random perturbations in the lowest kilometer trigger convection.

Important: After compute_reference_state!, we must use set!(model, T=...) rather than set!(model, θ=...). The compute_reference_state! call recomputes the reference pressure, which changes the Exner function used to convert θ → T. Setting θ directly would produce incorrect temperatures in the stratosphere.

δT = 1//2  # K perturbation amplitude
zδ = 1000  # m perturbation depth
δq = 1e-4  # moisture perturbation amplitude (kg/kg)

Tᵢ(x, y, z) = Tᵇᵍ(z) + δT * (2rand() - 1) * (z < zδ)
qᵗᵢ(x, y, z) = max(0, qᵇᵍ(z) + δq * (2rand() - 1) * (z < zδ))

set!(model, T = Tᵢ, qᵗ = qᵗᵢ)

Simulation

We run for 4 days, which is sufficient for moist TC genesis and intensification.

simulation = Simulation(model; Δt=1, stop_time=4days)
conjure_time_step_wizard!(simulation, cfl=0.7)

Output and progress

u, v, w = model.velocities
θ = liquid_ice_potential_temperature(model)
s = @at (Center, Center, Center) sqrt(u^2 + v^2)
s₀ = Field(s, indices = (:, :, 1))

ρqᵗ = model.moisture_density
ρe = static_energy_density(model)
ℒˡ = Breeze.Thermodynamics.liquid_latent_heat(T₀, constants)
𝒬ᵀ = BoundaryConditionOperation(ρe, :bottom, model)
Jᵛ = BoundaryConditionOperation(ρqᵗ, :bottom, model)
𝒬 = Field(𝒬ᵀ + ℒˡ * Jᵛ)

function progress(sim)
    compute!(s₀)
    compute!(𝒬)
    umax = maximum(abs, u)
    vmax = maximum(abs, v)
    wmax = maximum(abs, w)
    s₀max = maximum(s₀)
    𝒬max = maximum(𝒬)
    θmin, θmax = extrema(θ)
    msg = @sprintf("(%d) t = %s, Δt = %s",
                   iteration(sim), prettytime(sim, false), prettytime(sim.Δt, false))
    msg *= @sprintf(", s₀ = %.1f m/s, max(𝒬) = %.1f W/m², max|U| ≈ (%d, %d, %d) m/s, θ ∈ [%d, %d] K",
                    s₀max, 𝒬max, umax, vmax, wmax, floor(θmin), ceil(θmax))
    @info msg
    return nothing
end

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

Horizontally-averaged profiles.

qᵗ = model.specific_moisture
ℋ = RelativeHumidity(model)

avg_outputs = (θ = Average(θ, dims=(1, 2)),
               qᵗ = Average(qᵗ, dims=(1, 2)),
               ℋ = Average(ℋ, dims=(1, 2)),
               w² = Average(w^2, dims=(1, 2)),
               wθ = Average(w * θ, dims=(1, 2)),
               wqᵗ = Average(w * qᵗ, dims=(1, 2)))

function save_parameters(file, model)
    file["parameters/β"] = β
    file["parameters/T₀"] = T₀
    file["parameters/Tᵗˢ"] = Tᵗˢ
    file["parameters/Ṫ"] = Ṫ
    file["parameters/f₀"] = 3e-4
    file["parameters/Cᴰ"] = Cᴰ
    file["parameters/Nx"] = Nx
    file["parameters/Nz"] = Nz
end

simulation.output_writers[:profiles] = JLD2Writer(model, avg_outputs;
                                                  filename = "tc_world_profiles.jld2",
                                                  schedule = TimeInterval(1day),
                                                  init = save_parameters,
                                                  overwrite_existing = true)

JLD2Writer scheduled on TimeInterval(1 day):
├── filepath: tc_world_profiles.jld2
├── 6 outputs: (θ, qᵗ, ℋ, w², wθ, wqᵗ)
├── array_type: Array{Float32}
├── including: [:grid, :thermodynamic_constants]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

Surface fields for tracking TC development.

surface_outputs = (; s, 𝒬)
simulation.output_writers[:surface] = JLD2Writer(model, surface_outputs;
                                                 filename = "tc_world_surface.jld2",
                                                 indices = (:, :, 1),
                                                 schedule = TimeInterval(30minutes),
                                                 overwrite_existing = true)

JLD2Writer scheduled on TimeInterval(30 minutes):
├── filepath: tc_world_surface.jld2
├── 2 outputs: (s, 𝒬)
├── array_type: Array{Float32}
├── including: [:grid, :thermodynamic_constants]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

Run

run!(simulation)

[ Info: Initializing simulation...
[ Info: (0) t = 0 s, Δt = 1.100 s, s₀ = 0.0 m/s, max(𝒬) = 30.4 W/m², max|U| ≈ (0, 0, 0) m/s, θ ∈ [299, 706] K
[ Info:     ... simulation initialization complete (2.006 minutes)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (7.320 seconds).
[ Info: (1000) t = 45.705 m, Δt = 3.768 s, s₀ = 33.3 m/s, max(𝒬) = 2352.5 W/m², max|U| ≈ (28, 31, 16) m/s, θ ∈ [297, 709] K
[ Info: (2000) t = 4.152 hrs, Δt = 15.725 s, s₀ = 7.0 m/s, max(𝒬) = 462.8 W/m², max|U| ≈ (25, 20, 6) m/s, θ ∈ [299, 663] K
[ Info: (3000) t = 8.281 hrs, Δt = 14.572 s, s₀ = 7.2 m/s, max(𝒬) = 497.7 W/m², max|U| ≈ (18, 22, 7) m/s, θ ∈ [299, 640] K
[ Info: (4000) t = 12.628 hrs, Δt = 16.535 s, s₀ = 7.1 m/s, max(𝒬) = 477.2 W/m², max|U| ≈ (19, 16, 8) m/s, θ ∈ [299, 624] K
[ Info: (5000) t = 16.980 hrs, Δt = 15.259 s, s₀ = 5.5 m/s, max(𝒬) = 358.5 W/m², max|U| ≈ (17, 21, 9) m/s, θ ∈ [300, 615] K
[ Info: (6000) t = 21.454 hrs, Δt = 16.182 s, s₀ = 6.7 m/s, max(𝒬) = 355.9 W/m², max|U| ≈ (17, 19, 9) m/s, θ ∈ [300, 609] K
[ Info: (7000) t = 1.076 d, Δt = 13.167 s, s₀ = 7.5 m/s, max(𝒬) = 352.7 W/m², max|U| ≈ (17, 15, 11) m/s, θ ∈ [300, 605] K
[ Info: (8000) t = 1.261 d, Δt = 17.103 s, s₀ = 8.2 m/s, max(𝒬) = 351.4 W/m², max|U| ≈ (14, 15, 12) m/s, θ ∈ [299, 603] K
[ Info: (9000) t = 1.444 d, Δt = 13.620 s, s₀ = 10.8 m/s, max(𝒬) = 383.0 W/m², max|U| ≈ (11, 13, 12) m/s, θ ∈ [299, 601] K
[ Info: (10000) t = 1.604 d, Δt = 15.413 s, s₀ = 11.6 m/s, max(𝒬) = 459.1 W/m², max|U| ≈ (11, 12, 13) m/s, θ ∈ [299, 600] K
[ Info: (11000) t = 1.743 d, Δt = 11.621 s, s₀ = 12.9 m/s, max(𝒬) = 379.4 W/m², max|U| ≈ (20, 18, 32) m/s, θ ∈ [299, 600] K
[ Info: (12000) t = 1.884 d, Δt = 11.332 s, s₀ = 17.5 m/s, max(𝒬) = 450.8 W/m², max|U| ≈ (15, 15, 17) m/s, θ ∈ [298, 600] K
[ Info: (13000) t = 2.022 d, Δt = 10.279 s, s₀ = 17.2 m/s, max(𝒬) = 573.2 W/m², max|U| ≈ (19, 20, 22) m/s, θ ∈ [299, 600] K
[ Info: (14000) t = 2.150 d, Δt = 11.150 s, s₀ = 21.5 m/s, max(𝒬) = 570.6 W/m², max|U| ≈ (21, 20, 19) m/s, θ ∈ [298, 600] K
[ Info: (15000) t = 2.269 d, Δt = 11.842 s, s₀ = 19.5 m/s, max(𝒬) = 532.1 W/m², max|U| ≈ (19, 18, 17) m/s, θ ∈ [297, 600] K
[ Info: (16000) t = 2.388 d, Δt = 13.342 s, s₀ = 24.7 m/s, max(𝒬) = 556.0 W/m², max|U| ≈ (24, 22, 27) m/s, θ ∈ [296, 600] K
[ Info: (17000) t = 2.503 d, Δt = 10.780 s, s₀ = 23.1 m/s, max(𝒬) = 568.3 W/m², max|U| ≈ (21, 20, 22) m/s, θ ∈ [297, 601] K
[ Info: (18000) t = 2.620 d, Δt = 8.243 s, s₀ = 24.0 m/s, max(𝒬) = 558.3 W/m², max|U| ≈ (23, 25, 16) m/s, θ ∈ [297, 602] K
[ Info: (19000) t = 2.728 d, Δt = 11.086 s, s₀ = 25.0 m/s, max(𝒬) = 624.1 W/m², max|U| ≈ (25, 25, 20) m/s, θ ∈ [295, 602] K
[ Info: (20000) t = 2.842 d, Δt = 9.156 s, s₀ = 25.1 m/s, max(𝒬) = 534.8 W/m², max|U| ≈ (24, 22, 15) m/s, θ ∈ [299, 602] K
[ Info: (21000) t = 2.953 d, Δt = 9.097 s, s₀ = 27.1 m/s, max(𝒬) = 550.1 W/m², max|U| ≈ (27, 26, 23) m/s, θ ∈ [298, 603] K
[ Info: (22000) t = 3.072 d, Δt = 9.129 s, s₀ = 27.7 m/s, max(𝒬) = 574.2 W/m², max|U| ≈ (25, 25, 13) m/s, θ ∈ [298, 603] K
[ Info: (23000) t = 3.179 d, Δt = 11.388 s, s₀ = 25.1 m/s, max(𝒬) = 719.1 W/m², max|U| ≈ (26, 24, 13) m/s, θ ∈ [298, 603] K
[ Info: (24000) t = 3.287 d, Δt = 7.991 s, s₀ = 27.1 m/s, max(𝒬) = 613.6 W/m², max|U| ≈ (29, 25, 24) m/s, θ ∈ [295, 603] K
[ Info: (25000) t = 3.392 d, Δt = 9.640 s, s₀ = 29.3 m/s, max(𝒬) = 620.1 W/m², max|U| ≈ (27, 26, 24) m/s, θ ∈ [295, 602] K
[ Info: (26000) t = 3.498 d, Δt = 9.573 s, s₀ = 29.8 m/s, max(𝒬) = 529.0 W/m², max|U| ≈ (28, 28, 26) m/s, θ ∈ [298, 602] K
[ Info: (27000) t = 3.610 d, Δt = 10.523 s, s₀ = 33.9 m/s, max(𝒬) = 633.8 W/m², max|U| ≈ (31, 31, 15) m/s, θ ∈ [298, 602] K
[ Info: (28000) t = 3.720 d, Δt = 8.338 s, s₀ = 30.4 m/s, max(𝒬) = 598.8 W/m², max|U| ≈ (27, 29, 24) m/s, θ ∈ [297, 602] K
[ Info: (29000) t = 3.823 d, Δt = 7.437 s, s₀ = 31.0 m/s, max(𝒬) = 665.2 W/m², max|U| ≈ (29, 30, 14) m/s, θ ∈ [298, 602] K
[ Info: (30000) t = 3.925 d, Δt = 12.161 s, s₀ = 30.2 m/s, max(𝒬) = 604.3 W/m², max|U| ≈ (30, 30, 8) m/s, θ ∈ [298, 602] K
[ Info: Simulation is stopping after running for 21.357 minutes.
[ Info: Simulation time 4 days equals or exceeds stop time 4 days.

Results: mean profile evolution

Evolution of horizontally-averaged potential temperature, vertical velocity variance, and the vertical potential temperature flux.

θt = FieldTimeSeries("tc_world_profiles.jld2", "θ")
qᵗt = FieldTimeSeries("tc_world_profiles.jld2", "qᵗ")
ℋt = FieldTimeSeries("tc_world_profiles.jld2", "ℋ")
w²t = FieldTimeSeries("tc_world_profiles.jld2", "w²")
wθt = FieldTimeSeries("tc_world_profiles.jld2", "wθ")
wqᵗt = FieldTimeSeries("tc_world_profiles.jld2", "wqᵗ")

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

fig = Figure(size=(900, 400), fontsize=10)

axθ = Axis(fig[1, 1], xlabel="θ (K)", ylabel="z (m)")
axqᵗ = Axis(fig[1, 2], xlabel="qᵗ (kg/kg)")
axℋ = Axis(fig[1, 3], xlabel="ℋ")
axw² = Axis(fig[1, 4], xlabel="w² (m²/s²)")
axwθ = Axis(fig[1, 5], xlabel="wθ (m²/s² K)")
axwqᵗ = Axis(fig[1, 6], xlabel="wqᵗ (m²/s² kg/kg)", ylabel="z (m)", yaxisposition=:right)

default_colours = Makie.wong_colors()
colors = [default_colours[mod1(n, length(default_colours))] for n in 1:Nt]
linewidth = 3
alpha = 0.6

for n in 1:Nt
    label = n == 1 ? "initial" : "t = $(prettytime(times[n]))"
    lines!(axθ, θt[n], color=colors[n]; label, linewidth, alpha)
    lines!(axqᵗ, qᵗt[n], color=colors[n]; linewidth, alpha)
    lines!(axℋ, ℋt[n], color=colors[n]; linewidth, alpha)
    lines!(axw², w²t[n], color=colors[n]; linewidth, alpha)
    lines!(axwθ, wθt[n], color=colors[n]; linewidth, alpha)
    lines!(axwqᵗ, wqᵗt[n], color=colors[n]; linewidth, alpha)
end

for ax in (axqᵗ, axℋ, axw², axwθ)
    hideydecorations!(ax, grid=false)
    hidespines!(ax, :t, :r, :l)
end

hidespines!(axθ, :t, :r)
hidespines!(axwqᵗ, :t, :l)
xlims!(axℋ, -0.1, 1.1)

Legend(fig[2, :], axθ, labelsize=12, orientation=:horizontal)

fig[0, :] = Label(fig, "TC World (β = $β): mean profile evolution",
                  fontsize=16, tellwidth=false)

fig

Surface wind speed snapshots

Snapshots of the surface wind speed field at early, middle, and late times show the evolution of convective organization and TC formation.

s_ts = FieldTimeSeries("tc_world_surface.jld2", "s")
𝒬_ts = FieldTimeSeries("tc_world_surface.jld2", "𝒬")

times = s_ts.times
Nt = length(times)

smax = maximum(s_ts)
slim = smax / 2
𝒬lim = maximum(𝒬_ts) / 8

fig = Figure(size=(1200, 800), fontsize=12)

s_heatmaps = []
𝒬_heatmaps = []
indices = ceil.(Int, [Nt / 3, 2Nt / 3, Nt])

for (i, idx) in enumerate(indices)
    xlabel = i == 1 ? "x (m)" : ""
    ylabel = i == 1 ? "y (m)" : ""
    title = "t = $(prettytime(times[idx]))"
    axs = Axis(fig[1, i]; aspect = 1, xlabel, ylabel, title)
    ax𝒬 = Axis(fig[2, i]; aspect = 1, xlabel, ylabel, title)
    s_hm = heatmap!(axs, s_ts[idx]; colormap=:speed, colorrange=(0, slim))
    push!(s_heatmaps, s_hm)
    𝒬_hm = heatmap!(ax𝒬, 𝒬_ts[idx]; colormap=:magma, colorrange=(0, 𝒬lim))
    push!(𝒬_heatmaps, 𝒬_hm)
end

Colorbar(fig[1, length(indices) + 1], s_heatmaps[end]; label="Surface wind speed (m/s)")
Colorbar(fig[2, length(indices) + 1], 𝒬_heatmaps[end]; label="Surface moisture flux (W/m²)")

fig[0, :] = Label(fig, "TC World (β = $β): surface wind and heat flux",
                  fontsize=16, tellwidth=false)

fig

Animation of surface wind speed

fig = Figure(size=(600, 550), fontsize=14)
ax = Axis(fig[1, 1]; xlabel="x (m)", ylabel="y (m)", aspect=1)

n = Observable(1)
title = @lift "TC World (β = $β) — t = $(prettytime(times[$n]))"
sn = @lift s_ts[$n]

hm = heatmap!(ax, sn; colormap=:speed, colorrange=(0, slim))
Colorbar(fig[1, 2], hm; label="Surface wind speed (m/s)")
fig[0, :] = Label(fig, title, fontsize=16, tellwidth=false)

CairoMakie.record(fig, "tc_world.mp4", 1:Nt, framerate=16) do nn
    n[] = nn
end

Discussion

This example demonstrates spontaneous tropical cyclone genesis in a rotating radiative-convective equilibrium setup, following Cronin and Chavas (2019). The surface wetness parameter β controls moisture availability: β = 1 (default) produces robust moist TC genesis, while β = 0 yields dry TCs.

The radiative forcing is a piecewise temperature tendency: constant cooling at 1 K/day in the troposphere (T > Tᵗˢ) and Newtonian relaxation toward Tᵗˢ with timescale τᵣ = 20 days in the stratosphere. Surface fluxes use bulk formulas with drag coefficient Cᴰ = 1.5 × 10⁻³ and gustiness 1 m/s. The f-plane Coriolis parameter is f₀ = 3 × 10⁻⁴ s⁻¹.

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.4
Commit 01a2eadb047 (2026-01-06 16:56 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.4
  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.3.3 `.`
  [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.11
  [9e8cae18] Oceananigans v0.104.5
  [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.