Merge pull request #708 from CliMA/aj/evaporation_fix

Use condensate in evaporation/sublimation

Author: trontrytel

Project.toml

@@ -1,6 +1,6 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
-version = "0.32"
+version = "0.33"
 authors = ["Climate Modeling Alliance"]
 
 [deps]

docs/src/API.md

@@ -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

docs/src/MicrophysicsNonEq.md

@@ -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

docs/src/plots/NonEqCondEvapRate.jl

@@ -98,7 +98,6 @@ function generate_cond_evap_rate(::Range_qₗ_T)
     )
 end
 
-# --- Derivative line plots ---
 
 minussign(s) = replace(s, "-" => "−")
 
@@ -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()
 

src/BulkMicrophysicsTendencies.jl

@@ -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)
@@ -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)
@@ -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(ρ)