P3 ice nucleation parameterizations

CliMA · Feb 17, 2024 · 2a2d62f · 2a2d62f
1 parent 2c27e0f
commit 2a2d62f
Show file tree

Hide file tree

Showing 11 changed files with 462 additions and 13 deletions.
diff --git a/docs/bibliography.bib b/docs/bibliography.bib
@@ -46,6 +46,14 @@ @Article{Alpert2022
 doi = "10.1039/D1EA00077B",
 }
 
+@article{BarklieGokhale1959,
+  title = {The freezing of supercooled wter drops. Alberta hal, 1958, and related studies},
+  author = {Barklie, R. H. D. and Gokhale, N. R.},
+  journal = {McGill University Stormy Weather Group Sci. Rep. MW-30},
+  pages = {43-64},
+  year = {1959},
+}
+
 @Article{Baumgartner2022,
 author = {Baumgartner, M. and Rolf, C. and Groo{\ss}, J.-U. and Schneider, J. and Schorr, T. and M\"ohler, O. and Spichtinger, P. and Kr\"amer, M.},
 title = {New investigations on homogeneous ice nucleation: the effects of water activity and water saturation formulations},
@@ -69,6 +77,16 @@ @article{Beheng1994
   doi = {10.1016/0169-8095(94)90020-5}
 }
 
+@article{Bigg1953,
+  title = {The supercooling of water.},
+  author = {Bigg, E.K.},
+  journal = {Proc. Phys. Soc.},
+  volume = {66B},
+  pages = {688-694},
+  year = {1953},
+  doi = {10.1088/0370-1301/66/8/309}
+}
+
 @article{BrownFrancis1995,
   title = {Improved Measurements of the Ice Water Content in Cirrus Using a Total-Water Probe},
   author = {Philip R. A. Brown and Peter N. Francis},
@@ -599,6 +617,17 @@ @article{Spichtinger2023
   DOI = {10.5194/acp-23-2035-2023}
 }
 
+@article{Thompson2004,
+  author = {Gregory Thompson and Roy M. Rasmussen and Kevin Manning},
+  title = {Explicit Forecasts of Winter Precipitation Using an Improved Bulk Microphysics Scheme. Part I: Description and Sensitivity Analysis},
+  journal = {Monthly Weather Review},
+  year = {2004},
+  volume = {132},
+  number = {2},
+  doi = {10.1175/1520-0493(2004)132<0519:EFOWPU>2.0.CO;2},
+  pages = {519 - 542},
+}
+
 @article{TripoliCotton1980,
   title = {A Numerical Investigation of Several Factors Contributing to the Observed Variable Intensity of Deep Convection over South Florida},
   author = {Tripoli, G.J. and Cotton, W.R.},

diff --git a/docs/src/API.md b/docs/src/API.md
@@ -90,6 +90,8 @@ HetIceNucleation.dust_activated_number_fraction
 HetIceNucleation.MohlerDepositionRate
 HetIceNucleation.deposition_J
 HetIceNucleation.ABIFM_J
+HetIceNucleation.P3_deposition_N_i
+HetIceNucleation.P3_het_N_i
 HetIceNucleation.INP_concentration_frequency
 ```
 

diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md
@@ -63,7 +63,7 @@ Limited experimental values for the free parameters ``a`` and ``S_0`` can be fou
     freezing occurs in a different ice nucleation mode
     (either a second deposition or other condensation type mode).
 
-## Water Activity Based Deposition Nucleation
+### Water Activity Based Deposition Nucleation
 The water activity based deposition nucleation model is analagous to ABIFM
   for immersion freezing (see [ABIFM for Sulphuric Acid Containing Droplets](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/#ABIFM-for-Sulphuric-Acid-Containing-Droplets)
   section below). It calculates a nucleation rate coefficient, ``J``, which
@@ -94,7 +94,22 @@ include("plots/activity_based_deposition.jl")
 ![](water_activity_depo_nuc.svg)
 
 
-## ABIFM for Sulphuric Acid Containing Droplets
+### Cooper (1986) Ice Crystal Number
+The P3 scheme as described in [MorrisonMilbrandt2015](@cite) follows
+  the implementation of [Thompson2004](@cite) where a parameterization from
+  Cooper 1986 was used. Number of ice crystals formed is derived from
+  ambient temperature.
+```math
+\begin{equation}
+  N_i = 0.005 exp[0.304(T_0 - T)]
+\end{equation}
+```
+Where ``T_0 = 273.15K`` and ``T`` is the ambient temperature.
+Because the parameterization is prone to overpredict at low temperatures,
+  ``N_i = N_i(T = 233K)`` for all values of ``N_i`` at which ``T < 233K``.
+
+## Immersion Freezing
+### 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
   classical nucleation theory (CNT). More on CNT can be found in [Karthika2016](@cite).
@@ -169,7 +184,23 @@ It is also important to note that this plot is reflective of cirrus clouds
   and shows only a very small temperature range. The two curves are slightly
   off because of small differences in parameterizations for vapor pressures.
 
-## Homogeneous Freezing for Sulphuric Acid Containing Droplets
+### Bigg (1953) Volume and Time Dependent Heterogeneous Freezing
+Heterogeneous freezing in the P3 scheme as described in [MorrisonMilbrandt2015](@cite)
+  follows the parameterization from [Bigg1953](@cite) with parameters from 
+  [BarklieGokhale1959](@cite). The number of ice nucleated in a timestep via
+  heterogeneous freezing is determined by
+```math
+\begin{equation}
+  N_{ice} = N_{liq} \left[ 1 - exp(-B V_{l} \Delta t exp(aT_s)) \right]
+\end{equation}
+```
+where `a` and `B` are parameters taken from [BarklieGokhale1959](@cite) for
+  rainwater as 0.65 and 2e-4 respectively. `T_s` is the difference between
+  273.15K and ambient temperature. `V_l` is the volume of droplets to be
+  frozen and `\Delta t` is timestep in which freezing occurs.
+
+## Homogeneous Freezing
+### Homogeneous Freezing for Sulphuric Acid Containing Droplets
 Homogeneous freezing occurs when supercooled liquid droplets freeze on their own.
   Closly based off [Koop2000](@cite), this parameterization determines a homoegneous nucleation
   rate coefficient, ``J_{hom}``, using water activity. The change in water activity,
@@ -185,7 +216,7 @@ The nucleation rate coefficient is determined with the cubic function from [Koop
 ```
 This parameterization is valid only when ``0.26 < \Delta a_w < 0.36`` and ``185K < T < 235K``.
 
-### Homogeneous Ice Nucleation Rate Coefficient
+### Homogeneous Freezing Example Figures
 Here is a comparison of our parameterization of ``J_{hom}`` compared to Koop 2000 as
   plotted in figure 1 of [Spichtinger2023](@cite). Our parameterization differs in the calculation
   of ``\Delta a_w``. We define water activity to be a ratio of saturated vapor pressures whereas

diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md
@@ -315,3 +315,21 @@ The following plots show the parcel model running with homogeneous freezing and
 include("../../parcel/Jensen_et_al_2022.jl")
 ```
 ![](Jensen_et_al_2022.svg)
+
+## P3 Ice Nucleation Parameterizations
+The parcel also includes ice nucleation parameterizations used in
+  the P3 scheme as described in [MorrisonMilbrandt2015](@cite). 
+  Deposition nucleation is based on the ice crystal number parameterization
+  from Cooper (1986). The heterogeneous freezing parameterization, which
+  follows Bigg(1953) with parameters from Barklie aand Gokhale (1959), is
+  treated as immersion freezing in the parcel. Homogeneous freezing happens
+  instantaneously at 233.15K.
+Shown below are three separate parcel simulations for deposition nucleation,
+  immersion freezing, and homogeneous freezing. Note that initial temperature
+  varies for each run. The deposition nucleation run does not conserve
+  INP number, while the other two freezing modes do. Updraft velocity is
+  set to 0.5 m/s.
+```@example
+include("../../parcel/P3_ice_nuc.jl")
+```
+![](P3_ice_nuc.svg)
diff --git a/parcel/P3_ice_nuc.jl b/parcel/P3_ice_nuc.jl
@@ -0,0 +1,125 @@
+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)
+Nₗ = FT(2000)
+Nᵢ = FT(0)
+rₗ = FT(1.25e-6)
+p₀ = FT(20000)
+T₀_dep = FT(235)
+T₀_het = FT(235)
+T₀_hom = FT(233.2)
+qᵥ = FT(8.3e-4)
+qₗ = FT(Nₗ * 4 / 3 * π * rₗ^3 * ρₗ / 1.2)
+qᵢ = FT(0)
+x_sulph = FT(0)
+
+# Moisture dependent initial conditions
+q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+ts_dep = TD.PhaseNonEquil_pTq(tps, p₀, T₀_dep, q)
+ts_het = TD.PhaseNonEquil_pTq(tps, p₀, T₀_het, q)
+ts_hom = TD.PhaseNonEquil_pTq(tps, p₀, T₀_hom, q)
+ρₐ_dep = TD.air_density(tps, ts_dep)
+ρₐ_het = TD.air_density(tps, ts_het)
+ρₐ_hom = TD.air_density(tps, ts_hom)
+Rₐ = TD.gas_constant_air(tps, q)
+eₛ_dep = TD.saturation_vapor_pressure(tps, T₀_dep, TD.Liquid())
+eₛ_het = TD.saturation_vapor_pressure(tps, T₀_het, TD.Liquid())
+eₛ_hom = TD.saturation_vapor_pressure(tps, T₀_hom, TD.Liquid())
+e = qᵥ * p₀ * R_v / Rₐ
+
+Sₗ_dep = FT(e / eₛ_dep)
+Sₗ_het = FT(e / eₛ_het)
+Sₗ_hom = FT(e / eₛ_hom)
+IC_dep = [Sₗ_dep, p₀, T₀_dep, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC_het = [Sₗ_het, p₀, T₀_het, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC_hom = [Sₗ_hom, p₀, T₀_hom, 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(0.5)                             # updraft speed
+α_m = FT(0.5)                           # accomodation coefficient
+const_dt = FT(0.1)                      # model timestep
+t_max = FT(50)
+aerosol = []
+ice_nucleation_modes_list = [["P3_Deposition"], ["P3_het"], ["P3_hom"]]
+growth_modes = ["Deposition"]
+droplet_size_distribution_list = [["Monodisperse"]]
+
+# Plotting
+fig = MK.Figure(resolution = (900, 600))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]")
+ax7 = MK.Axis(fig[1, 3], ylabel = "Temperature [K]")
+ax3 = MK.Axis(fig[2, 1], ylabel = "Ice Saturation [-]")
+ax4 = MK.Axis(fig[2, 2], ylabel = "ICNC [cm^-3]")
+ax8 = MK.Axis(fig[2, 3], ylabel = "Temperature [K]")
+ax5 = MK.Axis(fig[3, 1], ylabel = "Ice Saturation [-]", xlabel = "Time [s]")
+ax6 = MK.Axis(fig[3, 2], ylabel = "ICNC [cm^-3]", xlabel = "Time [s]")
+ax9 = MK.Axis(fig[3, 3], ylabel = "Temperature [K]", xlabel = "Time [s]")
+
+for ice_nucleation_modes in ice_nucleation_modes_list
+    nuc_mode = ice_nucleation_modes[1]
+    droplet_size_distribution = droplet_size_distribution_list[1]
+    p = (;
+        wps,
+        aps,
+        tps,
+        ip,
+        const_dt,
+        r_nuc,
+        w,
+        α_m,
+        aerosol,
+        ice_nucleation_modes,
+        growth_modes,
+        droplet_size_distribution,
+    )
+    # solve ODE
+    #! format: off
+    if "P3_Deposition" in ice_nucleation_modes
+        sol = run_parcel(IC_dep, FT(0), t_max, p)
+        MK.lines!(ax1, sol.t, S_i.(sol[3, :], (sol[1, :])), label = nuc_mode) # saturation
+        MK.lines!(ax2, sol.t, sol[9, :] * 1e-6)                               # ICNC
+        MK.lines!(ax7, sol.t, sol[3, :])                                      # Temperature
+        MK.axislegend(ax1, framevisible = false, labelsize = 12, orientation = :horizontal, position = :rt)
+    elseif "P3_het" in ice_nucleation_modes
+        sol = run_parcel(IC_het, FT(0), t_max, p)
+        MK.lines!(ax3, sol.t, S_i.(sol[3, :], (sol[1, :])), label = nuc_mode, color = :red) # saturation
+        MK.lines!(ax4, sol.t, sol[9, :] * 1e-6, color = :red)                               # ICNC
+        MK.lines!(ax8, sol.t, sol[3, :], color = :red)                                                    # Temperature
+        MK.axislegend(ax3, framevisible = false, labelsize = 12, orientation = :horizontal, position = :lt)
+    elseif "P3_hom" in ice_nucleation_modes
+        sol = run_parcel(IC_hom, FT(0), t_max, p)
+        MK.lines!(ax5, sol.t, S_i.(sol[3, :], (sol[1, :])), label = nuc_mode, color = :orange) # saturation
+        MK.lines!(ax6, sol.t, sol[9, :] * 1e-6, color = :orange)                               # ICNC
+        MK.lines!(ax9, sol.t, sol[3, :], color = :orange)                                                       # Temperature
+        MK.axislegend(ax5, framevisible = false, labelsize = 12, orientation = :horizontal, position = :lt)
+    end
+    #! format: on
+end
+
+MK.save("P3_ice_nuc.svg", fig)
diff --git a/parcel/parcel.jl b/parcel/parcel.jl
@@ -107,6 +107,11 @@ function parcel_model(dY, Y, p, t)
             dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
         end
     end
+    if "P3_Deposition" in ice_nucleation_modes
+        Nᵢ_depo = CMI_het.P3_deposition_N_i(ip.p3, T)
+        dN_act_dt_depo = max(FT(0), (Nᵢ_depo - N_ice) / const_dt)
+        dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
+    end
 
     dN_act_dt_immersion = FT(0)
     dqi_dt_new_immers = FT(0)
@@ -126,6 +131,19 @@ function parcel_model(dY, Y, p, t)
         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
+    if "P3_het" in ice_nucleation_modes
+        if "Monodisperse" in droplet_size_distribution
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air)
+        elseif "Gamma" in droplet_size_distribution
+            λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air)
+            #A = N_liq* λ^2
+            r_l = 2 / λ
+        end
+        V_l = FT(4 / 3 * π * r_l^3)
+        Nᵢ_het = CMI_het.P3_het_N_i(ip.p3, T, N_liq, V_l, const_dt)
+        dN_act_dt_immersion = max(FT(0), (Nᵢ_het - N_ice) / const_dt)
+        dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
+    end
 
     dN_act_dt_hom = FT(0)
     dqi_dt_new_hom = FT(0)
@@ -159,6 +177,23 @@ function parcel_model(dY, Y, p, t)
             dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_aero^3 * ρ_ice / ρ_air
         end
     end
+    if "P3_hom" in ice_nucleation_modes
+        if "Monodisperse" in droplet_size_distribution
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air)
+        elseif "Gamma" in droplet_size_distribution
+            λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air)
+            #A = N_liq* λ^2
+            r_l = 2 / λ
+        elseif "Lognormal" in droplet_size_distribution
+            r_l = R_mode_liq * exp(σ)
+        end
+        if T < FT(233.15) && N_liq > FT(0)
+            dN_act_dt_hom = N_liq / const_dt
+        else
+            dN_act_dt_hom = FT(0)
+        end
+        dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
+    end
 
     dN_ice_dt = dN_act_dt_depo + dN_act_dt_immersion + dN_act_dt_hom
     if "Jensen_Example" in droplet_size_distribution
@@ -306,12 +341,25 @@ function run_parcel(IC, t_0, t_end, p)
         print("\n    Water activity based deposition nucleation ")
         print("(Note that this mode only runs for monodisperse size distributions.) ")
     end
+    if "P3_Deposition" in ice_nucleation_modes
+        print("\n    P3 deposition nucleation ")
+        timestepper = ODE.Euler()
+    end
     if "ImmersionFreezing" in ice_nucleation_modes
-        print("\n    Immersion freezing ")
+        print("\n    Immersion freezing (ABIFM) ")
+        print("(Note that this mode only runs for monodisperse and gamma size distributions.) ")
+    end
+    if "P3_het" in ice_nucleation_modes
+        print("\n    P3 heterogeneous freezing ")
         print("(Note that this mode only runs for monodisperse and gamma size distributions.) ")
+        timestepper = ODE.Euler()
     end
     if "HomogeneousFreezing" in ice_nucleation_modes
-        print("\n    Homogeneous freezing ")
+        print("\n    Water activity based homogeneous freezing ")
+    end
+    if "P3_hom" in ice_nucleation_modes
+        print("\n    P3 homogeneous freezing ")
+        timestepper = ODE.Euler()
     end
     print("\n")
 
@@ -338,8 +386,14 @@ function run_parcel(IC, t_0, t_end, p)
     if "Lognormal" in droplet_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.")
+    if "MohlerAF_Deposition" in ice_nucleation_modes ||
+        "MohlerRate_Deposition" in ice_nucleation_modes ||
+        "ActivityBasedDeposition" in ice_nucleation_modes ||
+        "P3_Deposition" in ice_nucleation_modes
+        print("\nWARNING: Multiple deposition nucleation parameterizations chosen to run in one simulation. Parcel can only properly handle one for now.")
+    end
+    if "P3_het" in ice_nucleation_modes
+        print("\nWARNING: Assuming P3_het parameterization is a parameterization for immeresion freezing.")
     end
     println(" ")