diff --git a/docs/src/API.md b/docs/src/API.md index 0711e8470f..52d7be55be 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -108,7 +108,8 @@ HetIceNucleation.INP_concentration_mean # Homogeneous ice nucleation ```@docs HomIceNucleation -HomIceNucleation.homogeneous_J +HomIceNucleation.homogeneous_J_cubic +HomIceNucleation.homogeneous_J_linear ``` # Common utility functions diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md index 9eeffd9b17..4dcec044e3 100644 --- a/docs/src/IceNucleation.md +++ b/docs/src/IceNucleation.md @@ -243,6 +243,14 @@ Multiple sulphuric acid concentrations, ``x``, values that are obtained from pure water droplets. Though CliMA lines look far from the Spichtinger 2023 line, the lines seem to move closer as x approaches 0. +### Linear fit homogeneous freezing +To avoid complications in extrapolating the Koop 2000 parameterization outside of the valide temperature range, a linear fit was created. + +```@example +include("plots/linear_HOM_J.jl") +``` +![](linear_HOM_J.svg) + ## Water Activity Based vs P3 Ice Nucleation Parameterizations The water activity based models are compared with P3 ice nucleation parameterizations as described in [MorrisonMilbrandt2015](@cite) using an adiabatic parcel model with depositional growth. The diff --git a/docs/src/plots/HomFreezingPlots.jl b/docs/src/plots/HomFreezingPlots.jl index e202631a6b..aaac07e471 100644 --- a/docs/src/plots/HomFreezingPlots.jl +++ b/docs/src/plots/HomFreezingPlots.jl @@ -31,16 +31,16 @@ x_sulph = Vector{FT}([0.03, 0.04, 0.06]) # wt% sulphuric acid in dropl T in T_range ] -J1 = @. CMI.homogeneous_J(ip.homogeneous, Δa1) -J2 = @. CMI.homogeneous_J(ip.homogeneous, Δa2) -J3 = @. CMI.homogeneous_J(ip.homogeneous, Δa3) +J1 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa1) +J2 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa2) +J3 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa3) log10J_1 = @. log10(J1) log10J_2 = @. log10(J2) log10J_3 = @. log10(J3) Δa_range = range(0.27, stop = 0.32, length = 50) -J_given_Δa = @. CMI.homogeneous_J(ip.homogeneous, Δa_range) +J_given_Δa = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa_range) #! format: off # Spichtinger et al 2023 diff --git a/docs/src/plots/linear_HOM_J.jl b/docs/src/plots/linear_HOM_J.jl index 3e425d2d1e..8db26074e7 100644 --- a/docs/src/plots/linear_HOM_J.jl +++ b/docs/src/plots/linear_HOM_J.jl @@ -1,64 +1,33 @@ -import OrdinaryDiffEq as ODE import CairoMakie as MK -import Thermodynamics as TD import CloudMicrophysics as CM -import ClimaParams as CP -import Distributions as DS import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.HomIceNucleation as CMH -const CMO = CM.Common -FT = Float32 -const tps = TD.Parameters.ThermodynamicsParameters(FT) +const ip = CMP.IceNucleationParameters(FT) -# From Atkinson, Murray, and O'Sullivan 2016 -# (DOI:10.1021/acs.jpca.6b03843) -# Using overall volume fit -function Atkinson_J(T) - lnJ = -4.0106 * T + 963.706 - return exp(lnJ) -end +# Koop2000 parameterization +Koop_Δa = collect(0.26:0.0025:0.34) +KoopJ = @. CMH.homogeneous_J_cubic(ip.homogeneous, Koop_Δa) # SI units +log10KoopJ = @. log10(KoopJ .* 1e-6) -# Initializing -T_range = range(FT(234.0), stop = 237, length = 50) # air temperature -Sₗ = FT(1.1) -eₛ = [TD.saturation_vapor_pressure(tps, T, TD.Liquid()) for T in T_range] -e_range = Sₗ .* eₛ - -Δa = zeros(FT, length(T_range)) - -# Solving for Δa and J values -J = [Atkinson_J(T) for T in T_range] # [cm^-3 s^-1] -log10J = @. log10(J) -for (i, T) in enumerate(T_range) - Δa[i] = CMO.a_w_eT(tps, e_range[i], T) - CMO.a_w_ice(tps, T) -end - -# Obtaining a linear parameterization for J in terms of Δa -slope = (log10J[end] - log10J[1]) / (Δa[end] - Δa[1]) -intercept = slope * (Δa[1]) + log10J[1] +# Linear fit on Koop 2000 +M = [ones(length(Koop_Δa)) Koop_Δa] +linear_coeffs = M \ log10KoopJ +new_log10J = [linear_coeffs[2] * Delta_a + linear_coeffs[1] for Delta_a in Koop_Δa] # Plotting J vs Δa -fig = MK.Figure(resolution = (800, 700)) +fig = MK.Figure(resolution = (800, 600), fontsize = 18) ax1 = MK.Axis( fig[1, 1], - ylabel = "log10(J), J = [cm^-3 s^-1]", - xlabel = "Temperature [K]", - title = "log10(J) vs T", -) -ax2 = MK.Axis( - fig[2, 1], xlabel = "Δa_w [-]", ylabel = "log10(J), J = [cm^-3 s^-1]", - title = "log10(J) vs Δa_w", + title = "Linear fit to Koop2000", ) -MK.lines!(ax1, T_range, log10J, label = "ln(J) = -4.0106 T + 963.706",) -MK.lines!(ax2, Δa, log10J, - label = "log10(J) = $slope Δa + $intercept", -) +MK.lines!(ax1, Koop_Δa, new_log10J, label = "log10(J) = $(linear_coeffs[2]) Δa + $(linear_coeffs[1])") +MK.lines!(ax1, Koop_Δa, log10KoopJ, label = "Koop 2000 Parameterization") -MK.axislegend(ax1, position = :rt, labelsize = 18.0f0) -MK.axislegend(ax2, position = :rb, labelsize = 18.0f0) +MK.axislegend(ax1, position = :rb) #! format: on diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index f2ad1a1924..6874fc2508 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -245,11 +245,11 @@ end Returns the homogeneous freezing nucleation rate coefficient, `J`, in m^-3 s^-1 for sulphuric acid solutions. -Parameterization based on Atkinson et al 2016, DOI:10.1021/acs.jpca.6b03847. +Parameterization derived from a linear fit of the Koop 2000 parameterization, DOI: 10.1038/35020537. """ -function homogeneous_J_linear(ip::CMP.Atkinson2016, Δa_w::FT) where {FT} +function homogeneous_J_linear(ip::CMP.Koop2000, Δa_w::FT) where {FT} - logJ::FT = ip.c₁ * Δa_w + ip.c₂ + logJ::FT = ip.linear_c₁ * Δa_w + ip.linear_c₂ return FT(10)^(logJ) * 1e6 end diff --git a/src/parameters/IceNucleation.jl b/src/parameters/IceNucleation.jl index 02a97572ad..cd37ab239f 100644 --- a/src/parameters/IceNucleation.jl +++ b/src/parameters/IceNucleation.jl @@ -49,6 +49,10 @@ Base.@kwdef struct Koop2000{FT} <: ParametersType{FT} c₃::FT "coefficient [-]" c₄::FT + "coefficient [-]" + linear_c₁::FT + "coefficient [-]" + linear_c₂::FT end function Koop2000(td::CP.AbstractTOMLDict) @@ -59,6 +63,8 @@ function Koop2000(td::CP.AbstractTOMLDict) :Koop2000_J_hom_coeff2 => :c₂, :Koop2000_J_hom_coeff3 => :c₃, :Koop2000_J_hom_coeff4 => :c₄, + :Linear_J_hom_coeff1 => :linear_c₁, + :Linear_J_hom_coeff1 => :linear_c₂, ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index ae2e392c8d..22735aaeeb 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -629,7 +629,7 @@ end end end -@kernel function IceNucleation_homogeneous_J_kernel!( +@kernel function IceNucleation_homogeneous_J_cubic_kernel!( ip, output::AbstractArray{FT}, Delta_a_w, @@ -638,7 +638,20 @@ end i = @index(Group, Linear) @inbounds begin - output[1] = CMI_hom.homogeneous_J(ip.homogeneous, Delta_a_w[1]) + output[1] = CMI_hom.homogeneous_J_cubic(ip.homogeneous, Delta_a_w[1]) + end +end + +@kernel function IceNucleation_homogeneous_J_linear_kernel!( + ip, + output::AbstractArray{FT}, + Delta_a_w, +) where {FT} + + i = @index(Group, Linear) + + @inbounds begin + output[1] = CMI_hom.homogeneous_J_linear(ip.homogeneous, Delta_a_w[1]) end end @@ -1091,10 +1104,23 @@ function test_gpu(FT) x_sulph = ArrayType([FT(0.15)]) Delta_a_w = ArrayType([FT(0.2907389666103033)]) - kernel! = IceNucleation_homogeneous_J_kernel!(backend, work_groups) + kernel! = IceNucleation_homogeneous_J_cubic_kernel!(backend, work_groups) + kernel!(ip, output, Delta_a_w; ndrange) + + # test homogeneous_J_cubic is callable and returns a reasonable value + @test Array(output)[1] ≈ FT(2.66194650334444e12) + + dims = (1, 1) + (; output, ndrange) = setup_output(dims, FT) + + T = ArrayType([FT(220)]) + x_sulph = ArrayType([FT(0.15)]) + Delta_a_w = ArrayType([FT(0.2907389666103033)]) + + kernel! = IceNucleation_homogeneous_J_linear_kernel!(backend, work_groups) kernel!(ip, output, Delta_a_w; ndrange) - # test homogeneous_J is callable and returns a reasonable value + # test homogeneous_J_linear is callable and returns a reasonable value @test Array(output)[1] ≈ FT(2.66194650334444e12) end diff --git a/test/homogeneous_ice_nucleation_tests.jl b/test/homogeneous_ice_nucleation_tests.jl index 8bb37a07e7..2b7b6bdfc9 100644 --- a/test/homogeneous_ice_nucleation_tests.jl +++ b/test/homogeneous_ice_nucleation_tests.jl @@ -47,8 +47,41 @@ function test_homogeneous_J_cubic(FT) end end +function test_homogeneous_J_linear(FT) + + tps = TD.Parameters.ThermodynamicsParameters(FT) + H2SO4_prs = CMP.H2SO4SolutionParameters(FT) + ip = CMP.IceNucleationParameters(FT) + + TT.@testset "Homogeneous J (linear parameterization)" begin + + T_warm = FT(229.2) + T_cold = FT(228.8) + x_sulph = FT(0.1) + + # higher nucleation rate at colder temperatures + TT.@test CMH.homogeneous_J_linear( + ip.homogeneous, + CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold) - + CO.a_w_ice(tps, T_cold), + ) > CMH.homogeneous_J_linear( + ip.homogeneous, + CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm) - + CO.a_w_ice(tps, T_warm), + ) + + end +end + + println("Testing Float64") test_homogeneous_J_cubic(Float64) println("Testing Float32") test_homogeneous_J_cubic(Float32) + +println("Testing Float64") +test_homogeneous_J_linear(Float64) + +println("Testing Float32") +test_homogeneous_J_linear(Float32) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index cbe34b7e33..8851a12c5b 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -159,7 +159,8 @@ function benchmark_test(FT) (ip_frostenberg, INPC, T_air_cold), 150, ) - bench_press(CMI_hom.homogeneous_J, (ip.homogeneous, Delta_a_w), 230) + bench_press(CMI_hom.homogeneous_J_cubic, (ip.homogeneous, Delta_a_w), 230) + bench_press(CMI_hom.homogeneous_J_lienar, (ip.homogeneous, Delta_a_w), 230) # non-equilibrium bench_press(CMN.τ_relax, (liquid,), 10)