Atmospheric convection over a slab ocean vs a hydrostatic ocean
This example demonstrates coupling a Breeze atmospheric large eddy simulation (LES) with two different ocean models using NumericalEarth's EarthSystemModel framework:
- Slab ocean (10m depth) — a well-mixed layer that responds uniformly to surface fluxes
- Full hydrostatic ocean (50m depth) — with CATKE turbulent mixing and stratification
The atmosphere drives convective turbulence over a warm ocean surface. The coupling framework computes turbulent surface fluxes (sensible heat, latent heat, and momentum) using Monin–Obukhov similarity theory. These fluxes cool the ocean and heat the atmosphere, creating a two-way feedback loop.
By comparing the two models, we can see how ocean vertical mixing and stratification affect the SST response to atmospheric forcing.
using NumericalEarth
using Breeze
using Oceananigans
using Oceananigans.Units
using Printf
using Statistics: mean[ Info: Oceananigans will use 8 threads
Grid setup
We use a 2D domain in the x-z plane: 20 km wide and 10 km tall with 128 × 128 grid points. The Periodic x-topology allows convective cells to wrap around, and Flat y-topology makes this a 2D simulation.
Nxᵃᵗ = 64 # Atmosphere horizontal resolution (shared with ocean)
Nzᵃᵗ = 64 # Atmosphere vertical resolution
Nzᵒᶜ = 20 # Hydrostatic ocean vertical resolution
grid = RectilinearGrid(size = (Nxᵃᵗ, Nzᵃᵗ), halo = (5, 5),
x = (-10kilometers, 10kilometers),
z = (0, 10kilometers),
topology = (Periodic, Flat, Bounded))64×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── Periodic x ∈ [-10000.0, 10000.0) regularly spaced with Δx=312.5
├── Flat y
└── Bounded z ∈ [0.0, 10000.0] regularly spaced with Δz=156.25Two independent atmospheres
Each coupled model needs its own Breeze atmosphere instance because the EarthSystemModel writes boundary conditions into the atmosphere. Both are initialized identically.
Tᵒᶜ = 290 # K
θᵃᵗ = 250 # K
U₀ = 10 # m/s
coriolis = FPlane(latitude=33)
slab_ocean_atmos = atmosphere_simulation(grid; potential_temperature=θᵃᵗ, coriolis)
full_ocean_atmos = atmosphere_simulation(grid; potential_temperature=θᵃᵗ, coriolis)AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=250.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme:
│ ├── momentum: WENO{5, Float64, Float32}(order=9)
│ ├── ρθ: WENO{3, Float64, Float32}(order=5)
│ └── ρqᵗ: WENO{3, Float64, Float32}(order=5)
├── forcing: @NamedTuple{ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵗ::Returns{Float64}, ρe::Returns{Float64}}
├── tracers: ()
├── coriolis: FPlane{Float64}(f=7.94314e-5)
└── microphysics: SaturationAdjustmentAtmospheric initial conditions
We initialize both atmospheres with the reference potential temperature profile plus small random perturbations below 500 m. These perturbations seed convective instability, which develops into turbulent convection driven by surface heat fluxes. A background zonal wind U₀ provides a nonzero wind speed for the similarity theory flux computation.
reference_state = slab_ocean_atmos.dynamics.reference_state
θᵢ(x, z) = reference_state.potential_temperature + 0.1 * randn() * (z < 500)
set!(slab_ocean_atmos, θ=θᵢ, u=U₀)
set!(full_ocean_atmos, θ=θᵢ, u=U₀)Slab ocean (10m depth)
The slab ocean represents a well-mixed ocean layer of fixed depth. Its temperature is in Kelvin (the coupling framework handles the conversion).
sst_grid = RectilinearGrid(grid.architecture,
size = grid.Nx,
halo = grid.Hx,
x = (-10kilometers, 10kilometers),
topology = (Periodic, Flat, Flat))
slab_ocean = SlabOcean(sst_grid, depth=10)
set!(slab_ocean, T=Tᵒᶜ)Full hydrostatic ocean (50m depth with CATKE mixing)
The full ocean uses a HydrostaticFreeSurfaceModel with the default TEOS-10 equation of state and CATKE vertical mixing parameterization. The grid has 20 vertical levels (2.5m vertical resolution). We disable advection since this is primarily a 1D vertical mixing problem.
Since TEOS-10 expects temperature in degrees Celsius, the ocean temperature is initialized accordingly. The coupling framework automatically converts from Celsius to Kelvin for the flux computation.
ocean_grid = RectilinearGrid(grid.architecture,
size = (grid.Nx, Nzᵒᶜ),
halo = (grid.Hx, 5),
x = (-10kilometers, 10kilometers),
z = (-50, 0),
topology = (Periodic, Flat, Bounded))
ocean = ocean_simulation(ocean_grid; coriolis,
closure = CATKEVerticalDiffusivity(), # note ocean_simulation default does not work.
momentum_advection = nothing,
tracer_advection = nothing,
Δt = 2,
warn = false)
celsius_to_kelvin = 273.15
T₀ = Tᵒᶜ - celsius_to_kelvin # surface temperature in °C
Tᵢ(x, z) = T₀ + (z + 10) / 50 * (z < -10) # linear stratification below 10m
set!(ocean.model, T=Tᵢ, S=35)Coupled models
We disable gustiness in the similarity theory flux computation so the surface wind speed is determined entirely by the resolved velocity field.
atmosphere_ocean_fluxes = SimilarityTheoryFluxes(gustiness_parameter=0, minimum_gustiness=0)
slab_interfaces = ComponentInterfaces(slab_ocean_atmos, slab_ocean; atmosphere_ocean_fluxes)
full_interfaces = ComponentInterfaces(full_ocean_atmos, ocean; atmosphere_ocean_fluxes)
slab_model = AtmosphereOceanModel(slab_ocean_atmos, slab_ocean; interfaces = slab_interfaces)
full_model = AtmosphereOceanModel(full_ocean_atmos, ocean; interfaces = full_interfaces)
Δt = 5seconds
stop_time = 4hours
slab_sim = Simulation(slab_model; Δt, stop_time)
full_sim = Simulation(full_model; Δt, stop_time)Simulation of EarthSystemModel{CPU}(time = 0 seconds, iteration = 0)
├── Next time step: 5 seconds
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 4 hours
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│ ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│ ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│ └── nan_checker => Callback of NaNChecker for T_ocean on IterationInterval(100)
└── output_writers: OrderedDict with no entriesProgress callbacks
function slab_progress(sim)
atmos = sim.model.atmosphere
u, v, w = atmos.velocities
umax = maximum(abs, u)
wmax = maximum(abs, w)
sst = sim.model.ocean.temperature
sst_min = minimum(sst)
sst_max = maximum(sst)
msg = @sprintf("[Slab] Iter: %d, t = %s, max|u|: %.2e, max|w|: %.2e, SST: (%.3f, %.3f)",
iteration(sim), prettytime(sim), umax, wmax, sst_min, sst_max)
@info msg
return nothing
end
function full_progress(sim)
atmos = sim.model.atmosphere
u, v, w = atmos.velocities
umax = maximum(abs, u)
wmax = maximum(abs, w)
T = sim.model.ocean.model.tracers.T
Nz = size(sim.model.ocean.model.grid, 3)
sst_surface = view(interior(T), :, 1, Nz)
sst_min = minimum(sst_surface)
sst_max = maximum(sst_surface)
msg = @sprintf("[Full] Iter: %d, t = %s, max|u|: %.2e, max|w|: %.2e, SST: (%.3f, %.3f)",
iteration(sim), prettytime(sim), umax, wmax, sst_min, sst_max)
@info msg
return nothing
end
add_callback!(slab_sim, slab_progress, IterationInterval(400))
add_callback!(full_sim, full_progress, IterationInterval(400))Output writers
- Slab simulation: atmospheric θ and u (for heatmaps and profiles), plus slab SST.
- Full simulation: atmospheric θ, u, cloud water, and w, plus full-depth ocean temperature T.
u_s, v_s, w_s = slab_ocean_atmos.velocities
s_s = @at (Center, Center, Center) √(u_s^2 + v_s^2 + w_s^2)
θ_s = liquid_ice_potential_temperature(slab_ocean_atmos)
slab_sim.output_writers[:atmos] = JLD2Writer(slab_model, (; θ=θ_s, u=u_s, s=s_s),
filename = "slab_ocean_atmos",
schedule = TimeInterval(1minute),
overwrite_existing = true)
slab_sim.output_writers[:sst] = JLD2Writer(slab_model, (; SST=slab_ocean.temperature),
filename = "sst_slab",
schedule = TimeInterval(1minute),
overwrite_existing = true)
u_f, v_f, w_f = full_ocean_atmos.velocities
s_f = @at (Center, Center, Center) √(u_f^2 + v_f^2 + w_f^2)
θ_f = liquid_ice_potential_temperature(full_ocean_atmos)
qˡ_f = full_ocean_atmos.microphysical_fields.qˡ
full_sim.output_writers[:atmos] = JLD2Writer(full_model, (; θ=θ_f, u=u_f, w=w_f, s=s_f, qˡ=qˡ_f),
filename = "full_ocean_atmos",
schedule = TimeInterval(1minute),
overwrite_existing = true)
full_sim.output_writers[:ocean] = JLD2Writer(full_model, (; T=ocean.model.tracers.T),
filename = "ocean_full",
schedule = TimeInterval(1minute),
overwrite_existing = true)JLD2Writer scheduled on TimeInterval(1 minute):
├── filepath: ocean_full.jld2
├── 1 outputs: T
├── array_type: Array{Float32}
├── including: ()
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)Run both simulations sequentially
@info "Running slab ocean coupled simulation..."
run!(slab_sim)
@info "Slab ocean simulation complete."
@info "Running full ocean coupled simulation..."
run!(full_sim)
@info "Full ocean simulation complete."[ Info: Running slab ocean coupled simulation...
[ Info: Initializing simulation...
[ Info: [Slab] Iter: 0, t = 0 seconds, max|u|: 1.00e+01, max|w|: 0.00e+00, SST: (290.000, 290.000)
[ Info: ... simulation initialization complete (24.100 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (12.553 seconds).
[ Info: [Slab] Iter: 400, t = 33.333 minutes, max|u|: 1.77e+01, max|w|: 1.16e+01, SST: (289.933, 289.936)
[ Info: [Slab] Iter: 800, t = 1.111 hours, max|u|: 2.34e+01, max|w|: 1.90e+01, SST: (289.868, 289.883)
[ Info: [Slab] Iter: 1200, t = 1.667 hours, max|u|: 3.24e+01, max|w|: 1.20e+01, SST: (289.801, 289.821)
[ Info: [Slab] Iter: 1600, t = 2.222 hours, max|u|: 3.49e+01, max|w|: 3.14e+01, SST: (289.729, 289.752)
[ Info: [Slab] Iter: 2000, t = 2.778 hours, max|u|: 3.71e+01, max|w|: 2.86e+01, SST: (289.640, 289.666)
[ Info: [Slab] Iter: 2400, t = 3.333 hours, max|u|: 3.40e+01, max|w|: 3.39e+01, SST: (289.539, 289.574)
[ Info: [Slab] Iter: 2800, t = 3.889 hours, max|u|: 3.86e+01, max|w|: 3.69e+01, SST: (289.430, 289.483)
[ Info: Simulation is stopping after running for 1.611 minutes.
[ Info: Simulation time 4 hours equals or exceeds stop time 4 hours.
[ Info: Slab ocean simulation complete.
[ Info: Running full ocean coupled simulation...
[ Info: Initializing simulation...
[ Info: [Full] Iter: 0, t = 0 seconds, max|u|: 1.00e+01, max|w|: 0.00e+00, SST: (16.850, 16.850)
[ Info: ... simulation initialization complete (6.408 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (11.367 seconds).
[ Info: [Full] Iter: 400, t = 33.333 minutes, max|u|: 1.77e+01, max|w|: 9.81e+00, SST: (16.782, 16.784)
[ Info: [Full] Iter: 800, t = 1.111 hours, max|u|: 2.50e+01, max|w|: 2.09e+01, SST: (16.734, 16.740)
[ Info: [Full] Iter: 1200, t = 1.667 hours, max|u|: 2.78e+01, max|w|: 2.35e+01, SST: (16.684, 16.696)
[ Info: [Full] Iter: 1600, t = 2.222 hours, max|u|: 2.79e+01, max|w|: 3.53e+01, SST: (16.640, 16.653)
[ Info: [Full] Iter: 2000, t = 2.778 hours, max|u|: 2.81e+01, max|w|: 2.84e+01, SST: (16.593, 16.618)
[ Info: [Full] Iter: 2400, t = 3.333 hours, max|u|: 3.85e+01, max|w|: 3.41e+01, SST: (16.543, 16.576)
[ Info: [Full] Iter: 2800, t = 3.889 hours, max|u|: 4.85e+01, max|w|: 3.39e+01, SST: (16.489, 16.537)
[ Info: Simulation is stopping after running for 1.974 minutes.
[ Info: Simulation time 4 hours equals or exceeds stop time 4 hours.
[ Info: Full ocean simulation complete.
Animation
The animation has three columns:
- Left: atmosphere over slab ocean (θ, u) and SST comparison line plot
- Middle: atmosphere over full ocean (cloud water, w) and ocean temperature cross-section
- Right (narrow): horizontal-mean vertical profiles from both simulations
The SST comparison converts the full ocean surface temperature from °C to K so both curves are on the same Kelvin scale as the slab ocean.
using CairoMakie
θ_slab_ts = FieldTimeSeries("slab_ocean_atmos.jld2", "θ"; grid)
u_slab_ts = FieldTimeSeries("slab_ocean_atmos.jld2", "u"; grid)
s_slab_ts = FieldTimeSeries("slab_ocean_atmos.jld2", "s"; grid)
sst_slab_ts = FieldTimeSeries("sst_slab.jld2", "SST"; grid=sst_grid)
θ_full_ts = FieldTimeSeries("full_ocean_atmos.jld2", "θ"; grid)
u_full_ts = FieldTimeSeries("full_ocean_atmos.jld2", "u"; grid)
s_full_ts = FieldTimeSeries("full_ocean_atmos.jld2", "s"; grid)
qˡ_full_ts = FieldTimeSeries("full_ocean_atmos.jld2", "qˡ"; grid)
T_ocean_ts = FieldTimeSeries("ocean_full.jld2", "T"; grid=ocean_grid)
times = θ_slab_ts.times
Nt = length(times)
Nzᵒᶜ = size(ocean_grid, 3)20Convert from °C to K
for n in 1:Nt
T_ocean_ts[n].data .+= celsius_to_kelvin
endCoordinate arrays for manual line plots.
x_ocean = xnodes(ocean_grid, Center())-9843.75:312.5:9843.75Figure layout
fig = Figure(size = (1600, 900), fontsize = 12)
ax_θ = Axis(fig[1, 1], title="θₗᵢ (K) — atmos (slab ocean)", ylabel="z (m)")
ax_ss = Axis(fig[2, 1], title="√(u² + w²) (m/s) — atmos (slab ocean)", ylabel="z (m)")
ax_sst = Axis(fig[3, 1], title="SST (K)", xlabel="x (m)", ylabel="SST (K)")
ax_qˡ = Axis(fig[1, 2], title="Cloud water (kg/kg) — atmos (full ocean)", ylabel="z (m)")
ax_sf = Axis(fig[2, 2], title="√(u² + w²) — atmos (full ocean)", ylabel="z (m)")
ax_oT = Axis(fig[3, 2], title="Ocean T (Κ)", xlabel="x (m)", ylabel="z (m)")
ax_θp = Axis(fig[1, 3], title="⟨θ⟩(z)", xlabel="θ (K)", ylabel="z (m)")
ax_up = Axis(fig[2, 3], title="⟨u⟩(z)", xlabel="u (m/s)", ylabel="z (m)")
ax_Tp = Axis(fig[3, 3], title="⟨T⟩(z)", xlabel="T (Κ)", ylabel="z (m)")
colsize!(fig.layout, 3, Relative(0.15))
for ax in (ax_θ, ax_ss, ax_qˡ, ax_sf)
hidexdecorations!(ax, ticks=false)
endObservables
n = Observable(1)Observable(1)
Left column
θn = @lift θ_slab_ts[$n]
un = @lift u_slab_ts[$n]
ssn = @lift s_slab_ts[$n]
sstn_slab = @lift sst_slab_ts[$n]
ocean_sst = @lift interior(T_ocean_ts[$n], :, 1, Nzᵒᶜ)Observable([290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697, 290.0000003814697])
Middle column
qˡn = @lift qˡ_full_ts[$n]
sfn = @lift s_full_ts[$n]
oTn = @lift T_ocean_ts[$n]Observable(64×1×20 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
├── grid: 64×1×20 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 5×0×5 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: Flux, immersed: Nothing
└── data: 74×1×30 OffsetArray(view(::Array{Float64, 4}, :, :, :, 1), -4:69, 1:1, -4:25) with eltype Float64 with indices -4:69×1:1×-4:25
└── max=290.0, min=289.225, mean=289.68)
Right column — horizontal-mean profiles
θ_avg_slab = @lift Field(Average(θ_slab_ts[$n], dims=1))
θ_avg_full = @lift Field(Average(θ_full_ts[$n], dims=1))
u_avg_slab = @lift Field(Average(u_slab_ts[$n], dims=1))
u_avg_full = @lift Field(Average(u_full_ts[$n], dims=1))
T_avg_ocean = @lift Field(Average(T_ocean_ts[$n], dims=1))
sst_avg = @lift fill(mean(sst_slab_ts[$n]), 2)Observable([290.0, 290.0])
Plot
heatmap!(ax_θ, θn; colormap=:thermal, colorrange=(θᵃᵗ - 1, θᵃᵗ + 3))
heatmap!(ax_ss, ssn; colormap=:speed, colorrange=(0, 30))
heatmap!(ax_qˡ, qˡn; colormap=Reverse(:Blues_4), colorrange=(0, 5e-4))
heatmap!(ax_sf, sfn; colormap=:speed, colorrange=(0, 30))
heatmap!(ax_oT, oTn; colormap=:thermal, colorrange=(Tᵒᶜ - 1.5, Tᵒᶜ + 0.5))
lines!(ax_sst, sstn_slab; color=:red, linewidth=2, label="Slab (10m)")
lines!(ax_sst, x_ocean, ocean_sst; color=:blue, linewidth=2, label="Full")
axislegend(ax_sst, position=:rb)
ylims!(ax_sst, Tᵒᶜ - 0.7, Tᵒᶜ + 0.2)
lines!(ax_θp, θ_avg_slab; color=:red, linewidth=1.5, label="Slab")
lines!(ax_θp, θ_avg_full; color=:blue, linewidth=1.5, label="Full")
xlims!(ax_θp, θᵃᵗ - 1, θᵃᵗ + 4)
axislegend(ax_θp, position=:rt)
lines!(ax_up, u_avg_slab; color=:red, linewidth=1.5)
lines!(ax_up, u_avg_full; color=:blue, linewidth=1.5)
xlims!(ax_up, -10, 25)
lines!(ax_Tp, T_avg_ocean; color=:blue, linewidth=1.5, label="Full")
lines!(ax_Tp, sst_avg, [-50.0, 0.0]; color=:red, linewidth=1.5, label="Slab")
xlims!(ax_Tp, Tᵒᶜ - 1, Tᵒᶜ + 0.5)
title = @lift "Atmosphere–ocean coupling comparison, t = " * prettytime(times[$n])
Label(fig[0, 1:3], title, fontsize=16)Makie.Label()Record
@info "Rendering animation..."
record(fig, "breeze_over_two_oceans.mp4", 1:Nt; framerate=12) do nn
n[] = nn
end
@info "Animation saved."[ Info: Rendering animation...
[ Info: Animation saved.
This page was generated using Literate.jl.