Merge pull request #366 from CliMA/al/linear_HOM

Linear HOM J option

docs/src/API.md

@@ -109,7 +109,8 @@ HetIceNucleation.INP_concentration_mean
 # Homogeneous ice nucleation
 ```@docs
 HomIceNucleation
-HomIceNucleation.homogeneous_J
+HomIceNucleation.homogeneous_J_cubic
+HomIceNucleation.homogeneous_J_linear
 ```
 
 # Common utility functions

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

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

docs/src/plots/linear_HOM_J.jl

@@ -0,0 +1,41 @@
+import CairoMakie as MK
+import CloudMicrophysics as CM
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.HomIceNucleation as CMH
+
+FT = Float32
+const ip = CMP.IceNucleationParameters(FT)
+
+# 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)
+
+# 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, 600), fontsize = 18)
+ax1 = MK.Axis(
+    fig[1, 1],
+    xlabel = "Δa_w [-]",
+    ylabel = "log10(J), J = [cm^-3 s^-1]",
+    title = "Linear fit to Koop2000",
+)
+
+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 = :rb)
+
+#! format: on
+
+MK.save("linear_HOM_J.svg", fig)

parcel/ParcelTendencies.jl

@@ -180,7 +180,7 @@ function homogeneous_freezing(params::ABHOM, PSD, state)
         Δa_w = ips.homogeneous.Δa_w_max
     end
 
-    J = CMI_hom.homogeneous_J(ips.homogeneous, Δa_w)
+    J = CMI_hom.homogeneous_J_cubic(ips.homogeneous, Δa_w)
 
     return min(max(FT(0), J * Nₗ * PSD.Vₗ), Nₗ / const_dt)
 end

src/IceNucleation.jl

@@ -214,10 +214,11 @@ module HomIceNucleation
 import ..Parameters as CMP
 import Thermodynamics as TD
 
-export homogeneous_J
+export homogeneous_J_cubic
+export homogeneous_J_linear
 
 """
-    homogeneous_J(ip, Δa_w)
+    homogeneous_J_cubic(ip, Δa_w)
 
  - `ip` - a struct with ice nucleation parameters
  - `Δa_w` - change in water activity [-].
@@ -226,7 +227,7 @@ 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.
 """
-function homogeneous_J(ip::CMP.Koop2000, Δa_w::FT) where {FT}
+function homogeneous_J_cubic(ip::CMP.Koop2000, Δa_w::FT) where {FT}
 
     @assert Δa_w >= ip.Δa_w_min
     @assert Δa_w <= ip.Δa_w_max
@@ -236,4 +237,21 @@ function homogeneous_J(ip::CMP.Koop2000, Δa_w::FT) where {FT}
     return FT(10)^(logJ) * 1e6
 end
 
+"""
+    homogeneous_J_linear(ip, Δa_w)
+
+ - `ip` - a struct with ice nucleation parameters
+ - `Δa_w` - change in water activity [-].
+
+Returns the homogeneous freezing nucleation rate coefficient,
+`J`, in m^-3 s^-1 for sulphuric acid solutions.
+Parameterization derived from a linear fit of the Koop 2000 parameterization, DOI: 10.1038/35020537.
+"""
+function homogeneous_J_linear(ip::CMP.Koop2000, Δa_w::FT) where {FT}
+
+    logJ::FT = ip.linear_c₂ * Δa_w + ip.linear_c₁
+
+    return FT(10)^(logJ) * 1e6
+end
+
 end # end module

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_coeff2 => :linear_c₂,
     )
     parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
     FT = CP.float_type(td)

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
 
@@ -1088,14 +1101,21 @@ function test_gpu(FT)
         (; 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_kernel!(backend, work_groups)
+        kernel! =
+            IceNucleation_homogeneous_J_cubic_kernel!(backend, work_groups)
         kernel!(ip, output, Delta_a_w; ndrange)
 
-        # test homogeneous_J is callable and returns a reasonable value
+        # test homogeneous_J_cubic is callable and returns a reasonable value
         @test Array(output)[1] ≈ FT(2.66194650334444e12)
+
+        kernel! =
+            IceNucleation_homogeneous_J_linear_kernel!(backend, work_groups)
+        kernel!(ip, output, Delta_a_w; ndrange)
+
+        # test homogeneous_J_linear is callable and returns a reasonable value
+        @test Array(output)[1] ≈ FT(7.156568123338207e11)
     end
 
     @testset "Modal nucleation kernels" begin

test/homogeneous_ice_nucleation_tests.jl

@@ -10,13 +10,13 @@ import CloudMicrophysics.HomIceNucleation as CMH
 
 @info "Homogeneous Ice Nucleation Tests"
 
-function test_homogeneous_J(FT)
+function test_homogeneous_J_cubic(FT)
 
     tps = TD.Parameters.ThermodynamicsParameters(FT)
     H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
     ip = CMP.IceNucleationParameters(FT)
 
-    TT.@testset "Homogeneous J" begin
+    TT.@testset "Homogeneous J (cubic parameterization)" begin
 
         T_warm = FT(229.2)
         T_cold = FT(228.8)
@@ -25,30 +25,63 @@ function test_homogeneous_J(FT)
         Δa_w_too_large = FT(0.35)
 
         # higher nucleation rate at colder temperatures
-        TT.@test CMH.homogeneous_J(
+        TT.@test CMH.homogeneous_J_cubic(
             ip.homogeneous,
             CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold) -
             CO.a_w_ice(tps, T_cold),
-        ) > CMH.homogeneous_J(
+        ) > CMH.homogeneous_J_cubic(
             ip.homogeneous,
             CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm) -
             CO.a_w_ice(tps, T_warm),
         )
 
         # If Δa_w out of range
-        TT.@test_throws AssertionError("Δa_w >= ip.Δa_w_min") CMH.homogeneous_J(
+        TT.@test_throws AssertionError("Δa_w >= ip.Δa_w_min") CMH.homogeneous_J_cubic(
             ip.homogeneous,
             Δa_w_too_small,
         )
-        TT.@test_throws AssertionError("Δa_w <= ip.Δa_w_max") CMH.homogeneous_J(
+        TT.@test_throws AssertionError("Δa_w <= ip.Δa_w_max") CMH.homogeneous_J_cubic(
             ip.homogeneous,
             Δa_w_too_large,
         )
     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(Float64)
+test_homogeneous_J_linear(Float64)
 
 println("Testing Float32")
-test_homogeneous_J(Float32)
+test_homogeneous_J_linear(Float32)

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_linear, (ip.homogeneous, Delta_a_w), 230)
 
     # non-equilibrium
     bench_press(CMN.τ_relax, (liquid,), 10)