feat(P3): add more P3 ice tendencies

new tendencies:
- ice self-collection
- ice sublimation / vapor deposition
- empirical ice nucleation (Frostenberg mean)

numerics:
- add numeric epsilon for rime volume
- squash a few bugs
- improve docs

Author: haakon-e

src/BulkMicrophysicsTendencies.jl

@@ -32,7 +32,9 @@ import ..Microphysics1M as CM1
 import ..Microphysics2M as CM2
 import ..MicrophysicsNonEq as CMNonEq
 import ..P3Scheme as CMP3
+import ..HetIceNucleation as CM_HetIce
 import ...ThermodynamicsInterface as TDI
+import ..Common as CO
 
 export MicrophysicsScheme,
     Microphysics0Moment,
@@ -801,9 +803,14 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
 """
 @inline function bulk_microphysics_tendencies(
     ::Microphysics2Moment, mp::CMP.Microphysics2MParams{WR, ICE}, tps,
-    ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai,
-    q_ice = zero(ρ), n_ice = zero(ρ), q_rim = zero(ρ), b_rim = zero(ρ), logλ = zero(ρ),
+    ρ, T, q_tot,
+    q_lcl, n_lcl, q_rai, n_rai,
+    q_ice, n_ice, q_rim, b_rim, logλ,
 ) where {WR, ICE <: CMP.P3IceParams}
+    FT = eltype(ρ)
+    ϵₘ = UT.ϵ_numerics_2M_M(FT)
+    ϵₙ = UT.ϵ_numerics_2M_M(FT)
+    ϵB = UT.ϵ_numerics_P3_B(FT)
     # Clamp negative inputs to zero (robustness against numerical errors)
     ρ = UT.clamp_to_nonneg(ρ)
     q_tot = UT.clamp_to_nonneg(q_tot)
@@ -816,13 +823,26 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
     q_rim = UT.clamp_to_nonneg(q_rim)
     b_rim = UT.clamp_to_nonneg(b_rim)
 
+    # Convert to volumetric quantities for P3 functions
+    L_lcl = q_lcl * ρ  # [kg lcl / m³ air]
+    L_rai = q_rai * ρ  # [kg rai / m³ air]
+    N_lcl = n_lcl * ρ  # [1 / m³ air]
+    N_rai = n_rai * ρ  # [1 / m³ air]
+    L_ice = q_ice * ρ  # [kg ice / m³ air]
+    N_ice = n_ice * ρ  # [1 / m³ air]
+    L_rim = q_rim * ρ  # [kg rim / m³ air]
+    B_rim = b_rim * ρ  # [m³ rim / m³ air]
+    # Compute rime fraction and density
+    F_rim = ifelse(q_ice > ϵₘ, q_rim / q_ice, zero(L_rim))
+    ρ_rim = ifelse(b_rim > ϵB, q_rim / b_rim, zero(L_rim))  # [kg rim / m³ rim]
+
     # Unpack warm rain parameters (always present)
     aps = mp.warm_rain.air_properties
+    subdep = mp.warm_rain.subdep
 
     # Initialize ice-related tendencies
     dq_ice_dt = zero(ρ)
-    # TODO: When ice number concentration becomes prognostic, add:
-    # dn_ice_dt = zero(ρ)  # Ice number tendency (changes due to melting, aggregation)
+    dn_ice_dt = zero(ρ)  # Ice number tendency (changes due to melting, aggregation)
     dq_rim_dt = zero(ρ)
     db_rim_dt = zero(ρ)
 
@@ -833,51 +853,24 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
     dq_rai_dt = warm.dq_rai_dt
     dn_rai_dt = warm.dn_rai_dt
 
-    # Convert to number densities for remaining functions
-    N_lcl = ρ * n_lcl
-    N_rai = ρ * n_rai
-
     # --- P3 Ice Processes ---
     p3 = mp.ice.scheme
     vel = mp.ice.terminal_velocity
     pdf_c = mp.ice.cloud_pdf
     pdf_r = mp.ice.rain_pdf
+    heterogeneous = mp.ice.heterogeneous
+    deposition_condfreeze = mp.ice.deposition_condfreeze
+
 
     # Only compute ice processes if there is ice mass/number present
-    if (q_ice > zero(q_ice) || n_ice > zero(n_ice))
-        # Convert to volumetric quantities for P3 functions
-        L_ice = q_ice * ρ  # [kg/m³]
-        N_ice = n_ice * ρ  # [1/m³]
-        L_lcl = q_lcl * ρ  # [kg/m³]
-        N_lcl = n_lcl * ρ  # [1/m³]
-        L_rai = q_rai * ρ  # [kg/m³]
-        N_rai = n_rai * ρ  # [1/m³]
-
-        # Compute rime fraction and density
-        F_rim = ifelse(q_ice > zero(q_ice), q_rim / q_ice, zero(q_rim))
-        ρ_rim = ifelse(b_rim > zero(b_rim), q_rim * ρ / (b_rim * ρ), ρ)  # [kg/m³]
+    if q_ice > ϵₘ && n_ice > ϵₙ
+        state = CMP3.P3State(p3, L_ice, N_ice, F_rim, ρ_rim)
 
         # Liquid-ice collision sources (core P3 process)
         # Only compute if there is ice present
         if L_ice > zero(L_ice) && N_ice > zero(N_ice)
             coll = CMP3.bulk_liquid_ice_collision_sources(
-                p3,
-                logλ,
-                L_ice,
-                N_ice,
-                F_rim,
-                ρ_rim,
-                pdf_c,
-                pdf_r,
-                L_lcl,
-                N_lcl,
-                L_rai,
-                N_rai,
-                aps,
-                tps,
-                vel,
-                ρ,
-                T,
+                state, logλ, pdf_c, pdf_r, L_lcl, N_lcl, L_rai, N_rai, aps, tps, vel, ρ, T,
             )
 
             # Add collision tendencies
@@ -895,7 +888,6 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
         # Ice melting (above freezing temperature)
         T_freeze = TDI.TD.Parameters.T_freeze(tps)
         if T > T_freeze && L_ice > zero(L_ice)
-            state = CMP3.P3State(p3, L_ice, N_ice, F_rim, ρ_rim)
             melt = CMP3.ice_melt(vel, aps, tps, T, ρ, state, logλ)
 
             # Melting converts ice to rain
@@ -907,9 +899,67 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
         end
     end
 
-    # TODO: When ice number concentration is tracked, add dn_ice_dt to return tuple:
-    # return (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt, dq_ice_dt, dn_ice_dt, dq_rim_dt, db_rim_dt)
-    return (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt, dq_ice_dt, dq_rim_dt, db_rim_dt)
+    # --- --------------- ---
+    # --- Ice Nucleation  ---
+    # --- --------------- ---
+    
+    # Note: Intuitively, you choose one of two parameterization pathways,
+    # either 1. process-based (e.g. using CM_HetIce.deposition_J, etc.) or
+    # 2. empirical (e.g. using CM_HetIce.INP_concentration_mean)
+    # 1. In the former, you try to differentiate between different aerosol types
+    # and their properties, homogeneous freezing, heterogeneous deposition
+    # freezing, heterogeneous immersion freezing, etc.
+    # 2. The latter, which we adopt here, empirically relates some environmental
+    # conditions (e.g. temperature) to the number of ice crystals nucleated.
+    # We do not track the aerosol population, and do not consider different
+    # nucleation pathways explicitly.
+
+    # --- Ice Nucleation (empirical) ---
+    # Assume mass loss is mean condensate mass
+    m_lcl = ifelse(n_lcl > ϵₙ, q_lcl / n_lcl, zero(q_lcl))  # mean liquid mass
+
+    # TODO: The parameterixation should return a rate, `∂N/∂t`, not number changes `ΔN`
+    inpc = CM_HetIce.INP_concentration_mean(heterogeneous, T) / ρ  # [particles / kg air]
+    n_het = max(0, inpc - n_ice)  # [particles / kg air]
+    q_het = n_het * m_lcl  # [kg ice / kg air]
+
+    dn_ice_dt += n_het
+    dq_ice_dt += q_het
+    dq_rim_dt += q_het
+    db_rim_dt += q_het / p3.ρ_i  # ρ_i = 916.7 kg m⁻³, the density of solid bulk ice
+    dn_lcl_dt -= n_het
+    dq_lcl_dt -= q_het
+
+    # --- Cloud Droplet Condensation Freezing ---
+    # ref: `homogeneous_freezing` in `parcel/ParcelTendencies.jl`
+    # get mean diameter of cloud droplets, then convert to volume
+
+    # --- Ice Sublimation / Deposition ---
+    # Deposition/sublimation of cloud ice
+    ∂ₜq_ice_dep = CMNonEq.conv_q_vap_to_q_lcl_icl_MM2015(subdep, tps, q_tot, q_lcl, q_ice, q_rai, zero(q_ice), ρ, T)
+    # No ice deposition above freezing (lack of INPs)
+    ∂ₜq_ice_dep = ifelse(T > tps.T_freeze, min(∂ₜq_ice_dep, zero(T)), ∂ₜq_ice_dep)
+    # During sublimation, the number of ice particles decreases in proportion to the mean ice mass
+    # During deposition, the number of ice particles remain unchanged
+    ∂ₜn_ice_dep = ifelse(∂ₜq_ice_dep < 0, (n_ice / q_ice) * ∂ₜq_ice_dep, zero(∂ₜq_ice_dep))
+    dq_ice_dt += ∂ₜq_ice_dep
+    dn_ice_dt += ∂ₜn_ice_dep
+
+    # --- Ice Self-collection (Aggregation) ---
+    if q_ice > ϵₘ && n_ice > ϵₙ
+        state = CMP3.P3State(p3, L_ice, N_ice, F_rim, ρ_rim)
+        S_ice_agg = CMP3.ice_self_collection(state, logλ, aps, tps, vel, ρ, T)
+        dn_ice_dt -= S_ice_agg.dNdt / ρ
+    end
+
+    # --- Rain Heterogeneous Freezing ---
+    # TODO: Implement heterogeneous freezing of rain
+    # This process is currently missing in P3_processes.jl
+    # S_rai_frz = ...
+    # dq_rai_dt -= S_rai_frz
+    # dq_ice_dt += S_rai_frz
+
+    return (; dq_lcl_dt, dn_lcl_dt, dq_rai_dt, dn_rai_dt, dq_ice_dt, dn_ice_dt, dq_rim_dt, db_rim_dt)
 end
 
 end # module BulkMicrophysicsTendencies

src/Common.jl

@@ -463,6 +463,22 @@ See e.g., Seifert and Beheng (2006), https://doi.org/10.1007/s00703-005-0112-4,
 
 # Returns
 - `F_v(D)`: ventilation factor function (dimensionless)
+
+# Extended help
+
+The ventilation factor is given by
+
+```math
+F_v(D) = a_v + b_v ⋅ ∛Sc ⋅ √Re(D)
+```
+
+where
+
+```math
+Sc = ν_{air} / D_{vapor}
+Re(D) = D * v_{term}(D) / ν_{air}
+```
+
 """
 @inline function ventilation_factor(vent, aps, v_term)
     (; aᵥ, bᵥ) = vent

src/P3_integral_properties.jl

@@ -90,7 +90,9 @@ Integrate the function `f` over each subinterval of the integration bounds, `bnd
 """
 function integrate(f, bnds...; quad = ChebyshevGauss(100))
     # compute integral over each subinterval (a, b), (b, c), (c, d), ...
-    return sum(integrate(f, a, b; quad) for (a, b) in zip(Base.front(bnds), Base.tail(bnds)))
+    return UU.unrolled_sum(
+        integrate(f, a, b; quad) for (a, b) in zip(Base.front(bnds), Base.tail(bnds))
+    )
 end
 
 """

src/P3_particle_properties.jl

@@ -74,10 +74,10 @@ Create a [`P3State`](@ref) from [`CMP.ParametersP3`](@ref) and rime state parame
 function get_state(params::CMP.ParametersP3; L_ice, N_ice, F_rim, ρ_rim)
     # rime mass fraction must always be non-negative ...
     # ... and there must always be some unrimed part
-    @assert 0 ≤ F_rim < 1 "Rime mass fraction, `F_rim`, must be between 0 and 1"
+     0 ≤ F_rim < 1 || throw(DomainError(F_rim, "Rime mass fraction, F_rim, must be 0 ≤ F_rim < 1"))
     # rime density must be positive ...
     # ... and as a bulk ice density can't exceed the density of water
-    @assert 0 < ρ_rim ≤ params.ρ_l "Rime density, `ρ_rim`, must be between 0 and ρ_l"
+    0 ≤ ρ_rim ≤ params.ρ_l || throw(DomainError(ρ_rim, "Rime density, ρ_rim, must be 0 ≤ ρ_rim ≤ ρ_l"))
     return P3State(; params, L_ice, N_ice, F_rim, ρ_rim)
 end
 
@@ -211,7 +211,7 @@ get_D_cr(mass::CMP.MassPowerLaw, F_rim, ρ_g) = _get_threshold(mass, ρ_g * (1 -
 # Returns
 - `(; D_th, D_gr, D_cr, ρ_g)`: The thresholds for the size distribution,
     and the density of total (deposition + rime) ice mass for graupel [kg/m³]
-
+    If `F_rim = 0`, then we set `D_gr = D_cr = Inf`, and `ρ_g = NaN` (should not be used).
 
 See [`get_D_th`](@ref), [`get_D_gr`](@ref), [`get_D_cr`](@ref), and [`get_ρ_g`](@ref) for more details.
 """
@@ -220,23 +220,27 @@ function get_thresholds_ρ_g(state::P3State)
     return get_thresholds_ρ_g(params, F_rim, ρ_rim)
 end
 function get_thresholds_ρ_g(params::CMP.ParametersP3, F_rim, ρ_rim)
+    FT = eltype(F_rim)
     (; mass, ρ_i) = params
+    D_th = get_D_th(mass, ρ_i)
+    
     ρ_d = get_ρ_d(mass, F_rim, ρ_rim)
     ρ_g = get_ρ_g(F_rim, ρ_rim, ρ_d)
-
-    D_th = get_D_th(mass, ρ_i)
     D_gr = get_D_gr(mass, ρ_g)
     D_cr = get_D_cr(mass, F_rim, ρ_g)
 
+    # If unrimed, set D_gr and D_cr to infinity
+    D_gr = ifelse(iszero(F_rim), FT(Inf), D_gr)
+    D_cr = ifelse(iszero(F_rim), FT(Inf), D_cr)
+
     return (; D_th, D_gr, D_cr, ρ_g)
 end
 
-function get_bounded_thresholds(state::P3State, D_min = 0, D_max = Inf)
-    FT = eltype(state)
+function get_bounded_thresholds(
+    state::P3State{FT}, D_min::FT = FT(0), D_max::FT = FT(Inf)
+) where {FT}
     (; D_th, D_gr, D_cr) = get_thresholds_ρ_g(state)
-    thresholds = clamp.((FT(D_min), D_th, D_gr, D_cr, FT(D_max)), FT(D_min), FT(D_max))
-    bounded_thresholds = replace(thresholds, NaN => FT(D_max))
-    return bounded_thresholds
+    return clamp.((D_min, D_th, D_gr, D_cr, D_max), D_min, D_max)
 end
 
 """
@@ -253,9 +257,9 @@ Return the segments of the size distribution as a tuple of intervals.
 For example, if the thresholds are `(D_th, D_gr, D_cr)`, then the segments are:
 - `(0, D_th)`, `(D_th, D_gr)`, `(D_gr, D_cr)`, `(D_cr, Inf)`
 """
-function get_segments(state::P3State, D_min = 0, D_max = Inf)
-    thresholds = get_bounded_thresholds(state, D_min, D_max)
-    segments = tuple.(Base.front(thresholds), Base.tail(thresholds))
+function get_segments(state::P3State{FT}, D_min::FT = FT(0), D_max::FT = FT(Inf)) where {FT}
+    (D_min, D_th, D_gr, D_cr, D_max) = get_bounded_thresholds(state, D_min, D_max)
+    segments = ((D_min, D_th), (D_th, D_gr), (D_gr, D_cr), (D_cr, D_max))
     return segments
 end