add cloud terminal velocity for 2M based on Stokes regime fall speed and gamma size dist assumption (#626)

Author: sajjadazimi

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.26.2"
+version = "0.26.3"
 
 [deps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"

docs/src/API.md

@@ -353,6 +353,7 @@ Parameters.NumberAdjustmentHorn2012
 Parameters.Blk1MVelType
 Parameters.Blk1MVelTypeRain
 Parameters.Blk1MVelTypeSnow
+Parameters.StokesRegimeVelType
 Parameters.SB2006VelType
 Parameters.Chen2022VelType
 Parameters.Chen2022VelTypeSmallIce

docs/src/Microphysics2M.md

@@ -72,7 +72,7 @@ The cloud droplet number distribution, as a function of mass ``x`` [kg], is assu
     f_c(x) = A_c x^{ν_c} \exp\left(-B_c x^{μ_c} \right), 
 \end{align}
 ```
-We assume that ``ν_c`` and ``μ_c`` are fixed constants. Our default choice for these parameters are ``ν_c = 2`` and ``μ_c = 1``.
+We assume that ``ν_c`` and ``μ_c`` are fixed constants. Our default choice for these parameters are ``ν_c = 1`` and ``μ_c = 1``.
 
 The free parameters for cloud droplets (``A_c``, ``B_c``) can be found analytically by integrating over the assumed 
 mass distribution to find the prognostic variables
@@ -323,7 +323,7 @@ The default free parameter values are:
 
 |   symbol   | default value |
 |------------|---------------|
-|``ν``       | ``2``         |
+|``ν``       | ``1``         |
 |``A``       | ``400``       |
 |``a``       | ``0.7``       |
 |``b``       | ``3``         |
@@ -566,7 +566,7 @@ The number- and mass-weighted terminal velocities for rain water are obtained by
   \overline{v}_{r,\, k} = \frac{1}{M_r^k} \int_0^\infty x^k f_r(x) v(x) dx,
 \end{equation}
 ```
-where the superscript ``k`` indicates the moment number,  with``k=0`` for number density and ``k=1`` for mass.
+where the superscript ``k`` indicates the moment number,  with ``k=0`` for number density and ``k=1`` for mass.
 
 The individual terminal velocity of particles is approximated by
 ```math
@@ -622,6 +622,29 @@ include("plots/TerminalVelocity2M.jl")
 ```
 ![](2M_terminal_velocity_comparisons.svg)
 
+For cloud liquid droplets in two-moment microphysics schemes, we use the analytical Stokes-regime expression for the terminal velocity of an individual spherical particle:
+```math
+v(r) = \frac{2}{9} \frac{(\rho_{water} - \rho_{air}) g}{\mu_{air}} r^2
+```
+with ``\mu_{air} = \rho_{air} \, \nu_{air}`` and assuming constant kinematic viscosity ``\nu_{air}`` in our parameterization. 
+When droplet sizes follow a gamma distribution (as in Seifert & Beheng 2006), integrating this equation over the size spectrum yields the number- and mass-weighted mean terminal velocities: 
+```math
+\begin{align*}
+  \overline{v}_{c,\, k} = \frac{2}{9}\, \left(\frac{3}{4\, \pi\, \rho_{water}}\right)^{2/3}\, \left(\frac{\rho_{water}}{\rho_{air}}-1\right)\, \frac{g}{\nu_{air}}\, \frac{M_{k+2/3}}{M_k}
+\end{align*}
+```
+where
+- ``M_i`` is the ``i``-th moment of the droplet size distribution.
+- ``k`` is the weighting moment: ``k=0`` for number-mean velocity, ``k=1`` for mass-mean velocity.
+
+For a gamma distribution with shape parameters ``\nu_c`` and ``\mu_c``, the mean terminal velocities can be written analytically as:
+```math
+\begin{align*}
+  \overline{v}_{c,\, k} = \frac{\Gamma(\frac{k+2/3+\nu_c+1}{\mu_c})}{\Gamma(\frac{k+\nu_c+1}{\mu_c})}\left[\frac{\Gamma(\frac{\nu_c+1}{\mu_c})}{\Gamma(\frac{\nu_c+2}{\mu_c})}\right]^{2/3} v_{c,\, k}(\overline{r})
+\end{align*}
+```
+where ``v_{c,\, k}(\overline{r})`` is the fall speed of a particle with mean radius ``\overline{r}``.
+
 For the Chen et al. [Chen2022](@cite) terminal velocity parameterization,
   the number- and mass-weighted group terminal velocities are:
 ```math

docs/src/TerminalVelocity.md

@@ -1,8 +1,8 @@
 # Terminal velocity
 
-`CloudMicrophysics.jl` offers three parameterizations of the relationship between
-particle size and terminal velocity:
+`CloudMicrophysics.jl` offers several parameterizations of the relationship between particle size and terminal velocity:
  - A simple power-law used in the 1-moment microphysics scheme,
+ - An analytical Stokes-regime formulation for cloud liquid droplets in two-moment microphysics,
  - The rain terminal velocity used in Seifert and Beheng 2006 [SeifertBeheng2006](@cite),
  - The rain and ice terminal velocities described in Chen et. al. 2022  [Chen2022](@cite).
 
@@ -15,8 +15,7 @@ The above terminal velocities need to be averaged over the assumed
    and use the power-law formulation when deriving process rates such as accretion.
    The Chen et al. [Chen2022](@cite) terminal velocity is available in 1-moment scheme
    for rain and snow, but without re-deriving other process rates.
- - The 2-moment scheme can be run with either the Seifert and Beheng [SeifertBeheng2006](@cite)
-   or the Chen et al. [Chen2022](@cite) terminal velocity.
+ - In the 2-moment scheme, for rain, the Seifert and Beheng [SeifertBeheng2006] (@cite) or the Chen et al. [Chen2022](@cite) parameterizations can be used. For cloud liquid droplets, we use the analytical Stokes-regime terminal velocity formulation.
  - The P3 scheme can only be run with the Chen et al. [Chen2022](@cite) terminal velocity
    and uses it when deriving the process rates.
 See the relevant sections in 1M, 2M, P3 and non-equilibrium microphysics documentation
@@ -70,6 +69,15 @@ The coefficients for the snow terminal velocity are empirical.
 |``v_0^{sno}``  | coefficient in ``v_{term}(r)`` for snow | ``\frac{m}{s}`` | ``2^{9/4} r_0^{1/4}`` | eq (6b) [Grabowski1998](@cite) |
 |``v_e^{sno}``  | exponent in ``v_{term}(r)`` for snow    | -               | ``0.25``              | eq (6b) [Grabowski1998](@cite) |
 
+## Stokes-flow terminal velocity
+In the Stokes regime (`Re < 1`), the analytical fall speed of a spherical particle is given by
+```math
+v_{term}(r) = \frac{2}{9} \frac{(\rho_{water} - \rho_{air}) g}{\mu_{air}} r^2
+```
+where ``\mu_{air} = \rho_{air} \, \nu_{air}`` is the dynamic viscosity of air, and ``\nu_{air}`` is the kinematic viscosity. In general, ``\nu_{air}`` depends on temperature and pressure, but in our parameterization it is treated as a constant for simplicity.
+
+In two-moment cloud microphysics parameterizations (e.g., Seifert & Beheng 2006), this expression can be integrated over a gamma droplet size distribution to obtain number-weighted and mass-weighted mean fall velocities of cloud droplets.
+
 ## Seifert and Beheng 2006
 
 Seifert and Beheng [SeifertBeheng2006](@cite) uses an empirical relationship between the rain drop size

docs/src/plots/TerminalVelocity.jl

@@ -20,6 +20,7 @@ const SB2006 = CMP.SB2006(FT)
 const SB2006_no_lim = CMP.SB2006(FT, false)
 
 const Chen2022 = CMP.Chen2022VelType(FT)
+const STVel = CMP.StokesRegimeVelType(FT)
 const SB2006Vel = CMP.SB2006VelType(FT)
 const Blk1MVel = CMP.Blk1MVelType(FT)
 
@@ -122,11 +123,14 @@ end
 D_r_range = range(1e-6, stop = 6e-3, length = 1000)
 D_r_range_small = range(0.1 * 1e-6, stop = 100 * 1e-6, length = 1000)
 q_range = range(0, stop = 5 * 1e-3, length = 100)
+N_liq = 500e6
 
 #! format: off
 # velocity values for cloud particle sizes
 v_term_rain = CMO.particle_terminal_velocity(Chen2022.rain, ρ_air)
 v_term_small_ice = CMO.particle_terminal_velocity(Chen2022.small_ice, ρ_air, ice.ρᵢ)
+v_term_stokes = CMO.particle_terminal_velocity(STVel, ρ_air)
+ST_cloud_small = v_term_stokes.(D_r_range_small)
 SB_rain_small = [rain_terminal_velocity_individual_SB(SB2006Vel, ρ_air, D_r)                 for D_r in D_r_range_small]
 M1_rain_small = [terminal_velocity_individual_1M(Blk1MVel.rain, ρ_air, D_r)                  for D_r in D_r_range_small]
 M1_snow_small = [terminal_velocity_individual_1M(Blk1MVel.snow, ρ_air, D_r)                  for D_r in D_r_range_small]
@@ -137,6 +141,7 @@ Ch_snow_small = [snow_terminal_velocity_individual_Chen(snow, Chen2022.large_ice
 Ch_snow_small_oblate = [snow_terminal_velocity_individual_Chen_oblate(snow, Chen2022.large_ice, ρ_air, D_r) for D_r in D_r_range_small]
 Ch_snow_small_prolate = [snow_terminal_velocity_individual_Chen_prolate(snow, Chen2022.large_ice, ρ_air, D_r) for D_r in D_r_range_small]
 # velocity values for precip particle sizes
+ST_cloud = v_term_stokes.(D_r_range)
 SB_rain = [rain_terminal_velocity_individual_SB(SB2006Vel, ρ_air, D_r)                 for D_r in D_r_range]
 M1_rain = [terminal_velocity_individual_1M(Blk1MVel.rain, ρ_air, D_r)                  for D_r in D_r_range]
 M1_snow = [terminal_velocity_individual_1M(Blk1MVel.snow, ρ_air, D_r)                  for D_r in D_r_range]
@@ -165,13 +170,16 @@ bCh_rain = [CM1.terminal_velocity(rain, Chen2022.rain,     ρ_air, q) for q in q
 bCh_snow = [CM1.terminal_velocity(snow, Chen2022.large_ice, ρ_air, q) for q in q_range]
 bCh_snow_oblate =  [CM1.terminal_velocity(snow, Chen2022.large_ice, ρ_air, q, oblate) for q in q_range]
 bCh_snow_prolate = [CM1.terminal_velocity(snow, Chen2022.large_ice, ρ_air, q, prolate) for q in q_range]
+bSt_N_liq = [CM2.cloud_terminal_velocity(SB2006.pdf_c, STVel, q, ρ_air, N_liq)[1] for q in q_range]
+bSt_liq = [CM2.cloud_terminal_velocity(SB2006.pdf_c, STVel, q, ρ_air, N_liq)[2] for q in q_range]
 bCh_liq = [CMNe.terminal_velocity(liquid, Chen2022.rain,     ρ_air, q) for q in q_range]
 bCh_ice = [CMNe.terminal_velocity(ice,    Chen2022.small_ice, ρ_air, q) for q in q_range]
 
 # individual particle comparison - cloud size range
 p1 = PL.scatter(D_Gunn_Kinzer_small * 1e6, u_Gunn_Kinzer_small * 1e2, ms = 3, color = 4, label = "Gunn and Kinzer (1949)")
 p1 = PL.plot!(D_r_range_small  * 1e6, M1_rain_small * 1e2, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [cm/s]", label = "Rain 1M",   color = :skyblue1, margin=10mm)
 p1 = PL.plot!(D_r_range_small  * 1e6, M1_snow_small * 1e2, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [cm/s]", label = "Snow 1M",   color = :plum)
+p1 = PL.plot!(D_r_range_small  * 1e6, ST_cloud_small * 1e2, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [cm/s]", label = "Stokes",   color = :cadetblue, style = :dash)
 p1 = PL.plot!(D_r_range_small  * 1e6, SB_rain_small * 1e2, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [cm/s]", label = "Rain SB",   color = :cadetblue)
 p1 = PL.plot!(D_r_range_small  * 1e6, Ch_rain_small * 1e2, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [cm/s]", label = "Rain Chen", color = :blue)
 p1 = PL.plot!(D_r_range_small  * 1e6, Ch_snow_small         * 1e2, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [cm/s]", label = "Snow Chen", color = :darkviolet)
@@ -182,18 +190,22 @@ p1 = PL.plot!(D_r_range_small  * 1e6, Ch_ice_small * 1e2,  linewidth = 3, xlabel
 p2 = PL.scatter(D_Gunn_Kinzer * 1e3, u_Gunn_Kinzer, ms = 3, color = 4, legend=false)
 p2 = PL.plot!(D_r_range * 1e3, M1_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", legend=false, color = :skyblue1)
 p2 = PL.plot!(D_r_range * 1e3, M1_snow, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", legend=false, color = :plum)
+p2 = PL.plot!(D_r_range * 1e3, ST_cloud, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [m/s]", legend=false, color = :cadetblue, style = :dash)
 p2 = PL.plot!(D_r_range * 1e3, SB_rain, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [m/s]", legend=false, color = :cadetblue)
 p2 = PL.plot!(D_r_range * 1e3, Ch_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", legend=false, color = :blue)
 p2 = PL.plot!(D_r_range * 1e3, Ch_snow,         linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", legend=false, color = :darkviolet)
 p2 = PL.plot!(D_r_range * 1e3, Ch_snow_oblate,  linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", legend=false, color = :darkviolet, style = :dash)
 p2 = PL.plot!(D_r_range * 1e3, Ch_snow_prolate, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", legend=false, color = :darkviolet, style = :dot)
 p2 = PL.plot!(D_r_range * 1e3, Ch_ice,  linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", legend=false, color = :orange)
+p2 = PL.plot!(ylim=(-1.0, 11.0))
 # snow aspect ratio
 p3 = PL.plot(D_r_range  * 1e3,  Aspect_Ratio_oblate,  linewidth = 3, xlabel = "D [mm]", ylabel = "aspect ratio", label = "Snow aspect ratio oblate",  color = :darkviolet, ylim = (0, 100), style = :dash)
 p3 = PL.plot!(D_r_range * 1e3,  Aspect_Ratio_prolate, linewidth = 3, xlabel = "D [mm]", ylabel = "aspect ratio", label = "Snow aspect ratio prolate", color = :darkviolet, style = :dot)
 # group velocity comparison
-p4 = PL.plot(q_range * 1e3,  bCh_liq * 100,  linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [cm/s]", color = :green,      label="Liq Chen", margin=10mm)
-p4 = PL.plot!(q_range * 1e3, bCh_ice * 100, linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [cm/s]", color = :orange,     label="Ice Chen")
+p4 = PL.plot(q_range * 1e3,  bSt_liq * 100,  linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [cm/s]", color = :cadetblue, label="Liq Stokes", margin=10mm)
+p4 = PL.plot!(q_range * 1e3,  bSt_N_liq * 100,  linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [cm/s]", color = :cadetblue, style = :dash, label="Liq Num Stokes")
+p4 = PL.plot!(q_range * 1e3,  bCh_liq * 100,  linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [cm/s]", color = :orange,      label="Liq Chen")
+p4 = PL.plot!(q_range * 1e3, bCh_ice * 100, linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [cm/s]", color = :blue,     label="Ice Chen")
 
 p5 = PL.plot(q_range  * 1e3, bM1_rain, linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [m/s]", color = :skyblue1,   label="Rain 1M")
 p5 = PL.plot!(q_range * 1e3, bM1_snow, linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [m/s]", color = :plum,       label="Snow 1M")

src/Common.jl

@@ -341,6 +341,22 @@ function particle_terminal_velocity(velocity_params::CMP.TerminalVelocityType, 
     return v_term
 end
 
+"""
+    particle_terminal_velocity(velocity_params::CMP.StokesRegimeVelType{FT}, ρ::FT)
+
+ - `velocity_params` - set with free parameters
+ - `ρ` - air density
+
+Returns a function `v_term(D)` that computes the analytical fall speed of a cloud droplet as a function of
+its size (diameter, `D`) in the Stokes regime (Re < 1)
+"""
+function particle_terminal_velocity(vel::CMP.StokesRegimeVelType, ρ)
+    (; ρw, grav, ν_air) = vel
+    FT = eltype(ρ)
+    terminal_velocity_prefactor = FT(1 / 18) * (ρw / ρ - 1) * grav / ν_air
+    v_term(D) = terminal_velocity_prefactor * D^2
+    return v_term
+end
 
 """
     volume_sphere_D(D)

src/Microphysics2M.jl

@@ -598,6 +598,49 @@ function rain_self_collection_and_breakup(
     return (; sc, br)
 end
 
+"""
+    cloud_terminal_velocity(pdf_c, vel_params, q_liq, ρₐ, N_liq)
+
+Compute the number-averaged and mass-averaged terminal velocities of cloud droplets
+assuming a gamma size distribution for droplet mass and the analytical Stokes-regime terminal
+velocity of spherical particles.
+
+# Arguments
+- `pdf_c`: Cloud droplet size distribution parameters, [`CMP.CloudParticlePDF_SB2006`](@ref).
+- `vel_params`: Terminal velocity parameters, [`CMP.StokesRegimeVelType`](@ref).
+- `q_liq`: Cloud liquid water specific content [kg kg⁻¹].
+- `ρₐ`: Air density [kg m⁻³].
+- `N_liq`: Cloud droplet number concentration [m⁻³].
+
+# Returns
+A tuple containing the number- and mass-weighted mean fall velocities of cloud droplets in [m/s].
+Individual droplet terminal velocities follow v_{term}(D) = (1/18) (ρw - ρₐ) g D^2 / μ_air with 
+μ_air = ρₐ * ν_air and assuming constant ν_air.
+"""
+function cloud_terminal_velocity(
+    pdf_c::CMP.CloudParticlePDF_SB2006{FT},
+    (; ρw, grav, ν_air)::CMP.StokesRegimeVelType{FT},
+    q_liq, ρₐ, N_liq,
+) where {FT}
+
+    (; νc, μc) = pdf_c
+    (; Bc) = pdf_cloud_parameters_mass(pdf_c, q_liq, ρₐ, N_liq)
+
+    terminal_velocity_prefactor = FT(1 / 18) * (6 / pi / ρw)^(2 / 3) * (ρw / ρₐ - 1) * grav / ν_air
+    vt0 =
+        N_liq < eps(FT) ? FT(0) :
+        terminal_velocity_prefactor * DT.generalized_gamma_Mⁿ(νc, μc, Bc, N_liq, FT(2 / 3)) / N_liq
+    vt1 =
+        q_liq < eps(FT) ? FT(0) :
+        terminal_velocity_prefactor * DT.generalized_gamma_Mⁿ(νc, μc, Bc, N_liq, FT(5 / 3)) / ρₐ / q_liq
+
+    if q_liq >= 5e-3
+        @show terminal_velocity_prefactor, vt0, vt1
+    end
+    return (vt0, vt1)
+
+end
+
 """
     rain_terminal_velocity(SB2006, vel, q_rai, ρ, N_rai)
 

src/parameters/TerminalVelocity.jl

@@ -1,4 +1,4 @@
-export Blk1MVelType, SB2006VelType, Chen2022VelType
+export Blk1MVelType, StokesRegimeVelType, SB2006VelType, Chen2022VelType
 
 """
     Blk1MVelTypeRain
@@ -101,12 +101,38 @@ function Blk1MVelType(toml_dict::CP.AbstractTOMLDict)
     return Blk1MVelType{FT, typeof(rain), typeof(snow)}(rain, snow)
 end
 
+"""
+    StokesRegimeVelType
+
+The type for precipitation terminal velocity in the Stokes regime (Re < 1)
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+Base.@kwdef struct StokesRegimeVelType{FT} <: TerminalVelocityType{FT}
+    ρw::FT
+    ν_air::FT
+    grav::FT
+end
+
+StokesRegimeVelType(::Type{FT}) where {FT <: AbstractFloat} =
+    StokesRegimeVelType(CP.create_toml_dict(FT))
+
+function StokesRegimeVelType(td::CP.AbstractTOMLDict)
+    name_map = (;
+        :density_liquid_water => :ρw,
+        :kinematic_viscosity_of_air => :ν_air,
+        :gravitational_acceleration => :grav,
+    )
+    parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
+    FT = CP.float_type(td)
+    return StokesRegimeVelType{FT}(; parameters...)
+end
 
 """
     SB2006VelType
 
 The type for precipitation terminal velocity from Seifert and Beheng 2006
-(Defined only for rain)
 
 # Fields
 $(DocStringExtensions.FIELDS)
@@ -116,6 +142,9 @@ Base.@kwdef struct SB2006VelType{FT} <: TerminalVelocityType{FT}
     aR::FT
     bR::FT
     cR::FT
+    ρw::FT
+    ν_air::FT
+    grav::FT
 end
 
 SB2006VelType(::Type{FT}) where {FT <: AbstractFloat} =
@@ -127,6 +156,9 @@ function SB2006VelType(td::CP.AbstractTOMLDict)
         :SB2006_raindrops_terminal_velocity_coeff_aR => :aR,
         :SB2006_raindrops_terminal_velocity_coeff_bR => :bR,
         :SB2006_raindrops_terminal_velocity_coeff_cR => :cR,
+        :density_liquid_water => :ρw,
+        :kinematic_viscosity_of_air => :ν_air,
+        :gravitational_acceleration => :grav,
     )
     parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
     FT = CP.float_type(td)

test/microphysics2M_tests.jl

@@ -35,6 +35,7 @@ function test_microphysics2M(FT)
     tps = TDI.TD.Parameters.ThermodynamicsParameters(FT)
 
     # Terminal velocity parameters
+    STVel = CMP.StokesRegimeVelType(FT)
     SB2006Vel = CMP.SB2006VelType(FT)
     Chen2022Vel = CMP.Chen2022VelTypeRain(FT)
 
@@ -383,6 +384,42 @@ function test_microphysics2M(FT)
         end
     end
 
+    TT.@testset "2M_microphysics - cloud terminal velocity" begin
+        #setup
+        ρ = FT(1.1)
+        q_liq = FT(1e-3)
+        N_liq = FT(1e7)
+
+        (; ρw, grav, ν_air) = STVel
+
+        #action
+        vt_liq = CM2.cloud_terminal_velocity(SB2006.pdf_c, STVel, q_liq, ρ, N_liq)
+
+        (; νc, μc, ρw) = SB2006.pdf_c
+        (; Bc) = CM2.pdf_cloud_parameters_mass(SB2006.pdf_c, q_liq, ρ, N_liq)
+        terminal_velocity_prefactor = FT(2 / 9) * (3 / 4 / pi / ρw)^(2 / 3) * (ρw / ρ - 1) * grav / ν_air
+        vt0 = terminal_velocity_prefactor * DT.generalized_gamma_Mⁿ(νc, μc, Bc, N_liq, FT(2 / 3)) / N_liq
+        vt1 = terminal_velocity_prefactor * DT.generalized_gamma_Mⁿ(νc, μc, Bc, N_liq, FT(5 / 3)) / ρ / q_liq
+
+        #test
+        TT.@test vt_liq isa Tuple
+        TT.@test vt_liq[1] ≈ vt0 rtol = 1e-6
+        TT.@test vt_liq[2] ≈ vt1 rtol = 1e-6
+
+        TT.@test CM2.cloud_terminal_velocity(
+            SB2006.pdf_c, STVel, q_liq, ρ, FT(0),
+        )[1] ≈ 0 atol = eps(FT)
+        TT.@test CM2.cloud_terminal_velocity(
+            SB2006.pdf_c, STVel, FT(0), ρ, N_liq,
+        )[2] ≈ 0 atol = eps(FT)
+        TT.@test CM2.cloud_terminal_velocity(
+            SB2006.pdf_c, STVel, FT(0), ρ, FT(0),
+        )[1] ≈ 0 atol = eps(FT)
+        TT.@test CM2.cloud_terminal_velocity(
+            SB2006.pdf_c, STVel, FT(0), ρ, FT(0),
+        )[2] ≈ 0 atol = eps(FT)
+    end
+
     TT.@testset "2M_microphysics - Seifert and Beheng 2006 rain terminal velocity with limiters" begin
         #setup
         ρ = FT(1.1)