Tropical Cyclone World (Cronin and 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)
if CUDA.functional()
    CUDA.seed!(2019)
end

Oceananigans.defaults.FloatType = 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 domain 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 = 1152kilometers
paper_Nx = 576
Lx = Ly = paper_Lx / 4
Nx = Ny = paper_Nx / 8 |> Int
H = 28kilometers

Δ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_0 = 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 ($T_{ref} ≈ 26$ K vs $T_{actual} = 210$ K at 28 km), producing catastrophic buoyancy forces.

T₀ = 300     # K
p₀ = 101325  # Pa
constants = ThermodynamicConstants()

reference_state = ReferenceState(grid, constants;
                                 surface_pressure = p₀,
                                 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
Π₀ = (p₀ / 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^D = 1.5 × 10^{-3}$ and gustiness $U^g = 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^{ts}$ (troposphere), and Newtonian relaxation toward $T^{ts}$ with timescale $τ_r = 20$ days for $T ≤ T^{ts}$ (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{south_and_north::Nothing, west_and_east::Nothing, bottom_and_top::KernelAbstractions.Kernel{CUDA.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.gpu__fill_bottom_and_top_halo!)}}, @NamedTuple{south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Nothing, Nothing}, bottom_and_top::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Value, Float32}, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition}}}, 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=26kilometers, width=2kilometers)
ρw_sponge = Relaxation(rate=1/30, mask=sponge_mask)

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

Model

We use WENO schemes for 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, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│   ├── ρθ: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   └── ρqᵉ: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(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^{ts} = 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 $θ \to T$. Setting $θ$ directly would produce incorrect temperatures in the stratosphere.

δT = 1//2  # perturbation amplitude (K)
zδ = 1000  # perturbation depth (m)
δ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)
Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)

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ᵛ = specific_humidity(model)
ℋ = 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₀"] = coriolis.f
    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: [: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: [: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(𝒬) = 28.4 W/m², max|U| ≈ (0, 0, 0) m/s, θ ∈ [299, 706] K
[ Info:     ... simulation initialization complete (1.910 minutes)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (10.528 seconds).
[ Info: (1000) t = 44.430 m, Δt = 3.019 s, s₀ = 38.4 m/s, max(𝒬) = 2400.5 W/m², max|U| ≈ (38, 31, 16) m/s, θ ∈ [296, 708] K
[ Info: (2000) t = 4.041 hrs, Δt = 15.912 s, s₀ = 6.6 m/s, max(𝒬) = 468.6 W/m², max|U| ≈ (21, 21, 6) m/s, θ ∈ [298, 663] K
[ Info: (3000) t = 8.537 hrs, Δt = 16.692 s, s₀ = 5.6 m/s, max(𝒬) = 344.8 W/m², max|U| ≈ (21, 20, 5) m/s, θ ∈ [299, 638] K
[ Info: (4000) t = 13.216 hrs, Δt = 18.012 s, s₀ = 5.6 m/s, max(𝒬) = 399.1 W/m², max|U| ≈ (17, 18, 6) m/s, θ ∈ [298, 621] K
[ Info: (5000) t = 18.200 hrs, Δt = 19.768 s, s₀ = 5.6 m/s, max(𝒬) = 336.8 W/m², max|U| ≈ (17, 18, 10) m/s, θ ∈ [299, 612] K
[ Info: (6000) t = 23.371 hrs, Δt = 19.665 s, s₀ = 6.6 m/s, max(𝒬) = 286.4 W/m², max|U| ≈ (19, 18, 9) m/s, θ ∈ [299, 606] K
[ Info: (7000) t = 1.183 d, Δt = 16.766 s, s₀ = 10.2 m/s, max(𝒬) = 441.6 W/m², max|U| ≈ (13, 14, 9) m/s, θ ∈ [299, 602] K
[ Info: (8000) t = 1.375 d, Δt = 18.102 s, s₀ = 9.0 m/s, max(𝒬) = 330.1 W/m², max|U| ≈ (11, 12, 13) m/s, θ ∈ [298, 601] K
[ Info: (9000) t = 1.546 d, Δt = 13.580 s, s₀ = 11.3 m/s, max(𝒬) = 359.7 W/m², max|U| ≈ (15, 15, 19) m/s, θ ∈ [298, 600] K
[ Info: (10000) t = 1.704 d, Δt = 11.863 s, s₀ = 12.8 m/s, max(𝒬) = 524.2 W/m², max|U| ≈ (14, 14, 17) m/s, θ ∈ [298, 600] K
[ Info: (11000) t = 1.846 d, Δt = 14.836 s, s₀ = 14.7 m/s, max(𝒬) = 338.7 W/m², max|U| ≈ (14, 14, 13) m/s, θ ∈ [298, 600] K
[ Info: (12000) t = 1.979 d, Δt = 10.762 s, s₀ = 14.7 m/s, max(𝒬) = 421.0 W/m², max|U| ≈ (16, 15, 18) m/s, θ ∈ [297, 602] K
[ Info: (13000) t = 2.116 d, Δt = 13.561 s, s₀ = 15.8 m/s, max(𝒬) = 437.7 W/m², max|U| ≈ (15, 16, 25) m/s, θ ∈ [298, 603] K
[ Info: (14000) t = 2.241 d, Δt = 13.662 s, s₀ = 18.1 m/s, max(𝒬) = 451.4 W/m², max|U| ≈ (17, 16, 17) m/s, θ ∈ [298, 605] K
[ Info: (15000) t = 2.369 d, Δt = 12.160 s, s₀ = 16.7 m/s, max(𝒬) = 470.8 W/m², max|U| ≈ (17, 18, 26) m/s, θ ∈ [298, 608] K
[ Info: (16000) t = 2.498 d, Δt = 8.603 s, s₀ = 18.3 m/s, max(𝒬) = 479.7 W/m², max|U| ≈ (19, 19, 20) m/s, θ ∈ [297, 610] K
[ Info: (17000) t = 2.614 d, Δt = 8.095 s, s₀ = 19.4 m/s, max(𝒬) = 675.0 W/m², max|U| ≈ (22, 22, 29) m/s, θ ∈ [297, 612] K
[ Info: (18000) t = 2.725 d, Δt = 9.119 s, s₀ = 24.5 m/s, max(𝒬) = 624.5 W/m², max|U| ≈ (28, 24, 19) m/s, θ ∈ [296, 614] K
[ Info: (19000) t = 2.835 d, Δt = 9.185 s, s₀ = 22.1 m/s, max(𝒬) = 621.3 W/m², max|U| ≈ (20, 20, 22) m/s, θ ∈ [295, 614] K
[ Info: (20000) t = 2.954 d, Δt = 9.216 s, s₀ = 25.1 m/s, max(𝒬) = 538.5 W/m², max|U| ≈ (24, 24, 21) m/s, θ ∈ [297, 617] K
[ Info: (21000) t = 3.067 d, Δt = 7.474 s, s₀ = 24.6 m/s, max(𝒬) = 523.7 W/m², max|U| ≈ (23, 22, 17) m/s, θ ∈ [298, 619] K
[ Info: (22000) t = 3.169 d, Δt = 8.848 s, s₀ = 22.5 m/s, max(𝒬) = 606.3 W/m², max|U| ≈ (21, 20, 16) m/s, θ ∈ [298, 619] K
[ Info: (23000) t = 3.283 d, Δt = 13.822 s, s₀ = 24.8 m/s, max(𝒬) = 564.7 W/m², max|U| ≈ (26, 24, 14) m/s, θ ∈ [298, 620] K
[ Info: (24000) t = 3.404 d, Δt = 9.707 s, s₀ = 23.8 m/s, max(𝒬) = 586.0 W/m², max|U| ≈ (23, 22, 16) m/s, θ ∈ [298, 622] K
[ Info: (25000) t = 3.508 d, Δt = 8.893 s, s₀ = 22.9 m/s, max(𝒬) = 599.4 W/m², max|U| ≈ (22, 22, 12) m/s, θ ∈ [298, 623] K
[ Info: (26000) t = 3.619 d, Δt = 10.553 s, s₀ = 21.8 m/s, max(𝒬) = 548.9 W/m², max|U| ≈ (21, 20, 37) m/s, θ ∈ [298, 623] K
[ Info: (27000) t = 3.732 d, Δt = 10.612 s, s₀ = 25.1 m/s, max(𝒬) = 653.2 W/m², max|U| ≈ (23, 24, 15) m/s, θ ∈ [297, 622] K
[ Info: (28000) t = 3.833 d, Δt = 9.344 s, s₀ = 23.7 m/s, max(𝒬) = 724.0 W/m², max|U| ≈ (22, 21, 23) m/s, θ ∈ [297, 626] K
[ Info: (29000) t = 3.940 d, Δt = 9.251 s, s₀ = 26.9 m/s, max(𝒬) = 804.0 W/m², max|U| ≈ (24, 25, 20) m/s, θ ∈ [297, 623] K
[ Info: Simulation is stopping after running for 32.167 minutes.
[ Info: Simulation time 4 days equals or exceeds stop time 4 days.

Results: mean profile evolution

Evolution of horizontally-averaged potential temperature, specific and relative humidity, vertical velocity variance, vertical potential temperature flux, and vertical specific humidity 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ᵛ (10⁻⁵ 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
    color = colors[n]
    label = n == 1 ? "initial" : "t = $(prettytime(times[n]))"
    lines!(axθ, θt[n]; color, linewidth, alpha, label)
    lines!(axqᵛ, qᵛt[n]; color, linewidth, alpha)
    lines!(axℋ, ℋt[n]; color, linewidth, alpha)
    lines!(axw², w²t[n]; color, linewidth, alpha)
    lines!(axwθ, wθt[n]; color, linewidth, alpha)
    lines!(axwqᵛ, 1e5 * wqᵛt[n]; color, 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.

st = FieldTimeSeries("tc_world_surface.jld2", "s")
𝒬t = FieldTimeSeries("tc_world_surface.jld2", "𝒬")

times = st.times
Nt = length(times)

smax = maximum(st)
slim = smax / 2
𝒬lim = maximum(𝒬t) / 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, st[idx]; colormap=:speed, colorrange=(0, slim))
    push!(s_heatmaps, s_hm)
    𝒬_hm = heatmap!(ax𝒬, 𝒬t[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 st[$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^{ts}$) and Newtonian relaxation toward $T^{ts}$ with timescale $τ_r = 20$ days in the stratosphere. Surface fluxes use bulk formulas with drag coefficient $C^D = 1.5 × 10^{-3}$ and gustiness 1 m/s. The $f$-plane Coriolis parameter is $f_0 = 3 × 10^{-4}$ s⁻¹.

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.8 `.`
⌃ [052768ef] CUDA v5.11.2
  [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.