Skip to content

Commit

Permalink
Refactor parcel
Browse files Browse the repository at this point in the history
  • Loading branch information
@trontrytel
trontrytel committed Feb 23, 2024
1 parent 7fbe344 commit b7d6097
Show file tree
Hide file tree
Showing 17 changed files with 971 additions and 916 deletions.
20 changes: 10 additions & 10 deletions docs/src/IceNucleationParcel0D.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,7 @@ There are multiple ways of running deposition nucleation in the parcel.
above this value will result in nucleation in a different mode.
`"ActivityBasedDeposition"` will trigger a water activity based approach
from [Alpert2022](@cite). In this approach, ice production rate ``P_{ice, depo}``
is calculated from
is calculated from
```math
\begin{equation}
P_{ice, depo} = \left[ \frac{dN_i}{dt} \right]_{depo} = J_{depo}\;A_{aero}\;N_{aero}
Expand Down Expand Up @@ -252,14 +252,14 @@ 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 using the `"MohlerAF_Deposition"` approach
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.
The pale blue line uses the `"MohlerRate_Deposition"` approach.
We only run it for the first GCM timestep because the rate approach requires
the change in ice saturation over time. With the discontinuous jump in saturation,
the parameterization is unable to determine a proper nucleation rate. When we force
Expand All @@ -268,7 +268,7 @@ The pale blue line uses the `"MohlerRate_Deposition"` approach.
in the `"MohlerAF_Deposition"` approach for the first GCM timestep.

```@example
include("../../parcel/Tully_et_al_2023.jl")
include("../../parcel/Example_Tully_et_al_2023.jl")
```
![](cirrus_box.svg)

Expand All @@ -280,7 +280,7 @@ The water activity based parameterization for deposition nucleation shows
is no common aerosol type between the two parameterizations.

```@example
include("../../parcel/Deposition_Nucleation.jl")
include("../../parcel/Example_Deposition_Nucleation.jl")
```
![](deposition_nucleation.svg)

Expand All @@ -289,7 +289,7 @@ In the plots below, the parcel model is ran with only condensation (no ice or fr
It is compared to [Rogers1975](@cite).

```@example
include("../../parcel/Liquid_only.jl")
include("../../parcel/Example_Liquid_only.jl")
```
![](liquid_only_parcel.svg)

Expand All @@ -299,7 +299,7 @@ The plots below are the results of the adiabatic parcel model
yet been validated against literature.

```@example
include("../../parcel/Immersion_Freezing.jl")
include("../../parcel/Example_Immersion_Freezing.jl")
```
![](immersion_freezing.svg)

Expand All @@ -312,13 +312,13 @@ The following plots show the parcel model running with homogeneous freezing and
for demonstrative purposes.

```@example
include("../../parcel/Jensen_et_al_2022.jl")
include("../../parcel/Example_Jensen_et_al_2022.jl")
```
![](Jensen_et_al_2022.svg)

## P3 Ice Nucleation Parameterizations
The parcel also includes ice nucleation parameterizations used in
the P3 scheme as described in [MorrisonMilbrandt2015](@cite).
the P3 scheme as described in [MorrisonMilbrandt2015](@cite).
Deposition nucleation is based on the ice crystal number parameterization
from Cooper (1986). The heterogeneous freezing parameterization, which
follows Bigg(1953) with parameters from Barklie aand Gokhale (1959), is
Expand All @@ -330,6 +330,6 @@ Shown below are three separate parcel simulations for deposition nucleation,
INP number, while the other two freezing modes do. Updraft velocity is
set to 0.5 m/s.
```@example
include("../../parcel/P3_ice_nuc.jl")
include("../../parcel/Example_P3_ice_nuc.jl")
```
![](P3_ice_nuc.svg)
92 changes: 32 additions & 60 deletions parcel/Deposition_Nucleation.jl → parcel/Example_Deposition_Nucleation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,18 +5,11 @@ import CloudMicrophysics as CM
import CLIMAParameters as CP

# definition of the ODE problem for parcel model
include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
FT = Float32

# get free parameters
tps = TD.Parameters.ThermodynamicsParameters(FT)
aps = CMP.AirProperties(FT)
wps = CMP.WaterProperties(FT)
ip = CMP.IceNucleationParameters(FT)

# Constants
ρₗ = wps.ρw
R_v = TD.Parameters.R_v(tps)
R_d = TD.Parameters.R_d(tps)

# Initial conditions
Nₐ = FT(2000 * 1e3)
Expand All @@ -34,60 +27,45 @@ q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
ρₐ = TD.air_density(tps, ts)
Rₐ = TD.gas_constant_air(tps, q)
Rᵥ = TD.Parameters.R_v(tps)
eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
e = qᵥ * p₀ * R_v / Rₐ
e = eᵥ(qᵥ, p₀, Rₐ, Rᵥ)

Sₗ = FT(e / eₛ)
IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]

ξ(T) =
TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
TD.saturation_vapor_pressure(tps, T, TD.Ice())
S_i(T, S_liq) = ξ(T) * S_liq

# Simulation parameters passed into ODE solver
r_nuc = FT(1.25e-6) # assumed size of nucleated particles
w = FT(3.5 * 1e-2) # updraft speed
α_m = FT(0.5) # accomodation coefficient
const_dt = FT(0.1) # model timestep
t_max = FT(100)
aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)] # aersol type for DustDeposition
aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)] # aersol type for DepositionNucleation
ice_nucleation_modes_list =
[["MohlerRate_Deposition"], ["ActivityBasedDeposition"]]
aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)] # aersol type for DustDeposition
aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)] # aersol type for DepositionNucleation
deposition_modes = ["MohlerRate", "ActivityBased"]
growth_modes = ["Deposition"]
droplet_size_distribution_list = [["Monodisperse"]]
size_distribution = "Monodisperse"

# Plotting
fig = MK.Figure(resolution = (800, 600))
fig = MK.Figure(size = (800, 600))
ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]")
ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]", xlabel = "time")
ax4 = MK.Axis(fig[2, 2], ylabel = "N_ice [m^-3]", xlabel = "time")

for ice_nucleation_modes in ice_nucleation_modes_list
nuc_mode = ice_nucleation_modes[1]
droplet_size_distribution = droplet_size_distribution_list[1]

if nuc_mode == "MohlerRate_Deposition"
for deposition in deposition_modes
if deposition == "MohlerRate"
for aerosol in aerosol_1
p = (;
wps,
aps,
tps,
ip,
const_dt,
r_nuc,
w,
α_m,
aerosol,
ice_nucleation_modes,
growth_modes,
droplet_size_distribution,
params = parcel_params{FT}(
const_dt = const_dt,
r_nuc = r_nuc,
w = w,
aerosol = aerosol,
deposition = deposition,
growth_modes = growth_modes,
size_distribution = size_distribution,
)

# solve ODE
sol = run_parcel(IC, FT(0), t_max, p)
sol = run_parcel(IC, FT(0), t_max, params)

# Plot results
if aerosol == CMP.DesertDust(FT)
Expand All @@ -98,33 +76,27 @@ for ice_nucleation_modes in ice_nucleation_modes_list
MK.lines!( # saturation
ax1,
sol.t,
S_i.(sol[3, :], sol[1, :]),
S_i.(tps, sol[3, :], sol[1, :]),
label = aero_label,
)
MK.lines!(ax2, sol.t, sol[3, :]) # temperature
MK.lines!(ax3, sol.t, sol[6, :] * 1e3) # q_ice
MK.lines!(ax4, sol.t, sol[9, :]) # N_ice
end

elseif nuc_mode == "ActivityBasedDeposition"
elseif deposition == "ActivityBased"
for aerosol in aerosol_2
p = (;
wps,
aps,
tps,
ip,
const_dt,
r_nuc,
w,
α_m,
aerosol,
ice_nucleation_modes,
growth_modes,
droplet_size_distribution,
params = parcel_params{FT}(
const_dt = const_dt,
r_nuc = r_nuc,
w = w,
aerosol = aerosol,
deposition = deposition,
growth_modes = growth_modes,
size_distribution = size_distribution,
)

# solve ODE
sol = run_parcel(IC, FT(0), t_max, p)
sol = run_parcel(IC, FT(0), t_max, params)

# Plot results
if aerosol == CMP.Feldspar(FT)
Expand All @@ -137,7 +109,7 @@ for ice_nucleation_modes in ice_nucleation_modes_list
MK.lines!( # saturation
ax1,
sol.t,
S_i.(sol[3, :], sol[1, :]),
S_i.(tps, sol[3, :], sol[1, :]),
label = aero_label,
linestyle = :dash,
)
Expand Down
66 changes: 18 additions & 48 deletions parcel/Immersion_Freezing.jl → parcel/Example_Immersion_Freezing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,20 +5,14 @@ import CloudMicrophysics as CM
import CLIMAParameters as CP

# definition of the ODE problem for parcel model
include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
FT = Float32
# get free parameters
tps = TD.Parameters.ThermodynamicsParameters(FT)
aps = CMP.AirProperties(FT)
wps = CMP.WaterProperties(FT)
ip = CMP.IceNucleationParameters(FT)

# Constants
ρₗ = wps.ρw
R_v = TD.Parameters.R_v(tps)
R_d = TD.Parameters.R_d(tps)

# Initial conditions
ρₗ = wps.ρw
Nₐ = FT(0)
Nₗ = FT(500 * 1e3)
Nᵢ = FT(0)
Expand All @@ -31,26 +25,23 @@ qᵢ = FT(0)
x_sulph = FT(0.01)

# Moisture dependent initial conditions
q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
ρₐ = TD.air_density(tps, ts)
q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
R_v = TD.Parameters.R_v(tps)
Rₐ = TD.gas_constant_air(tps, q)
eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
e = qᵥ * p₀ * R_v / Rₐ

e = eᵥ(qᵥ, p₀, Rₐ, R_v)
Sₗ = FT(e / eₛ)
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
α_m = FT(0.5) # accomodation coefficient
const_dt = FT(1) # model timestep
t_max = FT(1200)
aerosol = CMP.Illite(FT)
ice_nucleation_modes = ["ImmersionFreezing"]
heterogeneous = "ImmersionFreezing"
growth_modes = ["Condensation", "Deposition"]
droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]]
size_distribution_list = ["Monodisperse", "Gamma"]

# Plotting
fig = MK.Figure(resolution = (800, 600))
Expand All @@ -67,47 +58,26 @@ ax6 = MK.Axis(
)
MK.ylims!(ax6, 1e-6, 1)

for droplet_size_distribution in droplet_size_distribution_list
p = (;
wps,
aps,
tps,
ip,
const_dt,
r_nuc,
w,
α_m,
aerosol,
ice_nucleation_modes,
growth_modes,
droplet_size_distribution,
for DSD in size_distribution_list
params = parcel_params{FT}(
const_dt = const_dt,
r_nuc = r_nuc,
w = w,
aerosol = aerosol,
heterogeneous = heterogeneous,
growth_modes = growth_modes,
size_distribution = DSD,
)
# solve ODE
sol = run_parcel(IC, FT(0), t_max, p)

DSD = droplet_size_distribution[1]
sol = run_parcel(IC, FT(0), t_max, params)

# Plot results
ξ(T) =
TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
TD.saturation_vapor_pressure(tps, T, TD.Ice())
S_i(T, S_liq) = ξ(T) * S_liq - 1

MK.lines!(ax1, sol.t / 60, S_i.(sol[3, :], sol[1, :]), label = DSD)
MK.lines!(ax1, sol.t / 60, S_i.(tps, sol[3, :], sol[1, :]) .- 1, label = DSD)
MK.lines!(ax2, sol.t / 60, sol[3, :])
MK.lines!(ax3, sol.t / 60, sol[6, :] * 1e3)
MK.lines!(ax4, sol.t / 60, sol[5, :] * 1e3)

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
MK.lines!(ax5, sol.t / 60, sol_Nₗ)
MK.lines!(ax6, sol.t / 60, sol_Nᵢ ./ (sol_Nₗ .+ sol_Nᵢ))
end
Expand Down
Loading

0 comments on commit b7d6097

Please sign in to comment.