Merge pull request #561 from CliMA/al/evap_sub_sink

Sinks for evaporation and sublimation in parcel

Author: amylu00

papers/ice_nucleation_2024/AIDA_calibrations.jl

@@ -22,7 +22,7 @@ global edf_data_names = [
     "in07_01_aida.edf", "in07_19_aida.edf",
 ]
 data_file_names = [
-    ["in05_17_aida.edf", "in05_18_aida.edf", "TROPIC04"],
+    ["in05_17_aida.edf", "in05_18_aida.edf"],
     ["in07_01_aida.edf"],
     ["in07_19_aida.edf"],
     ["ACI04_22"],
@@ -33,7 +33,7 @@ data_file_names = [
 ]
 batch_names = ["HOM", "IN0701", "IN0719", "ACI04_22", "EXP19", "EXP45", "DEP", "IMM"]
 start_time_list = [     # freezing onset
-    [Int32(150), Int32(180), Int32(0)],
+    [Int32(150), Int32(180)],
     [Int32(50)],
     [Int32(35)],
     [Int32(0)],
@@ -43,7 +43,7 @@ start_time_list = [     # freezing onset
     [Int32(0), Int32(0)],
 ]
 end_time_list = [       # approximate time freezing stops
-    [Int32(295), Int32(290), Int32(700)],
+    [Int32(295), Int32(290)],
     [Int32(375)],
     [Int32(375)],
     [Int32(300)],
@@ -54,7 +54,7 @@ end_time_list = [       # approximate time freezing stops
 ]
 moving_average_n = 5    # average every length(data) / n points
 updrafts = [            # updrafts matching AIDA cooling rate
-    [FT(1.5), FT(1.4), FT(1.4)],
+    [FT(1.5), FT(1.4)],
     [FT(1.5)],
     [FT(1.5)],
     [FT(1.5)],
@@ -146,7 +146,7 @@ for (batch_index, batch_name) in enumerate(batch_names)
     ### Plot for individual experiments.
     for (exp_index, data_file) in enumerate(data_file_name_list)
         if batch_name == "HOM"
-            exp_names = ["HOM_Batch_IN0517", "HOM_Batch_IN0518", "HOM_Batch_TROPIC04"]
+            exp_names = ["HOM_Batch_IN0517", "HOM_Batch_IN0518"]
         elseif batch_name == "DEP"
             exp_names = ["DEP_Batch_IN0701", "DEP_Batch_IN0719", "DEP_Batch_EXP45"]
         elseif batch_name == "IMM"

papers/ice_nucleation_2024/calibration_setup.jl

@@ -164,7 +164,7 @@ function AIDA_IN05_params(FT, w, t_max, t_profile, T_profile, P_profile; P3 = fa
     IN_mode = "ABHOM"
     const_dt = FT(1)
     prescribed_thermodynamics = true
-    aerosol_act = "AeroAct"
+    aerosol_act = "None"
     aerosol = CMP.Sulfate(FT)
     dep_nucleation = P3 ? "P3_dep" : "ABDINM"
     heterogeneous = P3 ? "P3_het" : "ABIFM"
@@ -200,9 +200,9 @@ function AIDA_IN05_IC(FT, data_file)
     ϵₘ = R_d / R_v
 
     if data_file == "in05_17_aida.edf"
-        Nₗ = FT(320.876 * 1e6)
         Nᵢ = FT(0.05 * 1e6)
-        Nₐ = FT(360 * 1e6) - Nₗ - Nᵢ
+        Nₗ = FT(360 * 1e6) -  Nᵢ
+        Nₐ = FT(0)
         r₀ = FT(5.18056 / 2 * 1e-6)
         p₀ = FT(894.409 * 1e2)
         T₀ = FT(237.871)
@@ -217,9 +217,9 @@ function AIDA_IN05_IC(FT, data_file)
         eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
         Sₗ = FT(e / eₛ)
     elseif data_file == "in05_18_aida.edf"
-        Nₗ = FT(212.885 * 1e6)
         Nᵢ = FT(0.787 * 1e6)
-        Nₐ = FT(275 * 1e6) - Nₗ - Nᵢ
+        Nₗ = FT(275 * 1e6)- Nᵢ
+        Nₐ = FT(0)
         r₀ = FT(6.83925 / 2 * 1e-6)
         p₀ = FT(878.345 * 1e2)
         T₀ = FT(237.234)

papers/ice_nucleation_2024/plots.jl

@@ -114,9 +114,9 @@ function plot_calibrated_parcels(
     calibrated_parcel_fig = MK.Figure(size = (1500, 800), fontsize = 20)
     ax_parcel_1 =
         MK.Axis(calibrated_parcel_fig[1, 1], ylabel = "saturation [-]", xlabel = "time [s]", title = "$exp_name")
-    ax_parcel_2 = MK.Axis(calibrated_parcel_fig[2, 1], ylabel = "liq mixing ratio [g/kg]", xlabel = "time [s]")
+    ax_parcel_2 = MK.Axis(calibrated_parcel_fig[2, 1], ylabel = "qₗ [kg/kg]", xlabel = "time [s]")
     ax_parcel_3 = MK.Axis(calibrated_parcel_fig[1, 2], ylabel = "temperature [K]", xlabel = "time [s]")
-    ax_parcel_4 = MK.Axis(calibrated_parcel_fig[2, 2], ylabel = "qᵢ [g/kg]", xlabel = "time [s]")
+    ax_parcel_4 = MK.Axis(calibrated_parcel_fig[2, 2], ylabel = "qᵢ [kg/kg]", xlabel = "time [s]")
     ax_parcel_5 = MK.Axis(calibrated_parcel_fig[3, 1], ylabel = "Nₗ [m^-3]", xlabel = "time [s]")
     ax_parcel_6 = MK.Axis(calibrated_parcel_fig[3, 2], ylabel = "Nᵢ [m^-3]", xlabel = "time [s]")
     ax_parcel_7 = MK.Axis(calibrated_parcel_fig[1, 3], ylabel = "pressure [Pa]", xlabel = "time [s]")
@@ -309,12 +309,12 @@ function plot_ICNC_overview(overview_data)
     overview_fig = MK.Figure(size = (1000, 800), fontsize = 24)
 
     exp_name_list = [
-        "IN05_17", "IN05_18", "TROPIC04",
+        "IN05_17", "IN05_18",
         "IN07_01", "IN07_19", "EXP45",
         "ACI04_22", "EXP19",
     ]
     position = [
-        overview_fig[1, 1], overview_fig[2, 1], overview_fig[3, 1],
+        overview_fig[1, 1], overview_fig[2, 1],
         overview_fig[1, 2], overview_fig[2, 2], overview_fig[3, 2],
         overview_fig[1, 3], overview_fig[2, 3],
     ]

parcel/ParcelModel.jl

@@ -122,13 +122,15 @@ function parcel_model(dY, Y, p, t)
 
     # Condensation/evaporation
     dqₗ_dt_ce = condensation(ce_params, PSD_liq, state, ρ_air)
+    dNₗ_dt_ce = dqₗ_dt_ce < FT(0) && qₗ > FT(0) && qₗ < FT(1e-6) ? Nₗ * dqₗ_dt_ce / qₗ : FT(0)
     # Deposition/sublimation
     dqᵢ_dt_ds = deposition(ds_params, PSD_ice, state, ρ_air)
+    dNᵢ_dt_ds = dqᵢ_dt_ds < FT(0) && qᵢ > FT(0) && qᵢ < FT(1e-6) ? Nᵢ * dqᵢ_dt_ds / qᵢ : FT(0)
 
     # number concentration and ...
-    dNᵢ_dt = dNᵢ_dt_dep + dNᵢ_dt_imm + dNᵢ_dt_hom
+    dNᵢ_dt = dNᵢ_dt_dep + dNᵢ_dt_imm + dNᵢ_dt_hom + dNᵢ_dt_ds
     dNₐ_dt = -dNᵢ_dt_dep - dNₗ_dt_act
-    dNₗ_dt = dNₗ_dt_act - dNᵢ_dt_imm - dNᵢ_dt_hom
+    dNₗ_dt = dNₗ_dt_act - dNᵢ_dt_imm - dNᵢ_dt_hom + dNₗ_dt_ce
     # ... water mass budget
     dqₗ_dt_v2l = dqₗ_dt_ce + dqₗ_dt_act
     dqᵢ_dt_l2i = dqᵢ_dt_imm + dqᵢ_dt_hom