Skip to content

Commit

Permalink
Add P3 collision and melting tendencies
Browse files Browse the repository at this point in the history
  • Loading branch information
@anastasia-popova
anastasia-popova committed May 22, 2024
1 parent 10505a9 commit e3788cc
Show file tree
Hide file tree
Showing 12 changed files with 1,314 additions and 105 deletions.
2 changes: 2 additions & 0 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ version = "0.20.1"
ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
Expand All @@ -24,6 +25,7 @@ ClimaParams = "0.10.6"
DataFrames = "1.6"
DocStringExtensions = "0.8, 0.9"
ForwardDiff = "0.10"
HCubature = "1.6"
MLJ = "0.20"
QuadGK = "2.9"
RootSolvers = "0.3, 0.4"
Expand Down
6 changes: 6 additions & 0 deletions docs/src/API.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,12 @@ P3Scheme.thresholds
P3Scheme.distribution_parameter_solver
```

# Terminal Velocity
```@docs
TerminalVelocity
TerminalVelocity.velocity_chen
```

# Aerosol model

```@docs
Expand Down
89 changes: 89 additions & 0 deletions docs/src/P3Scheme.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,3 +188,92 @@ They can be compared with Figure 2 from [MorrisonMilbrandt2015](@cite).
include("plots/P3TerminalVelocityPlots.jl")
```
![](MorrisonandMilbrandtFig2.svg)

## Collisions

Collisions are calculated through the following equation:

```math
\frac{dq_c}{dt} = \int_{0}^{\infty} \! \int_{0}^{\infty} \! \frac{E_ci}{\rho} N'_i(D_p) N'_c(D_c) A(D_i, D_c) m_{collected}(D_{collected}) |V_i(D_i) - V_c(D_c)| \mathrm{d}D_i \mathrm{d}D_c
```
where the values are defined as follows:
- subscript ``i`` - ice particles (either cloud ice or precipitating ice)
- subscript ``c`` - corresponding colliding species
- subscript ``collected`` - whichever species is being collected by the collisions
- ``E_ci`` - collision efficiency between particles
- ``N'_x`` - distribution of particles of x type
- ``A`` - collision kernel between ice and colliding species
- ``m_x(D)`` - mass of given species x at diameter D
- ``V_x(D)`` - Chen terminal velocity of given species x at diameter D

In our parametrization we take ``E_ci = 1`` and ``A(D_i, D_c) = \pi \frac{(D_i + D_c)^2}{4}``. We also assume that if the temperature is above freezing then ice particles get collected, otherwise the colliding species is collected by the ice.

Furthermore, the distributions of non-ice species are taken to be of the following forms:
- cloud water: ``N_0 D^8 e^{-\lambda D^3}``
- rain : ``N_0 e^{-\lambda D}``

``N_0`` and ``\lambda`` can be solved for respectively in each scheme through setting:
- ``N_c = \int_{0}^{\infty} N'_c(D) \mathrm{d}D``
- ``q_c = \int_{0}^{\infty} m_c(D) N'_c(D) \mathrm{d}D``
where ``m_c(D)`` is the mass of particle of size D. Rain and cloud water particles are assumed to be spherical therefore, ``m_c(D) = \pi \rho_w \frac{D^3}{6}`` for both (``\rho_w`` = density of water).

Collision rate comparisons between the P3 Scheme and the 1M microphysics scheme are shown below:

```@example
include("plots/P3TendenciesSandbox.jl")
```
![](CollisionComparisons.svg)

## Melting

Melting rates are calculated through the following equation:

```math
\frac{dq}{dt} = \frac{1}{2 \rho_a} \int_{0}^{\infty} \! \frac{dm(D)}{dt} N'(D) \mathrm{d}D
```

where:
- ``\frac{dm(D)}{dt} = \frac{dm(D)}{dD} \frac{dD}{dt}``
- ``m(D)`` - mass of ice particle with maximum dimension D
- ``\frac{dD}{dt} = \frac{4}{D \rho_w} \frac{K_{thermo}}{L_f} (T - T_{freeze}) F(D)``
- ``F(D) = a + b N_{sc}^{\frac{1}{3}} N_{Re}(D)^{\frac{1}{2}}`` - ventilation factor
- ``N_{sc} = \frac{v_{air}}{D_{vapor}}`` - Schmidt number
- ``N_{Re}(D) = \frac{D V_{term}(D)}{v_{air}}`` - Reynolds number
- ``V_{term}(D)`` - terminal velocity of ice particle with maximum dimension D

with constant values:
- ``\rho_a`` - density of air (1.2 kg / m^3)
- ``\rho_w`` - density of water (1000 kg / m^3)
- ``K_{thermo}`` - thermal conductivity of air
- ``L_f`` - latent heat of freezing
- ``T_{freeze}`` - freezing temperature (273.15 K)
- ``v_{air}`` - kinematic viscosity of air
- ``D_{vapor}`` - diffusivity of water
- ``a`` - 0.78
- ``b`` - 0.308

Melting rate comparisons between the P3 Scheme and the 1M Microphysics scheme are shown below:

![](MeltRateComparisons.svg)

## Heterogeneous Freezing

Heterogeneous freezing rates are calculated using the ABIFM parametrization defined in the Ice Nucleation folder. Specifically, we use the following equation (taking into account the distribution of rain particles):

```math
\frac{dN}{dt} = \int_{0}^{\infty} \! J_{ABIFM} A(D) N'(D) \mathrm{d}D
\frac{dq}{dt} = \int_{0}^{\infty} \! J_{ABIFM} A(D) N'(D) m(D) \mathrm{d}D
```

where
- ``J_{ABIFM}`` - calculated heterogeneous nucleation rate
- ``A(D)`` - surface area of a particle with maximum dimension D (``\pi D^2`` in the case of spherical rain particles)
- ``N'(D)`` - number distribution of rain particles
- ``m(D)`` - corresponding mass of a rain particle with dimension D

The change in number concentration and mass due to heterogeneous freezing are shown below for 249K, 250K, and 251K.

![](FreezeRateComparisons.svg)

![](MassFreezeRateComparisons.svg)
116 changes: 12 additions & 104 deletions docs/src/plots/P3SchemePlots.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,99 +9,7 @@ FT = Float64
const PSP3 = CMP.ParametersP3
p3 = CMP.ParametersP3(FT)

"""
mass_(p3, D, ρ, F_r)

- p3 - a struct with P3 scheme parameters
- D - maximum particle dimension [m]
- ρ - bulk ice density (ρ_i for small ice, ρ_g for graupel) [kg/m3]
- F_r - rime mass fraction [q_rim/q_i]
Returns mass as a function of size for differen particle regimes [kg]
"""
# for spherical ice (small ice or completely rimed ice)
mass_s(D::FT, ρ::FT) where {FT <: Real} = FT(π) / 6 * ρ * D^3
# for large nonspherical ice (used for unrimed and dense types)
mass_nl(p3::PSP3, D::FT) where {FT <: Real} = P3.α_va_si(p3) * D^p3.β_va
# for partially rimed ice
mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} =
P3.α_va_si(p3) / (1 - F_r) * D^p3.β_va

"""
p3_mass(p3, D, F_r, th)
- p3 - a struct with P3 scheme parameters
- D - maximum particle dimension
- F_r - rime mass fraction (q_rim/q_i)
- th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d)
Returns mass(D) regime, used to create figures for the docs page.
"""
function p3_mass(
p3::PSP3,
D::FT,
F_r::FT,
th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)),
) where {FT <: Real}
D_th = P3.D_th_helper(p3)
if D_th > D
return mass_s(D, p3.ρ_i) # small spherical ice
elseif F_r == 0
return mass_nl(p3, D) # large nonspherical unrimed ice
elseif th.D_gr > D >= D_th
return mass_nl(p3, D) # dense nonspherical ice
elseif th.D_cr > D >= th.D_gr
return mass_s(D, th.ρ_g) # graupel
elseif D >= th.D_cr
return mass_r(p3, D, F_r) # partially rimed ice
end
end

"""
A_(p3, D)
- p3 - a struct with P3 scheme parameters
- D - maximum particle dimension
Returns particle projected area as a function of size for different particle regimes
"""
# for spherical particles
A_s(D::FT) = FT(π) / 4 * D^2
# for nonspherical particles
A_ns(p3::PSP3, D::FT) = p3.γ * D^p3.σ
# partially rimed ice
A_r(p3::PSP3, F_r::FT, D::FT) = F_r * A_s(D) + (1 - F_r) * A_ns(p3, D)

"""
area(p3, D, F_r, th)
- p3 - a struct with P3 scheme parameters
- D - maximum particle dimension
- F_r - rime mass fraction (q_rim/q_i)
- th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d)
Returns area(D), used to create figures for the documentation.
"""
function area(
p3::PSP3,
D::FT,
F_r::FT,
th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)),
) where {FT <: Real}
# Area regime:
if P3.D_th_helper(p3) > D
return A_s(D) # small spherical ice
elseif F_r == 0
return A_ns(p3, D) # large nonspherical unrimed ice
elseif th.D_gr > D >= P3.D_th_helper(p3)
return A_ns(p3, D) # dense nonspherical ice
elseif th.D_cr > D >= th.D_gr
return A_s(D) # graupel
elseif D >= th.D_cr
return A_r(p3, F_r, D) # partially rimed ice
end

end

function define_axis(fig, row_range, col_range, title, ylabel, yticks, aspect)
return CMK.Axis(
Expand Down Expand Up @@ -143,13 +51,13 @@ function p3_relations_plot()
sol4_5 = P3.thresholds(p3, 400.0, 0.5)
sol4_8 = P3.thresholds(p3, 400.0, 0.8)
# m(D)
fig1_0 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw)
fig1_5 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
fig1_8 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
fig1_0 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw)
fig1_5 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
fig1_8 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
# a(D)
fig2_0 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw)
fig2_5 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
fig2_8 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
fig2_0 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw)
fig2_5 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
fig2_8 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
# plot verical lines
for ax in [ax1, ax2]
global d_tha = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
Expand All @@ -164,13 +72,13 @@ function p3_relations_plot()
sol_4 = P3.thresholds(p3, 400.0, 0.95)
sol_8 = P3.thresholds(p3, 800.0, 0.95)
# m(D)
fig3_200 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
fig3_400 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
fig3_800 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
fig3_200 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
fig3_400 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
fig3_800 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
# a(D)
fig3_200 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw)
fig3_400 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw)
fig3_800 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw)
fig3_200 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw)
fig3_400 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw)
fig3_800 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw)
# plot verical lines
for ax in [ax3, ax4]
global d_thb = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
Expand Down
Loading

0 comments on commit e3788cc

Please sign in to comment.