Add terminal velocities to P3 scheme

Project.toml

@@ -7,6 +7,7 @@ version = "0.18.0"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
@@ -29,6 +30,7 @@ DocStringExtensions = "0.8, 0.9"
 EnsembleKalmanProcesses = "1.1.5"
 ForwardDiff = "0.10"
 MLJ = "0.20"
+QuadGK = "2.9.4"
 RootSolvers = "0.3, 0.4"
 SpecialFunctions = "1, 2"
 Thermodynamics = "0.12.4"

docs/Project.toml

@@ -9,6 +9,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

docs/src/API.md

@@ -70,7 +70,6 @@ MicrophysicsFlexible.weighted_vt
 ```@docs
 P3Scheme
 P3Scheme.thresholds
-P3Scheme.p3_mass
 P3Scheme.distribution_parameter_solver
 ```
 

docs/src/P3Scheme.md

@@ -159,3 +159,37 @@ Using this approach we get the following relative errors for ``\lambda``
 include("plots/P3LambdaErrorPlots.jl")
 ```
 ![](P3LambdaHeatmap.png)
+
+## Terminal Velocity 
+
+We calculated moment-weighted terminal velocity speeds using the [Chen2022](@cite) velocity parametrization. 
+The mass total-weighted fall speed (``V_m``) and the number total-weighted fall speed (``V_n``) were calculated as follows: 
+
+```math
+V_m = \frac{\int_{0}^{\infty} \! V(D) m(D) N'(D) \mathrn{d}D}{\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D}
+```
+```math 
+V_n = \frac{\int_{0}^{\infty} \! V(D) N'(D) \mathrn{d}D}{\int_{0}^{\infty} \! N'(D) \mathrm{d}D}
+```
+
+where ``V(D)`` is the velocity given by the Chen parametrization for a particle of diameter D, defined below. 
+
+```math
+V(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D}
+```
+
+and ``\phi = (16 \rho_{ice}^2 A(D)^3) / (9 \pi m(D)^2)``. 
+
+Finally, we also analyzed the mass-weighted mean particle size ``D_m`` which is given by the following: 
+
+```math 
+D_m = \frac{\int_{0}^{\infty} \! D m(D) N'(D) \mathrn{d}D}{\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D}
+```
+
+Plotting these relationships for small, medium, and large ``D_m`` gives us the following: 
+
+```@example
+include("plots/P3TerminalVelocityPlots.jl")
+```
+
+![](MorrisonandMilbrandtFig2.svg)

docs/src/plots/P3SchemePlots.jl

@@ -9,6 +9,54 @@ 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)
 
@@ -43,17 +91,13 @@ function area(
     # Area regime:
     if P3.D_th_helper(p3) > D
         return A_s(D)                      # small spherical ice
-    end
-    if F_r == 0
+    elseif F_r == 0
         return A_ns(p3, D)                 # large nonspherical unrimed ice
-    end
-    if th.D_gr > D >= P3.D_th_helper(p3)
+    elseif th.D_gr > D >= P3.D_th_helper(p3)
         return A_ns(p3, D)                 # dense nonspherical ice
-    end
-    if th.D_cr > D >= th.D_gr
+    elseif th.D_cr > D >= th.D_gr
         return A_s(D)                      # graupel
-    end
-    if D >= th.D_cr
+    elseif D >= th.D_cr
         return A_r(p3, F_r, D)             # partially rimed ice
     end
 
@@ -99,9 +143,9 @@ 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.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)
+    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)
     # 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)
@@ -120,9 +164,9 @@ 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.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)
+    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)
     # 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)

docs/src/plots/P3TerminalVelocityPlots.jl

@@ -0,0 +1,137 @@
+import ClimaParams
+import CloudMicrophysics as CM
+import CloudMicrophysics.P3Scheme as P3
+import CloudMicrophysics.Parameters as CMP
+import CairoMakie as Plt
+
+const PSP3 = CMP.ParametersP3
+
+FT = Float64
+
+p3 = CMP.ParametersP3(FT)
+
+function get_values(
+    p3::PSP3,
+    Chen2022::CMP.Chen2022VelTypeSnowIce,
+    q::FT,
+    N::FT,
+    ρ_a::FT,
+    x_resolution::Int,
+    y_resolution::Int,
+) where {FT}
+    F_rs = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution)
+    ρ_rs = range(FT(25), stop = FT(975), length = y_resolution)
+
+    V_m = zeros(x_resolution, y_resolution)
+    D_m = zeros(x_resolution, y_resolution)
+
+    for i in 1:x_resolution
+        for j in 1:y_resolution
+            F_r = F_rs[i]
+            ρ_r = ρ_rs[j]
+
+            V_m[i, j] =
+                P3.terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)[2]
+            # get D_m in mm for plots
+            D_m[i, j] = 1e3 * P3.D_m(p3, q, N, ρ_r, F_r)
+        end
+    end
+    return (; F_rs, ρ_rs, V_m, D_m)
+end
+
+function make_axis_top(fig, col, title)
+    return Plt.Axis(
+        fig[1, col],
+        xlabel = "F_r",
+        ylabel = "ρ_r",
+        title = title,
+        width = 350, height = 350,
+        limits = (0, 1.0, 25, 975),
+        xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0],
+        yticks = [200, 400, 600, 800],
+    )
+end
+function make_axis_bottom(fig, col, title)
+    return Plt.Axis(
+        fig[3, col],
+        height = 350,
+        width = 350,
+        xlabel = "F_r",
+        ylabel = "ρ_r",
+        title = title,
+    )
+end
+
+function figure_2()
+    Chen2022 = CMP.Chen2022VelType(FT)
+    # density of air in kg/m^3
+    ρ_a = FT(1.2) #FT(1.293)
+    # small D_m
+    q_s = FT(0.0008)
+    N_s = FT(1e6)
+    # medium D_m
+    q_m = FT(0.22)
+    N_m = FT(1e6)
+    # large D_m
+    q_l = FT(0.7)
+    N_l = FT(1e6)
+    # get V_m and D_m
+    xres = 100
+    yres = 100
+    (F_rs, ρ_rs, V_ms, D_ms) = get_values(p3, Chen2022.snow_ice, q_s, N_s, ρ_a, xres, yres)
+    (F_rm, ρ_rm, V_mm, D_mm) = get_values(p3, Chen2022.snow_ice, q_m, N_m, ρ_a, xres, yres)
+    (F_rl, ρ_rl, V_ml, D_ml) = get_values(p3, Chen2022.snow_ice, q_l, N_l, ρ_a, xres, yres)
+
+    fig = Plt.Figure()
+
+    # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015
+
+    args = (color = :black, labels = true, levels = 3, linewidth = 1.5, labelsize = 18)
+
+    ax1 = make_axis_top(fig, 1, "Small Dₘ")
+    hm = Plt.contourf!(ax1, F_rs, ρ_rs, V_ms, levels = 0.402:0.001:0.413)
+    Plt.contour!(ax1, F_rs, ρ_rs, D_ms; args...)
+    Plt.Colorbar(fig[2, 1], hm, vertical = false)
+
+    ax2 = make_axis_top(fig, 2, "Medium Dₘ")
+    hm = Plt.contourf!(ax2, F_rm, ρ_rm, V_mm, levels = 2:0.1:4)
+    Plt.contour!(ax2, F_rm, ρ_rm, D_mm; args...)
+    Plt.Colorbar(fig[2, 2], hm, vertical = false)
+
+    ax3 = make_axis_top(fig, 3, "Large Dₘ")
+    hm = Plt.contourf!(ax3, F_rl, ρ_rl, V_ml, levels = 2:0.2:5)
+    Plt.contour!(ax3, F_rl, ρ_rl, D_ml; args...)
+    Plt.Colorbar(fig[2, 3], hm, vertical = false)
+
+    Plt.linkxaxes!(ax1, ax2)
+    Plt.linkxaxes!(ax2, ax3)
+    Plt.linkyaxes!(ax1, ax2)
+    Plt.linkyaxes!(ax2, ax3)
+
+    ## Plot D_m as second row of comparisons
+
+    ax4 = make_axis_bottom(fig, 1, "Small Dₘ vs F_r and ρ_r")
+    hm = Plt.contourf!(ax4, F_rs, ρ_rs, D_ms)
+    Plt.Colorbar(fig[4, 1], hm, vertical = false)
+
+    ax5 = make_axis_bottom(fig, 2, "Medium Dₘ vs F_r and ρ_r")
+    hm = Plt.contourf!(ax5, F_rm, ρ_rm, D_mm)
+    Plt.Colorbar(fig[4, 2], hm, vertical = false)
+
+    ax6 = make_axis_bottom(fig, 3, "Large Dₘ vs F_r and ρ_r")
+    hm = Plt.contourf!(ax6, F_rl, ρ_rl, D_ml)
+    Plt.Colorbar(fig[4, 3], hm, vertical = false)
+
+    Plt.linkxaxes!(ax1, ax4)
+    Plt.linkxaxes!(ax4, ax5)
+    Plt.linkxaxes!(ax5, ax6)
+    Plt.linkyaxes!(ax1, ax4)
+    Plt.linkyaxes!(ax4, ax5)
+    Plt.linkyaxes!(ax5, ax6)
+
+    Plt.resize_to_layout!(fig)
+    Plt.save("MorrisonandMilbrandtFig2.svg", fig)
+end
+
+# Terminal Velocity figure
+figure_2()

src/Common.jl

@@ -18,6 +18,7 @@ export a_w_eT
 export a_w_ice
 export Chen2022_vel_add
 export Chen2022_vel_coeffs_small
+export Chen2022_vel_coeffs_large
 
 """
     G_func(air_props, tps, T, Liquid())