add supporting P3 collisions sub-parameterizations

dry freeze rate & local rime density

Author: haakon-e

docs/src/API.md

@@ -167,10 +167,27 @@ P3Scheme.D_m
 P3Scheme.ice_particle_terminal_velocity
 P3Scheme.ice_terminal_velocity_number_weighted
 P3Scheme.ice_terminal_velocity_mass_weighted
+```
+
+### Processes
+
+#### Heterogeneous ice nucleation
+```@docs
 P3Scheme.het_ice_nucleation
+```
+
+#### Ice melting
+```@docs
 P3Scheme.ice_melt
 ```
 
+#### Collisions with liquid droplets
+Supporting methods:
+```@docs
+P3Scheme.compute_max_freeze_rate
+P3Scheme.compute_local_rime_density
+```
+
 ### Supporting integral methods
 
 ```@docs

src/P3_processes.jl

@@ -99,3 +99,130 @@ function ice_melt(
     dLdt = min(dLdt, L_ice / dt)
     return (; dNdt, dLdt)
 end
+
+"""
+    compute_max_freeze_rate(state, aps, tps, velocity_params, ρₐ, Tₐ)
+
+Returns a function `max_freeze_rate(Dᵢ)` that returns the maximum possible freezing rate [kg/s] 
+    for an ice particle of diameter `Dᵢ` [m]. Evaluates to 0 if T ≥ T_freeze.
+
+# Arguments
+- `state`: [`P3State`](@ref)
+- `aps`: [`CMP.AirProperties`](@ref) 
+- `tps`: `TDP.ThermodynamicsParameters`
+- `velocity_params`: The velocity parameterization, e.g. [`CMP.Chen2022VelType`](@ref)
+- `ρₐ`: air density [kg/m³]
+- `Tₐ`: air temperature [K]
+
+This rate represents the thermodynamic upper limit to collisional freezing, 
+which occurs when the heat transfer from the ice particle to the environment is 
+balanced by the latent heat of fusion.
+
+From Eq (A7) in Musil (1970), [Musil1970](@cite).
+"""
+function compute_max_freeze_rate(state, aps, tps, velocity_params, ρₐ, Tₐ)
+    (; D_vapor, K_therm) = aps
+    cp_l = TDP.cp_l(tps)
+    T_frz = TDP.T_freeze(tps)
+    Lᵥ = TD.latent_heat_vapor(tps, Tₐ)
+    L_f = TD.latent_heat_fusion(tps, Tₐ)
+    Tₛ = T_frz  # the surface of the ice particle is assumed to be at the freezing temperature
+    ΔT = Tₛ - Tₐ  # temperature difference between the surface of the ice particle and the air
+    Δρᵥ_sat = ρₐ * (  # saturation vapor density difference between the surface of the ice particle and the air
+        TD.q_vap_saturation(tps, Tₛ, ρₐ, TD.PhaseNonEquil) - 
+        TD.q_vap_saturation(tps, Tₐ, ρₐ, TD.PhaseNonEquil)
+    )
+    v_term = ice_particle_terminal_velocity(state, velocity_params, ρₐ)
+    F_v = CO.ventilation_factor(state.params.vent, aps, v_term)
+    function max_freeze_rate(Dᵢ)
+        Tₐ ≥ T_frz && return zero(Dᵢ)  # No collisional freezing above the freezing temperature
+        return 2 * (π * Dᵢ) * F_v(Dᵢ) * (K_therm * ΔT + Lᵥ * D_vapor * Δρᵥ_sat) / (L_f - cp_l * ΔT)
+    end
+    return max_freeze_rate
+end
+
+"""
+    compute_local_rime_density(state, velocity_params, ρₐ, T)
+
+Provides a function `ρ′ᵣ(Dᵢ, Dₗ)` that computes the local rime density [kg/m³] 
+    for a given ice particle diameter `Dᵢ` [m] and liquid particle diameter `Dₗ` [m].
+
+From Cober and List (1993).
+
+# Arguments
+- `state`: a [`P3State`](@ref) object
+- `velocity_params`: The velocity parameterization, e.g. [`CMP.Chen2022VelType`](@ref)
+- `ρₐ`: air density [kg/m³]
+- `T`: temperature [K]
+
+# Returns
+A function that computes the local rime density [kg/m³] using the equation:
+
+```math
+ρ′_r = 51 + 114 R_i - 5.5 R_i^2
+```
+where
+```math
+R_i = ( D_{liq} ⋅ |v_{liq} - v_{ice}| ) / ( 2 T_{sfc} )
+```
+and ``T_{sfc}`` is the surface temperature [°C], ``D_{liq}`` is the liquid particle
+diameter [m], ``v_{liq/ice}`` is the particle terminal velocity [m/s].
+So the units of ``R_i`` are [m² s⁻¹ °C⁻¹]. The units of ``ρ′_r`` are [kg/m³].
+
+They assume for simplicity that `T_sfc` equals `T`, the ambient air temperature.
+For real graupel, `T_sfc` is slightly higher than `T` due to latent heat release 
+of freezing liquid particles onto the ice particle. MM13 found little sensitivity 
+to "realistic" increases in `T_sfc`.
+
+Implementation follows Cober and List (1993), Eq. 16 and 17.
+See also the P3 fortran code, microphy_p3.f90, Line 3315-3323,
+which extends the range of the calculation to Rᵢ ≤ 12, the upper limit of which
+then equals the solid bulk ice density, ``ρ^⭒ = 900 kg/m^3``.
+
+Note that MM15 only uses this parameterization for collisions with cloud droplets.
+For rain drops, they use the solid bulk ice density, ``ρ^* = 900 kg/m^3``.
+We do not consider this distinction, and use this parameterization for all liquid particles.
+"""
+function compute_local_rime_density(state, velocity_params, ρₐ, T)
+    (; T_freeze) = state.params
+    T°C = T - T_freeze  # Convert to °C
+
+    v_ice = ice_particle_terminal_velocity(state, velocity_params, ρₐ)
+    v_liq = CO.particle_terminal_velocity(velocity_params.rain, ρₐ)
+    function ρ′ᵣ(Dᵢ, Dₗ)
+        v_term = abs(v_ice(Dᵢ) - v_liq(Dₗ))
+        Rᵢ = (Dₗ * v_term) / (2 * T°C)  # Eq. 16 in Cober and List (1993). Note: no `-` due to absolute value in v_term
+        return ρ′ᵣ_P3(Rᵢ)
+    end
+    return ρ′ᵣ
+end
+
+"""
+    ρ′ᵣ_P3(Rᵢ)
+
+Returns the local rime density [kg/m³] for a given Rᵢ [m² s⁻¹ °C⁻¹].
+
+Based on Cober and List (1993), Eq. 16 and 17 (valid for 1 ≤ Rᵢ ≤ 8).
+For 8 < Rᵢ ≤ 12, linearly interpolate between ρ′ᵣ(8) ≡ 611 kg/m³ and ρ⭒ = 900 kg/m³.
+See also the P3 fortran code
+"""
+function ρ′ᵣ_P3(Rᵢ)
+    # TODO: Externalize these parameters
+    a, b, c = 51, 114, -11//2 # coeffs for Eq. 17 in Cober and List (1993), converted to [kg / m³]
+    ρ⭒ = 900  # ρ^⭒: density of solid bulk ice
+    
+    Rᵢ = clamp(Rᵢ, 1, 12)  # P3 fortran code, microphy_p3.f90, Line 3315 clamps to 1 ≤ Rᵢ ≤ 12
+    
+    ρ′ᵣ_CL93(Rᵢ) = a + b * Rᵢ + c * Rᵢ^2  # Eq. 17 in Cober and List (1993), in [kg / m³], valid for 1 ≤ Rᵢ ≤ 8
+    ρ′ᵣ = if Rᵢ ≤ 8
+        ρ′ᵣ_CL93(Rᵢ)
+    else
+        # following P3 fortran code, microphy_p3.f90, Line 3323
+        #   https://github.com/P3-microphysics/P3-microphysics/blob/main/src/microphy_p3.f90#L3323
+        # for 8 < Rᵢ ≤ 12, linearly interpolate between ρ′ᵣ(8) ≡ 611 kg/m³ and ρ⭒ = 900 kg/m³
+        ρ′ᵣ8 = ρ′ᵣ_CL93(8)
+        f_ρ⭒ = (Rᵢ - 8) / (12 - 8)
+        (1 - f_ρ⭒) * ρ′ᵣ8 + f_ρ⭒ * ρ⭒  # Linear interpolation beyond 8.
+    end
+    return float(ρ′ᵣ)
+end

test/p3_tests.jl

@@ -561,6 +561,58 @@ function test_p3_melting(FT)
     end
 end
 
+function test_p3_bulk_liquid_ice_collisions(FT)
+    params = CMP.ParametersP3(FT)
+    vel_params = CMP.Chen2022VelType(FT)
+    aps = CMP.AirProperties(FT)
+    tps = TD.Parameters.ThermodynamicsParameters(FT)
+
+    (; T_freeze) = params
+
+    ρₐ = FT(1.2)
+    qᵢ = FT(1e-4)
+    Lᵢ = qᵢ * ρₐ
+    Nᵢ = FT(2e5) * ρₐ
+    F_rim = FT(0.8)
+    ρ_rim = FT(800)
+
+    state = P3.get_state(params; F_rim, ρ_rim, L_ice = Lᵢ, N_ice = Nᵢ)
+    logλ = P3.get_distribution_logλ(state)
+    D̄ = exp(-logλ)
+
+    @testset "maximum dry freezing rate" begin
+        # Below freezing, max freeze rate is non-zero (check against reference value)
+        Tₐ = T_freeze - 1//10
+        max_rate = P3.compute_max_freeze_rate(state, aps, tps, vel_params, ρₐ, Tₐ)
+        @test max_rate(D̄) ≈ FT(8.819134467143614e-13) rtol = 1e-4
+
+        # At freezing, max freeze rate is zero
+        Tₐ = T_freeze
+        max_rate = P3.compute_max_freeze_rate(state, aps, tps, vel_params, ρₐ, Tₐ)
+        @test iszero(max_rate(D̄))
+
+        # Above freezing, max freeze rate is zero
+        Tₐ = T_freeze + 1//10
+        max_rate = P3.compute_max_freeze_rate(state, aps, tps, vel_params, ρₐ, Tₐ)
+        @test iszero(max_rate(D̄))
+    end
+
+    @testset "local rime density" begin
+        Tₐ = T_freeze - 1//10
+        ρ′ᵣ = P3.compute_local_rime_density(state, vel_params, ρₐ, Tₐ)
+        @test ρ′ᵣ(D̄, D̄) ≈ FT(159.5) rtol = 1e-6
+
+        a, b, c = 51, 114, -11//2 # coeffs for Eq. 17 in Cober and List (1993), converted to [kg / m³]
+        ρ′ᵣ_CL93(Rᵢ) = a + b * Rᵢ + c * Rᵢ^2  # Eq. 17 in Cober and List (1993), in [kg / m³], valid for 1 ≤ Rᵢ ≤ 8
+        ρ⭒ = FT(900)  # ρ^⭒: density of solid bulk ice
+
+        @test P3.ρ′ᵣ_P3(FT(1)) == ρ′ᵣ_CL93(FT(1))
+        @test P3.ρ′ᵣ_P3(FT(8)) == ρ′ᵣ_CL93(FT(8))
+        @test P3.ρ′ᵣ_P3(FT(12)) == ρ⭒
+    end
+
+end
+
 
 TT.@testset "P3 tests ($FT)" for FT in (Float64, Float32)
     # state creation
@@ -578,5 +630,8 @@ TT.@testset "P3 tests ($FT)" for FT in (Float64, Float32)
     # processes
     test_p3_het_freezing(FT)
     test_p3_melting(FT)
+
+    # bulk liquid-ice collisions and related processes
+    test_p3_bulk_liquid_ice_collisions(FT)
 end
 nothing