From 3ce2f0dfc91213660a5bbb94d76da22f5ae6b0e0 Mon Sep 17 00:00:00 2001 From: amylu00 Date: Fri, 26 Apr 2024 13:26:21 -0700 Subject: [PATCH] AIDA ice nucleation calibrations --- papers/Project.toml | 2 + .../ice_nucleation_2024/AIDA_calibrations.jl | 212 ++++++++++++++++++ papers/ice_nucleation_2024/calibration.jl | 81 +++++-- .../ice_nucleation_2024/calibration_setup.jl | 110 ++++++++- ...plot_ensemble_mean.jl => perfect_calib.jl} | 14 ++ parcel/ParcelModel.jl | 5 +- parcel/ParcelParameters.jl | 1 + parcel/ParcelTendencies.jl | 13 +- src/ArtifactCalling.jl | 2 + 9 files changed, 410 insertions(+), 30 deletions(-) create mode 100644 papers/ice_nucleation_2024/AIDA_calibrations.jl rename papers/ice_nucleation_2024/{plot_ensemble_mean.jl => perfect_calib.jl} (76%) diff --git a/papers/Project.toml b/papers/Project.toml index 3fca2a3da1..59b2eaba70 100644 --- a/papers/Project.toml +++ b/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" diff --git a/papers/ice_nucleation_2024/AIDA_calibrations.jl b/papers/ice_nucleation_2024/AIDA_calibrations.jl new file mode 100644 index 0000000000..88f3398d75 --- /dev/null +++ b/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 diff --git a/papers/ice_nucleation_2024/calibration.jl b/papers/ice_nucleation_2024/calibration.jl index f4cc46773d..6501543b09 100644 --- a/papers/ice_nucleation_2024/calibration.jl +++ b/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) diff --git a/papers/ice_nucleation_2024/calibration_setup.jl b/papers/ice_nucleation_2024/calibration_setup.jl index 49e146d2f0..6a42fc5b77 100644 --- a/papers/ice_nucleation_2024/calibration_setup.jl +++ b/papers/ice_nucleation_2024/calibration_setup.jl @@ -8,6 +8,11 @@ function perf_model_params(FT, IN_mode) if IN_mode == "ABDINM" const_dt = FT(1) w = FT(0.4) + t_max = FT(100) + aerosol_act = "None" + aerosol = "None" + aero_σ_g = FT(0) + r_nuc = FT(1e-7) dep_nucleation = "ABDINM" heterogeneous = "None" homogeneous = "None" @@ -18,6 +23,11 @@ function perf_model_params(FT, IN_mode) elseif IN_mode == "ABIFM" const_dt = FT(1) w = FT(0.4) + t_max = FT(100) + aerosol_act = "None" + aerosol = "None" + aero_σ_g = FT(0) + r_nuc = FT(1e-7) dep_nucleation = "None" heterogeneous = "ABIFM" homogeneous = "None" @@ -28,6 +38,11 @@ function perf_model_params(FT, IN_mode) elseif IN_mode == "ABHOM" const_dt = FT(1) w = FT(1) + t_max = FT(100) + aerosol_act = "None" + aerosol = "None" + aero_σ_g = FT(0) + r_nuc = FT(1e-7) dep_nucleation = "None" heterogeneous = "None" homogeneous = "ABHOM" @@ -36,7 +51,8 @@ function perf_model_params(FT, IN_mode) liq_size_distribution = "Monodisperse" ice_size_distribution = "Monodisperse" end - params = (; const_dt, w, + params = (; const_dt, w, t_max, + aerosol_act, aerosol, r_nuc, aero_σ_g, # aerosol activation condensation_growth, deposition_growth, # growth liq_size_distribution, ice_size_distribution, # size distribution dep_nucleation, heterogeneous, homogeneous, # nucleation @@ -128,4 +144,96 @@ function perf_model_pseudo_data(FT, IN_mode, params, IC) y_truth = G_truth .+ rand(noise_dist) return [y_truth, Γ, coeff_true] end + +function AIDA_IN05_params(FT, w, t_max) + IN_mode = "ABHOM" + const_dt = FT(1) + aerosol_act = "AeroAct" + aerosol = CMP.Sulfate(FT) + dep_nucleation = "None" + heterogeneous = "None" + homogeneous = "ABHOM" + condensation_growth = "Condensation" + deposition_growth = "Deposition" + liq_size_distribution = "Gamma" + ice_size_distribution = "Gamma" + aero_σ_g = FT(2.3) + r_nuc = FT(1e-7) #FT(3.057e-6) + + params = (; const_dt, w, t_max, + aerosol_act, aerosol, r_nuc, aero_σ_g, # aerosol activation + condensation_growth, deposition_growth, # growth + liq_size_distribution, ice_size_distribution, # size distribution + dep_nucleation, heterogeneous, homogeneous, # ice nucleation + ) + return params +end + +function AIDA_IN05_IC(FT, data_file) + tps = TD.Parameters.ThermodynamicsParameters(FT) + wps = CMP.WaterProperties(FT) + + ρₗ = wps.ρw + ρᵢ = wps.ρi + R_d = TD.Parameters.R_d(tps) + R_v = TD.Parameters.R_v(tps) + ϵₘ = R_d / R_v + + if data_file == "in05_17_aida.edf" + # starting at t = 205 s (to match moving average freezing onset) + Nₗ = FT(297.136 * 1e6) + Nᵢ = FT(1.49 * 1e6) + Nₐ = FT(360 * 1e6) - Nₗ - Nᵢ + r₀ = FT(6.17664e-6) # FT(1e-7) + p₀ = FT(865.179 * 1e2) + T₀ = FT(236.91) + qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ + qᵢ = FT(Nᵢ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) + m_l = Nₗ * ρₗ * 4 * π / 3 * r₀^3 + m_i = Nᵢ * ρᵢ * 4 * π / 3 * r₀^3 + e = FT(27.0341) + qᵥ = (e / R_v / T₀) / ((p₀ - e) / (R_d * T₀) + e / R_v / T₀ + m_l + m_i) + q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rₐ = TD.gas_constant_air(tps, q) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + Sₗ = FT(e / eₛ) + elseif data_file == "in05_18_aida.edf" + Nₗ = FT(209.46 * 1e6) + Nᵢ = FT(1.53 * 1e6) + Nₐ = FT(275 * 1e6) - Nₗ - Nᵢ + r₀ = FT(7.03467 * 1e-6) + p₀ = FT(873.212 * 1e2) + T₀ = FT(237.134) + qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ + qᵢ = FT(Nᵢ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) + m_l = Nₗ * ρₗ * 4 * π / 3 * r₀^3 + m_i = Nᵢ * ρᵢ * 4 * π / 3 * r₀^3 + e = FT(28.1324) + qᵥ = (e / R_v / T₀) / ((p₀ - e) / (R_d * T₀) + e / R_v / T₀ + m_l + m_i) + q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rₐ = TD.gas_constant_air(tps, q) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + Sₗ = FT(e / eₛ) + elseif data_file == "in05_19_aida.edf" + Nₐ = FT(0) + Nₗ = FT(180 * 1e6) + Nᵢ = FT(0.49 * 1e6) + r₀ = FT(6.5e-6) # !!missing in dataset!! + p₀ = FT(722.852 * 1e2) + T₀ = FT(237.521) + qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ + qᵢ = FT(Nᵢ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) + m_l = Nₗ * ρₗ * 4 * π / 3 * r₀^3 + m_i = Nᵢ * ρᵢ * 4 * π / 3 * r₀^3 + e = FT(29.5341) # !!missing in dataset!! + qᵥ = (e / R_v / T₀) / ((p₀ - e) / (R_d * T₀) + e / R_v / T₀ + m_l + m_i) + q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) + Rₐ = TD.gas_constant_air(tps, q) + eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) + Sₗ = FT(e / eₛ) + end + return [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, FT(0)] #remove the last 2 elements, its J & r_l +end + + #! format: on diff --git a/papers/ice_nucleation_2024/plot_ensemble_mean.jl b/papers/ice_nucleation_2024/perfect_calib.jl similarity index 76% rename from papers/ice_nucleation_2024/plot_ensemble_mean.jl rename to papers/ice_nucleation_2024/perfect_calib.jl index d11e01789a..d22d80076e 100644 --- a/papers/ice_nucleation_2024/plot_ensemble_mean.jl +++ b/papers/ice_nucleation_2024/perfect_calib.jl @@ -66,5 +66,19 @@ for IN_mode in IN_mode_list plot_name = "perfect_calibration_$mode_label.svg" MK.save(plot_name, fig) + # Plotting calibrated parcel's ICNC if ABHOM + if IN_mode == "ABHOM" + fig2 = MK.Figure(size = (800, 600)) + ax3 = MK.Axis(fig2[1, 1], ylabel = "ICNC [m^3]", xlabel = "time [s]", title = IN_mode) + + soln = run_calibrated_model(FT, IN_mode, calibrated_parameters, params, IC) + soln_dflt = run_calibrated_model(FT, IN_mode, [FT(255.927125), FT(-68.553283)], params, IC) + + MK.lines!(ax3, soln.t, soln[9,:], label = "calibrated") + MK.lines!(ax3, soln_dflt.t, soln_dflt[9,:], label = "default", color = :orange) + plot_name = "perfect_calibration_ICNC_$mode_label.svg" + MK.save(plot_name, fig2) + end + end #! format: on diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl index 42d20ced0e..ad53d36963 100644 --- a/parcel/ParcelModel.jl +++ b/parcel/ParcelModel.jl @@ -17,6 +17,7 @@ Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT} deposition_growth = "None" liq_size_distribution = "Monodisperse" ice_size_distribution = "Monodisperse" + σ_g = FT(0) aerosol = Empty{FT}() aero_σ_g = FT(0) wps = CMP.WaterProperties(FT) @@ -285,7 +286,7 @@ function run_parcel(IC, t_0, t_end, pp) if pp.condensation_growth == "None" ce_params = Empty{FT}() elseif pp.condensation_growth == "Condensation" - ce_params = CondParams{FT}(pp.aps, pp.tps) + ce_params = CondParams{FT}(pp.aps, pp.tps, pp.const_dt) else throw("Unrecognized condensation growth mode") end @@ -304,7 +305,7 @@ function run_parcel(IC, t_0, t_end, pp) p = ( liq_distr = liq_distr, ice_distr = ice_distr, - aero_act_params, + aero_act_params = aero_act_params, dep_params = dep_params, imm_params = imm_params, hom_params = hom_params, diff --git a/parcel/ParcelParameters.jl b/parcel/ParcelParameters.jl index 9722c6a211..42312308de 100644 --- a/parcel/ParcelParameters.jl +++ b/parcel/ParcelParameters.jl @@ -79,6 +79,7 @@ end struct CondParams{FT} <: CMP.ParametersType{FT} aps::CMP.ParametersType{FT} tps::TDP.ThermodynamicsParameters{FT} + const_dt::FT end struct DepParams{FT} <: CMP.ParametersType{FT} diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl index d1bbe6bd07..d02fe24996 100644 --- a/parcel/ParcelTendencies.jl +++ b/parcel/ParcelTendencies.jl @@ -190,12 +190,9 @@ end function homogeneous_freezing(params::ABHOM, PSD_liq, state) FT = eltype(state) (; tps, ips, const_dt) = params - (; T, p_air, qᵥ, qₗ, qᵢ, Nₗ) = state + (; T, p_air, qᵥ, qₗ, qᵢ, Nₗ, Sₗ) = state - q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) - Rᵥ = TD.Parameters.R_v(tps) - R_air = TD.gas_constant_air(tps, q) - e = eᵥ(qᵥ, p_air, R_air, Rᵥ) + e = TD.saturation_vapor_pressure(tps, T, TD.Liquid()) * Sₗ Δa_w = CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T) J = CMI_hom.homogeneous_J_linear(ips.homogeneous, Δa_w) @@ -218,10 +215,12 @@ end function condensation(params::CondParams, PSD_liq, state, ρ_air) FT = eltype(state) - (; aps, tps) = params + (; aps, tps, const_dt) = params (; Sₗ, T, Nₗ, qₗ) = state Gₗ = CMO.G_func(aps, tps, T, TD.Liquid()) - return 4 * FT(π) / ρ_air * (Sₗ - 1) * Gₗ * PSD_liq.r * Nₗ + dqₗ_dt = 4 * FT(π) / ρ_air * (Sₗ - 1) * Gₗ * PSD_liq.r * Nₗ + + return qₗ + dqₗ_dt * const_dt > 0 ? dqₗ_dt : -qₗ / const_dt end function deposition(::Empty, PSD_ice, state, ρ_air) diff --git a/src/ArtifactCalling.jl b/src/ArtifactCalling.jl index 91956b253c..db42cccf00 100644 --- a/src/ArtifactCalling.jl +++ b/src/ArtifactCalling.jl @@ -9,7 +9,9 @@ export AIDA_ice_nucleation """ AIDA_ice_nucleation(data_file_name) + - `data_file_name` - name of the data file on Caltech box. + Returns the filepath of the data file in Caltech box. """ function AIDA_ice_nucleation(data_file_name::String)