Skip to content

Commit

Permalink
mission accomplished
Browse files Browse the repository at this point in the history
  • Loading branch information
@amylu00
amylu00 committed Feb 6, 2024
1 parent 6ebb1c5 commit 973d9c4
Show file tree
Hide file tree
Showing 3 changed files with 136 additions and 91 deletions.
14 changes: 11 additions & 3 deletions docs/src/IceNucleationParcel0D.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,11 @@ where:
- ``N_{act}`` is the number of activated ice particles.

## Deposition Nucleation on dust particles
The nucleation rate due to deposition nucleation on dust particles is parameterized as described in the [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/#Water-Activity-Based-Deposition-Nucleation) section.
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
Expand Down Expand Up @@ -234,13 +238,17 @@ 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. 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.

```@example
include("../../parcel/Tully_et_al_2023.jl")
Expand Down
188 changes: 104 additions & 84 deletions parcel/Tully_et_al_2023.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,76 +74,11 @@ 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 = [["MohlerAF_Deposition"], ["MohlerRate_Deposition"]] # ice nucleation modes
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]")
Expand All @@ -160,23 +95,108 @@ 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"]
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

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
Expand Down
25 changes: 21 additions & 4 deletions parcel/parcel.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,17 @@ function parcel_model(dY, Y, p, t)

dN_act_dt_depo = FT(0)
dqi_dt_new_depo = FT(0)
if "DustDeposition" in ice_nucleation_modes
if "MohlerAF_Deposition" in ice_nucleation_modes
AF = CMI_het.dust_activated_number_fraction(
aerosol,
ip.deposition,
S_i,
T,
)
dN_act_dt_depo = max(FT(0), AF * N_aer - N_ice) / const_dt
dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
end
if "MohlerRate_Deposition" in ice_nucleation_modes
a = T > ip.deposition.T_thr ? aerosol.a_warm : aerosol.a_cold
dN_act_dt_depo = max(FT(0), N_aer * a * (dY[1] * ξ))
dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
Expand Down Expand Up @@ -263,14 +273,21 @@ function run_parcel(IC, t_0, t_end, p)
println(" ")
println("Running parcel model with: ")
print("Ice nucleation modes: ")
if "DustDeposition" in ice_nucleation_modes
print("Deposition on dust particles ")
if "MohlerAF_Deposition" in ice_nucleation_modes
print("Deposition on dust particles using AF ")
timestepper = ODE.Euler()
end
if "MohlerRate_Deposition" in ice_nucleation_modes
print("Deposition nucleation based on rate eqn in Mohler 2006 ")
timestepper = ODE.Tsit5()
end
if "ImmersionFreezing" in ice_nucleation_modes
print("Immersion freezing ")
timestepper = ODE.Tsit5()
end
if "HomogeneousFreezing" in ice_nucleation_modes
print("Homogeneous freezing ")
timestepper = ODE.Tsit5()
end
print("\n")
print("Growth modes: ")
Expand Down Expand Up @@ -301,7 +318,7 @@ function run_parcel(IC, t_0, t_end, p)
problem = ODE.ODEProblem(parcel_model, IC, (FT(t_0), FT(t_end)), p)
sol = ODE.solve(
problem,
ODE.Euler(),
timestepper,
dt = const_dt,
reltol = 10 * eps(FT),
abstol = 10 * eps(FT),
Expand Down

0 comments on commit 973d9c4

Please sign in to comment.