Add KA13 J vs T plots

CliMA · Aug 31, 2023 · baae1e0 · baae1e0
1 parent f9e9677
commit baae1e0
Show file tree

Hide file tree

Showing 3 changed files with 176 additions and 54 deletions.
diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md
@@ -58,7 +58,7 @@ 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}
-  log_{10}J_{het} = m \Delta a_w + c
+  log_{10}J_{ABIFM} = m \Delta a_w + c
 \end{equation}
 ```
 !!! note
@@ -90,10 +90,10 @@ 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.
 
-Once ``J_{het}`` is calculated, it can be used to determine the ice production rate, ``P_{ice}``, per second via immersion freezing.
+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_{het}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
@@ -116,57 +116,26 @@ 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``.
 
 ## ABIFM Example Figures
+The following plot shows ``J`` as a function of ``\Delta a_w`` as compared to figure 1 in Knopf & Alpert 2013. Solution droplets were assumed to contain a constant 10% wt. sulphuric acid. Changing the concentration will simply shift the line, following Knopf & Alpert's parameterization. As such, this plot is just to prove that ``J`` is correctly parameterized as a function of ``\Delta a_w``, with no implications of whether ``\Delta a_w`` is properly parameterized. For more on water activity, please see above section.
 ```@example
-import Plots
-
-import CloudMicrophysics
-import CLIMAParameters
-import Thermodynamics
-
-const PL = Plots
-const IN = CloudMicrophysics.HetIceNucleation
-const CMP = CloudMicrophysics.Parameters
-const CT = CloudMicrophysics.CommonTypes
-const CO = CloudMicrophysics.Common
-const CP =  CLIMAParameters
-const TD = Thermodynamics
-
-include(joinpath(pkgdir(CloudMicrophysics), "test", "create_parameters.jl"))
-FT = Float64
-toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
-const prs = cloud_microphysics_parameters(toml_dict)
-thermo_params = CMP.thermodynamics_params(prs)
-
-# Initializing
-temp = collect(210.0:2:232.0) # air temperature
-x = 0.1                     # wt% sulphuric acid in droplets
-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
-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(prs, dust_type, Delta_a[it])
-        global it += 1
-end
-log10J_converted = @. log10(J*1e-4)
-
-# data read from Fig 4 in Knopf & Alpert 2013
-# using https://automeris.io/WebPlotDigitizer/
-KA13_Delta_a_obs = [0.13641, 0.16205, 0.21538, 0.23897, 0.24513, 0.24718, 0.25026, 0.25128, 0.25231, 0.25333, 0.25538, 0.25744, 0.25846, 0.25949, 0.26051, 0.26051, 0.26462, 0.26462, 0.26872, 0.26974, 0.27077, 0.27077, 0.27179, 0.27385, 0.27692, 0.27795, 0.27795, 0.27795, 0.28308, 0.28410, 0.28410, 0.28615, 0.28718, 0.28718, 0.29128, 0.29128, 0.29231, 0.29333, 0.29744, 0.29744, 0.29744, 0.29949, 0.30359, 0.30462, 0.30564, 0.30667, 0.31077, 0.31077, 0.31077]
-KA13_log10J_obs = [-3.51880, -3.20301, 2.21053, 2.57143, 2.25564, 3.56391, 3.20301, 2.25564, 3.78947, 4.42105, 3.51880, 2.84211, 4.15038, 3.24812, 3.78947, 4.37594, 3.38346, 4.46617, 4.06015, 4.73684, 4.06015, 3.60902, 6.13534, 4.51128, 4.37594, 4.82707, 4.96241, 5.23308, 3.92481, 5.36842, 5.63910, 5.81955, 4.60150, 4.96241, 5.50376, 6.00000, 5.14286, 5.77444, 5.41353, 6.09023, 5.77444, 5.14286, 6.18045, 5.86466, 5.54887, 5.27820, 6.09023, 5.77444, 5.54887]
-
-KA13_Delta_a_param = [0.10256, 0.35692, 0.21949]
-KA13_log10J_param = [-4.91729, 8.97744, 1.44361]
-
-PL.plot(Delta_a, log10J_converted, label="CliMA", xlabel="Delta a_w [unitless]", ylabel="log10(J) [cm^-2 s^-1]")
-PL.scatter!(KA13_Delta_a_obs, KA13_log10J_obs, markercolor = :black, label="paper observations")
-PL.plot!(KA13_Delta_a_param, KA13_log10J_param, linecolor = :red, label="paper parameterization")
-
-PL.savefig("Knopf_Alpert_fig_1.svg")
+include("ice_nucleation_plots/KnopfAlpert2013_fig1.jl")
 ```
 ![](Knopf_Alpert_fig_1.svg)
+
+The following plot shows J as a function of temperature as compared to figure 5a in Knopf & Alpert 2013.
+
+```@example
+include("ice_nucleation_plots/KnopfAlpert2013_fig5.jl")
+```
+![](KnopfAlpert2013_fig5.svg)
+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}}
+\end{equation}
+```
+where `T_dew` is the dewpoint (in this example it is -45C).
+
+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.
diff --git a/docs/src/ice_nucleation_plots/KnopfAlpert2013_fig1.jl b/docs/src/ice_nucleation_plots/KnopfAlpert2013_fig1.jl
@@ -0,0 +1,73 @@
+import Plots as PL
+
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+import Thermodynamics as TD
+
+const IN = CM.HetIceNucleation
+const CMP = CM.Parameters
+const CT = CM.CommonTypes
+const CO = CM.Common
+
+include(joinpath(pkgdir(CloudMicrophysics), "test", "create_parameters.jl"))
+FT = Float64
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+const prs = cloud_microphysics_parameters(toml_dict)
+
+# Initializing
+T_range = range(210, stop = 232, length = 100)  # air temperature
+x = 0.1                                         # wt% sulphuric acid in droplets
+dust_type = CT.KaoliniteType()                  # dust type
+Δa_w = [CO.Delta_a_w(prs, x, T) for T in T_range]   # difference in solution and ice water activity
+J = @. IN.ABIFM_J(prs, (dust_type,), Δa_w)          # J in SI units
+log10J_converted = @. log10(J * 1e-4)           # converts J into cm^-2 s^-1 and takes log
+
+# Knopf and Alpert 2013 Figure 4A
+# https://doi.org/10.1039/C3FD00035D
+
+#! format: off
+# data read from Fig 4 in Knopf & Alpert 2013
+# using https://automeris.io/WebPlotDigitizer/
+KA13_Delta_a_obs = [
+    0.13641, 0.16205, 0.21538, 0.23897, 0.24513, 0.24718, 0.25026,
+    0.25128, 0.25231, 0.25333, 0.25538, 0.25744, 0.25846, 0.25949,
+    0.26051, 0.26051, 0.26462, 0.26462, 0.26872, 0.26974, 0.27077,
+    0.27077, 0.27179, 0.27385, 0.27692, 0.27795, 0.27795, 0.27795,
+    0.28308, 0.28410, 0.28410, 0.28615, 0.28718, 0.28718, 0.29128,
+    0.29128, 0.29231, 0.29333, 0.29744, 0.29744, 0.29744, 0.29949,
+    0.30359, 0.30462, 0.30564, 0.30667, 0.31077, 0.31077, 0.31077,
+]
+KA13_log10J_obs = [
+    -3.51880, -3.20301, 2.21053, 2.57143, 2.25564, 3.56391, 3.20301, 2.25564,
+    3.78947, 4.42105, 3.51880, 2.84211, 4.15038, 3.24812, 3.78947, 4.37594,
+    3.38346, 4.46617, 4.06015, 4.73684, 4.06015, 3.60902, 6.13534, 4.51128,
+    4.37594, 4.82707, 4.96241, 5.23308, 3.92481, 5.36842, 5.63910, 5.81955,
+    4.60150, 4.96241, 5.50376, 6.00000, 5.14286, 5.77444, 5.41353, 6.09023,
+    5.77444, 5.14286, 6.18045, 5.86466, 5.54887, 5.27820, 6.09023, 5.77444, 5.54887,
+]
+#! format: on
+
+KA13_Delta_a_param = [0.10256, 0.35692, 0.21949]
+KA13_log10J_param = [-4.91729, 8.97744, 1.44361]
+
+PL.plot(
+    Delta_a,
+    log10J_converted,
+    label = "CliMA",
+    xlabel = "Delta a_w [unitless]",
+    ylabel = "log10(J) [cm^-2 s^-1]",
+)
+PL.scatter!(
+    KA13_Delta_a_obs,
+    KA13_log10J_obs,
+    markercolor = :black,
+    label = "paper observations",
+)
+PL.plot!(
+    KA13_Delta_a_param,
+    KA13_log10J_param,
+    linecolor = :red,
+    label = "paper parameterization",
+)
+
+PL.savefig("Knopf_Alpert_fig_1.svg")
diff --git a/docs/src/ice_nucleation_plots/KnopfAlpert2013_fig5.jl b/docs/src/ice_nucleation_plots/KnopfAlpert2013_fig5.jl
@@ -0,0 +1,80 @@
+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)
+
+# Knopf and Alpert 2013 Figure 5A: Cirrus
+# https://doi.org/10.1039/C3FD00035D 
+temp = [
+    228.20357,
+    228.33571,
+    228.50357,
+    228.75000,
+    228.92143,
+    229.16786,
+    229.39643,
+    229.52143,
+]
+KA13_fig5A_J = [
+    70170382.86704,
+    12101528.74384,
+    1277935.32665,
+    52710.05764,
+    6040.39151,
+    293.39697,
+    18.97171,
+    4.18121,
+]
+
+# Our parameterization
+dust_type = CMT.IlliteType()    # dust type
+T_range = range(228.2, stop = 229.6, length = 100)  # air temperature
+T_dew = FT(228.0)               # dew point temperature
+a_sol = [                       # water activity of solution droplet at equilibrium
+    TD.saturation_vapor_pressure(thermo_params, T_dew, TD.Liquid()) /
+    TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) for
+    T in T_range
+]
+a_ice = [                       # water activity of ice
+    TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) /
+    TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) for
+    T in T_range
+]
+Δa_w = @. max(abs(a_sol - a_ice), FT(0.0))
+J_ABIFM = @. CMI.ABIFM_J(prs, (dust_type,), Δa_w) * 1e-4 # converted from SI units to cm^-2 s^-1
+
+# Plot results
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(
+    fig[1, 1],
+    ylabel = "J_het [cm^-2 s^-1]",
+    xlabel = "Temp [K]",
+    yscale = log10,
+)
+
+MK.ylims!(ax1, FT(1), FT(1e8))
+MK.lines!(
+    ax1,
+    temp,
+    KA13_fig5A_J,
+    label = "KA13 Figure 5A",
+    linestyle = :dash,
+    color = :green,
+)
+MK.lines!(ax1, temp, J_ABIFM, label = "CliMA", color = :green)
+fig[1, 2] = MK.Legend(fig, ax1)
+
+MK.save("KnopfAlpert2013_fig5.svg", fig)