Skip to content

Commit

Permalink
Add P3 collision tendencies
Browse files Browse the repository at this point in the history
  • Loading branch information
@anastasia-popova
anastasia-popova committed May 10, 2024
1 parent c156fb3 commit cfa34d0
Show file tree
Hide file tree
Showing 6 changed files with 720 additions and 104 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.0"
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
35 changes: 35 additions & 0 deletions docs/src/P3Scheme.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,3 +188,38 @@ 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)
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
122 changes: 122 additions & 0 deletions docs/src/plots/P3TendenciesSandbox.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,122 @@
import ClimaParams
import CloudMicrophysics as CM
import CloudMicrophysics.P3Scheme as P3
import CloudMicrophysics.Parameters as CMP
import CairoMakie as Plt
import CloudMicrophysics.Microphysics1M as CM1

import Plots as PL

const PSP3 = CMP.ParametersP3

FT = Float64

p3 = CMP.ParametersP3(FT)
Chen2022 = CMP.Chen2022VelType(FT)

const liquid = CMP.CloudLiquid(FT)
const ice = CMP.CloudIce(FT)
const rain = CMP.Rain(FT)
const snow = CMP.Snow(FT)

const Blk1MVel = CMP.Blk1MVelType(FT)
const ce = CMP.CollisionEff(FT)

# Get expected N from n_0 as provided in the 1M scheme and q
function getN(q, n_0)
ρ_w = FT(1000)
λ = (8 * n_0 * π * ρ_w / q)^(1 / 4)
return n_0 / λ
end

qᵢ = FT(0.001)
qᵣ = FT(0.1)
q_c = FT(0.0005)
Nᵢ = FT(1e8)
Nᵣ = FT(1e8)
N_c = FT(1e8)
ρ_r = FT(500)
F_r = FT(0.5)
ρ_a = FT(1.2)
T = FT(300)
ρ = ρ_a

q_rain_range = range(1e-8, stop = 0.005, length = 15)
q_snow_range = q_rain_range
n_0_rain = 16 * 1e6
n_0_ice = 2 * 1e7

PL.plot(
q_rain_range * 1e3,
[
P3.ice_collisions(
"rain",
p3,
Chen2022,
qᵢ,
Nᵢ,
q,
N_c,
ρ_a,
F_r,
ρ_r,
T,
) for q in q_rain_range
],
linewidth = 3,
xlabel = "q_precipitation [g/kg]",
ylabel = "collision/accretion rate [1/s]",
label = "P3 rain and ice",
linestyle = :dash,
color = :blue,
)

PL.plot!(
q_rain_range * 1e3,
[
CM1.accretion(ice, rain, Blk1MVel.rain, ce, qᵢ, q_rai, ρ_a) for
q_rai in q_rain_range
],
linewidth = 3,
label = "1 Moment rain and ice",
color = :blue,
)

PL.plot!(
q_snow_range * 1e3,
[
P3.ice_collisions(
"cloud",
p3,
Chen2022,
q,
N_c,
q_c,
N_c,
ρ_a,
F_r,
ρ_r,
T,
FT(0.1),
) for q in q_snow_range
],
label = "P3 cloudwater and ice",
xlabel = "q_precipitation [g/kg]",
ylabel = "collision/accretion rate [1/s]",
linewidth = 3,
linestyle = :dash,
color = :red,
)

PL.plot!(
q_snow_range * 1e3,
[
CM1.accretion(liquid, snow, Blk1MVel.snow, ce, q_c, q_sno, ρ_a) for
q_sno in q_snow_range
],
linewidth = 3,
label = "1 Moment liquid and snow",
color = :red,
)

PL.savefig("CollisionComparisons.svg")
Loading

0 comments on commit cfa34d0

Please sign in to comment.