P3 collisions

src/P3_particle_properties.jl

@@ -102,7 +102,7 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_rim::FT) where {FT}
     @assert F_rim >= FT(0)   # rime mass fraction must be positive ...
     @assert F_rim < FT(1)    # ... and there must always be some unrimed part
 
-    if F_rim == FT(0)
+    if F_rim == FT(0) || ρ_r == FT(0)
         return (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0))
     else
         @assert ρ_r > FT(0)   # rime density must be positive ...
@@ -304,6 +304,33 @@ function p3_area(
     end
 end
 
+"""
+    K(p3, D1, D1)
+
+ - p3 - a struct with P3 scheme parameters
+ - D_ice - maximum dimension of ice particle
+ - D_other - maximum dimension of second particle (can be cloud water or rain)
+ - F_rim - rime mass fraction (L_rim/ (L_ice - L_liq))
+ - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d)
+
+Computes the geometric collision kernel of an ice particle colliding
+with a spherical cloud water or rain particle.
+We assume a spherical particle of the same
+area if the particle is nonspherical. We also assume collision efficiency
+of 1, so the geometric collision kernel is the total collision kernel.
+"""
+function K(
+    p3::PSP3,
+    D_ice::FT,
+    D_other::FT,
+    F_rim::FT,
+    F_liq::FT,
+    th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)),
+) where {FT}
+    A_ice = p3_area(p3, D_ice, F_rim, F_liq, th)
+    return FT(π) * ((A_ice / FT(π))^0.5 + D_other)^2
+end
+
 """
     ϕᵢ(p3, D, F_rim, th)
 

src/P3_processes.jl

@@ -1,3 +1,24 @@
+"""
+A structure containing the rates of change of the specific humidities and number
+densities of liquid and rain water.
+"""
+Base.@kwdef struct P3Rates{FT}
+    "Rate of change of total ice mass content"
+    dLdt_p3_tot::FT = FT(0)
+    "Rate of change of rime mass content"
+    dLdt_rim::FT = FT(0)
+    "Rate of change of mass content of liquid on ice"
+    dLdt_liq::FT = FT(0)
+    "Rate of change of rime volume"
+    ddtB_rim::FT = FT(0)
+    "Rate of change of ice number concentration"
+    dNdt_ice::FT = FT(0)
+    "Rate of change of rain mass content"
+    dLdt_rai::FT = FT(0)
+    "Rate of change of rain number concentration"
+    dNdt_rai::FT = FT(0)
+end
+
 """
     het_ice_nucleation(pdf_c, p3, tps, q, N, T, ρₐ, p, aerosol)
 
@@ -114,3 +135,206 @@ function ice_melt(
     dLdt = min(dLdt, L_ice / dt)
     return (; dNdt, dLdt)
 end
+
+"""
+    ice_collisions(type, p3, Chen2022, ρ_a, F_r, qᵣ, qᵢ, Nᵣ, Nᵢ, ρ, ρ_r, E_ri)
+
+ - type - defines what is colliding with the ice ("cloud" or "rain")
+ - p3 - a struct with P3 scheme parameters 
+ - Chen2022 - a struct with terminal velocity parameters as in Chen (2022) 
+ - L_p3_tot - total p3 ice mass content [kg/m3]
+ - N_ice - number mixing ratio of ice 
+ - L_c - mass content of colliding species
+ - N_c - number mixing ratio of colliding species
+ - ρ_a - density of air
+ - F_rim - rime mass fraction (L_rim / L_ice)
+ - ρ_rim - rime density (L_rim / B_rim)
+ - F_liq - liquid fraction (L_liq / L_p3_tot)
+ - T - temperature (in K)
+ - E_ci - collision efficiency between ice and colliding species
+
+ Returns the rate of collisions between p3 ice and liquid. 
+ Equivalent to the measure of QRCOL in Morrison and Mildbrandt (2015)
+"""
+function ice_collisions(
+    pdf::Union{
+        CMP.RainParticlePDF_SB2006{FT},
+        CMP.RainParticlePDF_SB2006_limited{FT},
+        CMP.CloudParticlePDF_SB2006{FT},
+    },
+    p3::PSP3,
+    Chen2022::CMP.Chen2022VelType,
+    aps::CMP.AirProperties{FT},
+    tps::TDP.ThermodynamicsParameters{FT},
+    ts::TD.PhaseNonEquil{FT},
+    L_p3_tot::FT,
+    N_ice::FT,
+    L_c::FT,
+    N_c::FT,
+    ρ_a::FT,
+    F_rim::FT,
+    ρ_rim::FT,
+    F_liq::FT,
+    T::FT,
+    dt::FT,
+    E_ci = FT(1),
+) where {FT}
+    # initialize rates
+    rates = P3Rates{FT}()
+    dLdt = rates.dLdt_p3_tot
+    dNdt = rates.dNdt_rai
+    dLdt_liq = rates.dLdt_liq
+    dLdt_rim = rates.dLdt_rim
+    ddt_B_rim = rates.ddtB_rim
+    dLdt_rai = rates.dLdt_rai
+    th = thresholds(p3, ρ_rim, F_rim)
+    (λ, N_0) =
+        distribution_parameter_solver(p3, L_p3_tot, N_ice, ρ_rim, F_rim, F_liq)
+    colliding_bound =
+        CM2.get_size_distribution_bound(pdf, L_c / ρ_a, N_c, ρ_a, FT(1e-6))
+    ice_bound = get_ice_bound(p3, λ, FT(1e-6))
+    N(D) = N′ice(p3, D, λ, N_0)
+    PSD_c(D) = CM2.size_distribution(pdf, D, L_c / ρ_a, ρ_a, N_c)
+    if T > p3.T_freeze
+        # liquid mass collected is a source for L_liq
+        f_warm(D_c, Dᵢ) =
+            E_ci / ρ_a *
+            PSD_c(D_c) *
+            N(Dᵢ) *
+            K(p3, Dᵢ, D_c, F_rim, F_liq, th) *
+            mass_s(D_c, p3.ρ_l) *
+            velocity_difference(
+                pdf,
+                D_c,
+                Dᵢ,
+                p3,
+                Chen2022,
+                ρ_a,
+                F_rim,
+                F_liq,
+                th,
+            )
+        (dLdt, error) = HC.hcubature(
+            d -> f_warm(d[1], d[2]),
+            (FT(0), FT(0)),
+            (colliding_bound, ice_bound),
+        )
+        dLdt = min(dLdt, L_c / dt)
+        dLdt_liq = dLdt
+        dLdt_rai = dLdt
+
+        # calculate sink in rain number proportional to dLdt
+        dNdt = N_c / L_c * dLdt
+        dNdt = min(dNdt, N_c / dt)
+    else
+        # need to calculate dry and wet growth rates 
+        # P3 uses Musil (1970)
+        # but it seems like we already have a dry growth parameterization here
+        # from Anastasia, so maybe I'll try adding the wet growth rate from Musil
+
+        # (see Cholette 2019 p 578) for this:
+        # for dry growth, collected mass is a source for L_rim, B_rim (ρ_rim = 900)
+        # and for wet growth, collected mass is a source for L_liq
+
+        # (see MM 2015 p 307 for this:)
+        # the total rate is assumed to be
+        # dLdt = abs(dLdt_wet - dLdt_dry)
+        # when dLdt_wet > dLdt_dry, dLdt = dLdt_wet - dLdt_dry
+        # and dLdt is transferred to L_liq
+        # when dLdt_dry > dLdt_wet, dLdt = dLdt_dry
+        # and dLdt is transferred to L_rim, B_rim (ρ_rim = 900)
+        f_cold(D_c, Dᵢ) =
+            E_ci / ρ_a *
+            PSD_c(D_c) *
+            N(Dᵢ) *
+            K(p3, Dᵢ, D_c, F_rim, F_liq, th) *
+            mass_s(D_c, p3.ρ_l) *
+            velocity_difference(
+                pdf,
+                D_c,
+                Dᵢ,
+                p3,
+                Chen2022,
+                ρ_a,
+                F_rim,
+                F_liq,
+                th,
+            )
+        (dLdt_dry, error) = HC.hcubature(
+            d -> f_cold(d[1], d[2]),
+            (FT(0), FT(0)),
+            (colliding_bound, ice_bound),
+        )
+
+        # wet growth (Musil 1970 eq A7)
+        # K_term = thermal conductivity of air
+        # HV = latent heat of vaporization
+        # D_vapor = diffusivity of water vapor
+        # Δρ = "sat vapor density at temp of ice minus
+        #       that at temp of cloud air"
+        # HF = latent heat of fusion
+        # C = specific heat of liquid water?
+        (; ν_air, D_vapor, K_therm) = aps
+        HV = TD.latent_heat_vapor(tps, T)
+        Δρ =
+            TD.q_vap_saturation(tps, p3.T_freeze, ρ_a, TD.PhaseNonEquil{FT}) -
+            TD.q_vap_saturation(tps, T, ρ_a, TD.PhaseNonEquil{FT})
+        HF = TD.latent_heat_fusion(tps, T)
+        C = FT(4184)
+        N_sc = ν_air / D_vapor
+        # Ice particle terminal velocity
+        # Reynolds number
+        v(D) = p3_particle_terminal_velocity(
+            p3,
+            D,
+            Chen2022,
+            ρ_a,
+            F_rim,
+            F_liq,
+            th,
+        )
+        N_Re(D) = D * v(D) / ν_air
+        # Ventillation factor
+        a(D) = p3.vent_a + p3.vent_b * N_sc^(1 / 3) * N_Re(D)^(1 / 2)
+        wet_rate(Dᵢ) =
+            2 * FT(π) * Dᵢ * a(Dᵢ) * (-K_therm * T + HV * D_vapor * Δρ) /
+            (HF + C * T)
+
+        # integrate the wet rate over the PSDs
+        # definitely wrong, because the rate should depend on D_c
+        # but it doesn't...
+        f_wet(D_c, Dᵢ) = wet_rate(Dᵢ) * N(Dᵢ) * PSD_c(D_c)
+        (dLdt_wet, error) = HC.hcubature(
+            d -> f_wet(d[1], d[2]),
+            (FT(0), FT(0)),
+            (colliding_bound, ice_bound),
+        )
+
+        if dLdt_wet > dLdt_dry
+            dLdt = dLdt_wet - dLdt_dry
+            dLdt = min(dLdt, L_c / dt)
+            dLdt_liq = dLdt
+            dLdt_rai = dLdt
+            # calculate sink in rain number proportional to dLdt
+            dNdt = N_c / L_c * dLdt
+            dNdt = min(dNdt, N_c / dt)
+        else
+            dLdt = dLdt_dry
+            dLdt = min(dLdt, L_c / dt)
+            dLdt_rai = dLdt
+            dLdt_rim = dLdt
+            ddt_B_rim = dLdt / FT(900)
+            # calculate sink in rain number proportional to dLdt
+            dNdt = N_c / L_c * dLdt
+            dNdt = min(dNdt, N_c / dt)
+        end
+    end
+
+    return P3Rates{FT}(
+        dLdt_p3_tot = dLdt,
+        dLdt_rai = dLdt_rai,
+        dLdt_liq = dLdt_liq,
+        dLdt_rim = dLdt_rim,
+        ddtB_rim = ddt_B_rim,
+    )
+end

src/P3_terminal_velocity.jl

@@ -102,6 +102,7 @@ function velocity_difference(
     Chen2022::CMP.Chen2022VelType,
     ρₐ::FT,
     F_rim::FT,
+    F_liq::FT,
     th,
     use_aspect_ratio = true,
 ) where {FT}