fix deposition warnings in parcel

box/Alpert_Knopf_2016_forward.jl

@@ -20,7 +20,7 @@ end
 include(joinpath(pkgdir(CM), "box", "box.jl"))
 
 # initial conditions following KA16 figure 4 (Cr1 case)
-A_droplet = FT(1e-5) * 1e-4 # INP surface area, m^2
+A_aero = FT(1e-5) * 1e-4 # INP surface area, m^2
 σg = 10
 N₀ = 1000
 N_ice = 0
@@ -31,10 +31,10 @@ tps = TD.Parameters.ThermodynamicsParameters(FT) # thermodynamics free parameter
 t_0 = 0
 t_end = 3310
 dt = 10
-A_sum = N₀ * A_droplet # Total surface area m^2 when all droplets are equal
+A_sum = N₀ * A_aero # Total surface area m^2 when all droplets are equal
 
 # A vector with surface areas sampled from a lognormal distribution and sorted
-A_distr = RD.rand(DS.LogNormal(log(A_droplet), log(σg)), N₀)
+A_distr = RD.rand(DS.LogNormal(log(A_aero), log(σg)), N₀)
 Aj_unsorted = zeros(N₀)
 for it in range(start = 1, stop = N₀, step = 1)
     Aj_unsorted[it] = A_distr[it]
@@ -45,7 +45,7 @@ Aj_sorted = sort(Aj_unsorted, rev = true)
 # initial condition for the ODE problem
 IC = [T_initial, A_sum, FT(N₀), FT(N_ice)]
 # additional model parameters
-p_def = (; tps, A_droplet, aerosol, cooling_rate, N₀)
+p_def = (; tps, A_aero, aerosol, cooling_rate, N₀)
 
 # run the box model without and with distribution sampling
 sol_cst = run_box(IC, t_0, t_end, (p_def..., const_dt = dt, flag = false))
@@ -80,7 +80,7 @@ MK.lines!(ax2, sol_cst[1, :],          ff_cst,               color = :orange, li
 MK.lines!(ax2, sol_var[1, :],          ff_var,               color = :green3, linewidth = 3, label = "variable A")
 
 idx = (sol_cst[3, :] .> 0)
-MK.lines!(ax3, sol_cst[1, idx], sol_cst[3, idx] .* A_droplet .* 1e4, color = :orange, linewidth = 3, label = "const A")
+MK.lines!(ax3, sol_cst[1, idx], sol_cst[3, idx] .* A_aero .* 1e4, color = :orange, linewidth = 3, label = "const A")
 idx = (sol_var[2, :] .> 0)
 MK.lines!(ax3, sol_var[1, idx], sol_var[2, idx] .* 1e4,              color = :green3, linewidth = 3, label = "variable A")
 MK.ylims!(ax3, 0.75e-6, 2e-1)

box/box.jl

@@ -11,7 +11,7 @@ import CloudMicrophysics.HetIceNucleation as CMI_het
 function box_model(dY, Y, p, t)
 
     # Get simulation parameters
-    (; tps, A_droplet, aerosol, cooling_rate) = p
+    (; tps, A_aero, aerosol, cooling_rate) = p
     # Numerical precision used in the simulation
     FT = eltype(Y)
 
@@ -28,7 +28,7 @@ function box_model(dY, Y, p, t)
     dN_ice_dt = FT(0)
     dN_liq_dt = FT(0)
     if N_liq > 0
-        dN_ice_dt = J_immer * N_liq * A_droplet
+        dN_ice_dt = J_immer * N_liq * A_aero
         dN_liq_dt = -dN_ice_dt
     end
 
@@ -45,7 +45,7 @@ end
 function box_model_with_probability(dY, Y, p, t)
 
     # Get simulation parameters
-    (; tps, const_dt, A_droplet, aerosol, cooling_rate, Aj_sorted, N₀) = p
+    (; tps, const_dt, A_aero, aerosol, cooling_rate, Aj_sorted, N₀) = p
     # Numerical precision used in the simulation
     FT = eltype(Y)
 
@@ -116,7 +116,7 @@ The named tuple p should contain:
  - const_dt - simulation timestep,
  - aerosol - insoluble aerosol parameters
  - cooling_rate - cooling rate
- - A_droplet - assumed surface area for freezing (when assuming constant A)
+ - A_aero - assumed surface area for freezing (when assuming constant A)
  - Aj_sorted - a vector with available surface area (when variable)
  - N₀ - initial number of liquid droplets
 """

docs/src/IceNucleationParcel0D.md

@@ -296,7 +296,8 @@ include("../../parcel/Liquid_only.jl")
 The plots below are the results of the adiabatic parcel model
   with immersion freezing, condensation growth, and deposition growth for
   both a monodisperse and gamma size distribution. Note that this has not
-  yet been validated against literature.
+  yet been validated against literature. Results are very sensitive to 
+  initial conditions.
 
 ```@example
 include("../../parcel/Immersion_Freezing.jl")

parcel/Immersion_Freezing.jl

@@ -20,15 +20,15 @@ R_d = TD.Parameters.R_d(tps)
 
 # Initial conditions
 Nₐ = FT(0)
-Nₗ = FT(500 * 1e3)
+Nₗ = FT(2000)
 Nᵢ = FT(0)
 r₀ = FT(1e-6)       # radius of droplets
 p₀ = FT(800 * 1e2)
 T₀ = FT(251)
 qᵥ = FT(8.1e-4)
 qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ
 qᵢ = FT(0)
-x_sulph = FT(0.01)
+x_sulph = FT(0)
 
 # Moisture dependent initial conditions
 q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -43,7 +43,7 @@ IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
 
 # Simulation parameters passed into ODE solver
 r_nuc = FT(0.5 * 1.e-4 * 1e-6)             # assumed size of nucleated particles
-w = FT(0.7)                                # updraft speed
+w = FT(0.4)                                # updraft speed
 α_m = FT(0.5)                              # accomodation coefficient
 const_dt = FT(1)                           # model timestep
 t_max = FT(1200)
@@ -59,13 +59,9 @@ ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
 ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]")
 ax4 = MK.Axis(fig[2, 2], ylabel = "q_liq [g/kg]")
 ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_liq")
-ax6 = MK.Axis(
-    fig[3, 2],
-    xlabel = "Time [min]",
-    ylabel = "N_ice / N_tot",
-    yscale = log10,
-)
-MK.ylims!(ax6, 1e-6, 1)
+ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_ice / N_tot")
+# ax7 = MK.Axis(fig[4, 1], ylabel = "r_liq [m]", yscale = log10)
+# MK.ylims!(ax7, 1e-6, 0.01)
 
 for droplet_size_distribution in droplet_size_distribution_list
     p = (;
@@ -101,15 +97,16 @@ for droplet_size_distribution in droplet_size_distribution_list
     sol_Nₗ = sol[8, :]
     sol_Nᵢ = sol[9, :]
     sol_qₗ = sol[5, :]
-    if DSD == "Monodisperse"
-        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ)
-    end
-    if DSD == "Gamma"
-        λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
-        rₗ = 2 ./ λ
-    end
+    # if DSD == "Monodisperse"
+    #     rₗ = max.(cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ), FT(0))
+    # end
+    # if DSD == "Gamma"
+    #     λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
+    #     rₗ = max.(2 ./ λ, FT(0))
+    # end
     MK.lines!(ax5, sol.t / 60, sol_Nₗ)
     MK.lines!(ax6, sol.t / 60, sol_Nᵢ ./ (sol_Nₗ .+ sol_Nᵢ))
+    # MK.lines!(ax7, sol.t / 60, rₗ)
 end
 
 MK.axislegend(

parcel/parcel.jl

@@ -127,8 +127,8 @@ function parcel_model(dY, Y, p, t)
             #A = N_liq* λ^2
             r_l = 2 / λ
         end
-        A_aero = 4 * FT(π) * r_nuc^2
-        dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_nuc)
+        A_aer = 4 * FT(π) * r_nuc^2
+        dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer)
         dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
     end
     if "P3_het" in ice_nucleation_modes
@@ -322,6 +322,7 @@ function run_parcel(IC, t_0, t_end, p)
     FT = eltype(IC)
     (; const_dt, ice_nucleation_modes, growth_modes, droplet_size_distribution) = p
     timestepper = ODE.Tsit5()
+    N_nuc_modes = 0
 
     println(" ")
     println("Running parcel model with: ")
@@ -331,20 +332,25 @@ function run_parcel(IC, t_0, t_end, p)
         print("\n    Deposition on dust particles using AF ")
         print("(Note that this mode only runs for monodisperse size distributions.) ")
         timestepper = ODE.Euler()
+        N_nuc_modes += 1
     end
     if "MohlerRate_Deposition" in ice_nucleation_modes
         print("\n    Deposition nucleation based on rate eqn in Mohler 2006 ")
         print("(Note that this mode only runs for monodisperse size distributions.) ")
         timestepper = ODE.Euler()
+        N_nuc_modes += 1
     end
     if "ActivityBasedDeposition" in ice_nucleation_modes
         print("\n    Water activity based deposition nucleation ")
         print("(Note that this mode only runs for monodisperse size distributions.) ")
+        N_nuc_modes += 1
     end
     if "P3_Deposition" in ice_nucleation_modes
         print("\n    P3 deposition nucleation ")
         timestepper = ODE.Euler()
+        N_nuc_modes += 1
     end
+
     if "ImmersionFreezing" in ice_nucleation_modes
         print("\n    Immersion freezing (ABIFM) ")
         print("(Note that this mode only runs for monodisperse and gamma size distributions.) ")
@@ -386,10 +392,7 @@ 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 ||
-        "P3_Deposition" in ice_nucleation_modes
+    if N_nuc_modes > 1
         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