Water Activity Documentation/Plots

docs/Project.toml

@@ -6,6 +6,8 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
+
 [compat]
 CLIMAParameters = "0.7.12"

docs/bibliography.bib

@@ -24,6 +24,18 @@ @article{Abdul-RazzakandGhan2000
   year = {2000}
 }
 
+@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},
+journal = {Atmospheric Chemistry and Physics},
+volume = {22},
+year = {2022},
+number = {1},
+pages = {65--91},
+URL = {https://acp.copernicus.org/articles/22/65/2022/},
+DOI = {10.5194/acp-22-65-2022}
+}
+
 @article{Beheng1994,
   title = {A parameterization of warm cloud microphysical conversion processes},
   author = {Beheng, K.D.},

docs/make.jl

@@ -16,6 +16,7 @@ pages = Any[
     "Smooth transition at thresholds" => "ThresholdsTransition.md",
     "Aerosol activation" => "AerosolActivation.md",
     "ARG2000 activation example" => "ARG2000_example.md",
+    "Water Activity" => "WaterActivity.md",
     "Ice Nucleation" => "IceNucleation.md",
     "Box model ice nucleation example" => "IceNucleationBox.md",
     "Aerosol Nucleation" => "AerosolNucleation.md",

docs/src/IceNucleation.md

@@ -2,21 +2,21 @@
 
 The `IceNucleation.jl` module includes
   the parameterization of activation of dust aerosol particles into ice crystals
-  via deposition of water vapor and water activity based parameterization of immersion freezing.
-These are heterogeneous ice nucleation processes.
+  via deposition of water vapor, water activity based parameterization of immersion freezing,
+  and water activity based parameterization of homogeneous 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).
 
 !!! note
 
-    Future work includes adding parameterizations
-    for other nucleation paths such as
-    heterogeneous immersion freezing or homogeneous freezing
-    and modeling the competition between them.
+    Future work includes refining the homogeneous freezing
+    parameterization and modeling the competition between
+    freezing modes.
 
 ## Activated fraction for deposition freezing on dust
 The parameterization models the activated fraction
@@ -93,7 +93,7 @@ where ``x`` is the weight fraction of sulphuric acid in the droplets
 Once ``J_{ABIFM}`` is calculated, it can be used to determine the ice production rate, ``P_{ice}``, per second via immersion freezing.
 ```math
 \begin{equation}
-  P_{ice} = J_{ABIFM}A(N_{tot}-N_{ice})
+  P_{ice} = J_{ABIFM}A(N_{tot} - N_{ice})
 \end{equation}
 ```
 where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number
@@ -137,10 +137,17 @@ include("ice_nucleation_plots/KnopfAlpert2013_fig5.jl")
 Note that water activity of the droplet was assumed equal to relative humidity so that:
 ```math
 \begin{equation}
-  a_{w} = \frac{p_{sol}(x_{sulph} = 0, T = T_{dew})}{p_{sat}}
+  a_{w} = \frac{p_{sat}(T = T_{dew})}{p_{sat}(T)}
 \end{equation}
 ```
-where `T_dew` is the dewpoint (in this example it is -45C).
+where `T_dew` is the dewpoint (in this example, it is constant at -45C).
+
+!!! note
+
+    For the same figure in Knopf & Alpert,
+    the mixed-phase cloud uses T_{dew} = T.
+    We are unsure when to use constant T_{dew}
+    or set it equal to T.
 
 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

docs/src/WaterActivity.md

@@ -0,0 +1,96 @@
+# Water Activity
+The `Common.jl` module includes
+  a parameterization for difference in water activity between a H2SO4
+  solution droplet and ice. This can be used in immersion and homogeneous
+  freezing parameterizations of nucleation rate coefficient, ``J``.
+
+``\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)``. When the
+  droplet is in equilibrium with its surroundings, ``a_w`` is equivalent to relative
+  humidity. Otherwise, a parameterization can be found in the `Common.jl` file and
+  comes from [Koop2002](@cite),
+```math
+\begin{equation}
+  a_w = \frac{p_{sol}}{p_{sat}}
+\end{equation}
+```
+```math
+\begin{equation}
+  a_{w,ice} = \frac{p_{i,sat}}{p_{sat}}
+\end{equation}
+```
+where ``p_{sol}`` is saturated vapor pressure of water above solution, ``p_{sat}``
+  is saturated vapor pressure above pure liquid water, and ``p_{i,sat}`` is saturated
+  vapor pressure above ice. ``p_{sol}`` is determined in mbar using a parameterization
+  for supercooled, binary ``H_2SO_4/H_2O`` solution from [Luo1995](@cite) which is only valid for ``185K < T < 235K``:
+```math
+\begin{equation}
+  ln(p_{sol}) = 23.306 - 5.3465x + 12xw_h - 8.19xw_h^2 + \frac{1}{T}(-5814 + 928.9x - 1876.7xw_h)
+\end{equation}
+```
+where ``x`` is the weight fraction of sulphuric acid in the droplets
+  (i.e. if droplets are 10% sulphuric acid by mass, ``x = 0.1``), ``w_h = 1.4408x``,
+  and temperature is in Kelvins. Note that ``x`` should vary as the size of the droplet
+  changes.
+
+!!! note
+
+    There is a need to find a parameterization for p_{sol}
+    at temperatures warmer than 235K for mixed phase clouds.
+
+For now, the equation used to find water activity of a droplet at equilibrium at temperatures warmer than 235K is taken from [Baumgartner2022](@cite) equation 4:
+```math
+\begin{equation}
+  a_w = S_i \frac{p_{i,sat}(T)}{p_{sat}(T)}
+\end{equation}
+```
+``S_i`` is water vapor saturation raito with respect to ice and defined by
+```math
+\begin{equation}
+  S_i = \frac{p_{v}}{p_{i,sat}(T)}
+\end{equation}
+```
+where ``p_v`` is the ambient partial pressure of water vapor.
+
+!!! note
+
+    Only use the above parameterization for S_i in the application
+    that it is not available. i.e. for a parcel model, use the S_i
+    calculated in the parcel model.
+
+``p_{sat}`` is currently defined using the `Thermodynamics.jl` library.
+  An alternative is using the parameterization for ``p_{sol}`` and setting
+  `x = 0`. We choose to use the `Thermodynamics.jl` to keep saturated vapor
+  pressure of pure water consistent throughout `CloudMicrophysic.jl`.
+
+Similarly, ``p_{i,sat}`` is defined through the `Thermodynamics.jl` library. It is
+  possible to calculate it using chemical potential as done in [Koop2000](@cite), however,
+  we will use the current parameterization for consistency since they give similar results.
+
+To verify that our parameterizations for ``p_{i,sat}`` and ``p_{sat}`` from
+    `Thermodynamics.jl` are valid at cold temperatures, we can compare to various
+    other vapor pressure parameterizations as found in [Baumgartner2022](@cite).
+    Here, ``CM`` refers to `CloudMicrophysics.jl` and all dashed curves are different
+    implementations of ``p_{i,sat}`` and ``p_{sat}`` (``p_{sat}`` labelled as
+    ``p_{liq}`` to emphasize ice vs liquid phase of the pure water).
+```@example
+include("water_activity_plots/p_sol_parameterizations.jl")
+```
+![](vap_pressure_vs_T.svg)
+
+To verify that our parameterizations for water activty using `Thermodynamics.jl`
+    is ok, we plot critical water activity (water activtiy at which freezing occurs)
+    against the other variations found in [Baumgartner2022](@cite). Though not
+    matching exactly, our parameterization is relatively close to the other parameterizations.
+```@example
+include("water_activity_plots/Baumgartner2022_fig5.jl")
+```
+![](Baumgartner2022_fig5.svg)
+Shown in red is the water activity over ice using our parameterization. With these two lines plotted (critical water activity of the droplet and ice water activity), we create a phase diagram. Under the red line is liquid, above the critical water activity is ice, and between the two curves is supercooled liquid.
+
+Another plot to test if our parameterization is reasonable is plotting against other parameterizations of water activity (as opposed to critical water activity) as a function of temperature. Plotted in green are various ways to compute water activity over ice. ``using p(0,T)`` refers to how the denominator, ``a_w``, is calculated. By default, this is parameterized assuming a pure liquid droplet with `Thermodynamics.jl`. ``using p(0,T)`` implies that the parameterization of vapor pressure of a solution droplet is used but setting concentration of H2SO4 to zero. ``using \mu`` refers to the parameterization used in [Koop2000](@cite) where water activity is dependent on chemical potential.
+```@example
+include("water_activity_plots/T_vs_wateractivity.jl")
+```
+![](T_vs_wateractivity.svg)
+Taking the difference between any pair of blue and green lines will give a ``\Delta a_w(T)``. Since all the blue lines are similar and all the green lines are similar, we can assume that our parameterization of pure liquid and ice water activities are reasonable.

docs/src/water_activity_plots/Baumgartner2022_fig5.jl

@@ -0,0 +1,62 @@
+using Roots
+
+import CairoMakie as MK
+
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+const CMT = CM.CommonTypes
+const CMO = CM.Common
+const CMI = CM.HetIceNucleation
+const CMP = CM.Parameters
+
+include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
+# Boiler plate code to have access to model parameters and constants
+FT = Float64
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+prs = cloud_microphysics_parameters(toml_dict)
+thermo_params = CMP.thermodynamics_params(prs)
+
+# Baumgartner at al 2022 Figure 5
+# https://acp.copernicus.org/articles/22/65/2022/acp-22-65-2022.pdf
+#! format: off
+Baumgartner2022_fig5_temp = [182.903, 189.941, 199.9706, 208.4164, 216.3929, 223.2551, 227.5366]
+Baumgartner2022_fig5_water_activity = [0.7930, 0.8129, 0.8416, 0.8679, 0.89781, 0.928486, 0.95039]
+
+T_range = range(190, stop = 234, length = 100)
+
+p_sat_liq = [TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) for T in T_range]  # sat vap pressure over pure liq water using TD package
+p_sat_ice = [TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) for T in T_range]     # sat vap pressure over ice using TD package
+a_w_ice = p_sat_ice ./ p_sat_liq
+
+radius = 1.0e-4                             # cm
+J_crit = 1 / (4 / 3 * pi * radius^3) / 60.0 # cm^-3 s^-1
+fun(Delta_a)= - 906.7 + 8502 * Delta_a - 26924 * Delta_a^2 + 29180 * Delta_a^3 - log10(J_crit) # homogeneous J from Koop 2000
+initial_guess_Delta_a = 0.34
+Delta_a_crit = find_zero(fun, initial_guess_Delta_a)
+
+a_w = [Delta_a_crit + (TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()))/(TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid())) for T in T_range]
+
+# various ways Baumgartner calculates vapor pressures
+Baumgartner_p_ice(T) = exp(9.550426 - 5723.265 / T + 3.53068*log(T) - 0.00728332*T)
+MK_p_liq(T) = exp(54.842763 - 6763.22 / T - 4.21*log(T) + 0.000367*T + tanh(0.0415*(T - 218.8)) * (53.878 - 1331.22 / T - 9.44523*log(T) + 0.014025*T))
+Tab_p_liq(T) = exp(100 * (18.4524 - 3505.15788/T - 330918.55/(T^2) + 12725068.26/(T^3)))
+Nach_p_liq(T) = exp(74.8727 - 7167.405/T - 7.77107*log(T) + 0.00505*T)
+
+a_w_MK = [Delta_a_crit + (Baumgartner_p_ice(T))/(MK_p_liq(T)) for T in T_range]
+a_w_Tab = [Delta_a_crit + (Baumgartner_p_ice(T))/(Tab_p_liq(T)) for T in T_range]
+a_w_Nach = [Delta_a_crit + (Baumgartner_p_ice(T))/(Nach_p_liq(T)) for T in T_range]
+a_w_Luo = [Delta_a_crit + (Baumgartner_p_ice(T))/(CMO.H2SO4_soln_saturation_vapor_pressure(prs, 0.0, T)) for T in T_range]
+
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], title = "Water Activity vs Temperature", ylabel = "a_w [-]", xlabel = "T [K]")
+MK.lines!(ax1, Baumgartner2022_fig5_temp, Baumgartner2022_fig5_water_activity, label = "Baumgartner2022", linestyle = :dash, color = :blue)
+MK.lines!(ax1, T_range, a_w, label = "CloudMicrophysics", color = :blue)
+MK.lines!(ax1, T_range, a_w_MK, label = "Baum_MK", color = :green)
+MK.lines!(ax1, T_range, a_w_Nach, label = "Baum_Nach", color = :lightgreen)
+MK.lines!(ax1, T_range, a_w_Luo, label = "Baum_Luo", color = :darkgreen)
+MK.lines!(ax1, T_range, a_w_ice, label = "a_w_ice", color = :red)
+MK.axislegend(position = :lt)
+MK.save("Baumgartner2022_fig5.svg", fig)
+#! format: on

docs/src/water_activity_plots/T_vs_wateractivity.jl

@@ -0,0 +1,44 @@
+import CairoMakie as MK
+
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+const CMT = CM.CommonTypes
+const CMO = CM.Common
+const CMI = CM.HetIceNucleation
+const CMP = CM.Parameters
+
+include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
+# Boiler plate code to have access to model parameters and constants
+FT = Float64
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+prs = cloud_microphysics_parameters(toml_dict)
+thermo_params = CMP.thermodynamics_params(prs)
+
+T_range = range(190, stop = 234, length = 100)
+x = FT(0.1)
+#! format: off
+p_sol_1 = [CMO.H2SO4_soln_saturation_vapor_pressure(prs, x, T) for T in T_range]     # p_sol for concentration x
+p_sol_0 = [CMO.H2SO4_soln_saturation_vapor_pressure(prs, 0.0, T) for T in T_range]   # sat vap pressure over pure liq water using p_sol eqn
+p_sat_liq = [TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) for T in T_range]  # sat vap pressure over pure liq water using TD package
+p_sat_ice = [TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) for T in T_range]     # sat vap pressure over ice using TD package
+
+a_w = p_sol_1 ./ p_sat_liq                 # a_w current parameterization
+a_w_alternate = p_sol_1 ./ p_sol_0         # a_w if sat vapor pressure over pure liq water was calculated with p_sol eqn
+a_w_ice = p_sat_ice ./ p_sat_liq           # a_ice current parameterization
+a_w_ice_alternate = p_sat_ice ./ p_sol_0   # a_ice if sat vapor pressure over pure liq water was calculated with p_sol eqn
+a_w_ice_μ = [exp((210368 + 131.438*T - (3.32373e6 /T) - 41729.1*log(T))/(8.31441*T)) for T in T_range]  # a_ice using chemical potential parameterization
+
+# Plotting. NOTE: all ".* -1" is only to flip x axis
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], title = "Temperature vs Water Activity", ylabel = "T [K]", xlabel = "-a_w", limits = ((-1, -0.4), nothing))
+MK.lines!(ax1, a_w .* -1, T_range, label = "CM default a_w", color = :blue)
+MK.lines!(ax1, a_w_alternate .* -1, T_range, label = "a_w using p(0,T)", linestyle = :dash, color = :blue)
+MK.lines!(ax1, a_w_ice .* -1, T_range, label = "CM default a_w_ice", color = :green)
+MK.lines!(ax1, a_w_ice_alternate .* -1, T_range, label = "a_w_ice using p(0,T)", color = :green, linestyle = :dash)
+MK.lines!(ax1, a_w_ice_μ .* -1, T_range, label = "a_w_ice using μ", color = :lightgreen)
+MK.axislegend()
+
+MK.save("T_vs_wateractivity.svg", fig)
+#! format: on

docs/src/water_activity_plots/p_sol_parameterizations.jl

@@ -0,0 +1,48 @@
+using Roots
+
+import CairoMakie as MK
+
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+const CMT = CM.CommonTypes
+const CMO = CM.Common
+const CMI = CM.HetIceNucleation
+const CMP = CM.Parameters
+
+include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
+# Boiler plate code to have access to model parameters and constants
+FT = Float64
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+prs = cloud_microphysics_parameters(toml_dict)
+thermo_params = CMP.thermodynamics_params(prs)
+
+# Baumgartner at al 2022 Figure 5
+# https://acp.copernicus.org/articles/22/65/2022/acp-22-65-2022.pdf
+#! format: off
+Baumgartner2022_fig5_temp = [182.903, 189.941, 199.9706, 208.4164, 216.3929, 223.2551, 227.5366]
+Baumgartner2022_fig5_water_activity = [0.7930, 0.8129, 0.8416, 0.8679, 0.89781, 0.928486, 0.95039]
+
+T_range = range(190, stop = 234, length = 100)
+
+# our vapor pressure parameterizations
+p_sat_liq = [TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) for T in T_range]  # sat vap pressure over pure liq water using TD package
+p_sat_ice = [TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) for T in T_range]     # sat vap pressure over ice using TD package
+
+# vapor pressure parameterizations mentioned in Baumgartner 2022
+Baumgartner_p_ice(T) = exp(9.550426 - 5723.265 / T + 3.53068*log(T) - 0.00728332*T)
+MK_p_liq(T) = exp(54.842763 - 6763.22 / T - 4.21*log(T) + 0.000367*T + tanh(0.0415*(T - 218.8)) * (53.878 - 1331.22 / T - 9.44523*log(T) + 0.014025*T))
+Tab_p_liq(T) = exp(100 * (18.4524 - 3505.15788/T - 330918.55/(T^2) + 12725068.26/(T^3)))
+Nach_p_liq(T) = exp(74.8727 - 7167.405/T - 7.77107*log(T) + 0.00505*T)
+
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], title = "Saturated Vapor Pressure vs Temperature", ylabel = "Pressure [Pa]", xlabel = "T [K]")
+MK.lines!(ax1,T_range, [Baumgartner_p_ice(T) for T in T_range], label = "Baumgartner's p_ice", color = :blue, linestyle = :dash)
+MK.lines!(ax1,T_range, p_sat_ice, label = "CM's p_ice", color = :blue)
+MK.lines!(ax1,T_range, [MK_p_liq(T) for T in T_range], label = "MK_p_liq", color = :green, linestyle = :dash)
+MK.lines!(ax1,T_range, [Nach_p_liq(T) for T in T_range], label = "Nach_p_liq", color = :lightgreen, linestyle = :dash)
+MK.lines!(ax1, T_range, [TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) for T in T_range], label = "CM's p_liq", color = :green)
+MK.axislegend(position = :lt)
+MK.save("vap_pressure_vs_T.svg", fig)
+#! format: on