Skip to content

Commit

Permalink
Add Chen et al 2022 terminal velocities to 1-moment scheme
Browse files Browse the repository at this point in the history
  • Loading branch information
@apoorva-thanvantri
apoorva-thanvantri committed Aug 16, 2023
1 parent abec798 commit fcaa787
Show file tree
Hide file tree
Showing 12 changed files with 913 additions and 379 deletions.
5 changes: 5 additions & 0 deletions docs/src/API.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,8 @@ Common.logistic_function
Common.logistic_function_integral
Common.H2SO4_soln_saturation_vapor_pressure
Common.Delta_a_w
Common.Chen2022_vel_add
Common.Chen2022_vel_coeffs
```

# Common utility types
Expand All @@ -115,6 +117,9 @@ CommonTypes.B1994Type
CommonTypes.TC1980Type
CommonTypes.LD2004Type
CommonTypes.SB2006Type
CommonTypes.AbstractTerminalVelocityType
CommonTypes.Blk1MVelType
CommonTypes.SB2006VelType
CommonTypes.Chen2022Type
CommonTypes.AbstractAerosolType
CommonTypes.ArizonaTestDustType
Expand Down
285 changes: 263 additions & 22 deletions docs/src/Microphysics1M.md
View file

Large diffs are not rendered by default.

4 changes: 1 addition & 3 deletions docs/src/Microphysics2M.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -403,12 +403,10 @@ Below, we compare the individual terminal velocity formulas for [Chen2022](@cite
We also compare bulk number weighted [ND] and mass weighted [MD]
terminal velocities for both formulas integrated over the size distribution from [SeifertBeheng2006](@cite).
We also show the mass weighted terminal velocity from the 1-moment scheme.

```@example
include("TerminalVelocityComparisons.jl")
```
![](terminal_velocity_bulk_comparisons.svg)
![](terminal_velocity_individual_raindrop.svg)
![](2M_terminal_velocity_comparisons.svg)

### Accretion and Autoconversion

Expand Down
318 changes: 214 additions & 104 deletions docs/src/TerminalVelocityComparisons.jl
View file
Original file line number Diff line number Diff line change
@@ -1,37 +1,42 @@
import Plots as PL

FT = Float64

import CloudMicrophysics as CM
import CLIMAParameters as CP

FT = Float64

const CMT = CM.CommonTypes
const CM1 = CM.Microphysics1M
const CM2 = CM.Microphysics2M
const CMP = CM.Parameters

include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
const param_set = cloud_microphysics_parameters(toml_dict)
const CMO = CM.Common

const rain = CMT.RainType()
const liquid = CMT.LiquidType()
const ice = CMT.IceType()
const snow = CMT.SnowType()
const SB2006 = CMT.SB2006Type()
const Chen2022 = CMT.Chen2022Type()
const SB2006Vel = CMT.SB2006VelType()
const Blk1MVel = CMT.Blk1MVelType()
const APS = CMP.AbstractCloudMicrophysicsParameters

include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
const param_set = cloud_microphysics_parameters(toml_dict)

"""
rain_terminal_velocity_individual_C(ρ, scheme, D_r)
rain_terminal_velocity_individual_Chen(param_set, ρ, D_r)
- `ρ`` - air density
- `scheme` - type for parameterization
- `param_set` - set with free parameters
- `ρ` - air density
- `D_r` - diameter of the raindrops
Returns the fall velocity of a raindrop using multiple gamma-function terms
to increase accuracy for `scheme == Chen2022Type`
Returns the fall velocity of a raindrop from Chen et al 2022
"""
function rain_terminal_velocity_individual_Chen(
param_set,
ρ::FT,
scheme::CMT.Chen2022Type,
D_r::FT, #in m
) where {FT <: Real}

Expand Down Expand Up @@ -63,6 +68,15 @@ function rain_terminal_velocity_individual_Chen(
return v
end

"""
rain_terminal_velocity_individual_SB(param_set, ρ, D_r)
- `param_set` - set with free parameters
- `ρ` - air density
- `D_r` - diameter of the raindrops
Returns the fall velocity of a raindrop from Seifert and Beheng 2006
"""
function rain_terminal_velocity_individual_SB(
param_set,
ρ::FT,
Expand All @@ -78,95 +92,191 @@ function rain_terminal_velocity_individual_SB(
return v
end

q_liq_range = range(1e-8, stop = 1e-3, length = 1000)
q_rai_range = range(1e-8, stop = 1e-3, length = 1000)
D_r_range = range(50e-6, stop = 9e-3, length = 1000)
N_d_range = range(1e7, stop = 1e9, length = 1000)
q_liq = 5e-4
# q_rai = 5e-4
ρ_air = 1.0 # kg m^-3
N_rai = 1e4 #this value matters for the individual terminal velocity

PL.plot(
q_rai_range * 1e3,
[
CM2.rain_terminal_velocity(param_set, SB2006, q_rai, ρ_air, N_rai)[1]
for q_rai in q_rai_range
],
linewidth = 3,
xlabel = "q_rain [g/kg]",
ylabel = "terminal velocity [m/s]",
label = "Rain-SB2006 [ND]",
)
PL.plot!(
q_rai_range * 1e3,
[
CM2.rain_terminal_velocity(param_set, Chen2022, q_rai, ρ_air, N_rai)[1]
for q_rai in q_rai_range
],
linewidth = 3,
xlabel = "q_rain [g/kg]",
ylabel = "terminal velocity [m/s]",
label = "Rain-Chen2022 [ND]",
)
PL.plot!(
q_rai_range * 1e3,
[
CM2.rain_terminal_velocity(param_set, SB2006, q_rai, ρ_air, N_rai)[2]
for q_rai in q_rai_range
],
linewidth = 3,
xlabel = "q_rain [g/kg]",
ylabel = "terminal velocity [m/s]",
label = "Rain-SB2006 [M]",
)
PL.plot!(
q_rai_range * 1e3,
[
CM2.rain_terminal_velocity(param_set, Chen2022, q_rai, ρ_air, N_rai)[2]
for q_rai in q_rai_range
],
linewidth = 3,
xlabel = "q_rain [g/kg]",
ylabel = "terminal velocity [m/s]",
label = "Rain-Chen2022 [M]",
)
PL.plot!(
q_rai_range * 1e3,
[
CM1.terminal_velocity(param_set, rain, ρ_air, q_rai) for
q_rai in q_rai_range
],
linewidth = 3,
xlabel = "q_rain [g/kg]",
ylabel = "terminal velocity [m/s]",
label = "Rain-Ogura-1-moment",
)

PL.title!("Terminal Velocity vs. Specific Humidity of Rain", titlefontsize = 9)
PL.savefig("terminal_velocity_bulk_comparisons.svg")

PL.plot(
D_r_range,
[
rain_terminal_velocity_individual_SB(param_set, ρ_air, D_r) for
D_r in D_r_range
],
linewidth = 3,
xlabel = "D [m]",
ylabel = "terminal velocity [m/s]",
label = "Rain-SB2006",
)
PL.plot!(
D_r_range,
[
rain_terminal_velocity_individual_Chen(param_set, ρ_air, Chen2022, D_r)
for D_r in D_r_range
],
linewidth = 3,
xlabel = "D [m]",
ylabel = "terminal velocity [m/s]",
label = "Rain-Chen2022",
)
PL.title!("Terminal Velocity vs. Diameter", titlefontsize = 13)
PL.savefig("terminal_velocity_individual_raindrop.svg")
"""
terminal_velocity_individual_1M(prs, precip, ρ, D_r)
- `prs` - set with free parameters
- `precip` - precipitation type (rain or snow)
- `ρ` - air density
- `D_r` - particle diameter
Returns the fall velocity of a raindrop or snow from 1-moment scheme
"""
function terminal_velocity_individual_1M(
prs::APS,
precip::CMT.AbstractPrecipType,
ρ::FT,
D_r::FT,
) where {FT <: Real}
_χv::FT = CM1.χv(prs, precip)
_v0::FT = CM1.v0(prs, ρ, precip)
_ve::FT = CM1.ve(prs, precip)
_r0::FT = CM1.r0(prs, precip)
_Δv::FT = CM1.Δv(prs, precip)
vt = _χv * _v0 * (D_r / (2 * _r0))^(_Δv + _ve)
return vt
end

"""
ice_terminal_velocity_individual_Chen(param_set, ρ, D_r)
- `param_set` - set with free parameters
- `ρ` - air density
- `D_r` - diameter of the raindrops
Returns the fall velocity of an ice particle from Chen et al 2022
"""
function ice_terminal_velocity_individual_Chen(
prs::APS,
ρ::FT,
D_r::FT, #in m
) where {FT <: Real}

ρ_i::FT = CMP.ρ_cloud_ice(prs)
As = CMO.As(prs, ρ_i)
Bs = CMO.Bs(prs, ρ_i)
Cs = CMO.Cs(prs, ρ_i)
Es = CMO.Es(prs, ρ_i)
Fs = CMO.Fs(prs, ρ_i)
Gs = CMO.Gs(prs, ρ_i)
ai = [Es * ρ^(As), Fs * ρ^(As)]
bi = [Bs + ρ * (Cs), Bs + ρ * (Cs)]
ci = [0, Gs]

D_r = D_r * 1000 #D_r is in mm in the paper --> multiply D_r by 1000
v = 0
for i in 1:2
v += (ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i]))
end
return v
end

"""
snow_terminal_velocity_individual_Chen(prs, ρ, D_r)
- `prs` - set with free parameters
- `ρ` - air density
- `D_r` - particle diameter
Returns the fall velocity of snow from Chen et al 2022
"""
function snow_terminal_velocity_individual_Chen(
prs::APS,
ρ::FT,
D_r::FT,
) where {FT <: Real}
D_r = D_r * 1000
_m0::FT = CM1.m0(prs, snow)
_me::FT = CM1.me(prs, snow)
_Δm::FT = CM1.Δm(prs, snow)
_a0::FT = CM1.a0(prs, snow)
_ae::FT = CM1.ae(prs, snow)
_Δa::FT = CM1.Δa(prs, snow)
_χm::FT = CM1.χm(prs, snow)
_χa::FT = CM1.χa(prs, snow)
_r0::FT = CM1.r0(prs, snow)
ρ_i::FT = CMP.ρ_cloud_ice(prs)
m0_comb::FT = _m0 * _χm
a0_comb::FT = _a0 * _χa
me_comb::FT = _me + _Δm
ae_comb::FT = _ae + _Δa
α = FT(-1 / 3)

As = CMO.As(prs, ρ_i)
Bs = CMO.Bs(prs, ρ_i)
Cs = CMO.Cs(prs, ρ_i)
Es = CMO.Es(prs, ρ_i)
Fs = CMO.Fs(prs, ρ_i)
Gs = CMO.Gs(prs, ρ_i)
ai = [Es * ρ^(As), Fs * ρ^(As)]
bi = [Bs + ρ * (Cs), Bs + ρ * (Cs)]
ci = [0, Gs]

aspect_ratio =
((16 * (a0_comb)^3 * (ρ_i)^2) / (9 * π * (m0_comb)^2)) *
(D_r / (2 * _r0 * 1000))^(3 * ae_comb - 2 * me_comb)
aspect_ratio = aspect_ratio^(α)
vt = 0
for i in 1:2
vt = vt + aspect_ratio * ai[i] * ((D_r))^(bi[i]) * exp(-1 * ci[i] * D_r)
end
return vt
end

function aspect_ratio_snow_1M(prs::APS, D_r::FT) where {FT <: Real}
_m0::FT = CM1.m0(prs, snow)
_me::FT = CM1.me(prs, snow)
_Δm::FT = CM1.Δm(prs, snow)
_a0::FT = CM1.a0(prs, snow)
_ae::FT = CM1.ae(prs, snow)
_Δa::FT = CM1.Δa(prs, snow)
_χm::FT = CM1.χm(prs, snow)
_χa::FT = CM1.χa(prs, snow)
_r0::FT = CM1.r0(prs, snow)
ρ_i::FT = CMP.ρ_cloud_ice(prs)

m0_comb::FT = _m0 * _χm
a0_comb::FT = _a0 * _χa
me_comb::FT = _me + _Δm
ae_comb::FT = _ae + _Δa

aspect_ratio =
((16 * (a0_comb)^3 * (ρ_i)^2) / (9 * π * (m0_comb)^2)) *
(D_r / (2 * _r0))^(3 * ae_comb - 2 * me_comb)
return aspect_ratio
end

ρ_air, ρ_air_ground = 1.2, 1.22
q_liq, q_ice, q_tot = 5e-4, 5e-4, 20e-3
N_rai = 1e7 #this value matters for the individual terminal velocity
q_rain_range = range(1e-8, stop = 5e-3, length = 100)
D_r_range = range(1e-6, stop = 5e-3, length = 1000)
N_d_range = range(1e6, stop = 1e9, length = 1000)
D_r_ar_range = range(1e-6, stop = 6.25e-4, length = 1000)

#! format: off

SB_rain_bN = [CM2.rain_terminal_velocity(param_set, SB2006, SB2006Vel, q_rai, ρ_air, N_rai)[1] for q_rai in q_rain_range]
Ch_rain_bN = [CM2.rain_terminal_velocity(param_set, SB2006, Chen2022, q_rai, ρ_air, N_rai)[1] for q_rai in q_rain_range]
SB_rain_bM = [CM2.rain_terminal_velocity(param_set, SB2006, SB2006Vel, q_rai, ρ_air, N_rai)[2] for q_rai in q_rain_range]
Ch_rain_bM = [CM2.rain_terminal_velocity(param_set, SB2006, Chen2022, q_rai, ρ_air, N_rai)[2] for q_rai in q_rain_range]
M1_rain_bM = [CM1.terminal_velocity(param_set, rain, Blk1MVel, ρ_air, q_rai) for q_rai in q_rain_range]

SB_rain = [rain_terminal_velocity_individual_SB(param_set, ρ_air, D_r) for D_r in D_r_range]
Ch_rain = [rain_terminal_velocity_individual_Chen(param_set, ρ_air, D_r) for D_r in D_r_range]
M1_rain = [terminal_velocity_individual_1M(param_set, rain, ρ_air, D_r) for D_r in D_r_range]
Ch_snow = [snow_terminal_velocity_individual_Chen(param_set, ρ_air, D_r) for D_r in D_r_range]
M1_snow = [terminal_velocity_individual_1M(param_set, snow, ρ_air, D_r) for D_r in D_r_range]
Ch_ice = [ice_terminal_velocity_individual_Chen(param_set, ρ_air, D_r) for D_r in D_r_range]

Aspect_Ratio = [aspect_ratio_snow_1M(param_set, D_r) for D_r in D_r_ar_range]

# 2 Moment Scheme Figures

# single drop comparison
p1 = PL.plot(D_r_range *1e3, SB_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006", color = :blue)
p1 = PL.plot!(D_r_range * 1e3, Ch_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen2022", color = :red)
p1 = PL.plot!(D_r_range * 1e3, M1_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-1M", color = :green)

# group velocity comparison
p2 = PL.plot(q_rain_range * 1e3, SB_rain_bN, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [ND]", linestyle = :dot, color = :blue)
p2 = PL.plot!(q_rain_range * 1e3, Ch_rain_bN, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen2022 [ND]", linestyle = :dot, color = :red)
p2 = PL.plot!(q_rain_range * 1e3, SB_rain_bM, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [M]", color = :blue)
p2 = PL.plot!(q_rain_range * 1e3, Ch_rain_bM, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen2022 [M]", color = :red)
p2 = PL.plot!(q_rain_range * 1e3, M1_rain_bM, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-default-1M [M]", color = :green)
PL.plot(p1, p2, layout = (1, 2), size = (800, 500), dpi = 300)
PL.savefig("2M_terminal_velocity_comparisons.svg")

# 1 Moment Scheme Figures

# individual particle comparison
PL.plot(D_r_range * 1e3, Ch_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen", color = :red)
PL.plot!(D_r_range * 1e3, M1_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-default-1M", color = :green)
PL.plot!(D_r_range * 1e3, Ch_snow, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Snow-Chen", color = :red, linestyle =:dot)
PL.plot!(D_r_range * 1e3, M1_snow, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Snow-default-1M", color = :green, linestyle = :dot)
PL.plot!(D_r_range * 1e3, Ch_ice, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Ice-Chen", color = :blue, linestyle =:dot)
PL.savefig("1M_individual_terminal_velocity_comparisons.svg")

PL.plot(D_r_ar_range * 1e3, Aspect_Ratio, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Aspect Ratio", color = :red)
PL.savefig("aspect_ratio_1M_snow.svg")
#! format: on
Loading

0 comments on commit fcaa787

Please sign in to comment.