AIDA ice nucleation calibrations

papers/Project.toml

@@ -1,10 +1,12 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
+ClimaUtilities = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513"
 CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"

papers/ice_nucleation_2024/AIDA_calibrations.jl

@@ -0,0 +1,212 @@
+import ClimaParams as CP
+import CloudMicrophysics as CM
+import CloudMicrophysics.Parameters as CMP
+import Thermodynamics as TD
+import CairoMakie as MK
+
+# To grab data
+import DelimitedFiles
+using LazyArtifacts
+using ClimaUtilities.ClimaArtifacts
+
+FT = Float64
+include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration.jl"))
+ips = CMP.IceNucleationParameters(FT)
+
+function unpack_data(data_file_name)
+    file_path = CM.ArtifactCalling.AIDA_ice_nucleation(data_file_name)
+    return DelimitedFiles.readdlm(file_path, skipstart = 125)
+end
+moving_average(data, n) =
+    [sum(@view data[i:(i + n - 1)]) / n for i in 1:(length(data) - (n - 1))]
+
+data_file_names = ["in05_17_aida.edf", "in05_18_aida.edf", "in05_19_aida.edf"]
+plot_names = ["IN0517", "IN0518", "IN0519"]
+
+start_time_list = [195, 180, 80] .+ 100   # AIDA time starts at -100 seconds. Add freezing onset as needed.
+end_time_list = [290, 290, 170] .+ 100     # approximate time freezing stops
+moving_average_n = 20                      # average every n points
+
+updrafts = [FT(1.5), FT(1.4), FT(5)]
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+R_v = TD.Parameters.R_v(tps)
+R_d = TD.Parameters.R_d(tps)
+ϵₘ = R_d / R_v
+
+
+for (exp_index, data_file_name) in enumerate(data_file_names)
+    #! format: off
+    ### Unpacking experiment-specific variables
+    plot_name = plot_names[exp_index]
+    w = updrafts[exp_index]
+    start_time = start_time_list[exp_index]
+    end_time = end_time_list[exp_index]
+
+    ## Moving average to smooth data
+    moving_average_start = Int32(start_time + (moving_average_n / 2 - 1))
+    moving_average_end = Int32(end_time - (moving_average_n / 2))
+    t_max =
+        (end_time - 100 - (moving_average_n / 2)) -     # AIDA time starts at -100 seconds (t = 0 at index 101)
+        (start_time - 100 + (moving_average_n / 2 - 1)) # duration of simulation
+
+    params = AIDA_IN05_params(FT, w, t_max)
+    IC = AIDA_IN05_IC(FT, data_file_name)
+
+    ### Check for data file in AIDA_data folder.
+    chamber_data = unpack_data(data_file_name)
+    ICNC = chamber_data[start_time:end_time, 6] .* 1e6 # Nᵢ [m^-3]
+
+    frozen_frac = ICNC ./ (IC[7] + IC[8] + IC[9])
+    frozen_frac_moving_mean = moving_average(frozen_frac, moving_average_n)
+    ICNC_moving_avg = moving_average(ICNC, moving_average_n)
+    Γ = 0.15 * LinearAlgebra.I * (maximum(frozen_frac_moving_mean) - minimum(frozen_frac_moving_mean)) # coeff is an estimated of the noise
+
+    ### Plotting AIDA data.
+    AIDA_data_fig = MK.Figure(size = (800, 600), fontsize = 24)
+    ax1 = MK.Axis(AIDA_data_fig[1, 1], ylabel = "ICNC [m^-3]", xlabel = "time [s]", title = "AIDA data $plot_name")
+    MK.lines!(
+        ax1,
+        chamber_data[start_time:end_time, 1],
+        ICNC,
+        label = "Raw AIDA",
+        color =:blue,
+        linestyle =:dash,
+        linewidth = 2,
+    )
+    MK.lines!(
+        ax1,
+        chamber_data[moving_average_start:moving_average_end, 1],
+        ICNC_moving_avg,
+        label = "AIDA moving mean",
+        linewidth = 2.5,
+        color =:blue,
+    )
+    MK.axislegend(ax1, framevisible = true, labelsize = 12, position = :rc)
+    MK.save("$plot_name"*"_ICNC.svg", AIDA_data_fig)
+
+    ### Calibration.
+    output = calibrate_J_parameters(FT, "ABHOM", params, IC, frozen_frac_moving_mean, Γ)
+    n_iterations = size(output[3])[1]
+    n_ensembles = size(output[3][1])[2]
+    calibrated_parameters = [output[1], output[2]]
+    calibrated_ensemble_means = ensemble_means(output[3], n_iterations, n_ensembles)
+    # Calibrated parcel
+    parcel = run_calibrated_model(FT, "ABHOM", calibrated_parameters, params, IC)
+    #parcel_default = run_calibrated_model(FT, "ABHOM", [FT(255.927125), FT(-68.553283)], params, IC)
+
+    ### Plots.
+    ## Calibrated coefficients.
+    #  Did they converge?
+    calibrated_coeffs_fig = MK.Figure(size = (600, 600), fontsize = 24)
+    ax3 = MK.Axis(calibrated_coeffs_fig[1, 1], ylabel = "m coefficient [-]", title = "$plot_name")
+    ax4 = MK.Axis(calibrated_coeffs_fig[2, 1], ylabel = "c coefficient [-]", xlabel = "iteration #")
+    MK.lines!(ax3, collect(1:n_iterations), calibrated_ensemble_means[1], label = "ensemble mean", color = :orange, linewidth = 2.5)
+    MK.lines!(ax4, collect(1:n_iterations), calibrated_ensemble_means[2], label = "ensemble mean", color = :orange, linewidth = 2.5)
+    MK.save("$plot_name"*"_calibrated_coeffs_fig.svg", calibrated_coeffs_fig)
+
+    ## Calibrated parcel simulations.
+    #  Does the calibrated parcel give reasonable outputs?
+    tps = TD.Parameters.ThermodynamicsParameters(FT)
+    chamber_S_l = zeros(FT, size(chamber_data[start_time:end_time, 4], 1), size(chamber_data[start_time:end_time, 4], 2))
+    for (i, e) in enumerate(chamber_data[start_time:end_time, 4])
+        temp = chamber_data[start_time:end_time, 3][i]
+        eₛ = TD.saturation_vapor_pressure(tps, temp, TD.Liquid())
+        chamber_S_l[i] = e / eₛ
+    end
+    chamber_r_l = zeros(FT, size(chamber_data[start_time:end_time, 8], 1), size(chamber_data[start_time:end_time, 8], 2))
+    for (i, r_l) in enumerate(chamber_data[start_time:end_time, 8] .* 1e-6)
+        chamber_r_l[i] = r_l < 0 ? FT(0) : r_l
+    end
+
+    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 = "$plot_name")
+    ax_parcel_2 = MK.Axis(calibrated_parcel_fig[2, 1], ylabel = "liq mixing ratio [g/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_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]")
+    ax_parcel_8 = MK.Axis(calibrated_parcel_fig[2, 3], ylabel = "Nₐ [m^-3]", xlabel = "time [s]")
+    MK.lines!(ax_parcel_1, parcel.t .+ (moving_average_n / 2), parcel[1, :], label = "calibrated", color = :orange) # label = "liquid"
+    #MK.lines!(ax_parcel_1, parcel_default.t .+ (moving_average_n / 2), parcel_default[1, :], label = "default", color = :darkorange2)
+    MK.lines!(ax_parcel_1, parcel.t .+ (moving_average_n / 2), S_i.(tps, parcel[3, :], parcel[1, :]), label = "ice", color = :green)
+    # MK.lines!(
+    #     ax_parcel_1,
+    #     chamber_data[start_time:end_time, 1] .- (start_time - 100),
+    #     vec(chamber_S_l),
+    #     label = "chamber",
+    #     color = :blue
+    # )
+    MK.lines!(ax_parcel_2, parcel.t .+ (moving_average_n / 2), parcel[5, :], color = :orange)
+    #MK.lines!(ax_parcel_2, parcel_default.t .+ (moving_average_n / 2), parcel_default[5,:], color = :darkorange2)
+    MK.lines!(ax_parcel_3, parcel.t .+ (moving_average_n / 2), parcel[3, :], color = :orange)
+    #MK.lines!(ax_parcel_3, parcel_default.t .+ (moving_average_n / 2), parcel_default[3, :], color = :darkorange2)
+    MK.lines!(
+        ax_parcel_3,
+        chamber_data[start_time:end_time, 1] .- (start_time - 100),
+        chamber_data[start_time:end_time, 3],
+        color = :blue
+    )
+    MK.lines!(ax_parcel_4, parcel.t .+ (moving_average_n / 2), parcel[6, :], color = :orange)
+    #MK.lines!(ax_parcel_4, parcel_default.t .+ (moving_average_n / 2), parcel_default[6, :], color = :darkorange2)
+    MK.lines!(ax_parcel_5, parcel.t .+ (moving_average_n / 2), parcel[8, :], color = :orange)
+    #MK.lines!(ax_parcel_5, parcel_default.t .+ (moving_average_n / 2), parcel_default[8, :], color = :darkorange2)
+    MK.lines!(ax_parcel_6, parcel.t .+ (moving_average_n / 2), parcel[9, :], color = :orange)
+    # MK.lines!(ax_parcel_6, parcel_default.t .+ (moving_average_n / 2), parcel_default[9, :], color = :darkorange2)
+    MK.lines!(ax_parcel_6, 
+        chamber_data[start_time:end_time, 1] .- (start_time - 100),
+        ICNC,
+        color = :blue
+    )
+    MK.lines!(ax_parcel_7, parcel.t .+ (moving_average_n / 2), parcel[2, :], color = :orange)
+    MK.lines!(
+        ax_parcel_7,
+        chamber_data[start_time:end_time, 1] .- (start_time - 100),
+        chamber_data[start_time:end_time, 2] .* 1e2,
+        color = :blue
+    )
+    MK.lines!(ax_parcel_8, parcel.t .+ (moving_average_n / 2), parcel[7, :], color = :orange)
+    MK.save("$plot_name"*"_calibrated_parcel_fig.svg", calibrated_parcel_fig)
+
+    ## Comparing AIDA data and calibrated parcel.
+    #  Does calibrated parcel look like observations?
+    ICNC_comparison_fig = MK.Figure(size = (700, 600), fontsize = 24)
+    ax_compare = MK.Axis(ICNC_comparison_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "$plot_name")
+    # MK.lines!(
+    #     ax_compare,
+    #     parcel_default.t .+ (moving_average_n / 2),
+    #     parcel_default[9, :]./  (IC[7] + IC[8] + IC[9]),       # IC[8] if starting with droplets (no aerosol activation)
+    #     label = "CM.jl Parcel (Pre-Calib)",
+    #     linewidth = 2.5,
+    #     color =:orangered,
+    # )
+    MK.lines!(
+        ax_compare,
+        parcel.t .+ (moving_average_n / 2),
+        parcel[9, :]./ (IC[7] + IC[8] + IC[9]),               # IC[8] if starting with droplets (no aerosol activation)
+        label = "CM.jl Parcel (Calibrated)",
+        linewidth = 2.5,
+        color =:orange,
+    )
+    MK.lines!(
+        ax_compare,
+        chamber_data[start_time:end_time, 1] .- (start_time - 100),
+        frozen_frac,
+        label = "Raw AIDA",
+        linewidth = 2,
+        color =:blue,
+        linestyle =:dash,
+    )
+    MK.lines!(
+        ax_compare,
+        chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100),
+        frozen_frac_moving_mean,
+        label = "AIDA Moving Avg",
+        linewidth = 2.5,
+        color =:blue
+    )
+    MK.axislegend(ax_compare, framevisible = false, labelsize = 20, position = :lt)
+    MK.save("$plot_name"*"_ICNC_comparison_fig.svg", ICNC_comparison_fig)
+
+    #! format: on
+end

papers/ice_nucleation_2024/calibration.jl

@@ -19,14 +19,12 @@ include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration_setup
 function run_model(p, coefficients, IN_mode, FT, IC)
     # grabbing parameters
     m_calibrated, c_calibrated = coefficients
-    (; const_dt, w, deposition_growth) = p
-    (; liq_size_distribution, ice_size_distribution) = p
-
-    t_max = FT(100)
+    (; const_dt, w, t_max, aerosol_act, aerosol, r_nuc) = p
+    (; deposition_growth, condensation_growth) = p
+    (; dep_nucleation, heterogeneous, homogeneous) = p
+    (; liq_size_distribution, ice_size_distribution, aero_σ_g) = p
 
     if IN_mode == "ABDINM"
-        (; dep_nucleation) = p
-
         # overwriting
         override_file = Dict(
             "China2017_J_deposition_m_Kaolinite" =>
@@ -43,6 +41,7 @@ function run_model(p, coefficients, IN_mode, FT, IC)
             w = w,
             aerosol = overwrite,
             deposition = dep_nucleation,
+            condensation_growth = condensation_growth,
             deposition_growth = deposition_growth,
             liq_size_distribution = liq_size_distribution,
             ice_size_distribution = ice_size_distribution,
@@ -53,8 +52,6 @@ function run_model(p, coefficients, IN_mode, FT, IC)
         return sol[9, :] .* 1e10 # ICNC magnified
 
     elseif IN_mode == "ABIFM"
-        (; heterogeneous, condensation_growth) = p
-
         # overwriting
         override_file = Dict(
             "KnopfAlpert2013_J_ABIFM_m_Kaolinite" =>
@@ -82,8 +79,47 @@ function run_model(p, coefficients, IN_mode, FT, IC)
         return sol[9, :] # ICNC
 
     elseif IN_mode == "ABHOM"
-        (; homogeneous) = p
+        # overwriting
+        override_file = Dict(
+            "Linear_J_hom_coeff2" =>
+                Dict("value" => m_calibrated, "type" => "float"),
+            "Linear_J_hom_coeff1" =>
+                Dict("value" => c_calibrated, "type" => "float"),
+        )
+        ip_calibrated = CP.create_toml_dict(FT; override_file)
+        overwrite = CMP.IceNucleationParameters(ip_calibrated)
+
+        # run parcel with new coefficients
+        local params = parcel_params{FT}(
+            const_dt = const_dt,
+            w = w,
+            aerosol_act = aerosol_act,
+            aerosol = aerosol,
+            aero_σ_g = aero_σ_g,
+            homogeneous = homogeneous,
+            condensation_growth = condensation_growth,
+            deposition_growth = deposition_growth,
+            liq_size_distribution = liq_size_distribution,
+            ice_size_distribution = ice_size_distribution,
+            ips = overwrite,
+            r_nuc = r_nuc,
+        )
+
+        # solve ODE
+        local sol = run_parcel(IC, FT(0), t_max, params)
+        return sol[9, :] ./  (IC[7] + IC[8] + IC[9])  # frozen fraction
+    end
+end
 
+function run_calibrated_model(FT, IN_mode, coefficients, p, IC)
+    # grabbing parameters
+    m_calibrated, c_calibrated = coefficients
+    (; const_dt, w, t_max, aerosol_act, aerosol) = p
+    (; deposition_growth, condensation_growth) = p
+    (; homogeneous) = p
+    (; liq_size_distribution, ice_size_distribution, aero_σ_g) = p
+
+    if IN_mode == "ABHOM"
         # overwriting
         override_file = Dict(
             "Linear_J_hom_coeff2" =>
@@ -98,7 +134,11 @@ function run_model(p, coefficients, IN_mode, FT, IC)
         local params = parcel_params{FT}(
             const_dt = const_dt,
             w = w,
+            aerosol_act = aerosol_act,
+            aerosol = aerosol,
+            aero_σ_g = aero_σ_g,
             homogeneous = homogeneous,
+            condensation_growth = condensation_growth,
             deposition_growth = deposition_growth,
             liq_size_distribution = liq_size_distribution,
             ice_size_distribution = ice_size_distribution,
@@ -107,11 +147,10 @@ function run_model(p, coefficients, IN_mode, FT, IC)
 
         # solve ODE
         local sol = run_parcel(IC, FT(0), t_max, params)
-        return sol[9, :] # ICNC
+        return sol
     end
 end
 
-# Creating noisy pseudo-observations
 function create_prior(FT, IN_mode, ; perfect_model = false)
     # TODO - add perfect_model flag to plot_ensemble_mean.jl
     observation_data_names = ["m_coeff", "c_coeff"]
@@ -130,15 +169,17 @@ function create_prior(FT, IN_mode, ; perfect_model = false)
             c_stats = [FT(-70), FT(1), -Inf, Inf]
         end
     elseif perfect_model == false
-        if IN_mode == "ABDINM"      # TODO - dependent on dust type
-            m_stats = [FT(20), FT(1), FT(0), Inf]
-            c_stats = [FT(-1), FT(1), -Inf, Inf]
-        elseif IN_mode == "ABIFM"   # TODO - dependent on dust type
-            m_stats = [FT(50), FT(1), FT(0), Inf]
-            c_stats = [FT(-7), FT(1), -Inf, Inf]
+        if IN_mode == "ABDINM"
+            # m_stats = [FT(20), FT(1), FT(0), Inf]
+            # c_stats = [FT(-1), FT(1), -Inf, Inf]
+            println("Calibration for ABDINM with AIDA not yet implemented.")
+        elseif IN_mode == "ABIFM"
+            # m_stats = [FT(50), FT(1), FT(0), Inf]
+            # c_stats = [FT(-7), FT(1), -Inf, Inf]
+            println("Calibration for ABIFM with AIDA not yet implemented.")
         elseif IN_mode == "ABHOM"
-            m_stats = [FT(255.927125), FT(1), FT(0), Inf]
-            c_stats = [FT(-68.553283), FT(1), -Inf, Inf]
+            m_stats = [FT(260.927125), FT(25), FT(0), Inf]
+            c_stats = [FT(-68.553283), FT(10), -Inf, Inf]
         end
     end
 
@@ -167,7 +208,7 @@ function calibrate_J_parameters(FT, IN_mode, params, IC, y_truth, Γ,; perfect_m
 
     prior = create_prior(FT, IN_mode, perfect_model = perfect_model)
     N_ensemble = 10      # runs N_ensemble trials per iteration
-    N_iterations = 15    # number of iterations the inverse problem goes through
+    N_iterations = 150    # number of iterations the inverse problem goes through
 
     # Generate initial ensemble and set up EKI
     initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble)