feat(P3): add more P3 ice tendencies

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

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,
@@ -804,6 +806,9 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch.
     ρ, 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(ρ),
 ) where {WR, ICE <: CMP.P3IceParams}
+    ϵₘ = 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 +821,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 > ϵₘ, L_rim / L_ice, zero(L_rim))
+    ρ_rim = ifelse(b_rim > ϵB, L_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 +851,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 > ϵₙ
 
         # 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,
+                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,
             )
 
             # Add collision tendencies
@@ -907,9 +898,69 @@ 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 Modes ---
+    # --- -------------------- ---
+
+    # --- Ice Nucleation (Deposition) ---
+    # Assume nucleated particles have diameter 1μm --> nucleated mass (per particle)
+    D_nuc_ice = FT(1e-6) # 1μm
+    m_nuc = p3.ρ_i * CO.volume_sphere_D(D_nuc_ice)  # [kg]
+
+    # TODO: The parameterixation should return a rate, `∂N/∂t`, not number changes `ΔN`
+    n_dep = CM_HetIce.P3_deposition_N_i(deposition_condfreeze, T) / ρ  # [particles / kg air]
+    m_dep = n_dep * m_nuc  # [kg ice / kg air]
+    dn_ice_dt += n_dep
+    dq_ice_dt += m_dep
+    dq_rim_dt += m_dep
+    db_rim_dt += m_dep / 900
+
+    # --- Heterogeneous Ice Nucleation (Immersion freezing) ---
+    # 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 / 900  # = ρ^* TODO: make this a parameter?
+    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_processes.jl

@@ -382,8 +382,8 @@ function ∫liquid_ice_collisions(
     # Initialize integration buffers by evaluating a representative integral
     p = FT(0.00001)
     ice_bounds = integral_bounds(state, logλ; p)
-    bounds_c = CM2.get_size_distribution_bounds(psd_c, L_c, ρₐ, N_c, p)
-    bounds_r = CM2.get_size_distribution_bounds(psd_r, L_r, ρₐ, N_r, p)
+    bounds_c = CM2.get_size_distribution_bounds(psd_c, L_c / ρₐ, ρₐ, N_c, p)
+    bounds_r = CM2.get_size_distribution_bounds(psd_r, L_r / ρₐ, ρₐ, N_r, p)
 
     # Integrand components
     # NOTE: We assume collision efficiency, shape (spherical), and terminal velocity is the
@@ -488,3 +488,75 @@ function bulk_liquid_ice_collision_sources(
         (∂ₜq_c, ∂ₜq_r, ∂ₜN_c, ∂ₜN_r, ∂ₜL_rim, ∂ₜL_ice, ∂ₜB_rim)
     )
 end
+
+function collision_cross_section_ice_ice(state, D_1, D_2)
+    r_eff(D) = √(ice_area(state, D) / π)
+    return π * (r_eff(D_1) + r_eff(D_2))^2  # collision cross section
+end
+
+"""
+    volumetric_ice_ice_collision_rate_integrand(state, velocity_params, ρₐ)
+
+Returns a function that computes the volumetric collision rate integrand for ice-ice collisions [m³/s].
+
+# Arguments
+- `state`: [`P3State`](@ref)
+- `velocity_params`: velocity parameterization, e.g. [`CMP.Chen2022VelType`](@ref)
+- `ρₐ`: air density
+
+# Returns
+A function `(D_1, D_2) -> E * K * |vᵢ(D_1) - vᵢ(D_2)|` where:
+- `D_1` and `D_2` are the (maximum) diameters of the ice particles
+- `E` is the collision efficiency
+- `K` is the collision cross section
+- `vᵢ` is the terminal velocity of ice particles
+"""
+function volumetric_ice_ice_collision_rate_integrand(velocity_params, ρₐ, state)
+    v_ice = ice_particle_terminal_velocity(velocity_params, ρₐ, state)
+    function integrand(D_1::FT, D_2::FT) where {FT}
+        E = FT(1)  # Collision efficiency
+        K = collision_cross_section_ice_ice(state, D_1, D_2)
+        return E * K * abs(v_ice(D_1) - v_ice(D_2))
+    end
+    return integrand
+end
+
+"""
+    ice_self_collection(state, logλ, aps, tps, vel, ρₐ, T; ∫kwargs...)
+
+Computes the ice self-collection (aggregation) rate, which decreases the ice number concentration
+while leaving mass, rime mass, and rime volume unchanged.
+
+# Arguments
+- `state`: [`P3State`](@ref)
+- `logλ`: the log of the slope parameter [log(1/m)]
+- `aps`: [`CMP.AirProperties`](@ref)
+- `tps`: thermodynamics parameters
+- `vel`: the velocity parameterization, e.g. [`CMP.Chen2022VelType`](@ref)
+- `ρₐ`: air density [kg/m³]
+- `T`: temperature [K]
+
+# Returns
+A `NamedTuple` of `(; dNdt)`, where:
+1. `dNdt`: ice number concentration tendency due to self-collection [1/m³/s] (always positive or zero, represents a loss rate)
+"""
+function ice_self_collection(state, logλ, aps, tps, vel, ρₐ, T; ∫kwargs...)
+    n_i = DT.size_distribution(state, logλ)
+    ∂ₜV = volumetric_ice_ice_collision_rate_integrand(vel, ρₐ, state)
+
+    p = eps(one(ρₐ))
+    ice_bounds = integral_bounds(state, logλ; p)
+
+    function inner_integral(D_1)
+        (rate_at_D1,) = integrate(ice_bounds...; ∫kwargs...) do D_2
+            return SA.SVector(∂ₜV(D_1, D_2) * n_i(D_2))
+        end
+        return rate_at_D1 * n_i(D_1)
+    end
+
+    (total_rate,) = integrate(inner_integral, ice_bounds...; ∫kwargs...)
+
+    # The 0.5 factor accounts for double-counting in self-collection
+    dNdt = (1 // 2) * total_rate
+    return (; dNdt)
+end

src/Utilities.jl

@@ -6,7 +6,7 @@ Contains pure numerical operations with no physics dependencies.
 """
 module Utilities
 
-export clamp_to_nonneg, ϵ_numerics, ϵ_numerics_2M_M, ϵ_numerics_2M_N
+export clamp_to_nonneg, ϵ_numerics, ϵ_numerics_2M_M, ϵ_numerics_2M_N, ϵ_numerics_P3_B
 
 """
     clamp_to_nonneg(x)
@@ -44,4 +44,12 @@ Numerical epsilon for 2-moment number calculations.
 """
 @inline ϵ_numerics_2M_N(FT) = eps(FT)
 
+"""
+    ϵ_numerics_P3_B(FT)
+
+Numerical epsilon for P3 bulk microphysics mass calculations
+    relating to rim volume, B_rim.
+"""
+@inline ϵ_numerics_P3_B(FT) = eps(FT)
+
 end # module

src/parameters/Microphysics2MParams.jl

@@ -9,17 +9,19 @@ Parameters for 2-moment warm rain processes (Seifert-Beheng 2006).
 - `seifert_beheng::SB`: SB2006 — all warm rain parameters (autoconversion, accretion, etc.)
 - `air_properties::AP`: AirProperties — air properties for evaporation
 """
-@kwdef struct WarmRainParams2M{SB, AP, CE} <: ParametersType
+@kwdef struct WarmRainParams2M{SB, AP, CE, SD} <: ParametersType
     seifert_beheng::SB
     air_properties::AP
     condevap::CE
+    subdep::SD
 end
 # Construct WarmRainParams2M from a ClimaParams TOML dictionary
 WarmRainParams2M(toml_dict::CP.ParamDict; is_limited = true) =
     WarmRainParams2M(;
         seifert_beheng = SB2006(toml_dict; is_limited),
         air_properties = AirProperties(toml_dict),
         condevap = CondEvap2M(toml_dict),
+        subdep = SubDep2M(toml_dict),
     )
 
 Base.show(io::IO, mime::MIME"text/plain", x::WarmRainParams2M) =
@@ -28,19 +30,39 @@ Base.show(io::IO, mime::MIME"text/plain", x::WarmRainParams2M) =
 """
     P3IceParams
 
-Parameters for P3 ice-phase processes (optional).
+Parameters for P3 ice-phase processes.
 
 # Fields
-- `scheme::P3`: ParametersP3 — P3 scheme parameters
-- `terminal_velocity::VL`: Chen2022VelType — terminal velocity for ice
-- `cloud_pdf::PDc`: CloudParticlePDF_SB2006 — cloud droplet size distribution
-- `rain_pdf::PDr`: RainParticlePDF_SB2006 — rain drop size distribution
+$(DocStringExtensions.FIELDS)
+
+# Constructor
+
+The main constructor is
+```
+P3IceParams(toml_dict::CP.ParamDict; is_limited = true)
+```
+which constructs the parameterization with components:
+- `scheme` = [`ParametersP3`](@ref)
+- `terminal_velocity` = [`Chen2022VelType`](@ref)
+- `cloud_pdf` = [`CloudParticlePDF_SB2006`](@ref)
+- `rain_pdf` = [`RainParticlePDF_SB2006`](@ref)
+- `heterogeneous` = [`Frostenberg2023`](@ref)
+- `deposition_condfreeze` = [`MorrisonMilbrandt2014`](@ref)
+
 """
-@kwdef struct P3IceParams{P3, VL, PDc, PDr} <: ParametersType
+@kwdef struct P3IceParams{P3, VL, PDc, PDr, HET, DEP} <: ParametersType
+    "The core P3 scheme parameters"
     scheme::P3
+    "The terminal velocity parameterization"
     terminal_velocity::VL
+    "The cloud droplet size distribution"
     cloud_pdf::PDc
+    "The rain drop size distribution"
     rain_pdf::PDr
+    "The heterogeneous ice nucleation parameters"
+    heterogeneous::HET
+    "The deposition ice nucleation and condensation freezing parameters"
+    deposition_condfreeze::DEP  # TODO: Combine or split these structs?
 end
 Base.show(io::IO, mime::MIME"text/plain", x::P3IceParams) =
     ShowMethods.verbose_show_type_and_fields(io, mime, x)
@@ -51,6 +73,8 @@ P3IceParams(toml_dict::CP.ParamDict; is_limited = true) =
         terminal_velocity = Chen2022VelType(toml_dict),
         cloud_pdf = CloudParticlePDF_SB2006(toml_dict),
         rain_pdf = RainParticlePDF_SB2006(toml_dict; is_limited),
+        heterogeneous = Frostenberg2023(toml_dict),
+        deposition_condfreeze = MorrisonMilbrandt2014(toml_dict),
     )
 
 """