diff --git a/docs/src/Microphysics1M.md b/docs/src/Microphysics1M.md index ad1682ea7..bbd171d3f 100644 --- a/docs/src/Microphysics1M.md +++ b/docs/src/Microphysics1M.md @@ -109,13 +109,93 @@ where: - ``\rho`` is the density of air. !!! note + Assuming a constant drag coefficient is an approximation and it should + be size and flow dependent, see for example + [here](https://www1.grc.nasa.gov/beginners-guide-to-aeronautics/drag-of-a-sphere/). + We are in the process of updating the 1-moment microphysics scheme to formulae from [Chen2022](@cite). + Other possibilities: [Khvorostyanov2002](@cite) or [Karrer2020](@cite) + +[Chen2022](@cite) provides a terminal velocity parameterisation + based on an empirical fit to a high accuracy model. +The terminal velocity depends on particle shape, size and density, + consideres the deformation effects of large rain drops, + as well as size-specific air density dependence. +The fall speed of individual particles $v(D)$ is parameterized as: +```math +\begin{equation} + v_{term}(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D} +\end{equation} +``` +where: + - ``D`` is the particle diameter, + - ``a_i``, ``b_i``, ``c_i`` are the free parameters, + - ``\phi`` is the aspect ratio, and + - ``\kappa`` is a parameter that depends on the particle shape (``\kappa=0`` for spheres, ``\kappa=-1/3`` for oblate and ``\kappa=1/6`` for prolate spheroids). + +For ice and snow ``j=2`` and for rain ``j=3``, to account for deformation at larger sizes. +For rain and ice we assume ``\phi=1`` (spherical). +For snow we assume ``\kappa = -1/3`` and + find the aspect ratio that is consistent with the assumed ``m(r)`` and ``a(r)`` relationships. +The aspect ratio is defined as: +```math +\begin{equation} + \phi \equiv c/a +\end{equation} +``` +where: + - ``a`` is the basal plane axial half-length, and + - ``c`` is perpendicular to the basal plane. + +The volume of a spheroid can be represented as ``V_p = 4\pi/3 \; a^2 c`` + and the area can be represented as ``A_p = \pi a c``. +It follows that + ``c = (4A_p^2) / (3 \pi V_p)``, + ``a = (3V_p) / (4A_p)``, and + ``\phi = (16 A_p^3) / (9 \pi V_p^2)``. +The volume and area are defined by the assumed power-law size relations + ``V_p = m(r) / (\rho_{ice})``, ``A_p = a(r)``. +As a result the terminal velocity of individual snow particles as: +```math +\begin{equation} + v_{term}(r) = \left(\frac{16 \; \rho_{ice}^2 \; a_0^3 \; (r/r_0)^{3a_e}}{9 \pi \; m_0^2 \; (r/r_0)^{2 m_e}} \right)^{\kappa} \sum_{i=1}^{2} \; a_i (2r)^{b_i} e^{-2 c_i r} +\end{equation} +``` +where $r$ is the radius of the particle. + +The rain ``a_i``, ``b_i``, and ``c_i`` are listed in the table below. +The formula is applicable when ``D > 0.1 mm``, +$q$ refers to ``q = e^{0.115231 \; \rho_a}``, where ``\rho_a`` is air density [kg/m3]. +The units are: [v] = m/s, [D]=mm, [``a_i``] = ``mm^{-b_i} m/s``, [``b_i``] is dimensionless, [``c_i``] = 1/mm. + +| ``i`` | ``a_i`` | ``b_i`` | ``c_i`` | +|---------|-------------------------------------------|--------------------------------------|---------------------| +| 1 | `` 0.044612 \; q`` | ``2.2955 \; -0.038465 \; \rho_a`` | ``0`` | +| 2 | ``-0.263166 \; q`` | ``2.2955 \; -0.038465 \; \rho_a`` | ``0.184325`` | +| 3 | ``4.7178 \; q \; (\rho_a)^{-0.47335}`` | ``1.1451 \; -0.038465 \; \rho_a`` | ``0.184325`` | + +The ice and snow ``a_i``, ``b_i``, and ``c_i`` are listed in the table below. +The formula is applicable when ``D < 0.625 mm``. + +| ``i`` | ``a_i`` | ``b_i`` | ``c_i`` | +|-------|-----------------------------|-------------------------|-----------| +| 1 | ``E_s (\rho_a)^{A_s}`` | ``B_s + C_s \rho_a`` | ``0`` | +| 2 | ``F_s (\rho_a)^{A_s}`` | ``B_s + C_s \rho_a`` | ``G_s`` | + +| Coefficient | Formula | +|--------------|------------------------------------------------------------------------------------------------| +| ``A_s`` | ``0.00174079 \log{(\rho_{ice})}^2 − 0.0378769 \log{(\rho_{ice})} - 0.263503`` | +| ``B_s`` | ``(0.575231 + 0.0909307 \log{(\rho_{ice})} + 0.515579 / \sqrt{\rho_{ice}})^{-1}`` | +| ``C_s`` | ``-0.345387 + 0.177362 \, \exp{(-0.000427794 \rho_{ice})} + 0.00419647 \sqrt{\rho_{ice}}`` | +| ``E_s`` | ``-0.156593 - 0.0189334 \log{(\rho_{ice})}^2 + 0.1377817 \sqrt{\rho_{ice}}`` | +| ``F_s`` | ``- \exp{[-3.35641 - 0.0156199 \log{\rho_{ice}}^2 + 0.765337 \log{\rho_{ice}}]}`` | +| ``G_s`` | ``(-0.0309715 + 1.55054 / \log{(\rho_{ice})} - 0.518349 log{(\rho_{ice})} / \rho_{ice})^{-1}`` | + +The terminal velocity formulae from the current default 1-moment scheme and [Chen2022](@cite) are plotted below. +```@example +include("TerminalVelocityComparisons.jl") +``` +![](1M_terminal_velocity_comparisons.svg) - It would be great to replace the above simple power laws - with more accurate relationships. - For example: - [Khvorostyanov2002](@cite) - or - [Karrer2020](@cite) ## Assumed particle size distributions @@ -259,8 +339,7 @@ F(r) = a_{vent} + ## Terminal velocity -One way to define the mass weighted terminal velocity ``v_t`` (following -[Ogura1971](@cite)) is: +The mass weighted terminal velocity ``v_t`` (following [Ogura1971](@cite)) is: ```math \begin{equation} v_t = \frac{\int_0^\infty n(r) \, m(r) \, v_{term}(r) \, dr} @@ -268,8 +347,8 @@ One way to define the mass weighted terminal velocity ``v_t`` (following \label{eq:vt} \end{equation} ``` -Integrating over the assumed Marshall-Palmer distribution and using the - ``m(r)`` and ``v_{term}(r)`` relationships results in +Integrating the default 1-moment ``m(r)`` and ``v_{term}(r)`` relationships + over the assumed Marshall-Palmer distribution results in group terminal velocity: ```math \begin{equation} v_t = \chi_v \, v_0 \, \left(\frac{1}{r_0 \, \lambda}\right)^{v_e + \Delta_v} @@ -278,152 +357,30 @@ Integrating over the assumed Marshall-Palmer distribution and using the \end{equation} ``` -!!! note - - Assuming a constant drag coefficient is an approximation and it should - be size and flow dependent, see for example - [here](https://www1.grc.nasa.gov/beginners-guide-to-aeronautics/drag-of-a-sphere/). - In general we should implement these terminal velocity parameterizations: - [Khvorostyanov2002](@cite) - or - [Karrer2020](@cite) - -Another way to implement the mass weighted terminal velocity formulas is following [Chen2022](@cite). -[Chen2022](@cite) uses multiple gamma-function terms in their parameterizations to increase accuracy. -Focusing on the mass weighted terminal velocity of rain, the fall speed of individual raindrops $v(r)$ is parameterized as: -```math -\begin{equation} - v(r) = (\phi)^{\kappa} \Sigma_{i = 1, j} \; a_i (2r)^{b_i} e^{-2c_ir} -\end{equation} -``` -where $r$ is the radius of the particle. -The shape of larger raindrops can deform; to handle this, [Chen2022](@cite) calculates the fall speed using three gamma-function terms ($j = 3$). Coefficients $a_i$, $b_i$, and $c_i$ account for deformation at larger sizes - the aspect ratio $\phi$ itself is assumed to be 1 and $\kappa$ is 0 (corresponding to a spherical raindrop). $a_i$, $b_i$, and $c_i$ are listed in the table below and also account for air-density dependence and similarly important factors. - -Table: Coefficients for Terminal Velocity Calculations: These are applicable when $r > 50 \mu m$ (where $r$ is the particle radius). -$q$ refers to $q = e^{0.115231 \rho_a}$ where $\rho_a$ refers to air density. - -| $i$ | $a_i$ | $b_i$ | $c_i$ | -|----------------------------|-----------------------------|--------------------------|------------------------| -| 1 | $0.044612 q$ | $2.2955-0.038465\rho_a$ | $0$ | -| 2 | $-0.263166 q$ | $2.2955-0.038465\rho_a$ | $0.184325$ | -| 3 | $4.7178 q (\rho_a)^{-0.47335}$ | $1.1451-0.038465\rho_a$ | $0.184325$ | - -Integrating over the assumed Marshall-Palmer distribution, we can then find the mean terminal velocity to be: +Integrating [Chen2022](@cite) formulae for rain and ice + over the assumed Marshall-Palmer size distribution, + results in group terminal velocity: ```math \begin{equation} - \overline{v_t} = \frac{\int_0^\infty v(r) \, (2r)^k \, n(r) \, dr} - {\int_0^\infty (2r)^k \, n(r) \, dr} - = (\phi)^{\kappa} \Sigma_{i} \frac{a_i \lambda^{\delta} \Gamma(b_i + \delta)} - {(\lambda + c_i)^{b_i + \delta} \; \Gamma(\delta)} + v_t = \sum_{i=1}^{j} \frac{a_i \lambda^{\delta} \Gamma(b_i + \delta)}{(\lambda + c_i)^{b_i + 4} \; \Gamma(4)} \end{equation} ``` -where $\Gamma$ is the standard gamma function integral and $\delta = \mu + k + 1$. $\overline{v_t}$ corresponds to the mass-weighted mean terminal velocity when $k = 3$. - -For the mass-weighted terminal velocity of ice particles (assumed to be spherical), [Chen2022](@cite) first parameterizes the fall speed of individual ice crystals (similar to individual raindrops): +Finally, integrating [Chen2022](@cite) formulae for snow + over the assumed Marshall-Palmer distribution, + results in group terminal velocity: ```math \begin{equation} - v(D) = (\phi)^{\kappa} \Sigma_{i = 1, j} \; a_i (2r)^{b_i} e^{-2 c_i r} + v_t = \sum_{i=1}^{2} t_i \frac{\Gamma(3a_e \kappa - 2 m_e \kappa + b_i + k + 1)}{\Gamma(k+1)} \end{equation} ``` -where r is the radius of the particle. -[Chen2022](@cite) calculates the fall speed using two gamma-function terms ($j = 2$). Coefficients $a_i$, $b_i$, and $c_i$ are different from those used in calculating raindrop fall speed - they are listed in the table below and also account for air-density dependence and similarly important factors. The aspect ratio $\phi$ itself is assumed to be 1 and $\kappa$ is 0 (corresponding to a spherical ice crystal). - -Table: Coefficients for Terminal Velocity Calculations: These are applicable when $r < 0.3125 mm$ (where $r$ is the particle radius). -$\rho_a$ refers to air density (kg $m^{-3}$); coefficients $A_s$, $B_s$, $C_s$, $E_s$, $F_s$, and $G_s$ will be defined below this table. - -| $i$ | $a_i$ | $b_i$ | $c_i$ | -|----------------------------|-----------------------------|--------------------------|------------------------| -| 1 | $E_s (\rho_a)^{A_s}$ | $B_s + C_s \rho_a$ | $0$ | -| 2 | $F_s (\rho_a)^{A_s}$ | $B_s + C_s \rho_a$ | $G_s$ | - -Table: Details of coefficients $A_s$, $B_s$, $C_s$, $E_s$, $F_s$, and $G_s$. $\rho_i$ is the apparent density of ice particles in (kg $m^{-3}$). $\log$ refers to natural logarithm in the table below. - -| Coefficient | Formula | -|--------------|---------------------------------------------------------------| -| $A_s$ | $0.00174079 (\log{(\rho_i)})^2 − 0.0378769 \log{(\rho_i)} - 0.263503$ | -| $B_s$ | $(0.575231 + 0.0909307 \log{(\rho_i)} + 0.515579 / \sqrt{\rho_i})^{-1}$ | -| $C_s$ | $-0.345387 + 0.177362 \, \exp{(-0.000427794 \rho_i)} + 0.00419647 \sqrt{\rho_i}$ | -| $E_s$ | $-0.156593 - 0.0189334 (\log{(\rho_i)})^2 + 0.1377817 \sqrt{\rho_i}$ | -| $F_s$ | $- \exp{[-3.35641 - 0.0156199 (\log{\rho_i})^2 + 0.765337 \log{\rho_i}]}$ | -| $G_s$ | $(-0.0309715 + 1.55054 / (\log{(\rho_i)}) - 0.518349 (log{(\rho_i)}/\rho_i))^{-1}$ | - -Once again integrating over the assumed Marshall-Palmer distribution, we can then find the mean terminal velocity of spherical ice crystals to be: -```math -\begin{equation} - \overline{v_t} = \frac{\int_0^\infty v(r) \, (2r)^k \, n(r) \, dr} - {\int_0^\infty (2r)^k \, n(r) \, dr} - = (\phi)^{\kappa} \Sigma_{i} \frac{a_i \lambda^{\delta} \Gamma(b_i + \delta)} - {(\lambda + c_i)^{b_i + \delta} \; \Gamma(\delta)} -\end{equation} -``` -where $\Gamma$ is the standard gamma function integral and $\delta = \mu + k + 1$. $\overline{v_t}$ corresponds to the mass-weighted mean terminal velocity when $k = 3$. - -For the mass-weighted mean terminal velocity for snow, which is aspherical ice crystals, we must first consider the shape. [Chen2022](@cite) defines the aspect ratio to be $\phi \equiv c/a$ where $a$ is basal plane axial half-length and $c$ is perpendicular to the basal plane. -Since the volume of a spheroid can be represented as $V_p = (4 \pi a^2 c)/3$ and the area can be represented as $A_p = \pi a c$, we can define $c = (4 (A_p)^2)/(3 \pi V_p)$ and $a = (3 V_p)/(4 A_p)$. -So we can then say $\phi = c/a = ((4 (A_p)^2)/(3 \pi V_p))/((3 V_p)/(4 A_p)) = (16 (A_p)^3)/(9\pi(V_p)^2)$. -We also know -```math -\begin{equation} -m(r) = \chi_m \, m_0 \left(\frac{r}{r_0}\right)^{m_e + \Delta_m} -\end{equation} -``` - -```math -\begin{equation} -A_p = a(r)= \chi_a \, a_0 \left(\frac{r}{r_0}\right)^{a_e + \Delta_a} -\end{equation} -``` -and $V_p = m(r)/(\rho_i)$. Using this, we can then say - -```math -\begin{equation} -\phi = \frac{16 (\rho_i)^2 (a_0)^3 (r/r_0)^{3a_e}}{9 \pi (m_0)^2 (r/r_0)^{2 m_e}} -\end{equation} -``` - -We can then parameterize the terminal velocity of individual snow particles to be: -```math -\begin{equation} - v(r) = (\phi)^{\kappa} \Sigma_{i = 1, j} \; a_i (2r)^{b_i} e^{-2 c_i r} = \Sigma_{i = 1, j} \; (\frac{16 (\rho_i)^2 (a_0)^3 (r/r_0)^{3a_e}}{9 \pi (m_0)^2 (r/r_0)^{2 m_e}})^{\kappa} a_i (2r)^{b_i} e^{-2 c_i r} -\end{equation} -``` -where $r$ is the radius of the particle. -[Chen2022](@cite) calculates the fall speed using two gamma-function terms ($j = 2$). Coefficients $a_i$, $b_i$, and $c_i$ are the same as used above for spherical ice crystals. $\kappa = -1/3$ (corresponding to oblate ice crystals). - -Once again integrating over the assumed Marshall-Palmer distribution, we can then find the mean terminal velocity of snow to be: -```math -\begin{equation} - \overline{v_t} = \frac{\int_0^\infty v(r) \, (2r)^k \, n(r) \, dr} - {\int_0^\infty (2r)^k \, n(r) \, dr} - = \Sigma_{i} \frac{(16 (a_0)^3 (\rho_i)^2)^{\kappa} a_i (2)^{b_i} (2 c_i \lambda)^{-(3a_e\kappa - 2m_e\kappa + b_i + k)-1} \Gamma(3a_e\kappa - 2 m_e \kappa + b_i + k + 1)} - {(9 \pi (m_0)^2)^{\kappa} (r_0)^{3a_e\kappa - 2m_e\kappa} (\lambda)^{-k-1} \Gamma(k+1)} -\end{equation} -``` -To simplify this, we combine the constants into one: -```math -\begin{equation} -t_i = \frac{(16 (a_0)^3 (\rho_i)^2)^{\kappa} a_i (2)^{b_i} (2 c_i \lambda)^{-(3 a_e \kappa - 2 m_e \kappa + b_i + k)-1}} -{(9 \pi (m_0)^2)^{\kappa} (r_0)^{3 a_e \kappa - 2 m_e \kappa} (\lambda)^{-k-1}} -\end{equation} -``` - -So we get the bulk velocity to be +where: ```math \begin{equation} - \overline{v_t} = \frac{\int_0^\infty v(r) \, (2r)^k \, n(r) \, dr} - {\int_0^\infty (2r)^k \, n(r) \, dr} - = \Sigma_{i} \frac{t_i \Gamma(3a_e \kappa - 2 m_e \kappa + b_i + k + 1)} - {\Gamma(k+1)} +t_i = \frac{[16 a_0^3 \rho_{ice}^2]^{\kappa} \; a_i \; 2^{b_i} [2 c_i \lambda]^{-(3 a_e \kappa - 2 m_e \kappa + b_i + k)-1}} +{[9 \pi m_0^2]^{\kappa} \; r_0^{3 a_e \kappa - 2 m_e \kappa} \lambda^{-k-1}} \end{equation} ``` -where $\Gamma$ is the standard gamma function integral. $\overline{v_t}$ corresponds to the mass-weighted mean terminal velocity when $k = 3$. - -Below, we compare the terminal velocity formulas from [Ogura1971](@cite) and [Chen2022](@cite): -```@example -include("TerminalVelocityComparisons.jl") -``` -![](1M_terminal_velocity_comparisons.svg) - +and $k = 3$. ## Rain autoconversion diff --git a/docs/src/TerminalVelocityComparisons.jl b/docs/src/TerminalVelocityComparisons.jl index 7418b9767..05c5c85ba 100644 --- a/docs/src/TerminalVelocityComparisons.jl +++ b/docs/src/TerminalVelocityComparisons.jl @@ -202,18 +202,18 @@ p2 = PL.plot(q_rain_range * 1e3, SB_rain_bN, linewidth = 3, xlabel = "q_rain [g 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-Ogura-1M [M]", color = :green) +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)) +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") -PL.plot!(D_r_range * 1e3, M1_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-Ogura") -PL.plot!(D_r_range * 1e3, Ch_snow, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Snow-Chen",) -PL.plot!(D_r_range * 1e3, M1_snow, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Snow-Ogura",) +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) # TODO - add a plot with individual ice crystal terminal velocity PL.savefig("1M_terminal_velocity_comparisons.svg")