diff --git a/Project.toml b/Project.toml index 5064e7b155..02f1811ef7 100644 --- a/Project.toml +++ b/Project.toml @@ -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" @@ -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" diff --git a/docs/src/API.md b/docs/src/API.md index 8b25a07112..d04352f6e1 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -65,6 +65,12 @@ P3Scheme.thresholds P3Scheme.distribution_parameter_solver ``` +# Terminal Velocity +```@docs +TerminalVelocity +TerminalVelocity.velocity_chen +``` + # Aerosol model ```@docs diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index a197f080fb..e1d7ae3c04 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -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) \ No newline at end of file diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl index c106b8f991..b3bc2ecee4 100644 --- a/docs/src/plots/P3SchemePlots.jl +++ b/docs/src/plots/P3SchemePlots.jl @@ -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( @@ -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) @@ -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) diff --git a/docs/src/plots/P3TendenciesSandbox.jl b/docs/src/plots/P3TendenciesSandbox.jl new file mode 100644 index 0000000000..23f6d4bf0c --- /dev/null +++ b/docs/src/plots/P3TendenciesSandbox.jl @@ -0,0 +1,291 @@ +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 Thermodynamics as TD + +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.005) +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 = 100) +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ᵢ, + 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") + +const tps = TD.Parameters.ThermodynamicsParameters(FT) +const aps = CMP.AirProperties(FT) + +# 1 Moment Parameters to match graph +T, p = 273.15 + 15, 90000.0 +ϵ = 1.0 / TD.Parameters.molmass_ratio(tps) +p_sat = TD.saturation_vapor_pressure(tps, T, TD.Ice()) +q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0)) +q_tot = 15e-3 +q_vap = 0.15 * q_sat +q_liq = 0.0 +q_ice = q_tot - q_vap - q_liq +q = TD.PhasePartition(q_tot, q_liq, q_ice) +R = TD.gas_constant_air(tps, q) +ρ = p / R / T +T = 273.15 + +PL.plot( + q_snow_range * 1e3, + [ + CM1.snow_melt(snow, Blk1MVel.snow, aps, tps, q_sno, ρ, T + 2) for + q_sno in q_snow_range + ], + xlabel = "q_snow [g/kg]", + linewidth = 3, + ylabel = "snow melt rate [1/s]", + label = "1M: T=2C", + color = :blue, +) +PL.plot!( + q_snow_range * 1e3, + [ + CM1.snow_melt(snow, Blk1MVel.snow, aps, tps, q_sno, ρ, T + 4) for + q_sno in q_snow_range + ], + xlabel = "q_snow [g/kg]", + linewidth = 3, + label = "1M: T=4C", + color = :red, +) +PL.plot!( + q_snow_range * 1e3, + [ + CM1.snow_melt(snow, Blk1MVel.snow, aps, tps, q_sno, ρ, T + 6) for + q_sno in q_snow_range + ], + xlabel = "q_snow [g/kg]", + linewidth = 3, + label = "1M: T=6C", + color = :green, +) + +PL.plot!( + q_snow_range * 1e3, + [ + P3.p3_melt(p3, Chen2022, aps, tps, q_sno, Nᵢ, T + 2, ρ_a, F_r, ρ_r) for + q_sno in q_snow_range + ], + xlabel = "q_snow [g/kg]", + linewidth = 3, + label = "P3: T=2C", + linestyle = :dash, + color = :blue, +) + +PL.plot!( + q_snow_range * 1e3, + [ + P3.p3_melt(p3, Chen2022, aps, tps, q_sno, Nᵢ, T + 4, ρ_a, F_r, ρ_r) for + q_sno in q_snow_range + ], + xlabel = "q_snow [g/kg]", + linewidth = 3, + label = "P3: T=4C", + linestyle = :dash, + color = :red, +) + +PL.plot!( + q_snow_range * 1e3, + [ + P3.p3_melt(p3, Chen2022, aps, tps, q_sno, Nᵢ, T + 6, ρ_a, F_r, ρ_r) for + q_sno in q_snow_range + ], + xlabel = "q_snow [g/kg]", + linewidth = 3, + label = "P3: T=6C", + linestyle = :dash, + color = :green, +) + +PL.savefig("MeltRateComparisons.svg") + +T = FT(250) +qᵥ = FT(8.1e-4) +aero_type = CMP.Illite(FT) +mass = true + +PL.plot( + q_rain_range * 1e3, + [ + P3.p3_rain_het_freezing(mass, tps, q, N_c, T - 1, ρ_a, qᵥ, aero_type) + for q in q_rain_range + ], + xlabel = "q_rain [g/ kg]", + linewidth = 3, + label = "249K", + title = "Δ Mass / Second", + yscale = :log10, +) + +PL.plot!( + q_rain_range * 1e3, + [ + P3.p3_rain_het_freezing(mass, tps, q, N_c, T, ρ_a, qᵥ, aero_type) for + q in q_rain_range + ], + linewidth = 3, + label = "250K", +) + +PL.plot!( + q_rain_range * 1e3, + [ + P3.p3_rain_het_freezing(mass, tps, q, N_c, T + 1, ρ_a, qᵥ, aero_type) + for q in q_rain_range + ], + linewidth = 3, + label = "251K", +) + +PL.savefig("MassFreezeRateComparisons.svg") + +PL.plot( + q_rain_range * 1e3, + [ + P3.p3_rain_het_freezing(!mass, tps, q, N_c, T - 1, ρ_a, qᵥ, aero_type) + for q in q_rain_range + ], + xlabel = "q_rain [g/ kg]", + linewidth = 3, + label = "249K", + title = "Δ Number Concentration / Second", + yscale = :log10, +) + +PL.plot!( + q_rain_range * 1e3, + [ + P3.p3_rain_het_freezing(!mass, tps, q, N_c, T, ρ_a, qᵥ, aero_type) for + q in q_rain_range + ], + linewidth = 3, + label = "250K", +) + +PL.plot!( + q_rain_range * 1e3, + [ + P3.p3_rain_het_freezing(!mass, tps, q, N_c, T + 1, ρ_a, qᵥ, aero_type) + for q in q_rain_range + ], + linewidth = 3, + label = "251K", +) + +PL.savefig("FreezeRateComparisons.svg") diff --git a/src/CloudMicrophysics.jl b/src/CloudMicrophysics.jl index 16987404bc..527c59e852 100644 --- a/src/CloudMicrophysics.jl +++ b/src/CloudMicrophysics.jl @@ -14,6 +14,7 @@ include("Common.jl") include("Microphysics0M.jl") include("Microphysics1M.jl") include("Microphysics2M.jl") +include("TerminalVelocity.jl") include("P3Scheme.jl") include("MicrophysicsNonEq.jl") include("AerosolModel.jl") diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl index 6992b0b82f..192223ae0c 100644 --- a/src/P3Scheme.jl +++ b/src/P3Scheme.jl @@ -14,10 +14,16 @@ module P3Scheme import SpecialFunctions as SF import QuadGK as QGK import RootSolvers as RS +import HCubature as HC + +import Thermodynamics as TD +import Thermodynamics.Parameters as TDP import ClimaParams as CP import CloudMicrophysics.Parameters as CMP import CloudMicrophysics.Common as CO +import CloudMicrophysics as CM +import CloudMicrophysics.TerminalVelocity as TV const PSP3 = CMP.ParametersP3 @@ -156,6 +162,102 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT} end end +""" + 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} = α_va_si(p3) * D^p3.β_va +# for partially rimed ice +mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} = + α_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 = 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) where {FT <: Real} = FT(π) / 4 * D^2 +# for nonspherical particles +A_ns(p3::PSP3, D::FT) where {FT <: Real} = p3.γ * D^p3.σ +# partially rimed ice +A_r(p3::PSP3, F_r::FT, D::FT) where {FT <: Real} = + F_r * A_s(D) + (1 - F_r) * A_ns(p3, D) + +""" + p3_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 p3_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 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 >= 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 + else + throw("D not in range") + end +end + # Some wrappers to cast types from SF.gamma # (which returns Float64 even when the input is Float32) Γ(a::FT, z::FT) where {FT <: Real} = FT(SF.gamma(a, z)) @@ -475,7 +577,7 @@ function terminal_velocity( # Get the thresholds for different particles regimes (; D_cr, D_gr, ρ_g, ρ_d) = thresholds(p3, ρ_r, F_r) D_th = D_th_helper(p3) - D_ct = FT(0.000625) # TODO add to ClimaParams + D_ct = Chen2022.cutoff # TODO add to ClimaParams # Get the shape parameters of the particle size distribution (λ, N_0) = distribution_parameter_solver(p3, q, N, ρ_r, F_r) @@ -667,4 +769,488 @@ function D_m(p3::PSP3, q::FT, N::FT, ρ_r::FT, F_r::FT) where {FT} return n / q end +""" + get_rain_parameters(q, N) + + - q - mass mixing ratio of rain + - N - number mixing ratio of rain + + returns the parameters λ and N_0 where N' = N_0 * exp(-λ * D) +""" +function get_rain_parameters(q::FT, N::FT) where {FT} #MOVE TO 2M + ρ_l = FT(1000) + λ = (π * ρ_l * N / q)^(1 / 3) + N_0 = N * λ + return (λ, N_0) +end + +""" + get_cloud_parameters(q, N) + + - q - mass mixing ratio of rain + - N - number mixing ratio of rain + + returns the parameters λ and N_0 where N' = N_0 * D ^ 8 * exp(-λ ^3 * D ^ 3) +""" +function get_cloud_parameters(q::FT, N::FT) where {FT} #MOVE TO 2M + ρ_l = FT(1000) + ρ_a = FT(1.2) + λ = (π * ρ_l * N / (2 * q * ρ_a))^(1 / 3) + N_0 = 3 / 2 * N * λ^9 + return (λ, N_0) +end + +""" + N′_rain(D, q, N) + + - D - diameter of particle + - q - mass mixing ratio + - N - number mixing ratio + + Returns the distribution of rain particles (assumed to be of the form + N'(D) = N0 * exp(-λD)) at given D +""" +function N′rain(D::FT, q::FT, N::FT) where {FT} #MOVE TO 2M + (λ, N_0) = get_rain_parameters(q, N) + return N_0 * exp(-λ * D) +end + +""" + N′_cloud(D, q, N) + + - D - diameter of particle + - q - mass mixing ratio + - N - number mixing ratio + + Returns the distribution of cloud particles (assumed to be of the form + N'(D) = N0 * D ^ 8 * exp(-λ D ^ 3)) at given D +""" +function N′cloud(D::FT, q::FT, N::FT) where {FT} #MOVE TO 2M + (λ, N_0) = get_cloud_parameters(q, N) + return N_0 * D^8 * exp(-λ^3 * D^3) +end + +""" + N′_ice(D, p3, λ, N0) + + - D - diameter of particle + - p3 - a struct containing P3 scheme parameters + - λ - shape parameter of distribution + - N0 - shape parameter of distribution + + Returns the distribution of ice particles (assumed to be of the form + N'(D) = N0 * D ^ μ * exp(-λD)) at given D +""" +function N′ice(D::FT, p3::PSP3, λ::FT, N_0::FT) where {FT} + return N_0 * D^DSD_μ(p3, λ) * exp(-λ * D) +end + +""" + get_ice_bounds(p3, λ, tolerance) + + - p3 - a struct containing p3 parameters + - λ - shape parameters of ice distribution + - tolerance - tolerance for how much of the distribution we want to integrate over + + Returns the bound on the distribution that would guarantee that 1-tolerance + of the ice distribution is integrated over. This is calculated by setting + N_0(1 - tolerance) = ∫ N'(D) dD from 0 to bound and solving for bound. + This was further simplified to cancel out the N_0 from both sides. + The guess was calculated through a linear approximation extrapolated from + numerical solutions. +""" +function get_ice_bound(p3, λ::FT, tolerance::FT) where {FT} + ice_problem(x) = + tolerance - Γ(1 + DSD_μ(p3, λ), FT(exp(x)) * λ) / Γ(1 + DSD_μ(p3, λ)) + guess = log(19 / 6 * (DSD_μ(p3, λ) - 1) + 39) - log(λ) + log_ice_x = + RS.find_zero( + ice_problem, + RS.SecantMethod(guess - 1, guess), + RS.CompactSolution(), + RS.RelativeSolutionTolerance(eps(FT)), + 5, + ).root + + return exp(log_ice_x) +end + +""" + get_rain_bound(q, N, tolerance) + + - q - mass mixing ratio of rain + - N - number mixing ratio of rain + - tolerance - tolerance for integration error + + Returns the value such that N_0(1 - tolerance) = ∫ N'(D) dD from 0 to bound + solved analytically. +""" +function get_rain_bound(q::FT, N::FT, tolerance::FT) where {FT} + rain_λ, = get_rain_parameters(q, N) + return -1 / rain_λ * log(tolerance) +end + +""" + get_cloud_bound(q, N, tolerance) + +- q - mass mixing ratio of cloud water +- N - number mixing ratio of cloud_water +- tolerance - tolerance for integration error + +Returns the value such that 1- tolerance of the whole cloud distribution +is integrated over i.e, N_0(1 - tolerance) = ∫ N'(D) dD from 0 to bound +further simplified by cancelling out N_0 from both sides. +The guess was calculated through a linear approximation of the bounds from +numerical solutions. +""" +function get_cloud_bound(q::FT, N::FT, tolerance::FT) where {FT} + cloud_λ, = get_cloud_parameters(q, N) + cloud_problem(x) = + tolerance - + exp(-exp(x)^3 * cloud_λ^3) * + (1 + exp(x)^3 * cloud_λ^3 + 1 / 2 * exp(x)^6 * cloud_λ^6) + guess = + log(0.5) + + (log(0.00025) - log(0.5)) / (log(1e12) - log(1e2)) * + (log(cloud_λ^3) - log(10^2)) + log_cloud_x = + RS.find_zero( + cloud_problem, + RS.NewtonsMethodAD(guess), + RS.CompactSolution(), + RS.RelativeSolutionTolerance(eps(FT)), + 5, + ).root + return exp(log_cloud_x) +end + +""" + integration_bounds(p3, tolerance, λ, N_0, Nᵢ, qᵣ, Nᵣ) + + - p3 - a struct containing P3 Scheme parameters + - tolerance - tolerance to which distributions need to be evaluated to + - λ - shape parameter of ice distribution + - N_0 - intercept size distribution of ice + - Nᵢ - number mixing ratio of ice + - q - mass mizing ratio of colliding species + - N - number mixing ratio of colliding species + + Returns the bounds over which to integrate rain, cloud, and ice distributions to ensure + coverage or more than (1 - tolerance) of each distribution +""" +function integration_bounds( + type::String, + p3::PSP3, + tolerance::FT, + λ::FT, + N_0::FT, + Nᵢ::FT, + q::FT, + N::FT, +) where {FT} + if type == "rain" + colliding_x = get_rain_bound(q, N, tolerance) + elseif type == "cloud" + colliding_x = get_cloud_bound(q, N, tolerance) + end + + ice_bound = get_ice_bound(p3, λ, tolerance) + return (2 * colliding_x, 2 * ice_bound) +end + +""" + a(D1, D2) + + - D1 - maximum dimension of first particle + - D2 - maximum dimension of second particle + + Returns the collision kernel (assumed to be of the form π(r1 + r2)^2) for the + two colliding particles +""" +function a(D1::FT, D2::FT) where {FT} + # TODO make this more accurate for non-spherical particles + return π * (D1 / 2 + D2 / 2)^2 +end + +""" + N′colliding(type, D, q, N) + + - type - defines what is colliding with ice ("rain" or "cloud") + - D - maximum dimension + - q - mass mixing ratio + - N - number mixing ratio + + Returns the corresponding distribution of particles depending on type +""" +function N′colliding(type::String, D::FT, q::FT, N::FT) where {FT} + if type == "rain" + return N′rain(D, q, N) + elseif type == "cloud" + return N′cloud(D, q, N) + else + throw("Undefined Type of Collisions") + end +end + +""" + vel_diff(type, D, Dᵢ, p3, Chen2022, ρ_a, F_r, th) + + - type - defines what is colliding with ice ("rain" or "cloud") + - D - maximum dimension of colliding particle + - Dᵢ - maximum dimension of ice particle + - p3 - a struct containing P3 parameters + - Chen2022 - a struct containing Chen 2022 velocity parameters + - ρ_a - density of air + - F_r - rime mass fraction (q_rim/ q_i) + - th - thresholds as calculated by thresholds() + + Returns the corresponding velocity difference of colliding particles depending on type +""" +function vel_diff( + type::String, + D::FT, + Dᵢ::FT, + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + ρ_a::FT, + F_r::FT, + th, +) where {FT} + if type == "rain" + return abs( + TV.velocity_chen( + Dᵢ, + Chen2022.snow_ice, + ρ_a, + p3_mass(p3, D, F_r, th), + p3_area(p3, D, F_r, th), + p3.ρ_i, + ) - TV.velocity_chen(D, Chen2022.rain, ρ_a), + ) + elseif type == "cloud" + return abs( + TV.velocity_chen( + Dᵢ, + Chen2022.snow_ice, + ρ_a, + p3_mass(p3, D, F_r, th), + p3_area(p3, D, F_r, th), + p3.ρ_i, + ), + ) + end +end + +""" + ice_collisions(type, p3, Chen2022, ρ_a, F_r, qᵣ, qᵢ, Nᵣ, Nᵢ, ρ, ρ_r, E_ri) + + - type - defines what is colliding with the ice ("cloud" or "rain") + - p3 - a struct with P3 scheme parameters + - Chen2022 - a struct with terminal velocity parameters as in Chen (2022) + - qᵢ - mass mixing ratio of ice + - Nᵢ - number mixing ratio of ice + - q_c - mass mixing ratio of colliding species + - N_c - number mixing ratio of colliding species + - ρ_a - density of air + - F_r - rime mass fraction (q_rim/ q_i) + - ρ_r - rime density (q_rim/B_rim) + - T - temperature (in K) + - E_ci - collision efficiency between ice and colliding species + + Returns the rate of collisions between cloud ice and rain + Equivalent to the measure of QRCOL in Morrison and Mildbrandt (2015) +""" +function ice_collisions( + type::String, + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + qᵢ::FT, + Nᵢ::FT, + q_c::FT, + N_c::FT, + ρ_a::FT, + F_r::FT, + ρ_r::FT, + T::FT, + E_ci = FT(1), +) where {FT} + ρ_l = p3.ρ_l + th = thresholds(p3, ρ_r, F_r) + (λ, N_0) = distribution_parameter_solver(p3, qᵢ, Nᵢ, ρ_r, F_r) + (colliding_bound, ice_bound) = + integration_bounds(type, p3, eps(FT), λ, N_0, Nᵢ, q_c, N_c) + + if T > 273.15 + f_warm(D_c, Dᵢ) = + E_ci / ρ_a * + N′colliding(type, D_c, q_c, N_c) * + N′ice(Dᵢ, p3, λ, N_0) * + a(D_c, Dᵢ) * + p3_mass(p3, Dᵢ, F_r, th) * + vel_diff(type, D_c, Dᵢ, p3, Chen2022, ρ_a, F_r, th) + (dqdt, error) = HC.hcubature( + d -> f_warm(d[1], d[2]), + (FT(0), FT(0)), + (colliding_bound, ice_bound), + ) + else + f_cold(D_c, Dᵢ) = + E_ci / ρ_a * + N′colliding(type, D_c, q_c, N_c) * + N′ice(Dᵢ, p3, λ, N_0) * + a(D_c, Dᵢ) * + mass_s(D_c, ρ_l) * + vel_diff(type, D_c, Dᵢ, p3, Chen2022, ρ_a, F_r, th) + (dqdt, error) = HC.hcubature( + d -> f_cold(d[1], d[2]), + (FT(0), FT(0)), + (colliding_bound, ice_bound), + ) + end + + return dqdt +end + +""" + dmdt_mass(p3, D, F_r, th) + + - p3 - a struct containing p3 parameters + - D - maximum dimension of the particle + - F_r - rime mass fraction (q_rim/ q_i) + - th - thresholds as calculated by thresholds() + +Returns the value equivalent to dm(D)/dt * 4 / D for each P3 regime + 4 / D comes from dD/dt +""" +function dmdt_mass(p3, D, F_r, th) + D_th = D_th_helper(p3) + if D_th > D + return 2 * π * p3.ρ_i * D + elseif F_r == 0 + return 4 * α_va_si(p3) * p3.β_va * D^(p3.β_va - 2) + elseif th.D_gr > D >= D_th + return 4 * α_va_si(p3) * p3.β_va * D^(p3.β_va - 2) + elseif th.D_cr > D >= th.D_gr + return 2 * π * th.ρ_g * D + elseif D >= th.D_cr + return 4 * α_va_si(p3) / (1 - F_r) * p3.β_va * D^(p3.β_va - 2) + end +end + +""" + p3_melt(p3, Chen2022, aps, tps, q, N, ρ, T, ρ_a, F_r, ρ_r) + + - p3 - a struct containing p3 parameters + - Chen2022 - struct containing Chen 2022 velocity parameters + - aps - air properties + - tps - thermodynamics parameters + - q - mass mixing ratio of ice + - N - number mixing ratio of ice + - T - temperature (K) + - ρ_a - air density + - F_r - rime mass fraction (q_rim/ q_i) + - ρ_r - rime density (q_rim/B_rim) + + Returns the calculated melting rate of ice + Equivalent to the measure of QIMLT in Morrison and Mildbrandt (2015) +""" +function p3_melt( + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + aps::CMP.AirProperties{FT}, + tps::TDP.ThermodynamicsParameters{FT}, + q::FT, + N::FT, + T::FT, + ρ_a::FT, + F_r::FT, + ρ_r::FT, +) where {FT} + # Get constants + (; ν_air, D_vapor, K_therm) = aps + a = p3.vent_a + b = p3.vent_b + L_f = TD.latent_heat_fusion(tps, T) + T_freeze = p3.T_freeze + N_sc = ν_air / D_vapor + ρ_l = p3.ρ_l + + # Get distribution values + th = thresholds(p3, ρ_r, F_r) + (λ, N_0) = distribution_parameter_solver(p3, q, N, ρ_r, F_r) + + # Get bound + ice_bound = get_ice_bound(p3, λ, eps(FT)) + + # Define function pieces + N_re(D) = + D * TV.velocity_chen( + D, + Chen2022.snow_ice, + ρ_a, + p3_mass(p3, D, F_r, th), + p3_area(p3, D, F_r, th), + p3.ρ_i, + ) / ν_air + F(D) = a + b * N_sc^(1 / 3) * N_re(D)^(1 / 2) + dmdt(D) = + dmdt_mass(p3, D, F_r, th) / ρ_l * K_therm / L_f * (T - T_freeze) * F(D) + f(D) = 1 / (2 * ρ_a) * dmdt(D) * N′ice(D, p3, λ, N_0) + + # Integrate + (dqdt, error) = QGK.quadgk(d -> f(d), FT(0), 2 * ice_bound) + + return dqdt +end + +""" + p3_het_freezing(mass, tps, q, N, T, ρ_a, qᵥ, aero_type) + + - mass - true if calculating change in mass, false for change in number + - tps - thermodynamics parameters + - q - mass mixing ratio of rain + - N - number mixing ratio of rain + - T - temperature in K + - ρ_a - density of air + - qᵥ - mixing ratio of water vapor + - aero_type - type of aerosols present + + Returns the rate of hetergoeneous freezing within rain or cloud water + If mass false corresponds to NRHET in Morrison and Mildbrandt (2015) + If mass true corresponds to QRHET in Morrison and Mildbrandt (2015) +""" +function p3_rain_het_freezing( + mass::Bool, + tps::TDP.ThermodynamicsParameters{FT}, + q::FT, + N::FT, + T::FT, + ρ_a::FT, + qᵥ::FT, + aero_type, +) where {FT} + ρ_w = p3.ρ_l + + Rₐ = TD.gas_constant_air(tps, TD.PhasePartition(qᵥ)) + R_v = TD.Parameters.R_v(tps) + e = qᵥ * ρ_a * R_v / Rₐ + + a_w = CO.a_w_eT(tps, e, T) + a_w_ice = CO.a_w_ice(tps, T) + Δa_w = a_w - a_w_ice + J_immersion = CM.HetIceNucleation.ABIFM_J(aero_type, Δa_w) + + bound = get_rain_bound(q, N, eps(FT)) + if mass + f_rain_mass(D) = + J_immersion * N′rain(D, q, N) * mass_s(D, ρ_w) * π * D^2 + dqdt, = QGK.quadgk(d -> f_rain_mass(d), FT(0), 2 * bound) + return dqdt + else + f_rain(D) = J_immersion * N′rain(D, q, N) * π * D^2 + dNdt, = QGK.quadgk(d -> f_rain(d), FT(0), 2 * bound) + return dNdt + end +end + end diff --git a/src/TerminalVelocity.jl b/src/TerminalVelocity.jl new file mode 100644 index 0000000000..2b4d8ffa17 --- /dev/null +++ b/src/TerminalVelocity.jl @@ -0,0 +1,85 @@ +""" +Terminal velocity calculations including: +- Chen 2022 parametrization for ice and rain +""" +module TerminalVelocity + +import CloudMicrophysics.Common as CO +import CloudMicrophysics.Parameters as CMP + +export velocity_chen + +""" + ϕ_ice(mass, area, ρ_i) + + - mass - mass of particle + - area - area of particle + - ρ_i - density of ice + + Returns the aspect ratio (ϕ) for a particle with given mass, area, and ice density +""" +function ϕ_ice(mass::FT, area::FT, ρ_i::FT) where {FT} + return 16 * ρ_i^2 * area^3 / (9 * π * mass^2) +end + +""" + velocity_chen(p3, Chen2022, ρ_a, D, F_r, th) + + - D - maximum particle dimension + - Chen2022 - a struct with terminal velocity parameters as in Chen(2022) + - ρ_a - density of air + - mass - mass of particle + - area - area of particle + - ρ_i - density of ice + + Returns the terminal velocity of ice at given particle dimension using Chen 2022 parametrizations +""" +function velocity_chen( + D::FT, + Chen2022::CMP.Chen2022VelTypeSnowIce, + ρ_a::FT, + mass::FT, + area::FT, + ρ_i::FT, +) where {FT} + if D <= Chen2022.cutoff + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + else + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + end + + κ = FT(-1 / 6) + + v = 0 + for i in 1:2 + v += ai[i] * D^(bi[i]) * exp(-ci[i] * D) + end + + return ϕ_ice(mass, area, ρ_i)^κ * v +end + +""" + velocity_chen(p3, Chen2022, ρ_a, D) + + - D - maximum particle dimension + - Chen2022 - a struct with terminal velocity parameters as in Chen(2022) + - ρ_a - density of air + + Returns the terminal velocity of rain given Chen 2022 velocity parametrizations +""" +function velocity_chen( + D::FT, + Chen2022::CMP.Chen2022VelTypeRain, + ρ_a::FT, +) where {FT} + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + + v = 0 + for i in 1:3 + v += ai[i] * D^(bi[i]) * exp(-ci[i] * D) + end + + return v +end + +end diff --git a/src/parameters/MicrophysicsP3.jl b/src/parameters/MicrophysicsP3.jl index 773ee57db2..8e79cc7730 100644 --- a/src/parameters/MicrophysicsP3.jl +++ b/src/parameters/MicrophysicsP3.jl @@ -30,6 +30,12 @@ Base.@kwdef struct ParametersP3{FT} <: ParametersType{FT} c::FT "Limiter for shape parameter mu for ice. See eq 3 in Morrison and Milbrandt 2015. Units: [-]" μ_max::FT + "Water freeze temperature [K]" + T_freeze::FT + "ventillation factor a" + vent_a::FT + "ventillation factor b" + vent_b::FT end ParametersP3(::Type{FT}) where {FT <: AbstractFloat} = @@ -47,6 +53,9 @@ function ParametersP3(td::CP.AbstractTOMLDict) :Heymsfield_mu_coeff2 => :b, :Heymsfield_mu_coeff3 => :c, :Heymsfield_mu_cutoff => :μ_max, + :temperature_water_freeze => :T_freeze, + :p3_ventillation_a => :vent_a, + :p3_ventiallation_b => :vent_b, ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) diff --git a/src/parameters/TerminalVelocity.jl b/src/parameters/TerminalVelocity.jl index 8cb46fe373..3d60170590 100644 --- a/src/parameters/TerminalVelocity.jl +++ b/src/parameters/TerminalVelocity.jl @@ -157,6 +157,8 @@ struct Chen2022VelTypeSnowIce{FT} <: ParametersType{FT} Fl::FT Gl::FT Hl::FT + "cutoff for small vs large ice particle dimension [m]" + cutoff::FT "density of cloud ice [kg/m3]" ρᵢ::FT end @@ -177,6 +179,7 @@ function Chen2022VelTypeSnowIce(toml_dict::CP.AbstractTOMLDict) :Chen2022_table_B5_Fl => :Fl, :Chen2022_table_B5_Gl => :Gl, :Chen2022_table_B5_Hl => :Hl, + :Chen2022_ice_cutoff => :cutoff, :density_ice_water => :ρᵢ, ) p = CP.get_parameter_values(toml_dict, name_map, "CloudMicrophysics") @@ -197,6 +200,7 @@ function Chen2022VelTypeSnowIce(toml_dict::CP.AbstractTOMLDict) p.Gl[1] + p.Gl[2] * log(p.ρᵢ) * sqrt(p.ρᵢ) + p.Gl[3] / sqrt(p.ρᵢ) )^(-1), p.Hl[1] + p.Hl[2] * (p.ρᵢ)^(5 / 2) - p.Hl[3] * exp(-p.ρᵢ), + p.cutoff, p.ρᵢ, ) end diff --git a/test/p3_tests.jl b/test/p3_tests.jl index e2dab5a4f5..f7569fa79b 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -3,6 +3,10 @@ import ClimaParams import CloudMicrophysics as CM import CloudMicrophysics.P3Scheme as P3 import CloudMicrophysics.Parameters as CMP +import Thermodynamics as TD +import CloudMicrophysics.TerminalVelocity as TV + +import QuadGK as QGK @info "P3 Scheme Tests" @@ -205,6 +209,174 @@ function test_velocities(FT) end end +function test_tendencies(FT) + + p3 = CMP.ParametersP3(FT) + Chen2022 = CMP.Chen2022VelType(FT) + tps = TD.Parameters.ThermodynamicsParameters(FT) + aps = CMP.AirProperties(FT) + + TT.@testset "Collision Tendencies Smoke Test" begin + N = FT(1e8) + ρ_a = FT(1.2) + ρ_r = FT(500) + F_r = FT(0.5) + T_warm = FT(300) + T_cold = FT(200) + + qs = range(0.001, stop = 0.005, length = 5) + q_const = FT(0.05) + + cloud_expected_warm = + [6.341e-5, 0.0002099, 0.0004258, 0.0007047, 0.001042] + cloud_expected_cold = [0.002196, 0.003525, 0.004687, 0.005761, 0.006781] + rain_expected_warm = + [0.0003414, 0.0007199, 0.001117, 0.001528, 0.001951] + rain_expected_cold = [0.2919, 0.3002, 0.3059, 0.3104, 0.3140] + + for i in axes(qs, 1) + cloud_warm = P3.ice_collisions( + "cloud", + p3, + Chen2022, + qs[i], + N, + q_const, + N, + ρ_a, + F_r, + ρ_r, + T_warm, + ) + cloud_cold = P3.ice_collisions( + "cloud", + p3, + Chen2022, + qs[i], + N, + q_const, + N, + ρ_a, + F_r, + ρ_r, + T_cold, + ) + + TT.@test cloud_warm >= 0 + TT.@test cloud_warm ≈ cloud_expected_warm[i] rtol = 1e-3 + TT.@test cloud_cold >= 0 + TT.@test cloud_cold ≈ cloud_expected_cold[i] rtol = 1e-3 + + rain_warm = P3.ice_collisions( + "rain", + p3, + Chen2022, + qs[i], + N, + q_const, + N, + ρ_a, + F_r, + ρ_r, + T_warm, + ) + rain_cold = P3.ice_collisions( + "rain", + p3, + Chen2022, + qs[i], + N, + q_const, + N, + ρ_a, + F_r, + ρ_r, + T_cold, + ) + + TT.@test rain_warm >= 0 + TT.@test rain_warm ≈ rain_expected_warm[i] rtol = 1e-3 + TT.@test rain_cold >= 0 + TT.@test rain_cold ≈ rain_expected_cold[i] rtol = 1e-3 + end + end + + TT.@testset "Melting Tendencies Smoke Test" begin + N = FT(1e8) + ρ_a = FT(1.2) + ρ_r = FT(500) + F_r = FT(0.5) + T_freeze = FT(273.15) + + qs = range(0.001, stop = 0.005, length = 5) + + expected_melt = [0.0006982, 0.0009034, 0.001054, 0.001177, 0.001283] + + for i in axes(qs, 1) + rate = P3.p3_melt( + p3, + Chen2022, + aps, + tps, + qs[i], + N, + T_freeze + 2, + ρ_a, + F_r, + ρ_r, + ) + + TT.@test rate >= 0 + TT.@test rate ≈ expected_melt[i] rtol = 1e-3 + end + end + + TT.@testset "Heterogeneous Freezing Smoke Test" begin + T = FT(250) + N = FT(1e8) + ρ_a = FT(1.2) + qᵥ = FT(8.1e-4) + aero_type = CMP.Illite(FT) + + qs = range(0.001, stop = 0.005, length = 5) + + expected_freeze_q = + [1.109e-61, 3.519e-61, 6.918e-61, 1.117e-60, 1.621e-60] + expected_freeze_N = + [1.109e-51, 1.760e-51, 2.306e-51, 2.794e-51, 3.242e-51] + + for i in axes(qs, 1) + rate_mass = P3.p3_rain_het_freezing( + true, + tps, + qs[i], + N, + T, + ρ_a, + qᵥ, + aero_type, + ) + rate_num = P3.p3_rain_het_freezing( + false, + tps, + qs[i], + N, + T, + ρ_a, + qᵥ, + aero_type, + ) + + TT.@test rate_mass >= 0 + TT.@test rate_mass ≈ expected_freeze_q[i] rtol = 1e-3 + + TT.@test rate_num >= 0 + TT.@test rate_num ≈ expected_freeze_N[i] rtol = 1e-3 + end + end + +end + println("Testing Float32") test_p3_thresholds(Float32) #TODO - only works for Float64 now. We should switch the units inside the solver @@ -215,3 +387,4 @@ println("Testing Float64") test_p3_thresholds(Float64) test_p3_shape_solver(Float64) test_velocities(Float64) +test_tendencies(Float64) diff --git a/test/terminal_velocity_tests.jl b/test/terminal_velocity_tests.jl new file mode 100644 index 0000000000..ce1c1e15c8 --- /dev/null +++ b/test/terminal_velocity_tests.jl @@ -0,0 +1,55 @@ +import Test as TT +import ClimaParams +import CloudMicrophysics as CM +import CloudMicrophysics.P3Scheme as P3 +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.TerminalVelocity as TV + +@info "Terminal Velocity Tests" + +function test_chen_velocities(FT) + + p3 = CMP.ParametersP3(FT) + Chen2022 = CMP.Chen2022VelType(FT) + + ρ_a = FT(1.2) + + TT.@testset "Chen 2022 - Rain" begin + Ds = range(FT(1e-6), stop = FT(1e-5), length = 5) + expected = [0.002508, 0.009156, 0.01632, 0.02377, 0.03144] + for i in axes(Ds, 1) + vel = TV.velocity_chen(Ds[i], Chen2022.rain, ρ_a) + TT.@test vel >= 0 + TT.@test vel ≈ expected[i] rtol = 1e-3 + end + end + + TT.@testset "Chen 2022 - Ice" begin + F_r = FT(0.5) + ρ_r = FT(500) + th = P3.thresholds(p3, ρ_r, F_r) + # Allow for a D falling into every regime of the P3 Scheme + Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5) + expected = [0.08109, 0.3838, 0.5917, 0.8637, 1.148] + for i in axes(Ds, 1) + D = Ds[i] + vel = TV.velocity_chen( + D, + Chen2022.snow_ice, + ρ_a, + P3.p3_mass(p3, D, F_r, th), + P3.p3_area(p3, D, F_r, th), + p3.ρ_i, + ) + TT.@test vel >= 0 + TT.@test vel ≈ expected[i] rtol = 1e-3 + end + + end +end + +println("Testing Float32") +test_chen_velocities(Float32) + +println("Testing Float64") +test_chen_velocities(Float64)