diff --git a/Project.toml b/Project.toml index a52b5562c8..db52f01e2e 100644 --- a/Project.toml +++ b/Project.toml @@ -14,7 +14,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] -CLIMAParameters = ">=0.7.9" +CLIMAParameters = ">=0.7.12" CUDAKernels = "0.4" DocStringExtensions = "0.8, 0.9" ForwardDiff = "0.10" diff --git a/docs/Project.toml b/docs/Project.toml index 647a371981..5e2275fb09 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,3 +7,5 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" +[compat] +CLIMAParameters = ">=0.7.12" diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index 9d34a484ba..1e8462987e 100644 --- a/docs/src/IceNucleation.md +++ b/docs/src/IceNucleation.md @@ -8,7 +8,7 @@ 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) + The parameterization for immersion freezing is an implementation of [KnopfAlpert2013](@cite) and is valid for droplets containing sulphuric acid. !!! note @@ -45,16 +45,16 @@ Both parameters are dependent on aerosol properties and temperature. (either a second deposition or other condensation type mode). ## ABIFM for Sulphuric Acid Containing Droplets -Water Activity-Based Immersion Freezing Model (ABFIM) +Water Activity-Based Immersion Freezing Model (ABFIM) is a method of parameterizing immersion freezing inspired by the time-dependent - classical nucleation theory (CNT). More on CNT can be found in [Karthika2016](@cite). - The nucleation rate coefficient, ``J``, describes the number of ice nuclei formed per unit area + classical nucleation theory (CNT). More on CNT can be found in [Karthika2016](@cite). + The nucleation rate coefficient, ``J``, describes the number of ice nuclei formed per unit area per unit time and can be determined by the water activity, ``a_w``. This parameterization follows [KnopfAlpert2013](@cite), [Koop2002](@cite), [MurphyKoop2005](@cite), and [Luo1995](@cite). In this model, aerosols are assumed to contain an insoluble and soluble material. When immersed in water, the soluble material diffuses into the liquid water to create a sulphuric acid solution. -Using empirical coefficients, ``m`` and ``c``, from [KnopfAlpert2013](@cite), +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} @@ -63,10 +63,10 @@ Using empirical coefficients, ``m`` and ``c``, from [KnopfAlpert2013](@cite), ``` !!! note - Our source code for the nucleation rate coefficient returns + Our source code for the nucleation rate coefficient returns ``J`` in base SI units. -``\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)``. From [Koop2002](@cite), +``\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)``. From [Koop2002](@cite), ```math \begin{equation} a_w = \frac{p_{sol}}{p_{sat}} @@ -96,8 +96,8 @@ Once ``J_{het}`` is calculated, it can be used to determine the ice production r P_{ice} = J_{het}A(N_{tot}-N_{ice}) \end{equation} ``` -where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number - of ice nuclei, and ``N_{ice}`` is number of ice crystals already in the system. +where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number + of ice nuclei, and ``N_{ice}`` is number of ice crystals already in the system. ## Homogeneous Freezing for Sulphuric Acid Containing Droplets Homogeneous freezing occurs when supercooled liquid droplets freeze on their own. @@ -144,13 +144,13 @@ 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 +# 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(dust_type, Delta_a[it]) + J[it] = IN.ABIFM_J(prs, dust_type, Delta_a[it]) global it += 1 end log10J_converted = @. log10(J*1e-4) diff --git a/parcel/Project.toml b/parcel/Project.toml index ad7d47e792..d50d0d0ff8 100644 --- a/parcel/Project.toml +++ b/parcel/Project.toml @@ -8,3 +8,5 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" +[compat] +CLIMAParameters = ">=0.7.12" diff --git a/src/Common.jl b/src/Common.jl index b1a81f1617..25d51265eb 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -124,24 +124,41 @@ function logistic_function_integral(x::FT, x_0::FT, k::FT) where {FT <: Real} end """ - H2SO4_soln_saturation_vapor_pressure(x, T) + H2SO4_soln_saturation_vapor_pressure(prs, x, T) + - `prs` - a set with free parameters - `x` - wt percent sulphuric acid [unitless] - `T` - air temperature [K]. Returns the saturation vapor pressure above a sulphuric acid solution droplet in Pa. `x` is, for example, 0.1 if droplets are 10 percent sulphuric acid by weight """ -function H2SO4_soln_saturation_vapor_pressure(x::FT, T::FT) where {FT <: Real} +function H2SO4_soln_saturation_vapor_pressure( + prs::APS, + x::FT, + T::FT, +) where {FT <: Real} + + T_max::FT = CMP.H2SO4_sol_T_max(prs) + T_min::FT = CMP.H2SO4_sol_T_min(prs) + w_2::FT = CMP.H2SO4_sol_w_2(prs) + + c1::FT = CMP.H2SO4_sol_c1(prs) + c2::FT = CMP.H2SO4_sol_c2(prs) + c3::FT = CMP.H2SO4_sol_c3(prs) + c4::FT = CMP.H2SO4_sol_c4(prs) + c5::FT = CMP.H2SO4_sol_c5(prs) + c6::FT = CMP.H2SO4_sol_c6(prs) + c7::FT = CMP.H2SO4_sol_c7(prs) - @assert T < FT(235) - @assert T > FT(185) + @assert T < T_max + @assert T > T_min - w_h = 1.4408 * x + w_h = w_2 * x p_sol = exp( - 23.306 - 5.3465 * x + 12 * x * w_h - 8.19 * x * w_h^2 + - (-5814 + 928.9 * x - 1876.7 * x * w_h) / T, + c1 - c2 * x + c3 * x * w_h - c4 * x * w_h^2 + + (c5 + c6 * x - c7 * x * w_h) / T, ) * 100 # * 100 converts mbar --> Pa return p_sol end @@ -160,7 +177,7 @@ function Delta_a_w(prs::APS, x::FT, T::FT) where {FT <: Real} thermo_params = CMP.thermodynamics_params(prs) - p_sol = H2SO4_soln_saturation_vapor_pressure(x, T) + p_sol = H2SO4_soln_saturation_vapor_pressure(prs, x, T) p_sat = TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid()) p_ice = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index 9022dd0c5d..2759240177 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -22,14 +22,12 @@ S0_cold(prs::APS, ::CT.DesertDustType) = CMP.S0_cold_DD_Mohler2006(prs) a_warm(prs::APS, ::CT.DesertDustType) = CMP.a_warm_DD_Mohler2006(prs) a_cold(prs::APS, ::CT.DesertDustType) = CMP.a_cold_DD_Mohler2006(prs) -J_het_coeff_m(::CT.DesertDustType) = 22.62 -J_het_coeff_c(::CT.DesertDustType) = -1.35 - -J_het_coeff_m(::CT.KaoliniteType) = 54.58834 -J_het_coeff_c(::CT.KaoliniteType) = -10.54758 - -J_het_coeff_m(::CT.IlliteType) = 54.48075 -J_het_coeff_c(::CT.IlliteType) = -10.66873 +J_het_m(prs::APS, ::CT.DesertDustType) = CMP.J_ABIFM_m_KA2013_DesertDust(prs) +J_het_c(prs::APS, ::CT.DesertDustType) = CMP.J_ABIFM_c_KA2013_DesertDust(prs) +J_het_m(prs::APS, ::CT.KaoliniteType) = CMP.J_ABIFM_m_KA2013_Kaolinite(prs) +J_het_c(prs::APS, ::CT.KaoliniteType) = CMP.J_ABIFM_c_KA2013_Kaolinite(prs) +J_het_m(prs::APS, ::CT.IlliteType) = CMP.J_ABIFM_m_KA2013_Illite(prs) +J_het_c(prs::APS, ::CT.IlliteType) = CMP.J_ABIFM_c_KA2013_Illite(prs) """ dust_activated_number_fraction(prs, Si, T, dust_type) @@ -66,23 +64,31 @@ function dust_activated_number_fraction( end """ - ABIFM_J(dust_type, Δa_w) + ABIFM_J(prs, dust_type, Δa_w) - - `dust_type` - choosing aerosol type + - `prs` - set with free parameters + - `dust_type` - aerosol type - `Δa_w` - change in water activity [unitless]. -Returns the immersion freezing nucleation rate coefficient, `J`, in m^-2 s^-1 for sulphuric acid containing solutions. -For other solutions, p_sol should be adjusted accordingly. Delta_a_w can be found using the Delta_a_w function in Common.jl. -`m` and `c` constants are taken from Knopf & Alpert 2013. +Returns the immersion freezing nucleation rate coefficient, `J`, in m^-2 s^-1 +for sulphuric acid solutions. +For other solutions, p_sol should be adjusted accordingly. +Delta_a_w can be found using the Delta_a_w function in Common.jl. +The free parameters `m` and `c` are taken from Knopf & Alpert 2013 +see DOI: 10.1039/C3FD00035D """ -function ABIFM_J(dust_type::CT.AbstractAerosolType, Δa_w::FT) where {FT <: Real} +function ABIFM_J( + prs::APS, + dust_type::CT.AbstractAerosolType, + Δa_w::FT, +) where {FT <: Real} - m = J_het_coeff_m(dust_type) - c = J_het_coeff_c(dust_type) + m::FT = J_het_m(prs, dust_type) + c::FT = J_het_c(prs, dust_type) - logJ = m * Δa_w + c + logJ::FT = m * Δa_w + c - return max(0, 10^logJ * 100^2) # converts cm^-2 s^-1 to m^-2 s^-1 + return max(FT(0), FT(10)^logJ * FT(1e4)) # converts cm^-2 s^-1 to m^-2 s^-1 end end # end module @@ -101,23 +107,31 @@ const APS = CMP.AbstractCloudMicrophysicsParameters export homogeneous_J """ - homogeneous_J(Δa_w) + homogeneous_J(prs, Δa_w) - - `Δa_w` - change in water activity + - `prs` - a set with free parameters + - `Δa_w` - change in water activity [-]. -Returns the homogeneous freezing nucleation rate coefficient, `J`, in m^-3 s^-1 for sulphuric -acid containing solutions. Parameterization based off Koop 2000. +Returns the homogeneous freezing nucleation rate coefficient, +`J`, in m^-3 s^-1 for sulphuric acid solutions. +Parameterization based on Koop 2000, DOI: 10.1038/35020537. Delta_a_w can be found using the Delta_a_w function in Common.jl. """ -function homogeneous_J(Δa_w::FT) where {FT <: Real} +function homogeneous_J(prs::APS, Δa_w::FT) where {FT <: Real} + + Δa_w_min::FT = CMP.Koop2000_min_delta_aw(prs) + Δa_w_max::FT = CMP.Koop2000_max_delta_aw(prs) + c1::FT = CMP.Koop2000_J_hom_c1(prs) + c2::FT = CMP.Koop2000_J_hom_c2(prs) + c3::FT = CMP.Koop2000_J_hom_c3(prs) + c4::FT = CMP.Koop2000_J_hom_c4(prs) - @assert Δa_w > 0.26 - @assert Δa_w < 0.34 + @assert Δa_w > Δa_w_min + @assert Δa_w < Δa_w_max - logJ = -906.7 + 8502 * Δa_w - 26924 * Δa_w^2 + 29180 * Δa_w^3 - J = 10^(logJ) + logJ::FT = c1 + c2 * Δa_w - c3 * Δa_w^2 + c4 * Δa_w^3 - return J * 1e6 + return FT(10)^(logJ) * 1e6 end end # end module diff --git a/src/Parameters.jl b/src/Parameters.jl index f8d0d5bdd9..f6499177bf 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -164,6 +164,28 @@ Base.@kwdef struct CloudMicrophysicsParameters{FT, TP} <: S0_cold_DD_Mohler2006::FT a_warm_DD_Mohler2006::FT a_cold_DD_Mohler2006::FT + J_ABIFM_m_KA2013_DesertDust::FT + J_ABIFM_c_KA2013_DesertDust::FT + J_ABIFM_m_KA2013_Kaolinite::FT + J_ABIFM_c_KA2013_Kaolinite::FT + J_ABIFM_m_KA2013_Illite::FT + J_ABIFM_c_KA2013_Illite::FT + Koop2000_min_delta_aw::FT + Koop2000_max_delta_aw::FT + Koop2000_J_hom_c1::FT + Koop2000_J_hom_c2::FT + Koop2000_J_hom_c3::FT + Koop2000_J_hom_c4::FT + H2SO4_sol_T_max::FT + H2SO4_sol_T_min::FT + H2SO4_sol_w_2::FT + H2SO4_sol_c1::FT + H2SO4_sol_c2::FT + H2SO4_sol_c3::FT + H2SO4_sol_c4::FT + H2SO4_sol_c5::FT + H2SO4_sol_c6::FT + H2SO4_sol_c7::FT molmass_seasalt::FT rho_seasalt::FT osm_coeff_seasalt::FT diff --git a/test/common_functions_tests.jl b/test/common_functions_tests.jl index c78105e6a6..81678fd218 100644 --- a/test/common_functions_tests.jl +++ b/test/common_functions_tests.jl @@ -63,18 +63,20 @@ function test_H2SO4_soln_saturation_vapor_pressure(FT) x_sulph = FT(0.1) # If T out of range - TT.@test_throws AssertionError("T < FT(235)") CO.H2SO4_soln_saturation_vapor_pressure( + TT.@test_throws AssertionError("T < T_max") CO.H2SO4_soln_saturation_vapor_pressure( + prs, x_sulph, T_too_warm, ) - TT.@test_throws AssertionError("T > FT(185)") CO.H2SO4_soln_saturation_vapor_pressure( + TT.@test_throws AssertionError("T > T_min") CO.H2SO4_soln_saturation_vapor_pressure( + prs, x_sulph, T_too_cold, ) # p_sol should be higher at warmer temperatures - TT.@test CO.H2SO4_soln_saturation_vapor_pressure(x_sulph, T_warm) > - CO.H2SO4_soln_saturation_vapor_pressure(x_sulph, T_cold) + TT.@test CO.H2SO4_soln_saturation_vapor_pressure(prs, x_sulph, T_warm) > + CO.H2SO4_soln_saturation_vapor_pressure(prs, x_sulph, T_cold) end end @@ -97,15 +99,9 @@ end println("Testing Float64") test_H2SO4_soln_saturation_vapor_pressure(Float64) - - println("Testing Float32") test_H2SO4_soln_saturation_vapor_pressure(Float32) - - println("Testing Float64") test_Delta_a_w(Float64) - - println("Testing Float32") test_Delta_a_w(Float32) diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 7b206ee33c..6ca3215eb5 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -270,6 +270,7 @@ end end @kernel function test_Common_H2SO4_soln_saturation_vapor_pressure_kernel!( + prs, output::AbstractArray{FT}, x_sulph, T, @@ -278,7 +279,8 @@ end i = @index(Group, Linear) @inbounds begin - output[i] = CO.H2SO4_soln_saturation_vapor_pressure(x_sulph[i], T[i]) + output[i] = + CO.H2SO4_soln_saturation_vapor_pressure(prs, x_sulph[i], T[i]) end end @@ -297,6 +299,7 @@ end end @kernel function test_IceNucleation_ABIFM_J_kernel!( + prs, output::AbstractArray{FT}, Delta_a_w, ) where {FT} @@ -304,12 +307,13 @@ end i = @index(Group, Linear) @inbounds begin - output[1] = CMI_het.ABIFM_J(kaolinite, Delta_a_w[1]) - output[2] = CMI_het.ABIFM_J(illite, Delta_a_w[2]) + output[1] = CMI_het.ABIFM_J(prs, kaolinite, Delta_a_w[1]) + output[2] = CMI_het.ABIFM_J(prs, illite, Delta_a_w[2]) end end @kernel function test_IceNucleation_homogeneous_J_kernel!( + prs, output::AbstractArray{FT}, Delta_a_w, ) where {FT} @@ -317,7 +321,7 @@ end i = @index(Group, Linear) @inbounds begin - output[1] = CMI_hom.homogeneous_J(Delta_a_w[1]) + output[1] = CMI_hom.homogeneous_J(prs, Delta_a_w[1]) end end @@ -582,7 +586,7 @@ function test_gpu(FT) dev, work_groups, ) - event = kernel!(output, x_sulph, T, ndrange = ndrange) + event = kernel!(make_prs(FT), output, x_sulph, T, ndrange = ndrange) wait(dev, event) # test H2SO4_soln_saturation_vapor_pressure is callable and returns a reasonable value @@ -619,7 +623,7 @@ function test_gpu(FT) Delta_a_w = ArrayType([FT(0.16), FT(0.15)]) kernel! = test_IceNucleation_ABIFM_J_kernel!(dev, work_groups) - event = kernel!(output, Delta_a_w, ndrange = ndrange) + event = kernel!(make_prs(FT), output, Delta_a_w, ndrange = ndrange) wait(dev, event) # test if ABIFM_J is callable and returns reasonable values @@ -639,7 +643,7 @@ function test_gpu(FT) Delta_a_w = ArrayType([FT(0.2907389666103033)]) kernel! = test_IceNucleation_homogeneous_J_kernel!(dev, work_groups) - event = kernel!(output, Delta_a_w, ndrange = ndrange) + event = kernel!(make_prs(FT), output, Delta_a_w, ndrange = ndrange) wait(dev, event) # test homogeneous_J is callable and returns a reasonable value diff --git a/test/heterogeneous_ice_nucleation_tests.jl b/test/heterogeneous_ice_nucleation_tests.jl index e7567a36d6..2b6db419c1 100644 --- a/test/heterogeneous_ice_nucleation_tests.jl +++ b/test/heterogeneous_ice_nucleation_tests.jl @@ -87,8 +87,15 @@ function test_heterogeneous_ice_nucleation(FT) # higher nucleation rate at colder temperatures for dust in [CMT.IlliteType(), CMT.KaoliniteType(), CMT.DesertDustType()] - TT.@test CMI_het.ABIFM_J(dust, CO.Delta_a_w(prs, x_sulph, T_cold)) > - CMI_het.ABIFM_J(dust, CO.Delta_a_w(prs, x_sulph, T_warm)) + TT.@test CMI_het.ABIFM_J( + prs, + dust, + CO.Delta_a_w(prs, x_sulph, T_cold), + ) > CMI_het.ABIFM_J( + prs, + dust, + CO.Delta_a_w(prs, x_sulph, T_warm), + ) end end end @@ -96,6 +103,5 @@ end println("Testing Float64") test_heterogeneous_ice_nucleation(Float64) - println("Testing Float32") test_heterogeneous_ice_nucleation(Float32) diff --git a/test/homogeneous_ice_nucleation_tests.jl b/test/homogeneous_ice_nucleation_tests.jl index 676b7fd6f0..2fc48e5949 100644 --- a/test/homogeneous_ice_nucleation_tests.jl +++ b/test/homogeneous_ice_nucleation_tests.jl @@ -5,7 +5,7 @@ import CloudMicrophysics as CM import CLIMAParameters as CP const CO = CM.Common -const CMI_hom = CM.HomIceNucleation +const CMH = CM.HomIceNucleation include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) @@ -23,24 +23,25 @@ function test_homogeneous_J(FT) Δa_w_too_small = FT(0.25) Δa_w_too_large = FT(0.35) + # higher nucleation rate at colder temperatures - TT.@test CMI_hom.homogeneous_J(CO.Delta_a_w(prs, x_sulph, T_cold)) > - CMI_hom.homogeneous_J(CO.Delta_a_w(prs, x_sulph, T_warm)) + TT.@test CMH.homogeneous_J(prs, CO.Delta_a_w(prs, x_sulph, T_cold)) > + CMH.homogeneous_J(prs, CO.Delta_a_w(prs, x_sulph, T_warm)) # If Δa_w out of range - TT.@test_throws AssertionError("Δa_w > 0.26") CMI_hom.homogeneous_J( + TT.@test_throws AssertionError("Δa_w > Δa_w_min") CMH.homogeneous_J( + prs, Δa_w_too_small, ) - TT.@test_throws AssertionError("Δa_w < 0.34") CMI_hom.homogeneous_J( + TT.@test_throws AssertionError("Δa_w < Δa_w_max") CMH.homogeneous_J( + prs, Δa_w_too_large, ) end end - println("Testing Float64") test_homogeneous_J(Float64) - println("Testing Float32") test_homogeneous_J(Float32) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index b830407575..8d03bab51f 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -108,7 +108,7 @@ function benchmark_test(FT) # Common bench_press( CO.H2SO4_soln_saturation_vapor_pressure, - (x_sulph, T_air_cold), + (prs, x_sulph, T_air_cold), 50, ) bench_press(CO.Delta_a_w, (prs, x_sulph, T_air_cold), 230) @@ -119,8 +119,8 @@ function benchmark_test(FT) (prs, S_ice, T_air_2, dust), 50, ) - bench_press(CMI_het.ABIFM_J, (dust, Delta_a_w), 230) - bench_press(CMI_hom.homogeneous_J, (Delta_a_w), 230) + bench_press(CMI_het.ABIFM_J, (prs, dust, Delta_a_w), 230) + bench_press(CMI_hom.homogeneous_J, (prs, Delta_a_w), 230) # non-equilibrium bench_press(CMN.τ_relax, (prs, liquid), 10)