Skip to content

Use condensate in evaporation/sublimation #708

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
trontrytel merged 1 commit into from
Mar 25, 2026
Merged

Use condensate in evaporation/sublimation #708

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "CloudMicrophysics"
uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
version = "0.32"
version = "0.33"
authors = ["Climate Modeling Alliance"]

[deps]
Expand Down
3 changes: 0 additions & 3 deletions docs/src/API.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,10 @@ MicrophysicsNonEq
MicrophysicsNonEq.τ_relax
MicrophysicsNonEq.conv_q_vap_to_q_lcl_icl
MicrophysicsNonEq.conv_q_vap_to_q_lcl_icl_MM2015
MicrophysicsNonEq.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld
MicrophysicsNonEq.limit_MM2015_sinks
MicrophysicsNonEq.INP_limiter
MicrophysicsNonEq.terminal_velocity
MicrophysicsNonEq.dqcld_dT
MicrophysicsNonEq.gamma_helper
MicrophysicsNonEq.d2qcld_dT2
```

# 0-moment precipitation microphysics
Expand Down
30 changes: 0 additions & 30 deletions docs/src/MicrophysicsNonEq.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,38 +124,8 @@ include("plots/NonEqCondEvapRate.jl")
![](condensation_evaporation_ql_z.svg)
![](condensation_evaporation_ql_T.svg)

### Derivative of the tendency

For implicit time integration, we also compute the total derivative of the tendency
with respect to the cloud condensate, accounting for the implicit feedback
between the condensate and temperature through latent heat release
(``dT/dq = L/c_p``):

```math
\begin{equation}
\frac{d\dot{q}_{lcl}}{dq_{lcl}} = -\frac{1}{\tau_l}
- \frac{\dot{q}_{lcl}}{\Gamma_l} \left(\frac{L_v}{c_p}\right)^2 \frac{d^2 q_{sl}}{dT^2};
\;\;\;\;\;\;\;\;
\frac{d\dot{q}_{icl}}{dq_{icl}} = -\frac{1}{\tau_i}
- \frac{\dot{q}_{icl}}{\Gamma_i} \left(\frac{L_s}{c_p}\right)^2 \frac{d^2 q_{si}}{dT^2}
\end{equation}
```
where ``\dot{q}_{lcl}`` and ``\dot{q}_{icl}`` denote the condensation/evaporation
and deposition/sublimation tendencies respectively,
and the second derivative of saturation specific humidity is:
```math
\begin{equation}
\frac{d^2 q_{sl}}{dT^2} = q_{sl} \left[\left(\frac{L_v}{R_v T^2} - \frac{1}{T}\right)^2
+ \frac{1}{T^2} - \frac{2 L_v}{R_v T^3}\right]
\end{equation}
```

The leading-order term ``-1/\tau`` provides a simple approximation,
while the second term captures the higher-order correction from the
thermodynamic adjustment factor ``\Gamma``.

![](condensation_evaporation_deriv.svg)
![](condensation_evaporation_deriv_vs_qvap.svg)


## Cloud condensate sedimentation
Expand Down
195 changes: 0 additions & 195 deletions docs/src/plots/NonEqCondEvapRate.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,6 @@ function generate_cond_evap_rate(::Range_qₗ_T)
)
end

# --- Derivative line plots ---

minussign(s) = replace(s, "-" => "−")

Expand All @@ -122,200 +121,6 @@ function make_figure(case)
save(file_name, fig)
end

function make_deriv_line_plot()
# Temperature range
T_range = collect(220:0.5:310)

# Fixed parameters
qₜ = FT(10e-3) # 10 g/kg total water
qₗ = FT(1e-3) # 1 g/kg cloud liquid
qᵢ = FT(1e-3) # 1 g/kg cloud ice
qᵣ = FT(0)
qₛ = FT(0)
ρ = FT(1) # kg/m³

τ_relax = FT(1) # 1 s relaxation time
_liq = CM.Parameters.CloudLiquid(FT)
cm_liq = CM.Parameters.CloudLiquid(τ_relax, _liq.ρw, _liq.r_eff, _liq.N_0)
_ice = CM.Parameters.CloudIce(FT)
cm_ice = CM.Parameters.CloudIce(_ice.pdf, _ice.mass, _ice.r0, _ice.r_ice_snow, τ_relax, _ice.ρᵢ, _ice.r_eff, _ice.N_0)

# Compute tendency and derivatives for each temperature
S_liq = zeros(FT, length(T_range))
dS_liq_simple = zeros(FT, length(T_range))
dS_liq_full = zeros(FT, length(T_range))

S_ice = zeros(FT, length(T_range))
dS_ice_simple = zeros(FT, length(T_range))
dS_ice_full = zeros(FT, length(T_range))

for (i, T) in enumerate(T_range)
T_ft = FT(T)
s = CMNe.conv_q_vap_to_q_lcl_icl_MM2015(
cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft,
)
ds_s = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(
cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft,
)
ds_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(
cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft;
simplified = false,
)
S_liq[i] = s
dS_liq_simple[i] = ds_s
dS_liq_full[i] = ds_f

s = CMNe.conv_q_vap_to_q_lcl_icl_MM2015(
cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft,
)
ds_s = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(
cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft,
)
ds_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(
cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T_ft;
simplified = false,
)
S_ice[i] = s
dS_ice_simple[i] = ds_s
dS_ice_full[i] = ds_f
end

T_freeze = FT(TD.Parameters.T_freeze(thp))

# --- Figure: 2 panels ---
fig = Figure(size = (900, 450))

# Left panel: tendency vs T
ax1 = Axis(fig[1, 1];
xlabel = "Temperature (K)",
ylabel = "Tendency (1/s)",
title = "Condensation/Evaporation & Deposition/Sublimation\n(qₜ=10, qₗ=1, qᵢ=1 g/kg, ρ=1 kg/m³, τ=1 s)",
)
lines!(ax1, T_range, S_liq; color = :blue, label = "Liquid (cond/evap)")
lines!(ax1, T_range, S_ice; color = :magenta, label = "Ice (dep/subl)")
vlines!(ax1, [T_freeze]; color = :gray, linestyle = :dash, label = "T_freeze")
hlines!(ax1, [0]; color = :black, linewidth = 0.5)
text!(ax1, T_freeze + 2, minimum(S_ice) * 0.5;
text = "INP limiter\n(ice dep → 0)",
fontsize = 10, color = :gray,
)
axislegend(ax1; position = :lb, framevisible = false)

# Right panel: derivative vs T
ax2 = Axis(fig[1, 2];
xlabel = "Temperature (K)",
ylabel = "∂(tendency)/∂q (1/s)",
title = "Derivative of tendency w.r.t. cloud condensate",
)
lines!(ax2, T_range, dS_liq_full; color = :blue, label = "Liquid (full)")
lines!(ax2, T_range, dS_ice_full; color = :magenta, label = "Ice (full)")
# Simple derivative is -1/τ for both phases (identical when τ_liq == τ_ice)
lines!(ax2, T_range, dS_liq_simple; color = :black, linestyle = :dash, label = "Simple = −1/τ (both)")
vlines!(ax2, [T_freeze]; color = :gray, linestyle = :dash, label = "T_freeze")
text!(ax2, T_freeze + 2, minimum(dS_ice_full) * 0.5;
text = "INP limiter\n(∂ice → 0)",
fontsize = 10, color = :gray,
)
axislegend(ax2; position = :lt, framevisible = false)

save("condensation_evaporation_deriv.svg", fig)
end

function make_deriv_vs_qvap_plot()
# Available water vapor range (q_vap = q_tot - q_lcl - q_icl - q_rai - q_sno)
# We vary q_tot while keeping q_lcl, q_icl fixed, so q_vap changes
qₗ = FT(1e-3) # 1 g/kg cloud liquid (fixed)
qᵢ = FT(1e-3) # 1 g/kg cloud ice (fixed)
qᵣ = FT(0)
qₛ = FT(0)
ρ = FT(1) # kg/m³

# q_tot range: from 2 g/kg (all condensate, no vapor) to 20 g/kg
qₜ_range = collect(range(FT(2e-3), stop = FT(20e-3), length = 181))
qᵥ_range = qₜ_range .- qₗ .- qᵢ # available water vapor

τ_relax = FT(1)
_liq = CM.Parameters.CloudLiquid(FT)
cm_liq = CM.Parameters.CloudLiquid(τ_relax, _liq.ρw, _liq.r_eff, _liq.N_0)
_ice = CM.Parameters.CloudIce(FT)
cm_ice = CM.Parameters.CloudIce(_ice.pdf, _ice.mass, _ice.r0, _ice.r_ice_snow, τ_relax, _ice.ρᵢ, _ice.r_eff, _ice.N_0)

# Evaluate at two representative temperatures
T_vals = [FT(250), FT(280)]
T_labels = ["T=250K", "T=280K"]
T_colors_liq = [:royalblue, :cornflowerblue]
T_colors_ice = [:magenta, :hotpink]

fig = Figure(size = (900, 450))

ax1 = Axis(fig[1, 1];
xlabel = "Water vapor, qᵥ (g/kg)",
ylabel = "Tendency (1/s)",
title = "Tendency vs available water vapor\n(qₗ=1, qᵢ=1 g/kg, ρ=1 kg/m³, τ=1 s)",
)
ax2 = Axis(fig[1, 2];
xlabel = "Water vapor, qᵥ (g/kg)",
ylabel = "∂(tendency)/∂q (1/s)",
title = "Derivative vs available water vapor",
)

for (j, T) in enumerate(T_vals)
S_liq = zeros(FT, length(qₜ_range))
dS_liq_f = zeros(FT, length(qₜ_range))
dS_liq_s = zeros(FT, length(qₜ_range))
S_ice = zeros(FT, length(qₜ_range))
dS_ice_f = zeros(FT, length(qₜ_range))

for (i, qₜ) in enumerate(qₜ_range)
s = CMNe.conv_q_vap_to_q_lcl_icl_MM2015(
cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T,
)
ds_s = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(
cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T,
)
ds_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(
cm_liq, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T;
simplified = false,
)
S_liq[i] = s
dS_liq_f[i] = ds_f
dS_liq_s[i] = ds_s

s = CMNe.conv_q_vap_to_q_lcl_icl_MM2015(
cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T,
)
ds_f = CMNe.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(
cm_ice, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T;
simplified = false,
)
S_ice[i] = s
dS_ice_f[i] = ds_f
end

qᵥ_gkg = qᵥ_range * 1e3

# Left panel: tendency
lines!(ax1, qᵥ_gkg, S_liq; color = T_colors_liq[j], label = "Liquid $(T_labels[j])")
lines!(ax1, qᵥ_gkg, S_ice; color = T_colors_ice[j], label = "Ice $(T_labels[j])")

# Right panel: derivative
lines!(ax2, qᵥ_gkg, dS_liq_f; color = T_colors_liq[j], label = "Liquid $(T_labels[j])")
lines!(ax2, qᵥ_gkg, dS_ice_f; color = T_colors_ice[j], label = "Ice $(T_labels[j])")

if j == 1
lines!(ax2, qᵥ_gkg, dS_liq_s; color = :black, linestyle = :dash, label = "Simple = −1/τ")
end
end

hlines!(ax1, [0]; color = :black, linewidth = 0.5)
axislegend(ax1; position = :lt, framevisible = false)
axislegend(ax2; position = :lb, framevisible = false)

save("condensation_evaporation_deriv_vs_qvap.svg", fig)
end

make_figure(HydrostaticBalance_qₗ_z())
make_figure(Range_qₗ_T())
make_deriv_line_plot()
make_deriv_vs_qvap_plot()

6 changes: 3 additions & 3 deletions src/BulkMicrophysicsTendencies.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -434,7 +434,7 @@ and avoids quadratures over the derivatives.
vel = mp.terminal_velocity

# --- Cloud liquid: condensation + sink self-derivatives ---
∂tendency_∂q_lcl = CMNonEq.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(lcl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T)
∂tendency_∂q_lcl = -1 / CMNonEq.τ_relax(lcl)
# Autoconversion q_lcl → q_rai (sink)
S_acnv_lcl = CM1.conv_q_lcl_to_q_rai(rai.acnv1M, q_lcl, true)
# Accretion q_lcl + q_rai → q_rai (rate is exactly linear in q_lcl)
Expand All @@ -446,7 +446,7 @@ and avoids quadratures over the derivatives.
max(q_lcl, UT.ϵ_numerics(typeof(q_lcl)))

# --- Cloud ice: deposition + sink self-derivatives ---
∂tendency_∂q_icl = CMNonEq.∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld(icl, tps, q_tot, q_lcl, q_icl, q_rai, q_sno, ρ, T)
∂tendency_∂q_icl = -1 / CMNonEq.τ_relax(icl)
# Autoconversion q_icl → q_sno
S_acnv_icl = CM1.conv_q_icl_to_q_sno_no_supersat(sno.acnv1M, q_icl, true)
# Accretion q_icl + q_rai → q_sno (rate is exactly linear in q_icl)
Expand Down Expand Up @@ -519,7 +519,7 @@ derivatives; snow and cloud formation derivatives are zero for now.
N_rai = ρ * n_rai

# TODO: Cloud formation — 2M bulk_microphysics_tendencies does not call cloud formation yet;
# once it does, set these from MM2015 (same as 1M) via ∂conv_q_vap_to_q_lcl_icl_MM2015_∂q_cld.
# once it does, set these from MM2015 (same as 1M) via -1/τ_relax.
∂tendency_∂q_lcl = zero(ρ)
∂tendency_∂q_icl = zero(ρ)

Expand Down
Loading
Loading