Merge pull request #304 from CliMA/al/fix_mohler

replacing Mohler AF with nucleation rate
CliMA · Feb 10, 2024 · cb7302f · cb7302f
2 parents d2bc21f + 6d7a331
commit cb7302f
Show file tree

Hide file tree

Showing 9 changed files with 369 additions and 142 deletions.
diff --git a/docs/src/API.md b/docs/src/API.md
@@ -87,6 +87,7 @@ AerosolActivation.total_M_activated
 ```@docs
 HetIceNucleation
 HetIceNucleation.dust_activated_number_fraction
+HetIceNucleation.MohlerDepositionRate
 HetIceNucleation.deposition_J
 HetIceNucleation.ABIFM_J
 ```

diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md
@@ -22,28 +22,41 @@ The parameterization for deposition on dust particles is an implementation of
     freezing modes.
 
 ## Activated fraction for deposition freezing on dust
-The parameterization models the activated fraction
-  as an empirical function of ice saturation ratio,
-  see eq. (3) in [Mohler2006](@cite).
+There are 2 parameterizations from [Mohler2006](@cite) available: one
+  which calculates the activated fraction and one which calculates nucleation
+  rate. 
+The activated fraction parameterization follows eq. (3) in the paper.
 ```math
 \begin{equation}
 f_i(S_i) = exp[a(S_i - S_0)] - 1
 \end{equation}
 ```
-where:
+where 
   - ``f_i`` is the activated fraction
-      (the ratio of aerosol particles acting as ice nuclei to the total number of aerosol particles),
+     (the ratio of aerosol particles acting as ice nuclei to the total number of aerosol particles),
+  - ``a`` is a scaling parameter dependent on aerosol properties and temperature,
+  - ``S_i`` is the ice saturation ratio,
+  - ``S_0`` is an empirically derived threshold ice saturation ratio.
+The other parameterization models the nucleation rate of ice
+  as an empirical function of ice saturation ratio,
+  see eq. (5) in [Mohler2006](@cite).
+```math
+\begin{equation}
+\frac{dn_{ice}}{dt} = N_{aer} a \frac{dS_i}{dt}
+\end{equation}
+```
+where:
+  - ``N_{aer}`` is the number of available aerosol/ice nuclei,
+  - ``a`` is a scaling parameter dependent on aerosol properties and temperature,
   - ``S_i`` is the ice saturation ratio
-      (the ratio of water vapor partial pressure and the water vapor partial pressure at saturation over ice),
-  - ``S_0`` is the threshold ice saturation ratio,
-  - ``a`` is a scaling parameter dependent on aerosol properties and temperature.
+    (the ratio of water vapor partial pressure and the water vapor partial pressure at saturation over ice).
 
-Limited experimental values for the free parameters ``S_0`` and ``a`` can be found in [Mohler2006](@cite).
-Both parameters are dependent on aerosol properties and temperature.
+Limited experimental values for the free parameters ``a`` and ``S_0`` can be found in [Mohler2006](@cite). These
+  free parameters are strongly dependent on aerosol properties and temperature.
 
 !!! note
 
-    For a ``f_i`` values above 0.08 or ``S_i`` between 1.35 and 1.5,
+    For ``f_i`` values above 0.08 or ``S_i`` between 1.35 and 1.5,
     freezing occurs in a different ice nucleation mode
     (either a second deposition or other condensation type mode).
 
@@ -97,10 +110,10 @@ Once ``J_{ABIFM}`` is calculated, it can be used to determine the ice production
 per second via immersion freezing.
 ```math
 \begin{equation}
-  P_{ice} = J_{ABIFM}A(N_{tot} - N_{ice})
+  P_{ice} = J_{ABIFM}A(N_{aer} - N_{ice})
 \end{equation}
 ```
-where ``A`` is surface area of an individual ice nuclei, ``N_{tot}`` is total number
+where ``A`` is surface area of an individual ice nuclei, ``N_{aer}`` is total number
   of ice nuclei, and ``N_{ice}`` is number of ice crystals already in the system.
 
 ### ABIFM Example Figures

diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md
@@ -121,7 +121,7 @@ where ``\xi = \frac{e_{sl}}{e_{si}}`` is the ratio of saturation vapor pressure
 The crux of the problem is modeling the ``\frac{dq_l}{dt}`` and ``\frac{dq_i}{dt}``
   for different homogeneous and heterogeneous ice nucleation paths.
 
-## Condensation
+## Condensation growth
 
 The diffusional growth of individual cloud droplet is described by
 ([see discussion](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics1M/#Snow-autoconversion)),
@@ -165,7 +165,7 @@ For a gamma distribution of droplets ``n(r) = A \; r \; exp(-\lambda r)``,
 ``\bar{r} = \frac{2}{\lambda}``
 where ``\lambda = \left(\frac{32 \pi N_{tot} \rho_l}{q_l \rho_a}\right)^{1/3}``.
 
-## Deposition nucleation on dust particles
+## Deposition growth
 
 Similarly, for a case of a spherical ice particle growing through water vapor deposition
 ```math
@@ -196,15 +196,12 @@ It follows that
 where:
 - ``N_{act}`` is the number of activated ice particles.
 
-``N_{act}`` can be computed for example from
-  [activated fraction](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i``
-```math
-\begin{equation}
-  N_{act} = N_{aer} f_i
-\end{equation}
-```
-where:
-- ``N_{aer}`` is the number of available dust aerosol particles.
+## Deposition Nucleation on dust particles
+There are multiple ways of running deposition nucleation in the parcel.
+  `"MohlerAF_Deposition"` will trigger an activated fraction approach from [Mohler2006](@cite).
+  `"MohlerRate_Deposition"` will trigger a nucleation rate approach from [Mohler2006](@cite).
+The deposition nucleation methods are parameterized as described in
+  [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/).
 
 ## Immersion Freezing
 Following the water activity based immersion freezing model (ABIFM), the ABIFM derived
@@ -241,13 +238,20 @@ Here we show various example simulation results from the adiabatic parcel
   and homogeneous freezing with deposition growth.
 
 We start with deposition freezing on dust.
-The model is run three times for 30 minutes simulation time,
-  (shown by three different colors on the plot).
+The model is run three times using the `"MohlerAF_Deposition"` approach 
+  for 30 minutes simulation time, (shown by three different colors on the plot).
 Between each run the water vapor specific humidity is changed,
   while keeping all other state variables the same as at the last time step
   of the previous run.
 The prescribed vertical velocity is equal to 3.5 cm/s.
 Supersaturation is plotted for both liquid (solid lines) and ice (dashed lines).
+The pale blue line uses the `"MohlerRate_Deposition"`approach. 
+  We only run it for the first GCM timestep because the rate approach requires
+  the change in ice saturation over time. With the discontinuous jump in saturation,
+  the parameterization is unable to determine a proper nucleation rate. When we force
+  the initial ice crystal number concentration for this simulation to match
+  that in the `"MohlerAF_Deposition"` approach, we obtain the same results as
+  in the `"MohlerAF_Deposition"` approach for the first GCM timestep.
 
 ```@example
 include("../../parcel/Tully_et_al_2023.jl")

diff --git a/parcel/Tully_et_al_2023.jl b/parcel/Tully_et_al_2023.jl
@@ -74,76 +74,12 @@ function Tully_et_al_2023(FT)
     α_m = FT(0.5)                              # accomodation coefficient
     const_dt = 0.1                             # model timestep
     aerosol = CMP.DesertDust(FT)               # aerosol type
-    ice_nucleation_modes = ["DustDeposition"]  # switch on deposition on dust
+    ice_nucleation_modes_list =                # ice nucleation modes
+        [["MohlerAF_Deposition"], ["MohlerRate_Deposition"]]
     growth_modes = ["Deposition"]              # switch on deposition growth
     droplet_size_distribution = ["Monodisperse"]
-    p = (;
-        wps,
-        aps,
-        tps,
-        ip,
-        const_dt,
-        r_nuc,
-        w,
-        α_m,
-        aerosol,
-        ice_nucleation_modes,
-        growth_modes,
-        droplet_size_distribution,
-    )
-
-    # Simulation 1
-    IC1 = get_initial_condition(
-        tps,
-        p_0,
-        T_0,
-        q_vap_0,
-        q_liq_0,
-        q_ice_0,
-        N_aerosol,
-        N_droplets,
-        N_0,
-        x_sulph,
-    )
-    sol1 = run_parcel(IC1, 0, t_max, p)
-
-    # Simulation 2
-    # (alternatively set T and take q_vap from the previous simulation)
-    #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end])
-    IC2 = get_initial_condition(
-        tps,
-        sol1[2, end],
-        #sol1[3, end],
-        T2,
-        q_vap2,
-        q_liq_0,
-        sol1[6, end],
-        sol1[7, end],
-        sol1[8, end],
-        sol1[9, end],
-        x_sulph,
-    )
-    sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, p)
-
-    # Simulation 3
-    # (alternatively set T and take q_vap from the previous simulation)
-    #IC3 = get_initial_condition(sol2[2, end], sol2[3, end], T3, sol2[5, end], 0.0, sol2[6, end], sol2[7, end])
-    IC3 = get_initial_condition(
-        tps,
-        sol2[2, end],
-        #sol2[3, end],
-        T3,
-        q_vap3,
-        q_liq_0,
-        sol2[6, end],
-        sol2[7, end],
-        sol2[8, end],
-        sol2[9, end],
-        x_sulph,
-    )
-    sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, p)
-
-    # Plot results
+
+    # Plots
     fig = MK.Figure(resolution = (800, 600))
     ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [-]")
     ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
@@ -160,23 +96,110 @@ function Tully_et_al_2023(FT)
         TD.saturation_vapor_pressure(tps, T, TD.Ice())
     S_i(T, S_liq) = ξ(T) * S_liq - 1
 
-    sol = [sol1, sol2, sol3]
-    clr = ["blue", "orange", "green"]
-    for it in [1, 2, 3]
-        MK.lines!(ax1, sol[it].t * w, sol[it][1, :] .- 1, color = clr[it])
-        MK.lines!(ax2, sol[it].t * w, sol[it][3, :], color = clr[it])
-        MK.lines!(ax3, sol[it].t * w, sol[it][9, :] * 1e-3, color = clr[it])
-        MK.lines!(ax4, sol[it].t * w, sol[it][7, :] * 1e-3, color = clr[it])
-        MK.lines!(ax5, sol[it].t * w, sol[it][4, :] * 1e3, color = clr[it])
-        MK.lines!(ax6, sol[it].t * w, sol[it][6, :] * 1e3, color = clr[it])
-
-        MK.lines!(
-            ax1,
-            sol[it].t * w,
-            S_i.(sol[it][3, :], sol[it][1, :]),
-            linestyle = :dash,
-            color = clr[it],
+    for ice_nucleation_modes in ice_nucleation_modes_list
+        p = (;
+            wps,
+            aps,
+            tps,
+            ip,
+            const_dt,
+            r_nuc,
+            w,
+            α_m,
+            aerosol,
+            ice_nucleation_modes,
+            growth_modes,
+            droplet_size_distribution,
+        )
+
+        # Simulation 1
+        IC1 = get_initial_condition(
+            tps,
+            p_0,
+            T_0,
+            q_vap_0,
+            q_liq_0,
+            q_ice_0,
+            N_aerosol,
+            N_droplets,
+            N_0,
+            x_sulph,
         )
+        sol1 = run_parcel(IC1, 0, t_max, p)
+        if ice_nucleation_modes == ["MohlerAF_Deposition"]
+            # Simulation 2
+            # (alternatively set T and take q_vap from the previous simulation)
+            #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end])
+            IC2 = get_initial_condition(
+                tps,
+                sol1[2, end],
+                #sol1[3, end],
+                T2,
+                q_vap2,
+                q_liq_0,
+                sol1[6, end],
+                sol1[7, end],
+                sol1[8, end],
+                sol1[9, end],
+                x_sulph,
+            )
+            sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, p)
+
+            # Simulation 3
+            # (alternatively set T and take q_vap from the previous simulation)
+            #IC3 = get_initial_condition(sol2[2, end], sol2[3, end], T3, sol2[5, end], 0.0, sol2[6, end], sol2[7, end])
+            IC3 = get_initial_condition(
+                tps,
+                sol2[2, end],
+                #sol2[3, end],
+                T3,
+                q_vap3,
+                q_liq_0,
+                sol2[6, end],
+                sol2[7, end],
+                sol2[8, end],
+                sol2[9, end],
+                x_sulph,
+            )
+            sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, p)
+
+            # Plot results
+            sol = [sol1, sol2, sol3]
+            clr = ["blue", "orange", "green"]
+            #! format: off
+            for it in [1, 2, 3]
+                MK.lines!(ax1, sol[it].t * w, sol[it][1, :] .- 1, color = clr[it])
+                MK.lines!(ax2, sol[it].t * w, sol[it][3, :], color = clr[it])
+                MK.lines!(ax3, sol[it].t * w, sol[it][9, :] * 1e-3, color = clr[it])
+                MK.lines!(ax4, sol[it].t * w, sol[it][7, :] * 1e-3, color = clr[it])
+                MK.lines!(ax5, sol[it].t * w, sol[it][4, :] * 1e3, color = clr[it])
+                MK.lines!(ax6, sol[it].t * w, sol[it][6, :] * 1e3, color = clr[it])
+                MK.lines!(
+                    ax1,
+                    sol[it].t * w,
+                    S_i.(sol[it][3, :], sol[it][1, :]),
+                    linestyle = :dash,
+                    color = clr[it],
+                )
+            end
+            #! format: on
+
+        elseif ice_nucleation_modes == ["MohlerRate_Deposition"]
+            MK.lines!(ax1, sol1.t * w, sol1[1, :] .- 1, color = :lightblue)
+            MK.lines!(ax2, sol1.t * w, sol1[3, :], color = :lightblue)
+            MK.lines!(ax3, sol1.t * w, sol1[9, :] * 1e-3, color = :lightblue)
+            MK.lines!(ax4, sol1.t * w, sol1[7, :] * 1e-3, color = :lightblue)
+            MK.lines!(ax5, sol1.t * w, sol1[4, :] * 1e3, color = :lightblue)
+            MK.lines!(ax6, sol1.t * w, sol1[6, :] * 1e3, color = :lightblue)
+
+            MK.lines!(
+                ax1,
+                sol1.t * w,
+                S_i.(sol1[3, :], sol1[1, :]),
+                linestyle = :dash,
+                color = :lightblue,
+            )
+        end
     end
     MK.save("cirrus_box.svg", fig)
 end