Merge pull request #435 from CliMA/aj/het_freezing

Add P3 immersion freezing
CliMA · Aug 15, 2024 · 4549469 · 4549469
2 parents b24eedf + 2e18da2
commit 4549469
Show file tree

Hide file tree

Showing 6 changed files with 237 additions and 50 deletions.
diff --git a/docs/src/API.md b/docs/src/API.md
@@ -67,6 +67,7 @@ P3Scheme
 P3Scheme.thresholds
 P3Scheme.distribution_parameter_solver
 P3Scheme.ice_terminal_velocity
+P3Scheme.het_ice_nucleation
 ```
 
 # Aerosol model

diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md
@@ -215,6 +215,34 @@ include("plots/P3TerminalVelocityPlots.jl")
 ```
 ![](MorrisonandMilbrandtFig2.svg)
 
+## Heterogeneous Freezing
+
+Immersion freezing is parameterized based on water activity and follows the ABIFM
+  parameterization from [KnopfAlpert2013](@cite).
+See also the derivation notes about different
+  [ice nucleation parameterizations](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/).
+The immersion freezing nucleation rate is computed by numerically integrating
+  over the distribution of cloud droplets given by the 2-moment warm rain
+  microphysics scheme from [SeifertBeheng2006](@cite).
+The rate is limited by the available cloud droplet number concentration
+  and water content.
+```math
+\frac{dN}{dt} = \int_{0}^{D_{max}} \! J_{ABIFM} A_a(D) N'(D) \mathrm{d}D
+```
+```math
+\frac{dQ}{dt} = \int_{0}^{D_{max}} \! J_{ABIFM} A_a(D) N'(D) m(D) \mathrm{d}D
+```
+where
+- ``J_{ABIFM}`` - is the immersion freezing nucleation rate,
+- ``A_a(D)`` - is the assumed surface area of insoluble ice nucleating particles,
+- ``N'(D)`` - number distribution of cloud droplets,
+- ``m(D)`` - assumed mass of a cloud droplet as a function of its diameter.
+
+```@example
+include("plots/P3ImmersionFreezing.jl")
+```
+![](P3_het_ice_nucleation.svg)
+
 ## Acknowledgments
 
 Click on the P3 mascot duck to be taken to the repository

diff --git a/docs/src/plots/P3ImmersionFreezing.jl b/docs/src/plots/P3ImmersionFreezing.jl
@@ -0,0 +1,106 @@
+import CairoMakie as PL
+PL.activate!(type = "svg")
+
+import Thermodynamics as TD
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.P3Scheme as P3
+
+FT = Float64
+
+# thermodynamics parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+
+# helper functions
+function RH2qₜ(T, RH)
+    eᵥ_sat = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
+    eᵥ = RH * eᵥ_sat
+    qᵥ = 1 / (1 - tps.molmass_dryair / tps.molmass_water * (eᵥ - p) / eᵥ)
+    qₜ = qᵥ + qₗ + qᵢ
+    return qₜ
+end
+function p2ρ(T, RH)
+    return TD.air_density(tps, T, p, TD.PhasePartition(RH2qₜ(T, RH), qₗ, qᵢ))
+end
+
+# ambient conditions
+Nₗ = FT(500 * 1e6)
+qᵢ = FT(1 * 1e-3)
+qₗ = FT(1 * 1e-3)
+p = FT(800 * 1e2)
+
+# supported aerosol types
+dd = CMP.DesertDust(FT)
+il = CMP.Illite(FT)
+
+# model time step (for limiting)
+dt = FT(1)
+
+# plot data
+RH_range = range(0.8, stop = 1.2, length = 1000)
+T1 = FT(273.15 - 15)
+T2 = FT(273.15 - 35)
+
+#! format: off
+
+# limiters to not nucleate more mass and number than we have in liquid phase
+max_dLdt_T1 = [qₗ* p2ρ(T1, RH) / dt for RH in RH_range]
+max_dLdt_T2 = [qₗ* p2ρ(T2, RH) / dt for RH in RH_range]
+max_dNdt =    [Nₗ              / dt for RH in RH_range]
+
+dLdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, TD.PhasePartition(RH2qₜ(T1, RH), qₗ, qᵢ), Nₗ, RH, T1, p2ρ(T1, RH), dt).dLdt for RH in RH_range]
+dNdt_dd_T1 = [P3.het_ice_nucleation(dd, tps, TD.PhasePartition(RH2qₜ(T1, RH), qₗ, qᵢ), Nₗ, RH, T1, p2ρ(T1, RH), dt).dNdt for RH in RH_range]
+
+dLdt_il_T1 = [P3.het_ice_nucleation(il, tps, TD.PhasePartition(RH2qₜ(T1, RH), qₗ, qᵢ), Nₗ, RH, T1, p2ρ(T1, RH), dt).dLdt for RH in RH_range]
+dNdt_il_T1 = [P3.het_ice_nucleation(il, tps, TD.PhasePartition(RH2qₜ(T1, RH), qₗ, qᵢ), Nₗ, RH, T1, p2ρ(T1, RH), dt).dNdt for RH in RH_range]
+
+
+dLdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, TD.PhasePartition(RH2qₜ(T2, RH), qₗ, qᵢ), Nₗ, RH, T2, p2ρ(T2, RH), dt).dLdt for RH in RH_range]
+dNdt_dd_T2 = [P3.het_ice_nucleation(dd, tps, TD.PhasePartition(RH2qₜ(T2, RH), qₗ, qᵢ), Nₗ, RH, T2, p2ρ(T2, RH), dt).dNdt for RH in RH_range]
+
+dLdt_il_T2 = [P3.het_ice_nucleation(il, tps, TD.PhasePartition(RH2qₜ(T2, RH), qₗ, qᵢ), Nₗ, RH, T2, p2ρ(T2, RH), dt).dLdt for RH in RH_range]
+dNdt_il_T2 = [P3.het_ice_nucleation(il, tps, TD.PhasePartition(RH2qₜ(T2, RH), qₗ, qᵢ), Nₗ, RH, T2, p2ρ(T2, RH), dt).dNdt for RH in RH_range]
+
+# plotting
+fig = PL.Figure(size = (1500, 500), fontsize=22, linewidth=3)
+
+ax1 = PL.Axis(fig[1, 1]; yscale = log10)
+ax2 = PL.Axis(fig[1, 2]; yscale = log10)
+
+ax1.xlabel = "RH [%]"
+ax1.ylabel = "ice mass nucleation rate [g/m3/s]"
+ax2.xlabel = "RH [%]"
+ax2.ylabel = "ice number nucleation rate [1/cm3/s]"
+
+l_max_dLdt_T1 = PL.lines!(ax1, RH_range * 1e2,  max_dLdt_T1 * 1e3,  color = :thistle)
+l_max_dLdt_T2 = PL.lines!(ax1, RH_range * 1e2,  max_dLdt_T2 * 1e3,  color = :thistle)
+l_max_dNdt    = PL.lines!(ax2, RH_range * 1e2,  max_dNdt    * 1e-6, color = :thistle)
+
+
+l_dLdt_dd_T1 = PL.lines!(ax1, RH_range * 1e2,  dLdt_dd_T1 * 1e3,  color = :skyblue)
+l_dNdt_dd_T1 = PL.lines!(ax2, RH_range * 1e2,  dNdt_dd_T1 * 1e-6, color = :skyblue)
+
+l_dLdt_dd_T2 = PL.lines!(ax1, RH_range * 1e2,  dLdt_dd_T2 * 1e3,  color = :blue3)
+l_dNdt_dd_T2 = PL.lines!(ax2, RH_range * 1e2,  dNdt_dd_T2 * 1e-6, color = :blue3)
+
+
+l_dLdt_il_T1 = PL.lines!(ax1, RH_range * 1e2,  dLdt_il_T1 * 1e3,  color = :orchid)
+l_dNdt_il_T1 = PL.lines!(ax2, RH_range * 1e2,  dNdt_il_T1 * 1e-6, color = :orchid)
+
+l_dLdt_il_T2 = PL.lines!(ax1, RH_range * 1e2,  dLdt_il_T2 * 1e3,  color = :purple)
+l_dNdt_il_T2 = PL.lines!(ax2, RH_range * 1e2,  dNdt_il_T2 * 1e-6, color = :purple)
+
+PL.Legend(
+    fig[1, 3],
+    [l_max_dNdt,
+     l_dNdt_dd_T1, l_dNdt_dd_T2,
+     l_dNdt_il_T1, l_dNdt_il_T2],
+    [
+       "limit",
+       "T=-15C, desert dust",
+       "T=-35C, desert dust",
+       "T=-15C, illite",
+       "T=-35C, illite",
+    ],
+    framevisible = false,
+)
+PL.save("P3_het_ice_nucleation.svg", fig)
diff --git a/src/P3_processes.jl b/src/P3_processes.jl
@@ -1 +1,47 @@
-# TODO
+"""
+    het_ice_nucleation(pdf_c, p3, tps, q, N, T, ρₐ, p, aerosol)
+
+ - aerosol - aerosol parameters (supported types: desert dust, illite, kaolinite)
+ - tps - thermodynamics parameters
+ - qₚ - phase partition
+ - N_liq - cloud water number concentration
+ - RH - relative humidity
+ - T - temperature
+ - ρₐ - air density
+ - dt - model time step
+
+Returns a named tuple with ice number concentration and ice content
+hetergoeneous freezing rates from cloud droplets.
+"""
+function het_ice_nucleation(
+    aerosol::Union{CMP.DesertDust, CMP.Illite, CMP.Kaolinite},
+    tps::TDP.ThermodynamicsParameters{FT},
+    qₚ::TD.PhasePartition{FT},
+    N_liq::FT,
+    RH::FT,
+    T::FT,
+    ρₐ::FT,
+    dt::FT,
+) where {FT}
+    #TODO - Also consider rain freezing
+
+    # Immersion freezing nucleation rate coefficient
+    J = CM_HetIce.ABIFM_J(aerosol, RH - CO.a_w_ice(tps, T))
+
+    # Assumed erosol surface area
+    # TODO - Make it a parameter of ABIFM scheme
+    # We could consider making it a function of the droplet size distribution
+    A_aer = FT(1e-10)
+
+    dNdt = J * A_aer * N_liq
+    dLdt = J * A_aer * qₚ.liq * ρₐ
+
+    # nucleation rates are always positive definite...
+    dNdt = max(0, dNdt)
+    dLdt = max(0, dLdt)
+    # ... and dont exceed the available number and mass of water droplets
+    dNdt = min(dNdt, N_liq / dt)
+    dLdt = min(dLdt, qₚ.liq * ρₐ / dt)
+
+    return (; dNdt, dLdt)
+end
diff --git a/test/p3_tests.jl b/test/p3_tests.jl
@@ -459,54 +459,6 @@ end
 #        end
 #    end
 #
-#    TT.@testset "Heterogeneous Freezing Smoke Test" begin
-#        T = FT(250)
-#        N = FT(1e8)
-#        ρ_a = FT(1.2)
-#        qᵥ = FT(8.1e-4)
-#        aero_type = CMP.Illite(FT)
-#
-#        qs = range(0.001, stop = 0.005, length = 5)
-#
-#        expected_freeze_q =
-#            [2.036e-61, 6.463e-61, 1.270e-60, 2.052e-60, 2.976e-60]
-#        expected_freeze_N =
-#            [1.414e-51, 2.244e-51, 2.941e-51, 3.562e-51, 4.134e-51]
-#
-#        for i in axes(qs, 1)
-#            rate_mass = P3.p3_rain_het_freezing(
-#                true,
-#                pdf_r,
-#                p3,
-#                tps,
-#                qs[i],
-#                N,
-#                T,
-#                ρ_a,
-#                qᵥ,
-#                aero_type,
-#            )
-#            rate_num = P3.p3_rain_het_freezing(
-#                false,
-#                pdf_r,
-#                p3,
-#                tps,
-#                qs[i],
-#                N,
-#                T,
-#                ρ_a,
-#                qᵥ,
-#                aero_type,
-#            )
-#
-#            TT.@test rate_mass >= 0
-#            TT.@test rate_mass ≈ expected_freeze_q[i] rtol = 1e-3
-#
-#            TT.@test rate_num >= 0
-#            TT.@test rate_num ≈ expected_freeze_N[i] rtol = 1e-3
-#        end
-#    end
-#
 #end
 
 function test_integrals(FT)
@@ -574,6 +526,43 @@ function test_integrals(FT)
     end
 end
 
+function test_p3_het_freezing(FT)
+
+    TT.@testset "Heterogeneous Freezing Smoke Test" begin
+        tps = TD.Parameters.ThermodynamicsParameters(FT)
+        p3 = CMP.ParametersP3(FT)
+        SB2006 = CMP.SB2006(FT, false) # no limiters
+        pdf_c = SB2006.pdf_c
+        aerosol = CMP.Illite(FT)
+
+        N_liq = FT(1e8)
+        T = FT(244)
+        p = FT(500 * 1e2)
+
+        dt = FT(1)
+
+        expected_freeze_L =
+            [1.544e-22, 1.068e-6, 0.0001428, 0.0001428, 0.0001428, 0.0001428]
+        expected_freeze_N = [1.082e-10, 747647.5, N_liq, N_liq, N_liq, N_liq]
+        qᵥ_range = range(0.5e-3, stop = 1.5e-3, length = 6)
+
+        for it in range(1, 6)
+            qₚ = TD.PhasePartition(FT(qᵥ_range[it]), FT(2e-4), FT(0))
+            eᵥ_sat = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
+            ϵ = tps.molmass_water / tps.molmass_dryair
+            eᵥ = p * qᵥ_range[it] / (ϵ + qᵥ_range[it] * (1 - ϵ))
+            RH = eᵥ / eᵥ_sat
+            ρₐ = TD.air_density(tps, T, p, qₚ)
+            rate = P3.het_ice_nucleation(aerosol, tps, qₚ, N_liq, RH, T, ρₐ, dt)
+
+            TT.@test rate.dNdt >= 0
+            TT.@test rate.dLdt >= 0
+
+            TT.@test rate.dNdt ≈ expected_freeze_N[it] rtol = 1e-2
+            TT.@test rate.dLdt ≈ expected_freeze_L[it] rtol = 1e-2
+        end
+    end
+end
 
 println("Testing Float32")
 test_p3_thresholds(Float32)
@@ -588,5 +577,6 @@ test_p3_thresholds(Float64)
 test_p3_shape_solver(Float64)
 test_particle_terminal_velocities(Float64)
 test_bulk_terminal_velocities(Float64)
-#test_tendencies(Float64)
 test_integrals(Float64)
+test_p3_het_freezing(Float64)
+#test_tendencies(Float64)
diff --git a/test/performance_tests.jl b/test/performance_tests.jl
@@ -117,6 +117,7 @@ function benchmark_test(FT)
     N = FT(1e8)
 
     T_air_2 = FT(250)
+    RH_2 = FT(1.5)
     T_air_cold = FT(230)
     S_ice = FT(1.2)
     dSi_dt = FT(0.05)
@@ -171,6 +172,21 @@ function benchmark_test(FT)
         )
         bench_press(P3.D_m, (p3, q_ice, N, ρ_r, F_r), 1e5)
     end
+    # P3 ice nucleation
+    bench_press(
+        P3.het_ice_nucleation,
+        (
+            kaolinite,
+            tps,
+            TD.PhasePartition(q_tot, q_liq, q_ice),
+            N_liq,
+            RH_2,
+            T_air_2,
+            ρ_air,
+            Δt,
+        ),
+        150,
+    )
 
     # aerosol activation
     bench_press(