P3 ice nucleation parameterizations

docs/bibliography.bib

@@ -587,6 +587,17 @@ @article{Spichtinger2023
   DOI = {10.5194/acp-23-2035-2023}
 }
 
+@article{Thompson2004,
+  author = {Gregory Thompson and Roy M. Rasmussen and Kevin Manning},
+  title = {Explicit Forecasts of Winter Precipitation Using an Improved Bulk Microphysics Scheme. Part I: Description and Sensitivity Analysis},
+  journal = {Monthly Weather Review},
+  year = {2004},
+  volume = {132},
+  number = {2},
+  doi = {10.1175/1520-0493(2004)132<0519:EFOWPU>2.0.CO;2},
+  pages = {519 - 542},
+}
+
 @article{TripoliCotton1980,
   title = {A Numerical Investigation of Several Factors Contributing to the Observed Variable Intensity of Deep Convection over South Florida},
   author = {Tripoli, G.J. and Cotton, W.R.},

docs/src/API.md

@@ -89,6 +89,8 @@ HetIceNucleation
 HetIceNucleation.dust_activated_number_fraction
 HetIceNucleation.deposition_J
 HetIceNucleation.ABIFM_J
+HetIceNucleation.P3_deposition_N_i
+HetIceNucleation.P3_het_N_i
 ```
 
 # Homogeneous ice nucleation

docs/src/IceNucleation.md

@@ -21,7 +21,8 @@ The parameterization for deposition on dust particles is an implementation of
     parameterization and modeling the competition between
     freezing modes.
 
-## Activated fraction for deposition freezing on dust
+## Deposition Nucleation
+### Activated fraction for deposition freezing on dust
 The parameterization models the activated fraction
   as an empirical function of ice saturation ratio,
   see eq. (3) in [Mohler2006](@cite).
@@ -47,7 +48,7 @@ Both parameters are dependent on aerosol properties and temperature.
     freezing occurs in a different ice nucleation mode
     (either a second deposition or other condensation type mode).
 
-## Water Activity Based Deposition Nucleation
+### Water Activity Based Deposition Nucleation
 The water activity based deposition nucleation model is analagous to ABIFM
   for immersion freezing (see [ABIFM for Sulphuric Acid Containing Droplets](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/#ABIFM-for-Sulphuric-Acid-Containing-Droplets)
   section below). It calculates a nucleation rate coefficient, ``J``, which
@@ -66,7 +67,22 @@ where ``J`` is in units of ``cm^{-2}s^{-1}``. Note that our source code returns
   ``J`` in SI units. ``m`` and ``c`` are aerosol dependent coefficients. They will
   have different values than those for ABIFM.
 
-## ABIFM for Sulphuric Acid Containing Droplets
+### Cooper (1986) Ice Crystal Number
+The P3 scheme as described in [MorrisonMilbrandt2015](@cite) follows
+  the implementation of [Thompson2004](@cite) where a parameterization from
+  Cooper 1986 was used. Number of ice crystals formed is derived from
+  ambient temperature.
+```math
+\begin{equation}
+  N_i = 0.005 exp[0.304(T_0 - T)]
+\end{equation}
+```
+Where ``T_0 = 273.15K`` and ``T`` is the ambient temperature.
+Because the parameterization is prone to overpredict at low temperatures,
+  ``N_i = N_i(T = 233K)`` for all values of ``N_i`` at which ``T < 233K``.
+
+## Immersion Freezing
+### ABIFM for Sulphuric Acid Containing Droplets
 Water Activity-Based Immersion Freezing Model (ABFIM)
   is a method of parameterizing immersion freezing inspired by the time-dependent
   classical nucleation theory (CNT). More on CNT can be found in [Karthika2016](@cite).
@@ -141,7 +157,8 @@ It is also important to note that this plot is reflective of cirrus clouds
   and shows only a very small temperature range. The two curves are slightly
   off because of small differences in parameterizations for vapor pressures.
 
-## Homogeneous Freezing for Sulphuric Acid Containing Droplets
+## Homogeneous Freezing
+### Homogeneous Freezing for Sulphuric Acid Containing Droplets
 Homogeneous freezing occurs when supercooled liquid droplets freeze on their own.
   Closly based off [Koop2000](@cite), this parameterization determines a homoegneous nucleation
   rate coefficient, ``J_{hom}``, using water activity. The change in water activity,
@@ -157,7 +174,7 @@ The nucleation rate coefficient is determined with the cubic function from [Koop
 ```
 This parameterization is valid only when ``0.26 < \Delta a_w < 0.36`` and ``185K < T < 235K``.
 
-### Homogeneous Ice Nucleation Rate Coefficient
+### Homogeneous Freezing Example Figures
 Here is a comparison of our parameterization of ``J_{hom}`` compared to Koop 2000 as
   plotted in figure 1 of [Spichtinger2023](@cite). Our parameterization differs in the calculation
   of ``\Delta a_w``. We define water activity to be a ratio of saturated vapor pressures whereas

src/IceNucleation.jl

@@ -9,6 +9,8 @@ import Thermodynamics as TD
 export dust_activated_number_fraction
 export deposition_J
 export ABIFM_J
+export P3_deposition_N_i
+export P3_het_N_i
 
 """
     dust_activated_number_fraction(dust, ip, Si, T)
@@ -83,6 +85,49 @@ function ABIFM_J(
     return max(FT(0), FT(10)^logJ * FT(1e4)) # converts cm^-2 s^-1 to m^-2 s^-1
 end
 
+"""
+    P3_deposition_N_i(T)
+
+ - `T` - air temperature [K].
+
+Returns the number of ice nucleated via deposition nucleation with units of m^-3.
+From Thompson et al 2004 eqn 2 as used in Morrison & Milbrandt 2015.
+"""
+function P3_deposition_N_i(T::FT) where {FT}
+
+    T₀ = FT(273.15)
+    T_thres = FT(233)
+
+    Nᵢ = ifelse(
+        T < T_thres,
+        max(FT(0), FT(0.005 * exp(0.304 * (T₀ - T_thres)))),
+        max(FT(0), FT(0.005 * exp(0.304 * (T₀ - T)))),
+    )
+    return Nᵢ * 1e3  # converts L^-1 to m^-3
+end
+
+"""
+    P3_het_N_i(T)
+
+ - `T` - air temperature [K],
+ - `N_l` - number of droplets [m^-3],
+ - `B` - water-type dependent parameter,
+ - `V_l` - volume of droplets to be frozen,
+ - `a` - empirical parameter,
+ - `t` - time since start of simulation.
+
+Returns the number of ice nucleated via heterogeneous freezing with units of m^-3.
+From Pruppacher & Klett 1997 eqn (9-51) as used in Morrison & Milbrandt 2015.
+"""
+function P3_het_N_i(T::FT, N_l::FT, V_l::FT, t::FT) where {FT}
+
+    a = 0.65    # (celcius)^-1
+    B = 2e-4    # cm^-3 s^-1 for rain water
+    Tₛ = FT(273.15) - T
+
+    return N_l * exp(B * V_l * t * exp(a * Tₛ))
+end
+
 end # end module
 
 """

test/gpu_tests.jl

@@ -30,8 +30,8 @@ CUDA.allowscalar(false)
 const ArrayType = CuArray
 
 # For debugging on the CPU
-#const backend = CPU()
-#const ArrayType = Array
+# const backend = CPU()
+# const ArrayType = Array
 
 @info "GPU Tests" backend ArrayType
 
@@ -512,6 +512,33 @@ end
     end
 end
 
+@kernel function IceNucleation_P3_deposition_N_i_kernel!(
+    output::AbstractArray{FT},
+    T,
+) where {FT}
+
+    i = @index(Group, Linear)
+
+    @inbounds begin
+        output[1] = CMI_het.P3_deposition_N_i(T[1])
+    end
+end
+
+@kernel function IceNucleation_P3_het_N_i_kernel!(
+    output::AbstractArray{FT},
+    T,
+    N_l,
+    V_l,
+    t,
+) where {FT}
+
+    i = @index(Group, Linear)
+
+    @inbounds begin
+        output[1] = CMI_het.P3_het_N_i(T[1], N_l[1], V_l[1], t[1])
+    end
+end
+
 @kernel function IceNucleation_homogeneous_J_kernel!(
     ip,
     output::AbstractArray{FT},
@@ -868,6 +895,31 @@ function test_gpu(FT)
         dims = (1, 1)
         (; output, ndrange) = setup_output(dims, FT)
 
+        T = ArrayType([FT(240)])
+
+        kernel! = IceNucleation_P3_deposition_N_i_kernel!(backend, work_groups)
+        kernel!(output, T; ndrange)
+
+        # test if P3_deposition_N_i is callable and returns reasonable values
+        @test Array(output)[1] ≈ FT(119018.93920746)
+
+        dims = (1, 1)
+        (; output, ndrange) = setup_output(dims, FT)
+
+        T = ArrayType([FT(240)])
+        N_l = ArrayType([FT(2e4)])
+        V_l = ArrayType([FT(3e-6)])
+        t = ArrayType([FT(25)])
+
+        kernel! = IceNucleation_P3_het_N_i_kernel!(backend, work_groups)
+        kernel!(output, T, N_l, V_l, t; ndrange)
+
+        # test if P3_het_N_i is callable and returns reasonable values
+        @test Array(output)[1] ≈ FT(1.42814881631e19)
+
+        dims = (1, 1)
+        (; output, ndrange) = setup_output(dims, FT)
+
         T = ArrayType([FT(220)])
         x_sulph = ArrayType([FT(0.15)])
         Delta_a_w = ArrayType([FT(0.2907389666103033)])

test/heterogeneous_ice_nucleation_tests.jl

@@ -116,6 +116,22 @@ function test_heterogeneous_ice_nucleation(FT)
         end
     end
 
+    TT.@testset "P3 Deposition Nᵢ" begin
+
+        T_warm = FT(235)
+        T_cold = FT(234)
+
+        T_too_cold = FT(232)
+
+        # higher ice concentration at colder temperatures
+        TT.@test CMI_het.P3_deposition_N_i(T_cold) >
+                 CMI_het.P3_deposition_N_i(T_warm)
+
+        # if colder than threshold T, use threshold T
+        TT.@test CMI_het.P3_deposition_N_i(T_too_cold) ==
+                 CMI_het.P3_deposition_N_i(FT(233))
+    end
+
     TT.@testset "ABIFM J" begin
 
         T_warm_1 = FT(229.2)
@@ -148,6 +164,20 @@ function test_heterogeneous_ice_nucleation(FT)
             )
         end
     end
+
+    TT.@testset "P3 Heterogeneous Nᵢ" begin
+
+        T_warm = FT(235)
+        T_cold = FT(234)
+        N_liq = FT(2e5)
+        r_l = FT(2e-5)
+        V_l = FT(4 / 3 * FT(π) * r_l^3)
+        t = FT(25)
+
+        # higher ice concentration at colder temperatures
+        TT.@test CMI_het.P3_het_N_i(T_cold, N_liq, V_l, t) >
+                 CMI_het.P3_het_N_i(T_warm, N_liq, V_l, t)
+    end
 end
 
 println("Testing Float64")

test/performance_tests.jl

@@ -66,7 +66,7 @@ function benchmark_test(FT)
     mixed_nuc = CMP.MixedNucleationParameters(FT)
     h2so4_nuc = CMP.H2S04NucleationParameters(FT)
     organ_nuc = CMP.OrganicNucleationParameters(FT)
-    # air and thermodunamics parameters
+    # air and thermodynamics parameters
     aps = CMP.AirProperties(FT)
     tps = TD.Parameters.ThermodynamicsParameters(FT)
 
@@ -110,6 +110,9 @@ function benchmark_test(FT)
 
     x_sulph = FT(0.1)
     Delta_a_w = FT(0.27)
+    r_liq = FT(1e-6)
+    V_liq = FT(4 / 3 * π * r_liq^3)
+    time = FT(25)
 
     # P3 scheme
     bench_press(P3.thresholds, (p3, ρ_r, F_r), 12e6, 2048, 80)
@@ -138,7 +141,9 @@ function benchmark_test(FT)
         50,
     )
     bench_press(CMI_het.deposition_J, (kaolinite, Delta_a_w), 230)
+    bench_press(CMI_het.P3_deposition_N_i, (T_air_cold), 230)
     bench_press(CMI_het.ABIFM_J, (desert_dust, Delta_a_w), 230)
+    bench_press(CMI_het.P3_het_N_i, (T_air_cold, N_liq, V_liq, time), 230)
     bench_press(CMI_hom.homogeneous_J, (ip.homogeneous, Delta_a_w), 230)
 
     # non-equilibrium