Merge branch 'main' of https://github.com/CliMA/CloudMicrophysics.jl …

…into al/p3_ice_nuc

docs/bibliography.bib

@@ -143,6 +143,18 @@ @article{Dunne2016
   doi = {10.1126/science.aaf2649}
 }
 
+@Article{Frostenberg2023,
+author = {Frostenberg, H. C. and Welti, A. and Luhr, M. and Savre, J. and Thomson, E. S. and Ickes, L.},
+title = {The chance of freezing -- a conceptional study to parameterize temperature-dependent freezing by including randomness of ice-nucleating particle concentrations},
+journal = {Atmospheric Chemistry and Physics},
+volume = {23},
+year = {2023},
+number = {19},
+pages = {10883--10900},
+url = {https://acp.copernicus.org/articles/23/10883/2023/},
+doi = {10.5194/acp-23-10883-2023}
+}
+
 @article{Glassmeier2016,
   author = {Franziska Glassmeier and Ulrike Lohmann},
   title = {Constraining Precipitation Susceptibility of Warm-, Ice-, and Mixed-Phase Clouds with Microphysical Equations},

docs/src/API.md

@@ -92,6 +92,7 @@ HetIceNucleation.deposition_J
 HetIceNucleation.ABIFM_J
 HetIceNucleation.P3_deposition_N_i
 HetIceNucleation.P3_het_N_i
+HetIceNucleation.INP_concentration_frequency
 ```
 
 # Homogeneous ice nucleation
@@ -138,6 +139,7 @@ Parameters.Seasalt
 Parameters.Sulfate
 Parameters.AerosolActivationParameters
 Parameters.IceNucleationParameters
+Parameters.Frostenberg2023
 Parameters.H2SO4SolutionParameters
 Parameters.Mohler2006
 Parameters.Koop2000

docs/src/IceNucleation.md

@@ -6,14 +6,17 @@ The `IceNucleation.jl` module includes:
   - the parameterization of activation of dust aerosol particles into ice crystals
     via deposition of water vapor,
   - water activity based parameterization of immersion freezing,
-  - water activity based parameterization of homogeneous freezing.
+  - water activity based parameterization of homogeneous freezing,
+  - parametrization of temperature-dependent aerosol-independent ice nucleating particles concentration for immersion freezing.
+
 The parameterization for deposition on dust particles is an implementation of
   the empirical formulae from [Mohler2006](@cite)
   and is valid for two types of dust particles:
   Arizona Test Dust and desert dust from Sahara.
   The parameterization for immersion freezing is an implementation of [KnopfAlpert2013](@cite)
   and is valid for droplets containing sulphuric acid.
   The parameterization for homogeneous freezing is an implementation of [Koop2000](@cite).
+  The parametrization for ice nucleating particles concentration for immersion freezing is an implementation of [Frostenberg2023](@cite).
 
 !!! note
 
@@ -239,3 +242,37 @@ Multiple sulphuric acid concentrations, ``x``,
     The current parameterization in CloudMicrophysics.jl is not valid for \Delta a_w
     values that are obtained from pure water droplets. Though CliMA lines look far
     from the Spichtinger 2023 line, the lines seem to move closer as x approaches 0.
+
+
+## INP Concentration Frequency
+
+With this parametrization, the concentration of ice nucleating particles (INPs) is found based on the relative frequency distribution, which depends on the temperature, but does not depend on aerosol types. It is based on [Frostenberg2023](@cite) and is derived from measurements, for marine
+data sets. It is a lognormal distribution, described by
+
+```math
+\begin{equation}
+  D(\mu, \sigma^2) = \frac{1}{\sqrt{2 \pi} \sigma} exp \left(- \frac{[ln(a \cdot INPC) - \mu(T)]^2}{2 \sigma^2} \right),
+\end{equation}
+```
+where ``T`` is the temperature in degrees Celsius, ``INPC`` is the INP concentration in m``^{-3}`` and the temperature-dependent mean ``\mu(T)`` is given by
+```math
+\begin{equation}
+  \mu(T) = ln(-(b \cdot T)^9 \times 10^{-9})
+\end{equation}
+```
+ ``\sigma^2`` is the variance, ``a`` and ``b`` are coefficients. The parameters defined in [Frostenberg2023](@cite) for marine data sets are ``\sigma=1.37``, ``a=1`` m``^3``, ``b=1``/C.
+
+!!! note
+
+    Our implementation uses base SI units and takes ``T`` in Kelvin.
+
+The following plot shows the relative frequency distribution for INPCs, as a function of temperature (the same as figure 1 in [Frostenberg2023](@cite)). 
+
+```@example
+include("plots/Frostenberg_fig1.jl")
+```
+![](Frostenberg_fig1.svg)
+
+The following plot shows the relative frequency distribution for INPCs at the temperature ``T=-16 C``.
+
+![](Frostenberg_fig1_T16.svg)

docs/src/plots/Frostenberg_fig1.jl

@@ -0,0 +1,72 @@
+import Plots as PL
+
+import CloudMicrophysics as CM
+import CloudMicrophysics.HetIceNucleation as IN
+import CLIMAParameters as CP
+import CloudMicrophysics.Parameters as CMP
+
+FT = Float32
+ip = CMP.Frostenberg2023(FT)
+
+T_range_celsius = range(-40, stop = -2, length = 500)  # air temperature [°C]
+T_range_kelvin = T_range_celsius .+ 273.15  # air temperature [K]
+INPC_range = 10.0 .^ (range(-5, stop = 7, length = 500)) #ice nucleating particle concentration
+
+frequency = [
+    IN.INP_concentration_frequency(ip, INPC, T) > 0.015 ?
+    IN.INP_concentration_frequency(ip, INPC, T) : missing for
+    INPC in INPC_range, T in T_range_kelvin
+]
+mu = [exp(log(-(ip.b * T)^9 * 10^(-9))) for T in T_range_celsius]
+
+
+PL.contourf(
+    T_range_kelvin,
+    INPC_range,
+    frequency,
+    xlabel = "T [K]",
+    ylabel = "INPC [m⁻³]",
+    colorbar_title = "frequency",
+    yaxis = :log,
+    color = :lighttest,
+    gridlinewidth = 3,
+    ylims = (1e-5, 1e7),
+    lw = 0,
+)
+
+PL.plot!(
+    T_range_kelvin,
+    mu,
+    label = "median INPC",
+    legend = :topright,
+    color = :darkred,
+)
+
+PL.plot!(
+    repeat([257], 50),
+    10.0 .^ (range(-2, stop = 6, length = 50)),
+    label = "T = 257 K",
+    color = :darkgreen,
+    linestyle = :dash,
+)
+PL.savefig("Frostenberg_fig1.svg")
+
+
+#plotting the distribution for T=257 K
+
+T = FT(257.0) # [K]
+INPC_range = FT(10) .^ range(FT(-1), 4, length = 100)
+frequency = [IN.INP_concentration_frequency(ip, INPC, T) for INPC in INPC_range]
+
+PL.plot(
+    INPC_range,
+    frequency,
+    xlabel = "INPC [m⁻³]",
+    ylabel = "frequency",
+    xaxis = :log,
+    color = :darkgreen,
+    label = "T = 257 K",
+    legend = :topright,
+)
+
+PL.savefig("Frostenberg_fig1_T16.svg")

src/IceNucleation.jl

@@ -12,6 +12,8 @@ export deposition_J
 export ABIFM_J
 export P3_deposition_N_i
 export P3_het_N_i
+export INP_concentration_frequency
+
 
 """
     dust_activated_number_fraction(dust, ip, Si, T)
@@ -155,6 +157,31 @@ function P3_het_N_i(T::FT, N_l::FT, V_l::FT, Δt::FT) where {FT}
     return N_l * (1 - exp(-B * V_l_converted * Δt * exp(a * Tₛ)))
 end
 
+"""
+    INP_concentration_frequency(INPC,T)
+
+ - `params` - a struct with INPC(T) distribution parameters
+ - `INPC` - concentration of ice nucleating particles [m^-3]
+ - `T` - air temperature [K]
+
+Returns the relative frequency of a given INP concentration,
+depending on the temperature.
+Based on Frostenberg et al., 2023. See DOI: 10.5194/acp-23-10883-2023
+"""
+function INP_concentration_frequency(
+    params::CMP.Frostenberg2023,
+    INPC::FT,
+    T::FT,
+) where {FT}
+
+    T_celsius = T - 273.15
+    (; σ, a, b) = params
+
+    μ = log(-(b * T_celsius)^9 * 10^(-9))
+
+    return 1 / (sqrt(2 * FT(π)) * σ) * exp(-(log(a * INPC) - μ)^2 / (2 * σ^2))
+end
+
 end # end module
 
 """

src/parameters/IceNucleation.jl

@@ -1,4 +1,5 @@
 export IceNucleationParameters
+export Frostenberg2023
 
 """
     Mohler2006{FT}
@@ -87,3 +88,36 @@ function IceNucleationParameters(
     HOM = typeof(homogeneous)
     return IceNucleationParameters{FT, DEP, HOM}(deposition, homogeneous)
 end
+
+
+"""
+    Frostenberg2023{FT}
+
+Parameters for frequency distribution of INP concentration
+DOI: 10.5194/acp-23-10883-2023
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+Base.@kwdef struct Frostenberg2023{FT} <: ParametersType{FT}
+    "standard deviation"
+    σ::FT
+    "coefficient"
+    a::FT
+    "coefficient"
+    b::FT
+end
+
+Frostenberg2023(::Type{FT}) where {FT <: AbstractFloat} =
+    Frostenberg2023(CP.create_toml_dict(FT))
+
+function Frostenberg2023(td::CP.AbstractTOMLDict)
+    name_map = (;
+        :Frostenberg2023_standard_deviation => :σ,
+        :Frostenberg2023_a_coefficient => :a,
+        :Frostenberg2023_b_coefficient => :b,
+    )
+    parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
+    FT = CP.float_type(td)
+    return Frostenberg2023{FT}(; parameters...)
+end

test/gpu_tests.jl

@@ -595,6 +595,20 @@ end
     end
 end
 
+@kernel function IceNucleation_INPC_frequency_kernel!(
+    output::AbstractArray{FT},
+    ip_frostenberg,
+    INPC,
+    T,
+) where {FT}
+
+    i = @index(Group, Linear)
+
+    @inbounds begin
+        output[1] = CMI_het.INP_concentration_frequency(ip_frostenberg, INPC, T)
+    end
+end
+
 @kernel function IceNucleation_homogeneous_J_kernel!(
     ip,
     output::AbstractArray{FT},
@@ -681,6 +695,7 @@ function test_gpu(FT)
     feldspar = CMP.Feldspar(FT)
     ferrihydrite = CMP.Ferrihydrite(FT)
     ip = CMP.IceNucleationParameters(FT)
+    ip_frostenberg = CMP.Frostenberg2023(FT)
 
     @testset "Aerosol activation kernels" begin
         dims = (6, 2)
@@ -1026,6 +1041,18 @@ function test_gpu(FT)
         dims = (1, 1)
         (; output, ndrange) = setup_output(dims, FT)
 
+        T = FT(233)
+        INPC = FT(220000)
+
+        kernel! = IceNucleation_INPC_frequency_kernel!(backend, work_groups)
+        kernel!(output, ip_frostenberg, INPC, T; ndrange)
+
+        # test INPC_frequency is callable and returns a reasonable value
+        @test Array(output)[1] ≈ FT(0.26) rtol = 0.1
+
+        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

@@ -16,6 +16,7 @@ function test_heterogeneous_ice_nucleation(FT)
     tps = TD.Parameters.ThermodynamicsParameters(FT)
     H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
     ip = CMP.IceNucleationParameters(FT)
+    ip_frostenberg = CMP.Frostenberg2023(FT)
     # more parameters for aerosol properties
     ATD = CMP.ArizonaTestDust(FT)
     desert_dust = CMP.DesertDust(FT)
@@ -207,6 +208,21 @@ function test_heterogeneous_ice_nucleation(FT)
         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
+
+    TT.@testset "Frostenberg" begin
+
+        temperatures = FT.([233, 257])
+        INPCs = FT.([220000, 9])
+        frequencies = FT.([0.26, 0.08])
+
+        for (T, INPC, frequency) in zip(temperatures, INPCs, frequencies)
+            TT.@test CMI_het.INP_concentration_frequency(
+                ip_frostenberg,
+                INPC,
+                T,
+            ) ≈ frequency rtol = 0.1
+        end
+    end
 end
 
 println("Testing Float64")

test/performance_tests.jl

@@ -63,6 +63,7 @@ function benchmark_test(FT)
     kaolinite = CMP.Kaolinite(FT)
     ip = CMP.IceNucleationParameters(FT)
     H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
+    ip_frostenberg = CMP.Frostenberg2023(FT)
     # aerosol nucleation parameters
     mixed_nuc = CMP.MixedNucleationParameters(FT)
     h2so4_nuc = CMP.H2S04NucleationParameters(FT)
@@ -116,6 +117,8 @@ function benchmark_test(FT)
     V_liq = FT(4 / 3 * π * r_liq^3)
     Δt = FT(25)
 
+    INPC = FT(1e5)
+
     # P3 scheme
     bench_press(P3.thresholds, (p3, ρ_r, F_r), 12e6, 2048, 80)
 
@@ -151,6 +154,11 @@ function benchmark_test(FT)
     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, Δt), 230)
+    bench_press(
+        CMI_het.INP_concentration_frequency,
+        (ip_frostenberg, INPC, T_air_cold),
+        110,
+    )
     bench_press(CMI_hom.homogeneous_J, (ip.homogeneous, Delta_a_w), 230)
 
     # non-equilibrium