Merge #170

170: Move ice nucleation J params to CLIMAParams r=trontrytel a=trontrytel

This PR:
- moves ice nucleation parameters to CLIMAParameters
- specifies CLIMAParameters version in Compat

Co-authored-by: Anna Jaruga <ajaruga@caltech.edu>

Project.toml

@@ -14,7 +14,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [compat]
-CLIMAParameters = ">=0.7.9"
+CLIMAParameters = ">=0.7.12"
 CUDAKernels = "0.4"
 DocStringExtensions = "0.8, 0.9"
 ForwardDiff = "0.10"

docs/Project.toml

@@ -7,3 +7,5 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
+[compat]
+CLIMAParameters = ">=0.7.12"

docs/src/IceNucleation.md

@@ -8,7 +8,7 @@ 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) 
+  The parameterization for immersion freezing is an implementation of [KnopfAlpert2013](@cite)
   and is valid for droplets containing sulphuric acid.
 
 !!! note
@@ -45,16 +45,16 @@ Both parameters are dependent on aerosol properties and temperature.
     (either a second deposition or other condensation type mode).
 
 ## ABIFM for Sulphuric Acid Containing Droplets
-Water Activity-Based Immersion Freezing Model (ABFIM) 
+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). 
-  The nucleation rate coefficient, ``J``, describes the number of ice nuclei formed per unit area 
+  classical nucleation theory (CNT). More on CNT can be found in [Karthika2016](@cite).
+  The nucleation rate coefficient, ``J``, describes the number of ice nuclei formed per unit area
   per unit time and can be determined by the water activity, ``a_w``. This parameterization follows
   [KnopfAlpert2013](@cite), [Koop2002](@cite), [MurphyKoop2005](@cite), and [Luo1995](@cite). In this model,
   aerosols are assumed to contain an insoluble and soluble material. When immersed in water,
   the soluble material diffuses into the liquid water to create a sulphuric acid solution.
 
-Using empirical coefficients, ``m`` and ``c``, from [KnopfAlpert2013](@cite), 
+Using empirical coefficients, ``m`` and ``c``, from [KnopfAlpert2013](@cite),
   the heterogeneous nucleation rate coefficient in units of ``cm^{-2}s^{-1}`` can be determined by the linear equation
 ```math
 \begin{equation}
@@ -63,10 +63,10 @@ Using empirical coefficients, ``m`` and ``c``, from [KnopfAlpert2013](@cite),
 ```
 !!! note
 
-    Our source code for the nucleation rate coefficient returns 
+    Our source code for the nucleation rate coefficient returns
     ``J`` in base SI units.
 
-``\Delta a_w``is the difference between the water activity of the droplet, ``a_w``, and the water activity of ice at the same temperature, ``a_{w,ice}(T)``. From [Koop2002](@cite), 
+``\Delta a_w``is the difference between the water activity of the droplet, ``a_w``, and the water activity of ice at the same temperature, ``a_{w,ice}(T)``. From [Koop2002](@cite),
 ```math
 \begin{equation}
   a_w = \frac{p_{sol}}{p_{sat}}
@@ -96,8 +96,8 @@ Once ``J_{het}`` is calculated, it can be used to determine the ice production r
   P_{ice} = J_{het}A(N_{tot}-N_{ice})
 \end{equation}
 ```
-where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number 
-  of ice nuclei, and ``N_{ice}`` is number of ice crystals already in the system. 
+where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number
+  of ice nuclei, and ``N_{ice}`` is number of ice crystals already in the system.
 
 ## Homogeneous Freezing for Sulphuric Acid Containing Droplets
 Homogeneous freezing occurs when supercooled liquid droplets freeze on their own.
@@ -144,13 +144,13 @@ Delta_a = Vector{Float64}(undef, length(temp))
 J = Vector{Float64}(undef, length(temp))
 
 # Knopf and Alpert 2013 Figure 4A
-# https://doi.org/10.1039/C3FD00035D 
+# https://doi.org/10.1039/C3FD00035D
 dust_type = CT.KaoliniteType()
 
 it = 1
 for T in temp
         Delta_a[it] = CO.Delta_a_w(prs, x, T)
-        J[it] = IN.ABIFM_J(dust_type, Delta_a[it])
+        J[it] = IN.ABIFM_J(prs, dust_type, Delta_a[it])
         global it += 1
 end
 log10J_converted = @. log10(J*1e-4)

parcel/Project.toml

@@ -8,3 +8,5 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
+[compat]
+CLIMAParameters = ">=0.7.12"

src/Common.jl

@@ -124,24 +124,41 @@ function logistic_function_integral(x::FT, x_0::FT, k::FT) where {FT <: Real}
 end
 
 """
-    H2SO4_soln_saturation_vapor_pressure(x, T)
+    H2SO4_soln_saturation_vapor_pressure(prs, x, T)
 
+ - `prs` - a set with free parameters
  - `x` - wt percent sulphuric acid [unitless]
  - `T` - air temperature [K].
 
 Returns the saturation vapor pressure above a sulphuric acid solution droplet in Pa.
 `x` is, for example, 0.1 if droplets are 10 percent sulphuric acid by weight
 """
-function H2SO4_soln_saturation_vapor_pressure(x::FT, T::FT) where {FT <: Real}
+function H2SO4_soln_saturation_vapor_pressure(
+    prs::APS,
+    x::FT,
+    T::FT,
+) where {FT <: Real}
+
+    T_max::FT = CMP.H2SO4_sol_T_max(prs)
+    T_min::FT = CMP.H2SO4_sol_T_min(prs)
+    w_2::FT = CMP.H2SO4_sol_w_2(prs)
+
+    c1::FT = CMP.H2SO4_sol_c1(prs)
+    c2::FT = CMP.H2SO4_sol_c2(prs)
+    c3::FT = CMP.H2SO4_sol_c3(prs)
+    c4::FT = CMP.H2SO4_sol_c4(prs)
+    c5::FT = CMP.H2SO4_sol_c5(prs)
+    c6::FT = CMP.H2SO4_sol_c6(prs)
+    c7::FT = CMP.H2SO4_sol_c7(prs)
 
-    @assert T < FT(235)
-    @assert T > FT(185)
+    @assert T < T_max
+    @assert T > T_min
 
-    w_h = 1.4408 * x
+    w_h = w_2 * x
     p_sol =
         exp(
-            23.306 - 5.3465 * x + 12 * x * w_h - 8.19 * x * w_h^2 +
-            (-5814 + 928.9 * x - 1876.7 * x * w_h) / T,
+            c1 - c2 * x + c3 * x * w_h - c4 * x * w_h^2 +
+            (c5 + c6 * x - c7 * x * w_h) / T,
         ) * 100 # * 100 converts mbar --> Pa
     return p_sol
 end
@@ -160,7 +177,7 @@ function Delta_a_w(prs::APS, x::FT, T::FT) where {FT <: Real}
 
     thermo_params = CMP.thermodynamics_params(prs)
 
-    p_sol = H2SO4_soln_saturation_vapor_pressure(x, T)
+    p_sol = H2SO4_soln_saturation_vapor_pressure(prs, x, T)
     p_sat = TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid())
     p_ice = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice())
 

src/IceNucleation.jl

@@ -22,14 +22,12 @@ S0_cold(prs::APS, ::CT.DesertDustType) = CMP.S0_cold_DD_Mohler2006(prs)
 a_warm(prs::APS, ::CT.DesertDustType) = CMP.a_warm_DD_Mohler2006(prs)
 a_cold(prs::APS, ::CT.DesertDustType) = CMP.a_cold_DD_Mohler2006(prs)
 
-J_het_coeff_m(::CT.DesertDustType) = 22.62
-J_het_coeff_c(::CT.DesertDustType) = -1.35
-
-J_het_coeff_m(::CT.KaoliniteType) = 54.58834
-J_het_coeff_c(::CT.KaoliniteType) = -10.54758
-
-J_het_coeff_m(::CT.IlliteType) = 54.48075
-J_het_coeff_c(::CT.IlliteType) = -10.66873
+J_het_m(prs::APS, ::CT.DesertDustType) = CMP.J_ABIFM_m_KA2013_DesertDust(prs)
+J_het_c(prs::APS, ::CT.DesertDustType) = CMP.J_ABIFM_c_KA2013_DesertDust(prs)
+J_het_m(prs::APS, ::CT.KaoliniteType) = CMP.J_ABIFM_m_KA2013_Kaolinite(prs)
+J_het_c(prs::APS, ::CT.KaoliniteType) = CMP.J_ABIFM_c_KA2013_Kaolinite(prs)
+J_het_m(prs::APS, ::CT.IlliteType) = CMP.J_ABIFM_m_KA2013_Illite(prs)
+J_het_c(prs::APS, ::CT.IlliteType) = CMP.J_ABIFM_c_KA2013_Illite(prs)
 
 """
     dust_activated_number_fraction(prs, Si, T, dust_type)
@@ -66,23 +64,31 @@ function dust_activated_number_fraction(
 end
 
 """
-    ABIFM_J(dust_type, Δa_w)
+    ABIFM_J(prs, dust_type, Δa_w)
 
- - `dust_type` - choosing aerosol type
+ - `prs` - set with free parameters
+ - `dust_type` - aerosol type
  - `Δa_w` - change in water activity [unitless].
 
-Returns the immersion freezing nucleation rate coefficient, `J`, in m^-2 s^-1 for sulphuric acid containing solutions. 
-For other solutions, p_sol should be adjusted accordingly. Delta_a_w can be found using the Delta_a_w function in Common.jl. 
-`m` and `c` constants are taken from Knopf & Alpert 2013.
+Returns the immersion freezing nucleation rate coefficient, `J`, in m^-2 s^-1
+for sulphuric acid solutions.
+For other solutions, p_sol should be adjusted accordingly.
+Delta_a_w can be found using the Delta_a_w function in Common.jl.
+The free parameters `m` and `c` are taken from Knopf & Alpert 2013
+see DOI: 10.1039/C3FD00035D
 """
-function ABIFM_J(dust_type::CT.AbstractAerosolType, Δa_w::FT) where {FT <: Real}
+function ABIFM_J(
+    prs::APS,
+    dust_type::CT.AbstractAerosolType,
+    Δa_w::FT,
+) where {FT <: Real}
 
-    m = J_het_coeff_m(dust_type)
-    c = J_het_coeff_c(dust_type)
+    m::FT = J_het_m(prs, dust_type)
+    c::FT = J_het_c(prs, dust_type)
 
-    logJ = m * Δa_w + c
+    logJ::FT = m * Δa_w + c
 
-    return max(0, 10^logJ * 100^2) # converts cm^-2 s^-1 to m^-2 s^-1
+    return max(FT(0), FT(10)^logJ * FT(1e4)) # converts cm^-2 s^-1 to m^-2 s^-1
 end
 
 end # end module
@@ -101,23 +107,31 @@ const APS = CMP.AbstractCloudMicrophysicsParameters
 export homogeneous_J
 
 """
-    homogeneous_J(Δa_w)
+    homogeneous_J(prs, Δa_w)
 
- - `Δa_w` - change in water activity
+ - `prs` - a set with free parameters
+ - `Δa_w` - change in water activity [-].
 
-Returns the homogeneous freezing nucleation rate coefficient, `J`, in m^-3 s^-1 for sulphuric
-acid containing solutions. Parameterization based off Koop 2000.
+Returns the homogeneous freezing nucleation rate coefficient,
+`J`, in m^-3 s^-1 for sulphuric acid solutions.
+Parameterization based on Koop 2000, DOI: 10.1038/35020537.
 Delta_a_w can be found using the Delta_a_w function in Common.jl.
 """
-function homogeneous_J(Δa_w::FT) where {FT <: Real}
+function homogeneous_J(prs::APS, Δa_w::FT) where {FT <: Real}
+
+    Δa_w_min::FT = CMP.Koop2000_min_delta_aw(prs)
+    Δa_w_max::FT = CMP.Koop2000_max_delta_aw(prs)
+    c1::FT = CMP.Koop2000_J_hom_c1(prs)
+    c2::FT = CMP.Koop2000_J_hom_c2(prs)
+    c3::FT = CMP.Koop2000_J_hom_c3(prs)
+    c4::FT = CMP.Koop2000_J_hom_c4(prs)
 
-    @assert Δa_w > 0.26
-    @assert Δa_w < 0.34
+    @assert Δa_w > Δa_w_min
+    @assert Δa_w < Δa_w_max
 
-    logJ = -906.7 + 8502 * Δa_w - 26924 * Δa_w^2 + 29180 * Δa_w^3
-    J = 10^(logJ)
+    logJ::FT = c1 + c2 * Δa_w - c3 * Δa_w^2 + c4 * Δa_w^3
 
-    return J * 1e6
+    return FT(10)^(logJ) * 1e6
 end
 
 end # end module

src/Parameters.jl

@@ -164,6 +164,28 @@ Base.@kwdef struct CloudMicrophysicsParameters{FT, TP} <:
     S0_cold_DD_Mohler2006::FT
     a_warm_DD_Mohler2006::FT
     a_cold_DD_Mohler2006::FT
+    J_ABIFM_m_KA2013_DesertDust::FT
+    J_ABIFM_c_KA2013_DesertDust::FT
+    J_ABIFM_m_KA2013_Kaolinite::FT
+    J_ABIFM_c_KA2013_Kaolinite::FT
+    J_ABIFM_m_KA2013_Illite::FT
+    J_ABIFM_c_KA2013_Illite::FT
+    Koop2000_min_delta_aw::FT
+    Koop2000_max_delta_aw::FT
+    Koop2000_J_hom_c1::FT
+    Koop2000_J_hom_c2::FT
+    Koop2000_J_hom_c3::FT
+    Koop2000_J_hom_c4::FT
+    H2SO4_sol_T_max::FT
+    H2SO4_sol_T_min::FT
+    H2SO4_sol_w_2::FT
+    H2SO4_sol_c1::FT
+    H2SO4_sol_c2::FT
+    H2SO4_sol_c3::FT
+    H2SO4_sol_c4::FT
+    H2SO4_sol_c5::FT
+    H2SO4_sol_c6::FT
+    H2SO4_sol_c7::FT
     molmass_seasalt::FT
     rho_seasalt::FT
     osm_coeff_seasalt::FT

test/common_functions_tests.jl

@@ -63,18 +63,20 @@ function test_H2SO4_soln_saturation_vapor_pressure(FT)
         x_sulph = FT(0.1)
 
         # If T out of range
-        TT.@test_throws AssertionError("T < FT(235)") CO.H2SO4_soln_saturation_vapor_pressure(
+        TT.@test_throws AssertionError("T < T_max") CO.H2SO4_soln_saturation_vapor_pressure(
+            prs,
             x_sulph,
             T_too_warm,
         )
-        TT.@test_throws AssertionError("T > FT(185)") CO.H2SO4_soln_saturation_vapor_pressure(
+        TT.@test_throws AssertionError("T > T_min") CO.H2SO4_soln_saturation_vapor_pressure(
+            prs,
             x_sulph,
             T_too_cold,
         )
 
         # p_sol should be higher at warmer temperatures
-        TT.@test CO.H2SO4_soln_saturation_vapor_pressure(x_sulph, T_warm) >
-                 CO.H2SO4_soln_saturation_vapor_pressure(x_sulph, T_cold)
+        TT.@test CO.H2SO4_soln_saturation_vapor_pressure(prs, x_sulph, T_warm) >
+                 CO.H2SO4_soln_saturation_vapor_pressure(prs, x_sulph, T_cold)
     end
 end
 
@@ -97,15 +99,9 @@ end
 
 println("Testing Float64")
 test_H2SO4_soln_saturation_vapor_pressure(Float64)
-
-
 println("Testing Float32")
 test_H2SO4_soln_saturation_vapor_pressure(Float32)
-
-
 println("Testing Float64")
 test_Delta_a_w(Float64)
-
-
 println("Testing Float32")
 test_Delta_a_w(Float32)