CloudToppedMixedLayerModel
using DynamicalSystems, ConceptualClimateModels
import ConceptualClimateModels.CloudToppedMixedLayerModel as CTMLMFixed states of Stevens 2006
Here we want to recreate the same bulk boundary layer model defined in (Stevens, 2006), equations 31-33.
const cₚ = 1004.0 # heat capacity at constant pressure (J/K/kg)
eqs = [
CTMLM.mlm_dynamic(),
CTMLM.entrainment_velocity(:Stevens2006; use_augmentation = false),
## Figure 1 and Sec. 4.2 provide these values:
CTMLM.s₊ ~ 301200.0/cₚ,
CTMLM.q₊ ~ 1.56,
CTMLM.s₀ ~ CTMLM.s₊ - 12.5,
CTMLM.ρ₀ ~ 1,
CTMLM.q₀ ~ CTMLM.q_saturation(288.96 + 1.25),
CTMLM.ΔF_s ~ 40.0
]8-element Vector{Any}:
TimeDerivative[TimeDerivative(z_b(t), w_D(t) + w_e(t) + w_m(t), LiteralParameter{Float64}(1.1574074074074073e-5)), TimeDerivative(s_b(t), -(1//86400)*s_x(t) + ((SHF(t) - ΔF_s(t)) / (1004.0ρ₀(t)) + w_e(t)*Δ₊s(t)) / z_b(t), LiteralParameter{Float64}(1.1574074074074073e-5)), TimeDerivative(q_b(t), -(1//86400)*q_x(t) + ((LHF(t) - ΔF_q(t)) / (2530.0ρ₀(t)) + w_e(t)*Δ₊q(t)) / z_b(t), LiteralParameter{Float64}(1.1574074074074073e-5))]
w_e(t) ~ (e_e*ΔF_s(t)) / (1004.0Δ₊s(t)*ρ₀(t))
s₊(t) ~ 300.0
q₊(t) ~ 1.56
s₀(t) ~ -12.5 + s₊(t)
ρ₀(t) ~ 1
q₀(t) ~ 12.404970818808321
ΔF_s(t) ~ 40.0We now give these equations to the main function that creates a DynamicalSystem from the processes (and we also provide CTMLM to obtain default processes from)
ds = processes_to_coupledodes(eqs, CTMLM)3-dimensional CoupledODEs
deterministic: true
discrete time: false
in-place: false
dynamic rule: GeneratedFunctionWrapper
ODE solver: Tsit5
ODE kwargs: (abstol = 1.0e-6, reltol = 1.0e-6)
parameters: [0.9, 3.0e-6, 0.0012, 6.0, 0.0, 0.0, 0.0, 0.0, 11.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
time: 0.0
state: [1200.0, 290.0, 11.0]
We also make sure to use the same parameter values as in the paper:
set_parameter!(ds, :D, 4e-6)
set_parameter!(ds, :d_c, 0.0011)
set_parameter!(ds, :U, 0.008/0.0011) # U*d_c = 0.008
set_parameter!(ds, :e_e, 1.0) # effective emissivityAnd that's all. We can run the system to its steady state:
step!(ds, 100.0)3-dimensional CoupledODEs
deterministic: true
discrete time: false
in-place: false
dynamic rule: GeneratedFunctionWrapper
ODE solver: Tsit5
ODE kwargs: (abstol = 1.0e-6, reltol = 1.0e-6)
parameters: [1.0, 4.0e-6, 0.0011, 7.2727272727272725, 0.0, 0.0, 0.0, 0.0, 11.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
time: 100.11469236113098
state: [796.8127490039907, 287.4999999999997, 9.31524020812321]
leading to a steady state with height $z_b$ ($h_\infty$ in (Stevens, 2006)) of about 800 meters as in the paper. Extracting the variable $\sigma$ of Eq. 38 is also very easy:
Δs = CTMLM.s₊ - CTMLM.s₀ # symbolic variable not existing in the graph of the `ds`
σ = observe_state(ds, CTMLM.V * Δs / CTMLM.ΔF_s * cₚ)2.5100000000000002and we get the same value (note multiplication by cₚ, because s is in units of K).
If you want to see a list of all equations that compose the dynamical system then do
all_equations(ds)22-element Vector{Equation}:
Differential(t, 1)(z_b(t)) ~ (-w_D(t) - w_e(t)) / -1.1574074074074073e-5
Differential(t, 1)(s_b(t)) ~ (-((SHF(t) - ΔF_s(t)) / (1004.0ρ₀(t))) - w_e(t)*Δ₊s(t)) / (-1.1574074074074073e-5z_b(t))
Differential(t, 1)(q_b(t)) ~ (-(LHF(t) / (2530.0ρ₀(t))) - w_e(t)*Δ₊q(t)) / (-1.1574074074074073e-5z_b(t))
ρ₀(t) ~ 1
ΔF_s(t) ~ 40.0
s₊(t) ~ 300.0
q₊(t) ~ 1.56
q₀(t) ~ 12.404970818808321
V(t) ~ U*d_c
w_D(t) ~ -D*z_b(t)
⋮
ΔF_q(t) ~ 0
s₀(t) ~ -12.5 + s₊(t)
Δ₊s(t) ~ -s_b(t) + s₊(t)
Δ₊q(t) ~ -q_b(t) + q₊(t)
Δ₀q(t) ~ q_b(t) - q₀(t)
Δ₀s(t) ~ s_b(t) - s₀(t)
w_e(t) ~ (e_e*ΔF_s(t)) / (1004.0Δ₊s(t)*ρ₀(t))
LHF(t) ~ -2530.0V(t)*Δ₀q(t)*ρ₀(t)
SHF(t) ~ -1004.0V(t)*Δ₀s(t)*ρ₀(t)In the compiled documentation these render via LaTeX but running this in the REPL won't be as pretty :)
Adding clouds
We can add clouds to the mixed layer model using the same decoupling-based approach as in (Datseris, 2026) (while keeping SST a fixed boundary condition) just by including a couple more equations to the ones already defined, so that $C, \mathcal{D}$ have a process assigned to them. We will also augment ΔF to be partially proportional to C using a very simple ad-hoc approach.
eqs = [
CTMLM.mlm_dynamic(),
CTMLM.entrainment_velocity(:Stevens2006; use_augmentation = false),
## Cloud stuff
CTMLM.cf_dynamic(),
CTMLM.decoupling_variable(),
CTMLM.cloud_base_height(:Bolton1980),
CTMLM.CTRC ~ 10 + 40*CTMLM.C, # cloud top radiative cooling
CTMLM.ΔF_s ~ CTMLM.CTRC, # same equation
## rest the same
CTMLM.s₊ ~ 301200.0/cₚ,
CTMLM.q₊ ~ 1.56,
CTMLM.s₀ ~ CTMLM.s₊ - 12.5,
CTMLM.ρ₀ ~ 1,
CTMLM.q₀ ~ CTMLM.q_saturation(288.96 + 1.25),
]
ds = processes_to_coupledodes(eqs, CTMLM)4-dimensional CoupledODEs
deterministic: true
discrete time: false
in-place: false
dynamic rule: GeneratedFunctionWrapper
ODE solver: Tsit5
ODE kwargs: (abstol = 1.0e-6, reltol = 1.0e-6)
parameters: [0.9, 1.0, 0.2, 0.5, 1.0, 100.0, 2.0, 3.0e-6, 0.0012, 6.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
time: 0.0
state: [1200.0, 290.0, 11.0, 1.0]
And we run the model to a steady state:
set_parameter!(ds, :D, 4e-6)
set_parameter!(ds, :d_c, 0.0009)
set_parameter!(ds, :U, 6.8) # U*d_c = 0.008
set_parameter!(ds, :e_e, 1.0) # effective emissivity
step!(ds, 100.0)
observe_state.(ds, (:z_b, :q_b, :s_b, :C)) # get state variables by name(763.5014057710496, 8.794700233031584, 287.49999999999983, 0.7081938746298164)Climate change scenario: increasing SST
If you want to study multistability for alternate cloud states (Cumulus vs Stratocumulus), or perform continuations (like the climate change scenarios in (Datseris, 2026)), visit the documentation of Attractors.jl.
As a brief example we will perform a simple climate change scenario where SST increases with a constant rate. First, we modify the equations so that s₀, q₀ are derived from a prescribed SST:
@parameters SST = 290.0
eqs = [
CTMLM.mlm_dynamic(),
CTMLM.entrainment_velocity(:Stevens2006; use_augmentation = false),
## Cloud stuff
CTMLM.cf_dynamic(),
CTMLM.decoupling_variable(),
CTMLM.cloud_base_height(:Bolton1980),
CTMLM.CTRC ~ 10 + 40*CTMLM.C, # cloud top radiative cooling
CTMLM.ΔF_s ~ CTMLM.CTRC, # same equation
## Usage of SST
CTMLM.s₀ ~ SST,
CTMLM.q₀ ~ CTMLM.q_saturation(SST),
## rest the same
CTMLM.s₊ ~ 301200.0/cₚ,
CTMLM.q₊ ~ 1.56,
CTMLM.ρ₀ ~ 1,
]
ds = processes_to_coupledodes(eqs, CTMLM)
set_parameter!(ds, :D, 4e-6)
set_parameter!(ds, :d_c, 0.0009)
set_parameter!(ds, :U, 6.8) # U*d_c = 0.008
set_parameter!(ds, :e_e, 1.0) # effective emissivity
ds4-dimensional CoupledODEs
deterministic: true
discrete time: false
in-place: false
dynamic rule: GeneratedFunctionWrapper
ODE solver: Tsit5
ODE kwargs: (abstol = 1.0e-6, reltol = 1.0e-6)
parameters: [1.0, 290.0, 1.0, 0.2, 0.5, 1.0, 100.0, 2.0, 4.0e-6, 0.0009 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
time: 0.0
state: [1200.0, 290.0, 11.0, 1.0]
To perform global continuation we need to create an attractor mapper, which we create here by tessellating the state space. Have a look at the Attractors.jl tutorial if the next few lines of code are puzzling to you
grid = (
(1:1.0:25), # specific humidity
(270:5.0:330), # static energy
(0:100.0:3000), # height
(0:0.1:1), # cloud fraction
)
sampler, = statespace_sampler(grid)
mapper = AttractorsViaRecurrences(ds, grid)
ascm = AttractorSeedContinueMatch(mapper)
fs = basins_fractions(mapper, sampler; show_progress = false)
attractors = extract_attractors(mapper)Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()with the sampler, mapper, ascm data structures in order, we can easily now run a global continuation with changing SST:
prange = 290:1:310
pidx = :SST
fractions_cont, attractors_cont = global_continuation(ascm, prange, pidx, sampler; show_progress = false)
attractors_cont21-element Vector{Dict}:
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
⋮
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()
Dict{Int64, StateSpaceSet{4, Float64, SVector{4, Float64}}}()if you are not sure what the output means, no worries, just have a look at the Attractors.jl documentation. Here we visualize the cloud fraction:
using CairoMakie
# cloud fraction and height values of last point on the attractor
varvalues = [A -> A[end][4], A -> A[end][3]]
fig = plot_basins_attractors_curves(
fractions_cont, attractors_cont, varvalues, prange,
)
content(fig[2,1]).ylabel = "C"
content(fig[3,1]).ylabel = "z_b"
fig
we see for our ad hoc parameterizations, the dynamical system has no stable states for SST > 297. Before that cloud fraction decreases very fast from a Stratocumulus state to a Cumulus one.
The model diverges due to the definition of the entrainment velocity (which is proportional to ΔF_s/Δ₊s) in combination with a fixed s₊: the model reaches a state where Δ₊s is so small that height increases unnaturally well beyond the state space tessellation we defined.