Merge pull request #596 from CliMA/he/feat/p3-collisions-sub-parameterizations

add supporting P3 collisions sub-parameterizations

Author: haakon-e

docs/bibliography.bib

@@ -141,6 +141,19 @@ @article{Choletteetal2019
   pages= {561--582}
 }
 
+
+@article{CoberList1993,
+	title = {Measurements of the Heat and Mass Transfer Parameters Characterizing Conical Graupel Growth},
+	volume = {50},
+	doi = {10.1175/1520-0469(1993)050<1591:MOTHAM>2.0.CO;2},
+	pages = {1591--1609},
+	number = {11},
+	journal = {Journal of the Atmospheric Sciences},
+	author = {Cober, Stewart G. and List, Roland},
+  year = {1993},
+}
+
+
 @article{Desai2019,
   title = {Aerosol-Mediated Glaciation of Mixed-Phase Clouds: Steady-State Laboratory Measurements},
   author = {Desai, Neel and Chandrakar, KK and Kinney, G and Cantrell, W and Shaw, RA},
@@ -549,6 +562,19 @@ @article{Murray2012
 doi  = {10.1039/C2CS35200A}
 }
 
+
+@article{Musil1970,
+	title = {Computer Modeling of Hailstone Growth in Feeder Clouds},
+	volume = {27},
+	doi = {10.1175/1520-0469(1970)027<0474:CMOHGI>2.0.CO;2},
+	pages = {474--482},
+	number = {3},
+	journal = {Journal of the Atmospheric Sciences},
+	author = {Musil, Dennis J.},
+	year = {1970},
+}
+
+
 @article{Ogura1971,
   title = {Numerical simulation of the life cycle of a thunderstorm cell},
   author = {Ogura, Yoshimitsu and Takahashi, Tsutomu},

docs/src/API.md

@@ -109,6 +109,7 @@ CMP.SlopeLaw
 CMP.SlopePowerLaw{Float64}
 CMP.SlopeConstant{Float64}
 CMP.VentilationFactor
+CMP.LocalRimeDensity
 ```
 
 ## Obtain particle state
@@ -167,10 +168,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

docs/src/plots/P3Melting.jl

@@ -33,26 +33,26 @@ Fᵣ1 = FT(0.8)
 ρᵣ1 = FT(800)
 state = P3.get_state(params; F_rim = Fᵣ1, ρ_rim = ρᵣ1, L_ice = Lᵢ, N_ice = Nᵢ)
 logλ = P3.get_distribution_logλ(state)
-dLdt1 = [P3.ice_melt(state, logλ, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt).dLdt for ΔT in ΔT_range]
-dNdt1 = [P3.ice_melt(state, logλ, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt).dNdt for ΔT in ΔT_range]
+dLdt1 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
+dNdt1 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
 
 Fᵣ2 = FT(0.2)
 state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ1, L_ice = Lᵢ, N_ice = Nᵢ)
 logλ = P3.get_distribution_logλ(state)
-dLdt2 = [P3.ice_melt(state, logλ, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt).dLdt for ΔT in ΔT_range]
-dNdt2 = [P3.ice_melt(state, logλ, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt).dNdt for ΔT in ΔT_range]
+dLdt2 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
+dNdt2 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
 
 ρᵣ2 = FT(200)
 state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ2, L_ice = Lᵢ, N_ice = Nᵢ)
 logλ = P3.get_distribution_logλ(state)
-dLdt3 = [P3.ice_melt(state, logλ, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt).dLdt for ΔT in ΔT_range]
-dNdt3 = [P3.ice_melt(state, logλ, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt).dNdt for ΔT in ΔT_range]
+dLdt3 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dLdt for ΔT in ΔT_range]
+dNdt3 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ1, dt, state, logλ).dNdt for ΔT in ΔT_range]
 
 ρₐ2 = FT(0.5)
 state = P3.get_state(params; F_rim = Fᵣ2, ρ_rim = ρᵣ2, L_ice = Lᵢ, N_ice = Nᵢ)
 logλ = P3.get_distribution_logλ(state)
-dLdt4 = [P3.ice_melt(state, logλ, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt).dLdt for ΔT in ΔT_range]
-dNdt4 = [P3.ice_melt(state, logλ, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt).dNdt for ΔT in ΔT_range]
+dLdt4 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt, state, logλ).dLdt for ΔT in ΔT_range]
+dNdt4 = [P3.ice_melt(vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt, state, logλ).dNdt for ΔT in ΔT_range]
 
 # plotting
 fig = Makie.Figure(size = (1500, 500), fontsize=22, linewidth=3)

docs/src/plots/P3TerminalVelocityPlots.jl

@@ -32,8 +32,8 @@ function get_values(
             ρ_rim = ρ_rims[j]
             state = P3.get_state(params; F_rim, ρ_rim, L_ice = L, N_ice = N)
             logλ = P3.get_distribution_logλ(state)
-            V_m[i, j] = P3.ice_terminal_velocity_mass_weighted(state, logλ, Chen2022, ρ_a; use_aspect_ratio = false)
-            V_m_ϕ[i, j] = P3.ice_terminal_velocity_mass_weighted(state, logλ, Chen2022, ρ_a; use_aspect_ratio = true)
+            V_m[i, j] = P3.ice_terminal_velocity_mass_weighted(Chen2022, ρ_a, state, logλ; use_aspect_ratio = false)
+            V_m_ϕ[i, j] = P3.ice_terminal_velocity_mass_weighted(Chen2022, ρ_a, state, logλ; use_aspect_ratio = true)
             D_m[i, j] = P3.D_m(state, logλ)
             D_m_regimes[i, j] = D_m[i, j]
             ϕᵢ[i, j] = P3.ϕᵢ(state, D_m[i, j])

src/Common.jl

@@ -383,14 +383,15 @@ The ventilation factor parameterizes the increase in the mass and heat exchange
 # Returns
 - `F_v(D)`: The ventilation factor as a function of diameter, `D`
 
-See e.g. [SeifertBeheng2006](@cite) Eq. (24) for the definition of the ventilation factor.
+See e.g. [SeifertBeheng2006](@cite) Eq. (24), or [`CMP.VentilationFactor`](@ref),
+for the definition of the ventilation factor.
 """
 function ventilation_factor(vent, aps, v_term)
-    (; vent_a, vent_b) = vent
+    (; aᵥ, bᵥ) = vent
     (; ν_air, D_vapor) = aps
     N_sc = ν_air / D_vapor           # Schmidt number
     N_Re(D) = D * v_term(D) / ν_air  # Reynolds number
-    F_v(D) = vent_a + vent_b * ∛(N_sc) * √(N_Re(D))  # Ventilation factor
+    F_v(D) = aᵥ + bᵥ * ∛(N_sc) * √(N_Re(D))  # Ventilation factor
     return F_v
 end
 

src/P3_processes.jl

@@ -47,27 +47,25 @@ function het_ice_nucleation(
 end
 
 """
-    ice_melt(state, logλ, Chen2022, aps, tps, Tₐ, ρₐ, dt; ∫kwargs...)
+    ice_melt(velocity_params::CMP.Chen2022VelType, aps, tps, Tₐ, ρₐ, dt, state, logλ; ∫kwargs...)
 
 # Arguments
- - `state`: a [`P3State`](@ref) object
- - `logλ`: the log of the slope parameter [log(1/m)]
- - `Chen2022`: struct containing Chen 2022 velocity parameters
- - `aps`: air properties
+ - `velocity_params`: [`CMP.Chen2022VelType`](@ref)
+ - `aps`: [`CMP.AirProperties`](@ref)
  - `tps`: thermodynamics parameters
  - `Tₐ`: temperature (K)
  - `ρₐ`: air density
  - `dt`: model time step (for limiting the tendnecy)
+ - `state`: a [`P3State`](@ref) object
+ - `logλ`: the log of the slope parameter [log(1/m)]
 
 # Keyword arguments
  - `∫kwargs`: Named tuple of keyword arguments passed to [`∫fdD`](@ref)
 
 Returns the melting rate of ice (QIMLT in Morrison and Mildbrandt (2015)).
 """
 function ice_melt(
-    state::P3State, logλ, Chen2022::CMP.Chen2022VelType,
-    aps::CMP.AirProperties, tps::TDI.PS,
-    Tₐ, ρₐ, dt;
+    velocity_params::CMP.Chen2022VelType, aps::CMP.AirProperties, tps::TDI.PS, Tₐ, ρₐ, dt, state::P3State, logλ;
     ∫kwargs = (;),
 )
     # Note: process not dependent on `F_liq`
@@ -79,7 +77,7 @@ function ice_melt(
     (; L_ice, N_ice) = state
     (; T_freeze, vent) = state.params
 
-    v_term = ice_particle_terminal_velocity(state, Chen2022, ρₐ)
+    v_term = ice_particle_terminal_velocity(velocity_params, ρₐ, state)
     F_v = CO.ventilation_factor(vent, aps, v_term)
     N′ = size_distribution(state, logλ)
 
@@ -99,3 +97,105 @@ function ice_melt(
     dLdt = min(dLdt, L_ice / dt)
     return (; dNdt, dLdt)
 end
+
+"""
+    compute_max_freeze_rate(aps, tps, velocity_params, ρₐ, Tₐ, state)
+
+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
+- `aps`: [`CMP.AirProperties`](@ref) 
+- `tps`: `TDP.ThermodynamicsParameters`
+- `velocity_params`: velocity parameterization, e.g. [`CMP.Chen2022VelType`](@ref)
+- `ρₐ`: air density [kg/m³]
+- `Tₐ`: air temperature [K]
+- `state`: [`P3State`](@ref)
+
+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(aps, tps, velocity_params, ρₐ, Tₐ, state)
+    (; D_vapor, K_therm) = aps
+    cp_l = TDI.TD.Parameters.cp_l(tps)
+    T_frz = TDI.TD.Parameters.T_freeze(tps)
+    Lᵥ = TDI.Lᵥ(tps, Tₐ)
+    L_f = TDI.Lf(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
+            TDI.p2q(tps, Tₛ, ρₐ, TDI.saturation_vapor_pressure_over_ice(tps, Tₛ)) -
+            TDI.p2q(tps, Tₐ, ρₐ, TDI.saturation_vapor_pressure_over_ice(tps, Tₐ))
+        )
+    v_term = ice_particle_terminal_velocity(velocity_params, ρₐ, state)
+    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(velocity_params, ρₐ, T, state)
+
+Provides a function `ρ′_rim(Dᵢ, Dₗ)` that computes the local rime density [kg/m³] 
+    for a given ice particle diameter `Dᵢ` [m] and liquid particle diameter `Dₗ` [m].
+
+# Arguments
+- `velocity_params`: velocity parameterization, e.g. [`CMP.Chen2022VelType`](@ref)
+- `ρₐ`: air density [kg/m³]
+- `T`: temperature [K]
+- `state`: [`P3State`](@ref)
+
+# Returns
+A function that computes the local rime density [kg/m³] using the equation:
+
+```math
+ρ'_{rim} = a + b R_i + c R_i^2
+```
+where
+```math
+R_i = \\frac{ 10^6 ⋅ 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 ``ρ'_{rim}`` are [kg/m³].
+
+We 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. Morrison & Milbrandt (2013) 
+found little sensitivity to "realistic" increases in ``T_{sfc}``.
+
+See also [`LocalRimeDensity`](@ref CloudMicrophysics.Parameters.LocalRimeDensity).
+
+# Extended help
+
+ 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_i ≤ 12``, the upper limit of which
+ then equals the solid bulk ice density, ``ρ_ice = 916.7 kg/m^3``.
+
+ Note that Morrison & Milbrandt (2015) [MorrisonMilbrandt2015](@cite) only uses this 
+ parameterization for collisions with cloud droplets.
+ For rain drops, they use a value near 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(velocity_params, ρₐ, T, state)
+    (; T_freeze, ρ_rim_local) = state.params
+    T°C = T - T_freeze  # Convert to °C
+    μm = 1_000_000  # Note: m to μm factor, c.f. units of rₘ in Eq. 16 in Cober and List (1993)
+
+    v_ice = ice_particle_terminal_velocity(velocity_params, ρₐ, state)
+    v_liq = CO.particle_terminal_velocity(velocity_params.rain, ρₐ)
+    function ρ′_rim(Dᵢ, Dₗ)
+        v_term = abs(v_ice(Dᵢ) - v_liq(Dₗ))
+        Rᵢ = (Dₗ * μm * v_term) / (2 * T°C)  # Eq. 16 in Cober and List (1993). Note: no `-` due to absolute value in v_term
+        return ρ_rim_local(Rᵢ)
+    end
+    return ρ′_rim
+end

src/P3_terminal_velocity.jl

@@ -1,6 +1,7 @@
 
 """
-    ice_particle_terminal_velocity(state, velocity_params, ρₐ; [use_aspect_ratio])
+    ice_particle_terminal_velocity(velocity_params, ρₐ, state::P3State; [use_aspect_ratio])
+    ice_particle_terminal_velocity(velocity_params, ρₐ, params::CMP.ParametersP3, F_rim, ρ_rim; [use_aspect_ratio])
 
 Returns the terminal velocity of a single ice particle as a function of its size 
     (maximum dimension, `D`) using the Chen 2022 parametrization.
@@ -15,31 +16,34 @@ Returns the terminal velocity of a single ice particle as a function of its size
     of particle aspect ratio on its terminal velocity (default: `true`)
 """
 function ice_particle_terminal_velocity(
-    state::P3State, velocity_params::CMP.Chen2022VelType, ρₐ; use_aspect_ratio = true,
+    velocity_params::CMP.Chen2022VelType, ρₐ, state_components...; use_aspect_ratio = true,
 )
     function v_term(D::FT) where {FT}
         (; small_ice, large_ice) = velocity_params
         D_cutoff = small_ice.cutoff # Diameter cutoff between small and large ice regimes
         ρᵢ = FT(916.7) # ρᵢ = p3_density(p3, D, F_rim, th) # TODO: tmp
         v_term_small = CO.particle_terminal_velocity(small_ice, ρₐ, ρᵢ)
         v_term_large = CO.particle_terminal_velocity(large_ice, ρₐ, ρᵢ)
-        ϕ_factor = use_aspect_ratio ? cbrt(ϕᵢ(state, D)) : FT(1)
+        # `state_components...` is either `state` or `(params, F_rim, ρ_rim)`
+        ϕ_factor = use_aspect_ratio ? cbrt(ϕᵢ(state_components..., D)) : FT(1)
         vₜ = D <= D_cutoff ? v_term_small(D) : v_term_large(D)
         return ϕ_factor * vₜ
     end
     return v_term
 end
-
 """
-    ice_terminal_velocity_number_weighted(state, logλ velocity_params, ρₐ; [use_aspect_ratio], [∫kwargs...])
+    ice_terminal_velocity_number_weighted(
+        velocity_params::CMP.Chen2022VelType, ρₐ, state::P3State, logλ;
+        [use_aspect_ratio], [∫kwargs...],
+    )
 
 Return the terminal velocity of the number-weighted mean ice particle size.
 
 # Arguments
-- `state`: A [`P3State`](@ref)
-- `logλ`: The log of the slope parameter [log(1/m)]
 - `velocity_params`: A [`CMP.Chen2022VelType`](@ref) with terminal velocity parameters
 - `ρₐ`: Air density [kg/m³]
+- `state`: A [`P3State`](@ref)
+- `logλ`: The log of the slope parameter [log(1/m)]
 
 # Keyword arguments
  - `use_aspect_ratio`: Bool flag set to `true` if we want to consider the effects
@@ -49,15 +53,15 @@ Return the terminal velocity of the number-weighted mean ice particle size.
 See also [`ice_terminal_velocity_mass_weighted`](@ref)
 """
 function ice_terminal_velocity_number_weighted(
-    state::P3State, logλ, velocity_params::CMP.Chen2022VelType, ρₐ;
+    velocity_params::CMP.Chen2022VelType, ρₐ, state::P3State, logλ;
     use_aspect_ratio = true, ∫kwargs...,
 )
     (; N_ice, L_ice) = state
     if N_ice < eps(one(N_ice)) || L_ice < eps(one(L_ice))
         return zero(N_ice)
     end
 
-    v_term = ice_particle_terminal_velocity(state, velocity_params, ρₐ; use_aspect_ratio)
+    v_term = ice_particle_terminal_velocity(velocity_params, ρₐ, state; use_aspect_ratio)
     n = DT.size_distribution(state, logλ)
 
     # ∫n(D) v(D) dD
@@ -67,15 +71,15 @@ function ice_terminal_velocity_number_weighted(
 end
 
 """
-    ice_terminal_velocity_mass_weighted(state, logλ, velocity_params, ρₐ; [use_aspect_ratio], [∫kwargs...])
+    ice_terminal_velocity_mass_weighted(velocity_params::CMP.Chen2022VelType, ρₐ, state::P3State, logλ; [use_aspect_ratio], [∫kwargs...])
 
 Return the terminal velocity of the mass-weighted mean ice particle size.
 
 # Arguments
-- `state`: A [`P3State`](@ref)
-- `logλ`: The log of the slope parameter [log(1/m)]
 - `velocity_params`: A [`CMP.Chen2022VelType`](@ref) with terminal velocity parameters
 - `ρₐ`: Air density [kg/m³]
+- `state`: A [`P3State`](@ref)
+- `logλ`: The log of the slope parameter [log(1/m)]
 
 # Keyword arguments
  - `use_aspect_ratio`: Bool flag set to `true` if we want to consider the effects
@@ -85,15 +89,15 @@ Return the terminal velocity of the mass-weighted mean ice particle size.
 See also [`ice_terminal_velocity_number_weighted`](@ref)
 """
 function ice_terminal_velocity_mass_weighted(
-    state::P3State, logλ, velocity_params::CMP.Chen2022VelType, ρₐ;
+    velocity_params::CMP.Chen2022VelType, ρₐ, state::P3State, logλ;
     use_aspect_ratio = true, ∫kwargs...,
 )
     (; N_ice, L_ice) = state
     if N_ice < eps(one(N_ice)) || L_ice < eps(one(L_ice))
         return zero(L_ice)
     end
 
-    v_term = ice_particle_terminal_velocity(state, velocity_params, ρₐ; use_aspect_ratio)
+    v_term = ice_particle_terminal_velocity(velocity_params, ρₐ, state; use_aspect_ratio)
     n = DT.size_distribution(state, logλ)  # Number concentration at diameter D
 
     # ∫n(D) m(D) v(D) dD