From baae1e0f90c035514d7f9e7cbe218807d4c17002 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Fri, 18 Aug 2023 09:49:25 -0700 Subject: [PATCH] Add KA13 J vs T plots --- docs/src/IceNucleation.md | 77 ++++++------------ .../KnopfAlpert2013_fig1.jl | 73 +++++++++++++++++ .../KnopfAlpert2013_fig5.jl | 80 +++++++++++++++++++ 3 files changed, 176 insertions(+), 54 deletions(-) create mode 100644 docs/src/ice_nucleation_plots/KnopfAlpert2013_fig1.jl create mode 100644 docs/src/ice_nucleation_plots/KnopfAlpert2013_fig5.jl diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index 1e8462987e..e5e6e25001 100644 --- 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 new file mode 100644 index 0000000000..962b04e694 --- /dev/null +++ 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 new file mode 100644 index 0000000000..2943aa9d65 --- /dev/null +++ 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)