From 2aa178d3e34863b7ac37c9d2d880fa46d113e082 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Tue, 13 Feb 2024 01:17:25 -0800 Subject: [PATCH 1/3] Update to CLIMAParameters v0.9 and Thermodynamics v0.12 --- Project.toml | 4 +- docs/Project.toml | 2 +- parcel/Project.toml | 2 +- src/Common.jl | 8 ++-- src/parameters/Microphysics2M.jl | 22 +++++----- src/parameters/TerminalVelocity.jl | 67 ++++++++++-------------------- test/Project.toml | 1 + test/aerosol_activation_tests.jl | 1 + test/common_functions_tests.jl | 1 + test/microphysics0M_tests.jl | 1 + test/microphysics_noneq_tests.jl | 1 + test/p3_tests.jl | 1 + test/performance_tests.jl | 1 + 13 files changed, 48 insertions(+), 64 deletions(-) diff --git a/Project.toml b/Project.toml index 9a3a182899..20800d8b6c 100644 --- a/Project.toml +++ b/Project.toml @@ -12,10 +12,10 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] -CLIMAParameters = "0.8.6" +CLIMAParameters = "0.9" DocStringExtensions = "0.8, 0.9" ForwardDiff = "0.10" RootSolvers = "0.3, 0.4" SpecialFunctions = "1, 2" -Thermodynamics = "0.11" +Thermodynamics = "0.12" julia = "1.6" diff --git a/docs/Project.toml b/docs/Project.toml index b5b59f1221..20772a3894 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -16,6 +16,6 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] -CLIMAParameters = "0.8" +CLIMAParameters = "0.9" Documenter = "1.1" DocumenterCitations = "1.2" diff --git a/parcel/Project.toml b/parcel/Project.toml index 69532b37ad..842e4b8d4c 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -12,4 +12,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] -CLIMAParameters = "0.8" +CLIMAParameters = "0.9" diff --git a/src/Common.jl b/src/Common.jl index 056ade16fb..a712542bc8 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -218,12 +218,12 @@ function Chen2022_vel_coeffs( ρ::FT, ) where {FT} - (; ρ0, a1, a2, a3, a3_pow, b1, b2, b3, b_ρ, c1, c2, c3) = velo_scheme + (; ρ0, a, a3_pow, b, b_ρ, c) = velo_scheme q = exp(ρ0 * ρ) - ai = (a1 * q, a2 * q, a3 * q * ρ^a3_pow) - bi = (b1 - b_ρ * ρ, b2 - b_ρ * ρ, b3 - b_ρ * ρ) - ci = (c1, c2, c3) + ai = (a[1] * q, a[2] * q, a[3] * q * ρ^a3_pow) + bi = (b[1] - b_ρ * ρ, b[2] - b_ρ * ρ, b[3] - b_ρ * ρ) + ci = (c[1], c[2], c[3]) # unit conversions aiu = ai .* 1000 .^ bi diff --git a/src/parameters/Microphysics2M.jl b/src/parameters/Microphysics2M.jl index a74d45b8fa..042038da36 100644 --- a/src/parameters/Microphysics2M.jl +++ b/src/parameters/Microphysics2M.jl @@ -21,10 +21,10 @@ end function AcnvKK2000(td::CP.AbstractTOMLDict) name_map = (; - :KK2000_auctoconversion_coeff_A => :A, - :KK2000_auctoconversion_coeff_a => :a, - :KK2000_auctoconversion_coeff_b => :b, - :KK2000_auctoconversion_coeff_c => :c, + :KK2000_autoconversion_coeff_A => :A, + :KK2000_autoconversion_coeff_a => :a, + :KK2000_autoconversion_coeff_b => :b, + :KK2000_autoconversion_coeff_c => :c, ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) @@ -115,13 +115,13 @@ end function AcnvB1994(td::CP.AbstractTOMLDict) name_map = (; - :B1994_auctoconversion_coeff_C => :C, - :B1994_auctoconversion_coeff_a => :a, - :B1994_auctoconversion_coeff_b => :b, - :B1994_auctoconversion_coeff_c => :c, - :B1994_auctoconversion_coeff_N_0 => :N_0, - :B1994_auctoconversion_coeff_d_low => :d_low, - :B1994_auctoconversion_coeff_d_high => :d_high, + :B1994_autoconversion_coeff_C => :C, + :B1994_autoconversion_coeff_a => :a, + :B1994_autoconversion_coeff_b => :b, + :B1994_autoconversion_coeff_c => :c, + :B1994_autoconversion_coeff_N_0 => :N_0, + :B1994_autoconversion_coeff_d_low => :d_low, + :B1994_autoconversion_coeff_d_high => :d_high, :threshold_smooth_transition_steepness => :k, ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") diff --git a/src/parameters/TerminalVelocity.jl b/src/parameters/TerminalVelocity.jl index dbf82e7c11..9aa90bd105 100644 --- a/src/parameters/TerminalVelocity.jl +++ b/src/parameters/TerminalVelocity.jl @@ -156,36 +156,23 @@ end function Chen2022VelTypeSnowIce(toml_dict::CP.AbstractTOMLDict) # TODO: These should be array parameters. name_map = (; - :Chen2022_table_B3_As_coeff_1 => :A1, - :Chen2022_table_B3_As_coeff_2 => :A2, - :Chen2022_table_B3_As_coeff_3 => :A3, - :Chen2022_table_B3_Bs_coeff_1 => :B1, - :Chen2022_table_B3_Bs_coeff_2 => :B2, - :Chen2022_table_B3_Bs_coeff_3 => :B3, - :Chen2022_table_B3_Cs_coeff_1 => :C1, - :Chen2022_table_B3_Cs_coeff_2 => :C2, - :Chen2022_table_B3_Cs_coeff_3 => :C3, - :Chen2022_table_B3_Cs_coeff_4 => :C4, - :Chen2022_table_B3_Es_coeff_1 => :E1, - :Chen2022_table_B3_Es_coeff_2 => :E2, - :Chen2022_table_B3_Es_coeff_3 => :E3, - :Chen2022_table_B3_Fs_coeff_1 => :F1, - :Chen2022_table_B3_Fs_coeff_2 => :F2, - :Chen2022_table_B3_Fs_coeff_3 => :F3, - :Chen2022_table_B3_Gs_coeff_1 => :G1, - :Chen2022_table_B3_Gs_coeff_2 => :G2, - :Chen2022_table_B3_Gs_coeff_3 => :G3, + :Chen2022_table_B3_As => :As, + :Chen2022_table_B3_Bs => :Bs, + :Chen2022_table_B3_Cs => :Cs, + :Chen2022_table_B3_Es => :Es, + :Chen2022_table_B3_Fs => :Fs, + :Chen2022_table_B3_Gs => :Gs, :density_ice_water => :ρᵢ, ) p = CP.get_parameter_values(toml_dict, name_map, "CloudMicrophysics") FT = CP.float_type(toml_dict) return Chen2022VelTypeSnowIce{FT}( - p.A1 * (log(p.ρᵢ))^2 − p.A2 * log(p.ρᵢ) - p.A3, - FT(1) / (p.B1 + p.B2 * log(p.ρᵢ) + p.B3 / sqrt(p.ρᵢ)), - p.C1 + p.C2 * exp(p.C3 * p.ρᵢ) + p.C4 * sqrt(p.ρᵢ), - p.E1 - p.E2 * (log(p.ρᵢ))^2 + p.E3 * sqrt(p.ρᵢ), - -exp(p.F1 - p.F2 * (log(p.ρᵢ))^2 + p.F3 * log(p.ρᵢ)), - FT(1) / (p.G1 + p.G2 / (log(p.ρᵢ)) - p.G3 * log(p.ρᵢ) / p.ρᵢ), + p.As[2] * (log(p.ρᵢ))^2 − p.As[3] * log(p.ρᵢ) + p.As[1], + FT(1) / (p.Bs[1] + p.Bs[2] * log(p.ρᵢ) + p.Bs[3] / sqrt(p.ρᵢ)), + p.Cs[1] + p.Cs[2] * exp(p.Cs[3] * p.ρᵢ) + p.Cs[4] * sqrt(p.ρᵢ), + p.Es[1] - p.Es[2] * (log(p.ρᵢ))^2 + p.Es[3] * sqrt(p.ρᵢ), + -exp(p.Fs[1] - p.Fs[2] * (log(p.ρᵢ))^2 + p.Fs[3] * log(p.ρᵢ)), + FT(1) / (p.Gs[1] + p.Gs[2] / (log(p.ρᵢ)) - p.Gs[3] * log(p.ρᵢ) / p.ρᵢ), p.ρᵢ, ) end @@ -200,19 +187,13 @@ DOI: 10.1016/j.atmosres.2022.106171 # Fields $(DocStringExtensions.FIELDS) """ -Base.@kwdef struct Chen2022VelTypeRain{FT} <: ParametersType{FT} +Base.@kwdef struct Chen2022VelTypeRain{FT, N} <: ParametersType{FT} ρ0::FT - a1::FT - a2::FT - a3::FT + a::NTuple{N, FT} a3_pow::FT - b1::FT - b2::FT - b3::FT + b::NTuple{N, FT} b_ρ::FT - c1::FT - c2::FT - c3::FT + c::NTuple{N, FT} end Chen2022VelTypeRain(::Type{FT}) where {FT <: AbstractFloat} = @@ -221,21 +202,17 @@ Chen2022VelTypeRain(::Type{FT}) where {FT <: AbstractFloat} = function Chen2022VelTypeRain(td::CP.AbstractTOMLDict) name_map = (; :Chen2022_table_B1_q_coeff => :ρ0, - :Chen2022_table_B1_a1_coeff => :a1, - :Chen2022_table_B1_a2_coeff => :a2, - :Chen2022_table_B1_a3_coeff => :a3, + :Chen2022_table_B1_ai => :a, :Chen2022_table_B1_a3_pow_coeff => :a3_pow, - :Chen2022_table_B1_b1_coeff => :b1, - :Chen2022_table_B1_b2_coeff => :b2, - :Chen2022_table_B1_b3_coeff => :b3, + :Chen2022_table_B1_bi => :b, :Chen2022_table_B1_b_rho_coeff => :b_ρ, - :Chen2022_table_B1_c1_coeff => :c1, - :Chen2022_table_B1_c2_coeff => :c2, - :Chen2022_table_B1_c3_coeff => :c3, + :Chen2022_table_B1_ci => :c, ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") + # hack! + parameters = map(p -> p isa Vector ? Tuple(p) : p, parameters) FT = CP.float_type(td) - return Chen2022VelTypeRain{FT}(; parameters...) + return Chen2022VelTypeRain{FT, 3}(; parameters...) end """ diff --git a/test/Project.toml b/test/Project.toml index 541f70cf1f..63a215dc53 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,4 +12,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] +CLIMAParameters = "0.9" KernelAbstractions = "0.9" diff --git a/test/aerosol_activation_tests.jl b/test/aerosol_activation_tests.jl index 07f9ac5572..0e2534afe7 100644 --- a/test/aerosol_activation_tests.jl +++ b/test/aerosol_activation_tests.jl @@ -1,5 +1,6 @@ import Test as TT +import CLIMAParameters import CloudMicrophysics as CM import Thermodynamics as TD diff --git a/test/common_functions_tests.jl b/test/common_functions_tests.jl index 2b6632ea32..b36f665861 100644 --- a/test/common_functions_tests.jl +++ b/test/common_functions_tests.jl @@ -1,5 +1,6 @@ import Test as TT +import CLIMAParameters import Thermodynamics as TD import CloudMicrophysics as CM import CloudMicrophysics.Common as CO diff --git a/test/microphysics0M_tests.jl b/test/microphysics0M_tests.jl index 35c80cb832..a2cc4af08b 100644 --- a/test/microphysics0M_tests.jl +++ b/test/microphysics0M_tests.jl @@ -1,5 +1,6 @@ import Test as TT +import CLIMAParameters import Thermodynamics as TD import CloudMicrophysics.Parameters as CMP diff --git a/test/microphysics_noneq_tests.jl b/test/microphysics_noneq_tests.jl index b3c8af16d7..126fd326ba 100644 --- a/test/microphysics_noneq_tests.jl +++ b/test/microphysics_noneq_tests.jl @@ -1,5 +1,6 @@ import Test as TT +import CLIMAParameters import Thermodynamics as TD import CloudMicrophysics.Parameters as CMP diff --git a/test/p3_tests.jl b/test/p3_tests.jl index c5a0ecd68e..b5ddd2aa70 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -1,4 +1,5 @@ import Test as TT +import CLIMAParameters import CloudMicrophysics as CM import CloudMicrophysics.P3Scheme as P3 import CloudMicrophysics.Parameters as CMP diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 842eda4a9f..c892cef6d1 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -1,6 +1,7 @@ import Test as TT import BenchmarkTools as BT +import CLIMAParameters import Thermodynamics as TD import CloudMicrophysics.Common as CO From af494823d2cb4086f51fef0e15d3dad36ac9bc60 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Tue, 13 Feb 2024 16:53:37 -0800 Subject: [PATCH 2/3] Breaking release --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 20800d8b6c..cac9a0d805 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "CloudMicrophysics" uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" authors = ["Climate Modeling Alliance"] -version = "0.15.2" +version = "0.16.0" [deps] CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" From 3b316724d88f9dfedba47f57eac78512fefb89eb Mon Sep 17 00:00:00 2001 From: amylu00 Date: Tue, 30 Jan 2024 15:39:01 -0800 Subject: [PATCH 3/3] adding activity based dep to parcel --- docs/src/IceNucleation.md | 12 ++ docs/src/IceNucleationParcel0D.md | 38 ++++- docs/src/plots/activity_based_deposition.jl | 87 +++++++++++ parcel/Deposition_Nucleation.jl | 160 ++++++++++++++++++++ parcel/parcel.jl | 54 +++++-- 5 files changed, 329 insertions(+), 22 deletions(-) create mode 100644 docs/src/plots/activity_based_deposition.jl create mode 100644 parcel/Deposition_Nucleation.jl diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index ea3c24ae95..177aaa43f6 100644 --- a/docs/src/IceNucleation.md +++ b/docs/src/IceNucleation.md @@ -79,6 +79,18 @@ 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. +### Water activity based deposition nucleation plot +The following plot shows ``J`` as a function of ``\Delta a_w`` as compared to + figure 6 in Alpert et al 2013 and figure S5 in supplementary material of China et al 2017. Intent of this + plot is 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 +include("plots/activity_based_deposition.jl") +``` +![](water_activity_depo_nuc.svg) + + ## 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 diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md index 6be24cb69a..bf70908a40 100644 --- a/docs/src/IceNucleationParcel0D.md +++ b/docs/src/IceNucleationParcel0D.md @@ -198,18 +198,32 @@ where: ## Deposition Nucleation on dust particles There are multiple ways of running deposition nucleation in the parcel. - `"MohlerAF_Deposition"` will trigger an activated fraction approach from [Mohler2006](@cite). - `"MohlerRate_Deposition"` will trigger a nucleation rate approach from [Mohler2006](@cite). + `"MohlerAF_Deposition"` will trigger an activated fraction approach + from [Mohler2006](@cite). `"MohlerRate_Deposition"` will trigger a + nucleation rate approach from [Mohler2006](@cite). For both approaches, + there is no nucleation if saturation over ice exceeds 1.35 as conditions + above this value will result in nucleation in a different mode. + `"ActivityBasedDeposition"` will trigger a water activity based approach + from [Alpert2022](@cite). In this approach, ice production rate ``P_{ice, depo}`` + is calculated from +```math +\begin{equation} + P_{ice, depo} = \left[ \frac{dN_i}{dt} \right]_{depo} = J_{depo}\;A_{aero}\;N_{aero} + \label{eq:ActivityBasedDeposition_P_ice} +\end{equation} +``` +where ``N_{areo}`` is total number of unactiviated ice nucleating particles and + ``A_{aero}`` is surface area of those INP. The deposition nucleation methods are parameterized as described in [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/). ## Immersion Freezing Following the water activity based immersion freezing model (ABIFM), the ABIFM derived - nucleation rate coefficient, ``J_{immer}``, can be determined. The ice production rate,``P_{ice}``, + nucleation rate coefficient, ``J_{immer}``, can be determined. The ice production rate,``P_{ice, immer}``, per second via immersion freezing can then be calculating using ```math \begin{equation} - P_{ice} = \left[ \frac{dN_i}{dt} \right]_{immer} = J_{immer}A(N_{liq}) + P_{ice, immer} = \left[ \frac{dN_i}{dt} \right]_{immer} = J_{immer}A(N_{liq}) \label{eq:ABIFM_P_ice} \end{equation} ``` @@ -223,7 +237,7 @@ Homogeneous freezing follows the water-activity based model described in the The ice production rate from homogeneous freezing can then be determined: ```math \begin{equation} - P_{ice} = \left[ \frac{dN_i}{dt} \right]_{hom} = J_{hom}V(N_{liq}) + P_{ice, hom} = \left[ \frac{dN_i}{dt} \right]_{hom} = J_{hom}V(N_{liq}) \label{eq:hom_P_ice} \end{equation} ``` @@ -245,7 +259,7 @@ Between each run the water vapor specific humidity is changed, of the previous run. The prescribed vertical velocity is equal to 3.5 cm/s. Supersaturation is plotted for both liquid (solid lines) and ice (dashed lines). -The pale blue line uses the `"MohlerRate_Deposition"`approach. +The pale blue line uses the `"MohlerRate_Deposition"` approach. We only run it for the first GCM timestep because the rate approach requires the change in ice saturation over time. With the discontinuous jump in saturation, the parameterization is unable to determine a proper nucleation rate. When we force @@ -258,6 +272,18 @@ include("../../parcel/Tully_et_al_2023.jl") ``` ![](cirrus_box.svg) +The water activity based parameterization for deposition nucleation shows + similar outcomes when compared to the `"MohlerRate_Deposition"` approach. + Here, we run the parcel for 100 secs for all available aerosol types. The + solid lines correspond to the `"MohlerRate_Deposition"` approach while the + dashed lines correspond to `"ActivityBasedDeposition"`. Note that there + is no common aerosol type between the two parameterizations. + +```@example +include("../../parcel/Deposition_Nucleation.jl") +``` +![](deposition_nucleation.svg) + In the plots below, the parcel model is ran with only condensation (no ice or freezing) assuming either a monodisperse or a gamma distribution of droplets. It is compared to [Rogers1975](@cite). diff --git a/docs/src/plots/activity_based_deposition.jl b/docs/src/plots/activity_based_deposition.jl new file mode 100644 index 0000000000..551e16467a --- /dev/null +++ b/docs/src/plots/activity_based_deposition.jl @@ -0,0 +1,87 @@ +import Plots as PL + +import Thermodynamics as TD +import CloudMicrophysics as CM +import CloudMicrophysics.Common as CO +import CloudMicrophysics.HetIceNucleation as IN +import CloudMicrophysics.Parameters as CMP + +FT = Float32 +const tps = TD.Parameters.ThermodynamicsParameters(FT) +const feldspar = CMP.Feldspar(FT) +const ferrihydrite = CMP.Ferrihydrite(FT) +const kaolinite = CMP.Kaolinite(FT) + +# Initializing +Δa_w = range(FT(0), stop = FT(0.32), length = 50) # difference in solution and ice water activity +J_feldspar = @. IN.deposition_J(feldspar, Δa_w) # J in SI units +log10J_converted_feldspar = @. log10(J_feldspar * 1e-4) # converts J into cm^-2 s^-1 and takes log +J_ferrihydrite = @. IN.deposition_J(ferrihydrite, Δa_w) +log10J_converted_ferrihydrite = @. log10(J_ferrihydrite * 1e-4) +J_kaolinite = @. IN.deposition_J(kaolinite, Δa_w) +log10J_converted_kaolinite = @. log10(J_kaolinite * 1e-4) + +# Alpert et al 2022 Figure 6 +# https://doi.org/10.1039/d1ea00077b +# China et al 2017 Supplementary Figure S5 +# https://doi.org/10.1002/2016JD025817 + +# data read from Fig 6 in Alpert et al 2022 and +# China et al 2017 figure S5 +# using https://automeris.io/WebPlotDigitizer/ + +#! format: off +Alpert2022_Feldspar_Delta_a = [0.019459, 0.256216] +Alpert2022_Feldspar_log10J = [1.039735, 4.165563] + +Alpert2022_Ferrihydrite_Delta_a = [0.0989189, 0.336486] +Alpert2022_Ferrihydrite_log10J = [1.2781457, 4.21854] + +China2017_Delta_a = [0.01918, 0.02398, 0.02518, 0.03537, 0.07314, 0.07794, 0.10252, 0.10492, 0.1217, 0.1307, 0.15588, 0.16787, 0.20264, 0.23981, 0.256, 0.27638] +China2017_log10J = [-4.36923, -5.07692, -1.38462, -0.64615, 1.2, 1.35385, -0.58462, 1.90769, 2.06154, 2.24615, 2.52308, 0, 2.46154, 3.32308, 4, 4.43077, 5.26154] +#! format: on + +PL.plot( + Δa_w, + log10J_converted_feldspar, + label = "CM.jl Feldspar", + xlabel = "Delta a_w [unitless]", + ylabel = "log10(J) [cm^-2 s^-1]", + linestyle = :dash, + linecolor = :cyan, +) +PL.plot!( + Δa_w, + log10J_converted_ferrihydrite, + label = "CM.jl Ferrihydrite", + linestyle = :dash, + linecolor = :pink, +) +PL.plot!( + Δa_w, + log10J_converted_kaolinite, + label = "CM.jl Kaolinite", + linestyle = :dash, + linecolor = :orange, +) +PL.plot!( + Alpert2022_Feldspar_Delta_a, + Alpert2022_Feldspar_log10J, + linecolor = :blue, + label = "Alpert2022 Feldspar", +) +PL.plot!( + Alpert2022_Ferrihydrite_Delta_a, + Alpert2022_Ferrihydrite_log10J, + linecolor = :red, + label = "Alpert2022 Ferrihydrite", +) +PL.scatter!( + China2017_Delta_a, + China2017_log10J, + markercolor = :yellow, + markersize = 2, + label = "China2017 Kaolinite", +) + +PL.savefig("water_activity_depo_nuc.svg") diff --git a/parcel/Deposition_Nucleation.jl b/parcel/Deposition_Nucleation.jl new file mode 100644 index 0000000000..11cb530bae --- /dev/null +++ b/parcel/Deposition_Nucleation.jl @@ -0,0 +1,160 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK +import Thermodynamics as TD +import CloudMicrophysics as CM +import CLIMAParameters as CP + +# definition of the ODE problem for parcel model +include(joinpath(pkgdir(CM), "parcel", "parcel.jl")) +FT = Float32 +# get free parameters +tps = TD.Parameters.ThermodynamicsParameters(FT) +aps = CMP.AirProperties(FT) +wps = CMP.WaterProperties(FT) +ip = CMP.IceNucleationParameters(FT) + +# Constants +ρₗ = wps.ρw +R_v = TD.Parameters.R_v(tps) +R_d = TD.Parameters.R_d(tps) + +# Initial conditions +Nₐ = FT(2000 * 1e3) +Nₗ = FT(0) +Nᵢ = FT(0) +p₀ = FT(20000) +T₀ = FT(230) +qᵥ = FT(3.3e-4) +qₗ = FT(0) +qᵢ = FT(0) +x_sulph = FT(0) + +# Moisture dependent initial conditions +q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q) +ρₐ = TD.air_density(tps, ts) +Rₐ = TD.gas_constant_air(tps, q) +eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +e = qᵥ * p₀ * R_v / Rₐ + +Sₗ = FT(e / eₛ) +IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] + +ξ(T) = + TD.saturation_vapor_pressure(tps, T, TD.Liquid()) / + TD.saturation_vapor_pressure(tps, T, TD.Ice()) +S_i(T, S_liq) = ξ(T) * S_liq + +# Simulation parameters passed into ODE solver +r_nuc = FT(1.25e-6) # assumed size of nucleated particles +w = FT(3.5 * 1e-2) # updraft speed +α_m = FT(0.5) # accomodation coefficient +const_dt = FT(0.1) # model timestep +t_max = FT(100) +aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)] # aersol type for DustDeposition +aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)] # aersol type for DepositionNucleation +ice_nucleation_modes_list = + [["MohlerRate_Deposition"], ["ActivityBasedDeposition"]] +growth_modes = ["Deposition"] +droplet_size_distribution_list = [["Monodisperse"]] + +# Plotting +fig = MK.Figure(resolution = (800, 600)) +ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]") +ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") +ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]", xlabel = "time") +ax4 = MK.Axis(fig[2, 2], ylabel = "N_ice [m^-3]", xlabel = "time") + +for ice_nucleation_modes in ice_nucleation_modes_list + nuc_mode = ice_nucleation_modes[1] + droplet_size_distribution = droplet_size_distribution_list[1] + + if nuc_mode == "MohlerRate_Deposition" + for aerosol in aerosol_1 + p = (; + wps, + aps, + tps, + ip, + const_dt, + r_nuc, + w, + α_m, + aerosol, + ice_nucleation_modes, + growth_modes, + droplet_size_distribution, + ) + + # solve ODE + sol = run_parcel(IC, FT(0), t_max, p) + + # Plot results + if aerosol == CMP.DesertDust(FT) + aero_label = "DesertDust" + elseif aerosol == CMP.ArizonaTestDust(FT) + aero_label = "ArizonaTestDust" + end + MK.lines!( # saturation + ax1, + sol.t, + S_i.(sol[3, :], sol[1, :]), + label = aero_label, + ) + MK.lines!(ax2, sol.t, sol[3, :]) # temperature + MK.lines!(ax3, sol.t, sol[6, :] * 1e3) # q_ice + MK.lines!(ax4, sol.t, sol[9, :]) # N_ice + end + + elseif nuc_mode == "ActivityBasedDeposition" + for aerosol in aerosol_2 + p = (; + wps, + aps, + tps, + ip, + const_dt, + r_nuc, + w, + α_m, + aerosol, + ice_nucleation_modes, + growth_modes, + droplet_size_distribution, + ) + + # solve ODE + sol = run_parcel(IC, FT(0), t_max, p) + + # Plot results + if aerosol == CMP.Feldspar(FT) + aero_label = "Feldspar" + elseif aerosol == CMP.Ferrihydrite(FT) + aero_label = "Ferrihydrite" + elseif aerosol == CMP.Kaolinite(FT) + aero_label = "Kaolinite" + end + MK.lines!( # saturation + ax1, + sol.t, + S_i.(sol[3, :], sol[1, :]), + label = aero_label, + linestyle = :dash, + ) + MK.lines!(ax2, sol.t, sol[3, :], linestyle = :dash) # temperature + MK.lines!(ax3, sol.t, sol[6, :] * 1e3, linestyle = :dash) # q_ice + MK.lines!(ax4, sol.t, sol[9, :], linestyle = :dash) # N_ice + end + end +end + +MK.axislegend( + ax1, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 3, + position = :lt, +) + +MK.save("deposition_nucleation.svg", fig) diff --git a/parcel/parcel.jl b/parcel/parcel.jl index c11b2b97d0..afd017719e 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -76,7 +76,7 @@ function parcel_model(dY, Y, p, t) dN_act_dt_depo = FT(0) dqi_dt_new_depo = FT(0) if "MohlerAF_Deposition" in ice_nucleation_modes - if S_i > ip.deposition.Sᵢ_max + if S_i >= ip.deposition.Sᵢ_max println("Supersaturation exceeds the allowed value. No dust will be activated.") AF = FT(0) else @@ -89,9 +89,8 @@ function parcel_model(dY, Y, p, t) end dN_act_dt_depo = max(FT(0), AF * N_aer - N_ice) / const_dt dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air - end - if "MohlerRate_Deposition" in ice_nucleation_modes - if S_i > ip.deposition.Sᵢ_max + elseif "MohlerRate_Deposition" in ice_nucleation_modes + if S_i >= ip.deposition.Sᵢ_max println("Supersaturation exceeds the allowed value. No dust will be activated.") dN_act_dt_depo = FT(0) else @@ -99,6 +98,15 @@ function parcel_model(dY, Y, p, t) end dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air end + if "ActivityBasedDeposition" in ice_nucleation_modes + Δa_w = CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T) + J_deposition = CMI_het.deposition_J(aerosol, Δa_w) + if "Monodisperse" in droplet_size_distribution + A_aer = 4 * FT(π) * r_nuc^2 + dN_act_dt_depo = max(FT(0), J_deposition * N_aer * A_aer) + dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air + end + end dN_act_dt_immersion = FT(0) dqi_dt_new_immers = FT(0) @@ -114,8 +122,8 @@ function parcel_model(dY, Y, p, t) #A = N_liq* λ^2 r_l = 2 / λ end - A_aer = 4 * FT(π) * r_l^2 - dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer) + A_l = 4 * FT(π) * r_l^2 + dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_l) dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air end @@ -282,42 +290,56 @@ function run_parcel(IC, t_0, t_end, p) println(" ") println("Running parcel model with: ") + print("Ice nucleation modes: ") if "MohlerAF_Deposition" in ice_nucleation_modes - print("Deposition on dust particles using AF ") + print("\n Deposition on dust particles using AF ") + print("(Note that this mode only runs for monodisperse size distributions.) ") timestepper = ODE.Euler() end if "MohlerRate_Deposition" in ice_nucleation_modes - print("Deposition nucleation based on rate eqn in Mohler 2006 ") + print("\n Deposition nucleation based on rate eqn in Mohler 2006 ") + print("(Note that this mode only runs for monodisperse size distributions.) ") timestepper = ODE.Euler() end + if "ActivityBasedDeposition" in ice_nucleation_modes + print("\n Water activity based deposition nucleation ") + print("(Note that this mode only runs for monodisperse size distributions.) ") + end if "ImmersionFreezing" in ice_nucleation_modes - print("Immersion freezing ") + print("\n Immersion freezing ") + print("(Note that this mode only runs for monodisperse and gamma size distributions.) ") end if "HomogeneousFreezing" in ice_nucleation_modes - print("Homogeneous freezing ") + print("\n Homogeneous freezing ") end print("\n") + print("Growth modes: ") if "Condensation" in growth_modes - print("Condensation on liquid droplets",) + print("\n Condensation on liquid droplets",) + print("(Note that this mode only runs for monodisperse and gamma size distributions.) ") end if "Condensation_DSD" in growth_modes - print("Condensation on liquid droplets") + print("\n Condensation on liquid droplets") end if "Deposition" in growth_modes - print("Deposition on ice crystals") + print("\n Deposition on ice crystals") end print("\n") + print("Size distribution modes: ") if "Monodisperse" in droplet_size_distribution - print("Monodisperse size distribution") + print("\n Monodisperse size distribution") end if "Gamma" in droplet_size_distribution - print("Gamma size distribution") + print("\n Gamma size distribution") end if "Lognormal" in droplet_size_distribution - print("Lognormal size distribution") + print("\n Lognormal size distribution") + end + if "MohlerAF_Deposition" in ice_nucleation_modes || "MohlerRate_Deposition" in ice_nucleation_modes || "ActivityBasedDeposition" in ice_nucleation_modes + print("\nWARNING: Multiple deposition nucleation parameterizations chosen to run in one simulation. Parcel will only run MohlerAF_Deposition for now.") end println(" ")