F_r = 0 exception fixed, mass tests updated

docs/Project.toml

@@ -14,5 +14,6 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [compat]
+CLIMAParameters = "0.8"
 Documenter = "1.1"
 DocumenterCitations = "1.2"

docs/bibliography.bib

@@ -34,6 +34,18 @@ @article{Alpert2016
   year = {2016}
 }
 
+@Article{Alpert2022,
+author = "Alpert, Peter A. and Boucly, Anthony and Yang, Shuo and Yang, Huanyu and Kilchhofer, Kevin and Luo, Zhaochu and Padeste, Celestino and Finizio, Simone and Ammann, Markus and Watts, Benjamin",
+title = "Ice nucleation imaged with X-ray spectro-microscopy",
+journal = " Environ. Sci.: Atmos.",
+year = "2022",
+volume = "2",
+issue = "3",
+pages = "335-351",
+publisher = "RSC",
+doi = "10.1039/D1EA00077B",
+}
+
 @Article{Baumgartner2022,
 author = {Baumgartner, M. and Rolf, C. and Groo{\ss}, J.-U. and Schneider, J. and Schorr, T. and M\"ohler, O. and Spichtinger, P. and Kr\"amer, M.},
 title = {New investigations on homogeneous ice nucleation: the effects of water activity and water saturation formulations},
@@ -74,10 +86,22 @@ @article{Chen2022
   journal={Atmospheric Research},
   volume={293},
   number={106171},
-  doi = {https://doi.org/10.1016/j.atmosres.2022.106171},
+  doi = {10.1016/j.atmosres.2022.106171},
   year = {2022}
 }
 
+@article{China2017,
+  author = {China, Swarup and Alpert, Peter A. and Zhang, Bo and Schum, Simeon and Dzepina, Katja and Wright, Kendra and Owen, R. Chris and Fialho, Paulo and Mazzoleni, Lynn R. and Mazzoleni, Claudio and Knopf, Daniel A.},
+  title = {Ice cloud formation potential by free tropospheric particles from long-range transport over the Northern Atlantic Ocean},
+  journal = {Journal of Geophysical Research: Atmospheres},
+  volume = {122},
+  number = {5},
+  pages = {3065-3079},
+  keywords = {ice cloud, free troposphere, long-range transport, ice nucleation, atmospheric aging},
+  doi = {10.1002/2016JD025817},
+  year = {2017}
+}
+
 @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},
@@ -110,7 +134,7 @@ @article{Glassmeier2016
   address = {Boston MA, USA},
   volume = {73},
   number = {12},
-  doi = {https://doi.org/10.1175/JAS-D-16-0008.1},
+  doi = {10.1175/JAS-D-16-0008.1},
   pages = {5003 - 5023}
 }
 
@@ -164,7 +188,7 @@ @article{Karcher2006
   journal = {Journal of Geophysical Research: Atmospheres},
   volume = {111},
   number = {D1},
-  doi = {https://doi.org/10.1029/2005JD006219},
+  doi = {10.1029/2005JD006219},
   year = {2006}
 }
 
@@ -335,9 +359,7 @@ @article{Luo1995
 volume = {22},
 number = {3},
 pages = {247-250},
-doi = {https://doi.org/10.1029/94GL02988},
-url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/94GL02988},
-eprint = {https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/94GL02988},
+doi = {10.1029/94GL02988},
 year = {1995}
 }
 

docs/src/API.md

@@ -86,6 +86,7 @@ AerosolActivation.total_M_activated
 ```@docs
 HetIceNucleation
 HetIceNucleation.dust_activated_number_fraction
+HetIceNucleation.deposition_J
 HetIceNucleation.ABIFM_J
 ```
 
@@ -127,6 +128,8 @@ Parameters.ArizonaTestDust
 Parameters.DesertDust
 Parameters.Illite
 Parameters.Kaolinite
+Parameters.Feldspar
+Parameters.Ferrihydrite
 Parameters.Seasalt
 Parameters.Sulfate
 Parameters.AerosolActivationParameters

docs/src/IceNucleation.md

@@ -47,6 +47,25 @@ 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
+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
+  describes the number of ice nuclei formed per unit area of INP per unit time
+  dependent on the water activity criterion, ``\Delta  a_w``, and aerosol type.
+  The form of this empirical parameterization is taken from [KnopfAlpert2013](@cite).
+  Currently, we have parameters for kaolinite, feldspar, and ferrihydrite derived
+  from [China2017](@cite) and [Alpert2022](@cite).
+
+```math
+\begin{equation}
+  log_{10}J_{deposition} = m \Delta a_w + c
+\end{equation}
+```
+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
 Water Activity-Based Immersion Freezing Model (ABFIM)
   is a method of parameterizing immersion freezing inspired by the time-dependent
@@ -66,7 +85,8 @@ Using empirical coefficients, ``m`` and ``c``, from [KnopfAlpert2013](@cite),
 \end{equation}
 ```
 A parameterization for ``\Delta a_w`` can be found in `Common.jl`. More information on
-  it can be found in the `Water Activity` section.
+  it can be found in the `Water Activity` section. ``m`` and ``c`` here are different
+  from the ``m`` and ``c`` parameters for deposition nucleation.
 
 !!! note
 

docs/src/IceNucleationParcel0D.md

@@ -55,13 +55,15 @@ where:
 From the moist adiabatic assumption
 ```math
 \begin{equation}
-\frac{dT}{dt} = \frac{R_a T}{c_p p} \frac{dp}{dt} + \frac{L_v}{c_p} \frac{dq_l}{dt} + \frac{L_s}{c_p} \frac{d q_i}{dt}
+\frac{dT}{dt} = \frac{R_a T}{c_p p} \frac{dp}{dt} + \frac{L_v}{c_p} \frac{dq_{l,vap}}{dt} + \frac{L_s}{c_p} \frac{d q_{i,subl}}{dt} + \frac{L_f}{c_p} \frac{d q_{i,fus}}{dt}
 \end{equation}
 ```
 where:
-- ``q_l`` is the cloud liquid water specific humidity,
-- ``q_i`` is the cloud ice specific humidity,
-- ``L_s`` is the latent heat of sublimation.
+- ``q_{l,vap}`` is the cloud liquid water specific humidity from vaporization/condensation,
+- ``q_{i,subl}`` is the cloud ice specific humidity from sublimation/deposition,
+- ``q_{i,fus}`` is the cloud ice specific humidity from melting/freezing,
+- ``L_s`` is the latent heat of sublimation,
+- ``L_f`` is the latent heat of fusion.
 
 From hydrostatic balance and assuming constant vertical velocity:
 ```math
@@ -73,11 +75,11 @@ where:
 - ``g`` is the gravitational acceleration
 - ``w`` is the constant vertical velocity.
 
-Accounting for conservation of water, i.e. ``\frac{dq_v}{dt} + \frac{dq_l}{dt} + \frac{dq_i}{dt} = 0``,
+Accounting for conservation of water, i.e. ``\frac{dq_v}{dt} = - \frac{dq_{l,vap}}{dt} - \frac{dq_{i,subl}}{dt}``,
 and rearranging the terms
 ```math
 \begin{equation}
-\frac{dS_l}{dt} = a_1 w S_l - \left(a_2 + a_3 \right) S_l \frac{dq_l}{dt} - \left(a_2 + a_4\right) S_l \frac{dq_i}{dt}
+\frac{dS_l}{dt} = a_1 w S_l - \left(a_2 + a_3 \right)S_l  \frac{dq_{l,vap}}{dt} - \left(a_2 + a_4 \right) S_l \frac{dq_{i,subl}}{dt} - a_5 S_l \frac{dq_{i,fus}}{dt}
 \end{equation}
 ```
 where:
@@ -101,6 +103,11 @@ a_3 = \frac{L_v^2}{R_v T^2 c_p}
 a_4 = \frac{L_v L_s}{R_v T^2 c_p}
 \end{equation}
 ```
+```math
+\begin{equation}
+a_5 = \frac{L_v L_f}{R_v T^2 c_p}
+\end{equation}
+```
 
 Saturation ratio over ice can then be related to ``S_l`` by the relation
 ```math

parcel/parcel.jl

@@ -53,6 +53,7 @@ function parcel_model(dY, Y, p, t)
     R_a = TD.gas_constant_air(tps, q)
     cp_a = TD.cp_m(tps, q)
     L_subl = TD.latent_heat_sublim(tps, T)
+    L_fus = TD.latent_heat_fusion(tps, T)
     L_vap = TD.latent_heat_vapor(tps, T)
     ρ_air = TD.air_density(tps, ts)
     e = q_vap * p_a * R_v / R_a
@@ -62,6 +63,7 @@ function parcel_model(dY, Y, p, t)
     a2 = 1 / q_vap
     a3 = L_vap^2 / R_v / T^2 / cp_a
     a4 = L_vap * L_subl / R_v / T^2 / cp_a
+    a5 = L_vap * L_fus / R_v / cp_a / (T^2)
 
     # TODO - we should zero out all tendencies and augemnt them
     # TODO - add immersion, homogeneous, ...
@@ -135,12 +137,16 @@ function parcel_model(dY, Y, p, t)
         end
     end
 
+    dq_liq_dt_vap_to_liq = dql_dt_cond        # from liq-vap transitions
+    dq_ice_dt_vap_to_ice = dqi_dt_new_depo + dqi_dt_depo        # from ice-vap transitions
+    dq_ice_dt_liq_to_ice = dqi_dt_new_immers   # from ice-liq transitions
+
     # Update the tendecies
-    dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_immers
-    dq_liq_dt = dql_dt_cond - dqi_dt_new_immers
-    dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt - (a2 + a4) * S_liq * dq_ice_dt
+    dq_ice_dt = dq_ice_dt_vap_to_ice + dq_ice_dt_liq_to_ice
+    dq_liq_dt = dq_liq_dt_vap_to_liq - dq_ice_dt_liq_to_ice
+    dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt_vap_to_liq - (a2 + a4) * S_liq * dq_ice_dt_vap_to_ice - a5 * S_liq * dq_ice_dt_liq_to_ice
     dp_a_dt = -p_a * grav / R_a / T * w
-    dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt + L_subl / cp_a * dq_ice_dt
+    dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_vap_to_liq + L_fus / cp_a * dq_ice_dt_liq_to_ice + L_subl / cp_a * dq_ice_dt_vap_to_ice
     dq_vap_dt = -dq_ice_dt - dq_liq_dt
 
     # Set tendencies

src/IceNucleation.jl

@@ -7,6 +7,7 @@ import ..Parameters as CMP
 import Thermodynamics as TD
 
 export dust_activated_number_fraction
+export deposition_J
 export ABIFM_J
 
 """
@@ -41,7 +42,28 @@ function dust_activated_number_fraction(
 end
 
 """
-    ABIFM_J(dust, ip, Δa_w)
+    deposition_J(dust, Δa_w)
+
+ - `dust` - a struct with dust parameters
+ - `Δa_w` - change in water activity [unitless].
+
+Returns the deposition nucleation rate coefficient, `J`, in m^-2 s^-1
+for different minerals in liquid droplets.
+The free parameters `m` and `c` are derived from China et al (2017)
+see DOI: 10.1002/2016JD025817
+"""
+function deposition_J(
+    dust::Union{CMP.Ferrihydrite, CMP.Feldspar, CMP.Kaolinite},
+    Δa_w::FT,
+) where {FT}
+
+    logJ::FT = dust.deposition_m * Δa_w + dust.deposition_c
+
+    return max(FT(0), FT(10)^logJ * FT(1e4)) # converts cm^-2 s^-1 to m^-2 s^-1
+end
+
+"""
+    ABIFM_J(dust, Δa_w)
 
  - `dust` - a struct with dust parameters
  - `Δa_w` - change in water activity [unitless].
@@ -56,7 +78,7 @@ function ABIFM_J(
     Δa_w::FT,
 ) where {FT}
 
-    logJ::FT = dust.m * Δa_w + dust.c
+    logJ::FT = dust.ABIFM_m * Δa_w + dust.ABIFM_c
 
     return max(FT(0), FT(10)^logJ * FT(1e4)) # converts cm^-2 s^-1 to m^-2 s^-1
 end

src/P3Scheme.jl

@@ -120,32 +120,39 @@ Returns a named tuple containing:
 """
 function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT}
 
-    @assert ρ_r > FT(0)   # rime density must be positive ...
-    @assert ρ_r <= p3.ρ_l # ... and as a bulk ice density can't exceed the density of water
-    @assert F_r > FT(0)   # rime mass fraction must be positive ...
+    @assert F_r >= FT(0)   # rime mass fraction must be positive ...
     @assert F_r < FT(1)   # ... and there must always be some unrimed part
 
-    P3_problem(ρ_d) =
-        ρ_d - ρ_d_helper(
-            p3,
-            D_cr_helper(p3, F_r, ρ_g_helper(ρ_r, F_r, ρ_d)),
-            D_gr_helper(p3, ρ_g_helper(ρ_r, F_r, ρ_d)),
-        )
+    if F_r == 0
 
-    ρ_d =
-        RS.find_zero(
-            P3_problem,
-            RS.SecantMethod(FT(0), FT(1000)),
-            RS.CompactSolution(),
-        ).root
-    ρ_g = ρ_g_helper(ρ_r, F_r, ρ_d)
-
-    return (;
-        D_cr = D_cr_helper(p3, F_r, ρ_g),
-        D_gr = D_gr_helper(p3, ρ_g),
-        ρ_g,
-        ρ_d,
-    )
+        return (; D_cr = 0, D_gr = 0, ρ_g = 0, ρ_d = 0)
+
+    else
+        @assert ρ_r > FT(0)   # rime density must be positive ...
+        @assert ρ_r <= p3.ρ_l # ... and as a bulk ice density can't exceed the density of water
+
+        P3_problem(ρ_d) =
+            ρ_d - ρ_d_helper(
+                p3,
+                D_cr_helper(p3, F_r, ρ_g_helper(ρ_r, F_r, ρ_d)),
+                D_gr_helper(p3, ρ_g_helper(ρ_r, F_r, ρ_d)),
+            )
+
+        ρ_d =
+            RS.find_zero(
+                P3_problem,
+                RS.SecantMethod(FT(0), FT(1000)),
+                RS.CompactSolution(),
+            ).root
+        ρ_g = ρ_g_helper(ρ_r, F_r, ρ_d)
+
+        return (;
+            D_cr = D_cr_helper(p3, F_r, ρ_g),
+            D_gr = D_gr_helper(p3, ρ_g),
+            ρ_g,
+            ρ_d,
+        )
+    end
 end
 
 """

src/parameters/AerosolDesertDust.jl

@@ -20,9 +20,9 @@ struct DesertDust{FT} <: AerosolType{FT}
     "a for T < T_thr [-]"
     a_cold::FT
     "m coefficient for immersion freezing J [-]"
-    m::FT
+    ABIFM_m::FT
     "c coefficient for immersion freezing J [-]"
-    c::FT
+    ABIFM_c::FT
 end
 
 function DesertDust(

src/parameters/AerosolFeldspar.jl

@@ -0,0 +1,28 @@
+export Feldspar
+
+"""
+    Feldspar{FT}
+
+Parameters for Feldspar from Alpert et al 2022
+DOI: 10.1039/D1EA00077B
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+struct Feldspar{FT} <: AerosolType{FT}
+    "m coefficient for deposition nucleation J [-]"
+    deposition_m::FT
+    "c coefficient for deposition nucleation J [-]"
+    deposition_c::FT
+end
+
+function Feldspar(
+    ::Type{FT},
+    toml_dict::CP.AbstractTOMLDict = CP.create_toml_dict(FT),
+) where {FT}
+    (; data) = toml_dict
+    return Feldspar(
+        FT(data["Alpert2022_J_deposition_m_Feldspar"]["value"]),
+        FT(data["Alpert2022_J_deposition_c_Feldspar"]["value"]),
+    )
+end

src/parameters/AerosolFerrihydrite.jl

@@ -0,0 +1,28 @@
+export Ferrihydrite
+
+"""
+    Ferrihydrite{FT}
+
+Parameters for Ferrihydrite from Alpert et al 2022
+DOI: 10.1039/D1EA00077B
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+struct Ferrihydrite{FT} <: AerosolType{FT}
+    "m coefficient for deposition nucleation J [-]"
+    deposition_m::FT
+    "c coefficient for deposition nucleation J [-]"
+    deposition_c::FT
+end
+
+function Ferrihydrite(
+    ::Type{FT},
+    toml_dict::CP.AbstractTOMLDict = CP.create_toml_dict(FT),
+) where {FT}
+    (; data) = toml_dict
+    return Ferrihydrite(
+        FT(data["Alpert2022_J_deposition_m_Ferrihydrite"]["value"]),
+        FT(data["Alpert2022_J_deposition_c_Ferrihydrite"]["value"]),
+    )
+end