Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

changes after updating 2-moment distribution parameters #394

Merged
merged 1 commit into from
May 10, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7"
EmulatorModelsExt = ["DataFrames", "MLJ"]

[compat]
ClimaParams = "0.10.5"
ClimaParams = "0.10.6"
DataFrames = "1.6"
DocStringExtensions = "0.8, 0.9"
ForwardDiff = "0.10"
Expand Down
3 changes: 2 additions & 1 deletion docs/src/API.md
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,8 @@ Parameters.AccrTC1980
Parameters.LD2004
Parameters.VarTimescaleAcnv
Parameters.SB2006
Parameters.ParticlePDF_SB2006
Parameters.RainParticlePDF_SB2006
Parameters.CloudParticlePDF_SB2006
Parameters.AcnvSB2006
Parameters.AccrSB2006
Parameters.SelfColSB2006
Expand Down
6 changes: 3 additions & 3 deletions docs/src/Microphysics2M.md
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ where ``D`` is the drop diameter which is proportional to ``x^{1/3}``.
Through this proportionality relation, we can express the raindrops exponential distribution as a gamma distribution function
```math
\begin{align}
f_r(x)=A_rx^\nu_r e^{-B_rx^{\mu_r}},
f_r(x)=A_rx^{-\nu_r}_r e^{-B_rx^{\mu_r}},
\end{align}
```
where ``\nu_r=-\frac{2}{3}``, and ``\mu_r=\frac{1}{3}``.
Expand Down Expand Up @@ -103,7 +103,7 @@ where:
The function ``\phi_{acnv}(\tau)`` is used to correct the autoconversion rate for the undeveloped cloud droplet spectrum and the early stage rain evolution assumptions. This is a universal function which is obtained by fitting to numerical results of the SCE:
```math
\begin{equation}
\phi_{acnv}(\tau) = A \tau^b(1-\tau^b)^c,
\phi_{acnv}(\tau) = A \tau^a(1-\tau^a)^b,
\end{equation}
```
where
Expand All @@ -116,7 +116,7 @@ The default free parameter values are:
|``\nu`` | ``2`` |
|``A`` | ``400`` |
|``a`` | ``0.7`` |
|``c`` | ``3`` |
|``b`` | ``3`` |
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we still need to define c that is used in autoconversion? (line 106?)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There was an error in the documentation. Following the SB2006 paper and ClimaParams, c=4 is the accretion parameter, while in the autoconversion formula we have two parameters a and b, which are equal to 0.7 and 3. In the formula in line 106, b and c must be replaced with a and b respectively. I just change it in the last commit.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you!


The rate of change of raindrops number density is
``` math
Expand Down
6 changes: 3 additions & 3 deletions docs/src/plots/RainEvapoartionSB2006.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,11 @@ function rain_evaporation_CPU(SB2006, aps, tps, q, q_rai, ρ, N_rai, T)

(; ν_air, D_vapor) = aps
(; av, bv, α, β, ρ0) = SB2006.evap
ρw = SB2006.pdf.ρw
x_star = SB2006.pdf.xr_min
ρw = SB2006.pdf_r.ρw
x_star = SB2006.pdf_r.xr_min
G = CO.G_func(aps, tps, T, TD.Liquid())

xr = CM2.raindrops_limited_vars(SB2006.pdf, q_rai, ρ, N_rai).xr
xr = CM2.raindrops_limited_vars(SB2006.pdf_r, q_rai, ρ, N_rai).xr
Dr = (FT(6) / FT(π) / ρw)^FT(1 / 3) * xr^FT(1 / 3)

t_star = (FT(6) * x_star / xr)^FT(1 / 3)
Expand Down
97 changes: 43 additions & 54 deletions src/Microphysics2M.jl
Original file line number Diff line number Diff line change
Expand Up @@ -46,9 +46,9 @@ end
# mean terminal velocity of raindrops, and rain evaporation rates from Seifert and Beheng 2001

"""
raindrops_limited_vars(pdf, q_rai, ρ, N_rai)
raindrops_limited_vars(pdf_r, q_rai, ρ, N_rai)

- `pdf` - a struct with SB2006 size distribution parameters
- `pdf_r` - a struct with SB2006 raindrops size distribution parameters
- `q_rai` - rain water specific humidity
- `ρ` - air density
- `N_rai` raindrops number density
Expand All @@ -57,12 +57,12 @@ Returns a named tupple containing the mean mass of raindrops, xr, and the rate p
size distribution of raindrops (based on drops diameter), λr, limited within prescribed ranges
"""
function raindrops_limited_vars(
pdf::CMP.ParticlePDF_SB2006{FT},
pdf_r::CMP.RainParticlePDF_SB2006{FT},
q_rai::FT,
ρ::FT,
N_rai::FT,
) where {FT}
(; xr_min, xr_max, N0_min, N0_max, λ_min, λ_max, ρw) = pdf
(; xr_min, xr_max, N0_min, N0_max, λ_min, λ_max, ρw) = pdf_r

L_rai = ρ * q_rai
xr_0 = L_rai / N_rai
Expand Down Expand Up @@ -231,7 +231,7 @@ Returns the raindrops number density tendency due to collisions of raindrops
that produce larger raindrops (self-collection) for `scheme == SB2006Type`
"""
function rain_self_collection(
pdf::CMP.ParticlePDF_SB2006{FT},
pdf::CMP.RainParticlePDF_SB2006{FT},
self::CMP.SelfColSB2006{FT},
q_rai::FT,
ρ::FT,
Expand Down Expand Up @@ -268,7 +268,7 @@ Returns the raindrops number density tendency due to breakup of raindrops
that produce smaller raindrops for `scheme == SB2006Type`
"""
function rain_breakup(
pdf::CMP.ParticlePDF_SB2006{FT},
pdf::CMP.RainParticlePDF_SB2006{FT},
brek::CMP.BreakupSB2006{FT},
q_rai::FT,
ρ::FT,
Expand All @@ -294,8 +294,8 @@ end
"""
rain_self_collection_and_breakup(SB2006, q_rai, ρ, N_rai)

- `SB2006` - a struct with SB2006 parameters for size distribution
self collection and breakup
- `SB2006` - a struct with SB2006 parameters for raindrops size
distribution, self collection, and breakup
- `q_rai` - rain water specific humidity
- `ρ` - air density
- `N_rai` - raindrops number density
Expand All @@ -304,14 +304,14 @@ Returns a named tupple containing the raindrops self-collection and breakup rate
for `scheme == SB2006Type`
"""
function rain_self_collection_and_breakup(
(; pdf, self, brek)::CMP.SB2006{FT},
(; pdf_r, self, brek)::CMP.SB2006{FT},
q_rai::FT,
ρ::FT,
N_rai::FT,
) where {FT}

sc = rain_self_collection(pdf, self, q_rai, ρ, N_rai)
br = rain_breakup(pdf, brek, q_rai, ρ, N_rai, sc)
sc = rain_self_collection(pdf_r, self, q_rai, ρ, N_rai)
br = rain_breakup(pdf_r, brek, q_rai, ρ, N_rai, sc)

return (; sc, br)
end
Expand All @@ -332,7 +332,7 @@ Fall velocity of individual rain drops is parameterized:
- following Chen et. al 2022, DOI: 10.1016/j.atmosres.2022.106171 for `velo_scheme == Chen2022Type`
"""
function rain_terminal_velocity(
(; pdf)::CMP.SB2006{FT},
(; pdf_r)::CMP.SB2006{FT},
(; ρ0, aR, bR, cR)::CMP.SB2006VelType{FT},
q_rai::FT,
ρ::FT,
Expand All @@ -342,13 +342,13 @@ function rain_terminal_velocity(
if q_rai < eps(FT)
return (FT(0), FT(0))
end
λr = raindrops_limited_vars(pdf, q_rai, ρ, N_rai).λr
λr = raindrops_limited_vars(pdf_r, q_rai, ρ, N_rai).λr
vt0 = max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)))
vt1 = max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)^FT(4)))
return (vt0, vt1)
end
function rain_terminal_velocity(
(; pdf)::CMP.SB2006{FT},
(; pdf_r)::CMP.SB2006{FT},
vel::CMP.Chen2022VelTypeRain{FT},
q_rai::FT,
ρ::FT,
Expand All @@ -360,7 +360,7 @@ function rain_terminal_velocity(
# coefficients from Table B1 from Chen et. al. 2022
aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(vel, ρ)
# size distribution parameter
λ = raindrops_limited_vars(pdf, q_rai, ρ, N_rai).λr
λ = raindrops_limited_vars(pdf_r, q_rai, ρ, N_rai).λr

# eq 20 from Chen et al 2022
vt0 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 0))
Expand Down Expand Up @@ -401,7 +401,7 @@ specific humidity due to rain rain_evaporation, assuming a power law velocity re
fall velocity of individual drops and an exponential size distribution, for `scheme == SB2006Type`
"""
function rain_evaporation(
(; pdf, evap)::CMP.SB2006{FT},
(; pdf_r, evap)::CMP.SB2006{FT},
aps::CMP.AirProperties{FT},
tps::TDP.ThermodynamicsParameters{FT},
q::TD.PhasePartition{FT},
Expand All @@ -419,11 +419,11 @@ function rain_evaporation(

(; ν_air, D_vapor) = aps
(; av, bv, α, β, ρ0) = evap
x_star = pdf.xr_min
ρw = pdf.ρw
x_star = pdf_r.xr_min
ρw = pdf_r.ρw
G = CO.G_func(aps, tps, T, TD.Liquid())

xr = raindrops_limited_vars(pdf, q_rai, ρ, N_rai).xr
xr = raindrops_limited_vars(pdf_r, q_rai, ρ, N_rai).xr
Dr = (FT(6) / FT(π) / ρw)^FT(1 / 3) * xr^FT(1 / 3)

t_star = (FT(6) * x_star / xr)^FT(1 / 3)
Expand All @@ -447,11 +447,12 @@ function rain_evaporation(
return (; evap_rate_0, evap_rate_1)
end

"""
radar_reflectivity(struct, q_liq, q_rai, N_liq, N_rai, ρ_air, ρ_w, τ_q, τ_N)
"""
radar_reflectivity(structs, q_liq, q_rai, N_liq, N_rai, ρ_air, ρ_w, τ_q, τ_N)

- `struct` - type for 2-moment rain autoconversion parameterization
- `q_liq` - cloud water specific humidity
- `structs` - structs with SB2006 cloud droplets and raindrops
size distributions parameters
- `q_liq` - cloud water specific humidity
- `q_rai` - rain water specific humidity
- `N_liq` - cloud droplet number density
- `N_rai` - rain droplet number density
Expand All @@ -465,7 +466,7 @@ size distribuions normalized by the reflectivty of 1 millimiter drop in a volume
one meter cube
"""
function radar_reflectivity(
acnv::CMP.AcnvSB2006{FT},
(; pdf_c, pdf_r)::CMP.SB2006{FT},
q_liq::FT,
q_rai::FT,
N_liq::FT,
Expand All @@ -476,15 +477,8 @@ function radar_reflectivity(
τ_N::FT,
) where {FT}

# we assume a cloud droplets gamma distribution,
# where the parameters are νc=2 and μc=1
(; νc) = acnv

# we assume a raindrops exponential distribution
# which can be seen as a gamma distribution with
# parameters νr=-2/3 and μr=1/3
νr = FT(-2 / 3)
μr = FT(1 / 3)
(; νc, μc) = pdf_c
(; νr, μr) = pdf_r

# change of units for N_liq and N_rai
# from m^-3 to mm^-3
Expand All @@ -501,12 +495,12 @@ function radar_reflectivity(

Bc =
(xc == 0) ? FT(0) :
(SF.gamma(FT(νc + 1)) * xc / SF.gamma(FT(νc + 2)))^FT(-1)
(SF.gamma(FT(νc + 1) / μc) * xc / SF.gamma(FT(νc + 2) / μc))^(-μc)
Br =
(xr == 0) ? FT(0) :
(SF.gamma(FT(νr + 1) / μr) * xr / SF.gamma(FT(νr + 2) / μr))^(-μr)
Ac = (N_liq * Bc^(FT(νc + 1))) / SF.gamma(FT(νc + 1))
Ar = (μr * N_rai * Br^(FT(νr + 1) / μr)) / SF.gamma(FT(νr + 1) / μr)
Ac = μc * N_liq * Bc^(FT(νc + 1) / μc) / SF.gamma(FT(νc + 1) / μc)
Ar = μr * N_rai * Br^(FT(νr + 1) / μr) / SF.gamma(FT(νr + 1) / μr)

Zc = (Bc == 0) ? FT(0) : (FT(24) * Ac / (Bc^FT(5) * C^FT(2)))
Zr = (Br == 0) ? FT(0) : (FT(2160) * Ar / (Br^FT(7) * C^FT(2)))
Expand All @@ -515,11 +509,12 @@ function radar_reflectivity(
FT(10) * (log10((Zc + Zr) / Z₀) + log10(FT(1e-9)))
end

"""
effective_radius(struct, q_liq, q_rai, N_liq, N_rai, ρ_air, ρ_w, τ_q, τ_N)
"""
effective_radius(structs, q_liq, q_rai, N_liq, N_rai, ρ_air, ρ_w, τ_q, τ_N)

- `struct` - type for 2-moment rain autoconversion parameterization
- `q_liq` - cloud water specific humidity
- `structs` - structs with SB2006 cloud droplets and raindrops
size distribution parameters
- `q_liq` - cloud water specific humidity
- `q_rai` - rain water specific humidity
- `N_liq` - cloud droplet number density
- `N_rai` - rain droplet number density
Expand All @@ -532,7 +527,7 @@ Returns effective radius using the 2-moment scheme
cloud and rain particle size distributions
"""
function effective_radius(
acnv::CMP.AcnvSB2006{FT},
(; pdf_c, pdf_r)::CMP.SB2006{FT},
q_liq::FT,
q_rai::FT,
N_liq::FT,
Expand All @@ -543,15 +538,8 @@ function effective_radius(
τ_N::FT,
) where {FT}

# we assume a cloud droplets gamma distribution,
# where the parameters are νc=2 and μc=1
(; νc) = acnv

# we assume a raindrops exponential distribution
# which can be seen as a gamma distribution with
# parameters νr=-2/3 and μr=1/3
νr = FT(-2 / 3)
μr = FT(1 / 3)
(; νc, μc) = pdf_c
(; νr, μr) = pdf_r

# change of units for N_liq and N_rai
# from m^-3 to mm^-3
Expand All @@ -563,16 +551,17 @@ function effective_radius(

xc = (N_liq < τ_N) ? FT(0) : ((q_liq * ρ_air) / N_liq)
xr = (N_rai < τ_N) ? FT(0) : ((q_rai * ρ_air) / N_rai)

C = FT((4 / 3) * π * ρ_w)

Bc =
(xc == 0) ? FT(0) :
(SF.gamma(FT(νc + 1)) * xc / SF.gamma(FT(νc + 2)))^(-1)
(SF.gamma(FT(νc + 1) / μc) * xc / SF.gamma(FT(νc + 2) / μc))^(-μc)
Br =
(xr == 0) ? FT(0) :
((SF.gamma(FT(νr + 1) / μr) * xr / SF.gamma(FT(νr + 2) / μr))^(-μr))
Ac = (N_liq * Bc^(FT(νc + 1))) / SF.gamma(FT(νc + 1))
Ar = (μr * N_rai * Br^(FT(νr + 1) / μr)) / SF.gamma(FT(νr + 1) / μr)
(SF.gamma(FT(νr + 1) / μr) * xr / SF.gamma(FT(νr + 2) / μr))^(-μr)
Ac = μc * N_liq * Bc^(FT(νc + 1) / μc) / SF.gamma(FT(νc + 1) / μc)
Ar = μr * N_rai * Br^(FT(νr + 1) / μr) / SF.gamma(FT(νr + 1) / μr)

cloud_3moment = (Bc == 0) ? FT(0) : (FT(6) * Ac * C^3 / (Bc * C)^FT(4))
cloud_2moment =
Expand Down
Loading
Loading