diff --git a/docs/Project.toml b/docs/Project.toml index 76d874ac0..cb504ac03 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,6 +7,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 1350e5299..e7e23c03e 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -463,7 +463,7 @@ @article{MitchellHeymsfield2005 year = {2005}, volume = {62}, number = {5}, - doi = {https://doi.org/10.1175/JAS3413.1}, + doi = {10.1175/JAS3413.1}, pages= {1637--1644} } @@ -732,7 +732,7 @@ @article{Liu1997 volume = {123}, number = {542}, pages = {1789-1795}, - doi = {https://doi.org/10.1002/qj.49712354220}, + doi = {10.1002/qj.49712354220}, year = {1997} } diff --git a/docs/make.jl b/docs/make.jl index f937a1a4d..c3a366455 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -38,6 +38,7 @@ Parameterizations = [ "1-moment precipitation microphysics" => "Microphysics1M.md", "2-moment precipitation microphysics" => "Microphysics2M.md", "P3 Scheme" => "P3Scheme.md", + "Terminal Velocity" => "TerminalVelocity.md", "Non-equilibrium cloud formation" => "MicrophysicsNonEq.md", "Smooth transition at thresholds" => "ThresholdsTransition.md", "Aerosol activation" => "AerosolActivation.md", diff --git a/docs/src/API.md b/docs/src/API.md index 507fbbd7a..511556d68 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -10,6 +10,7 @@ MicrophysicsNonEq MicrophysicsNonEq.τ_relax MicrophysicsNonEq.conv_q_vap_to_q_liq_ice MicrophysicsNonEq.conv_q_vap_to_q_liq_ice_MM2015 +MicrophysicsNonEq.terminal_velocity ``` # 0-moment precipitation microphysics @@ -131,9 +132,11 @@ Common.H2SO4_soln_saturation_vapor_pressure Common.a_w_xT Common.a_w_eT Common.a_w_ice -Common.Chen2022_vel_add -Common.Chen2022_vel_coeffs_small -Common.Chen2022_vel_coeffs_large +Common.Chen2022_monodisperse_pdf +Common.Chen2022_exponential_pdf +Common.Chen2022_vel_coeffs_B1 +Common.Chen2022_vel_coeffs_B2 +Common.Chen2022_vel_coeffs_B4 ``` # Parameters diff --git a/docs/src/Microphysics1M.md b/docs/src/Microphysics1M.md index 0658dbf87..5a116c826 100644 --- a/docs/src/Microphysics1M.md +++ b/docs/src/Microphysics1M.md @@ -20,6 +20,7 @@ The cloud microphysics variables are expressed as specific humidities: Particles are assumed to follow power-law relationships involving the mass(radius), denoted by ``m(r)``, the cross section(radius), denoted by ``a(r)``, and the terminal velocity(radius), denoted by ``v_{term}(r)``, respectively. +See terminal velocity section for more details on the available terminal velocity options. The coefficients are defined in the [ClimaParams.jl](https://github.com/CliMA/ClimaParams.jl) package and are shown in the table below. @@ -34,15 +35,12 @@ m(r) = \chi_m \, m_0 \left(\frac{r}{r_0}\right)^{m_e + \Delta_m} ```math a(r) = \chi_a \, a_0 \left(\frac{r}{r_0}\right)^{a_e + \Delta_a} ``` -```math -v_{term}(r) = \chi_v \, v_0 \left(\frac{r}{r_0}\right)^{v_e + \Delta_v} -``` where: - ``r`` is the particle radius, - ``r_0`` is the typical particle radius used to nondimensionalize, - - ``m_0, \, m_e, \, a_0, \, a_e, \, v_0, \, v_e \,`` are the default + - ``m_0, \, m_e, \, a_0, \, a_e`` are the default coefficients, - - ``\chi_m``, ``\Delta_m``, ``\chi_a``, ``\Delta_a``, ``\chi_v``, ``\Delta_v`` + - ``\chi_m``, ``\Delta_m``, ``\chi_a``, ``\Delta_a`` are the coefficients that can be used during model calibration to expand around the default values. Without calibration all ``\chi`` parameters are set to 1 @@ -72,7 +70,6 @@ With that said, the assumption about the shape of the particles is used three |``m_e^{rai}`` | exponent in ``m(r)`` for rain | - | ``3`` | | |``a_0^{rai}`` | coefficient in ``a(r)`` for rain | ``m^2`` | ``\pi \, r_0^2`` | | |``a_e^{rai}`` | exponent in ``a(r)`` for rain | - | ``2`` | | -|``v_e^{rai}`` | exponent in ``v_{term}(r)`` for rain | - | ``0.5`` | | | | | | | | |``r_0^{ice}`` | typical ice crystal radius | ``m`` | ``10^{-5} `` | | |``m_0^{ice}`` | coefficient in ``m(r)`` for ice | ``kg`` | ``\frac{4}{3} \, \pi \, \rho_{ice} \, r_0^3`` | | @@ -83,122 +80,11 @@ With that said, the assumption about the shape of the particles is used three |``m_e^{sno}`` | exponent in ``m(r)`` for snow | - | ``2`` | eq (6b) [Grabowski1998](@cite) | |``a_0^{sno}`` | coefficient in ``a(r)`` for snow | ``m^2`` | ``0.3 \pi \, r_0^2`` | ``\alpha`` in eq(16b) [Grabowski1998](@cite).| |``a_e^{sno}`` | exponent in ``a(r)`` for snow | - | ``2`` | | -|``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) | where: - ``\rho_{water}`` is the density of water, - ``\rho_{ice}`` is the density of ice. -The terminal velocity of an individual rain drop is defined by the balance - between the gravitational acceleration (taking into account - the density difference between water and air) and the drag force. -Therefore the ``v_0^{rai}`` is defined as -```math -\begin{equation} -v_0^{rai} = - \left( - \frac{8}{3 \, C_{drag}} \left( \frac{\rho_{water}}{\rho} -1 \right) - \right)^{1/2} (g r_0^{rai})^{1/2} -\label{eq:vdrop} -\end{equation} -``` -where: - - ``g`` is the gravitational acceleration, - - ``C_{drag}`` is the drag coefficient, - - ``\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. - -Here we plot the terminal velocity formulae from the current default 1-moment scheme and [Chen2022](@cite). -We also show the aspect ratio of snow particles. - -```@example -include("plots/TerminalVelocityComparisons.jl") -``` -![](1M_individual_terminal_velocity_comparisons.svg) - -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}`` | - - ## Assumed particle size distributions The particle size distributions are assumed to follow @@ -287,7 +173,6 @@ They consist of: | symbol | definition | units | default value | reference | |----------------------------|-----------------------------------------------------------|--------------------------|------------------------|-----------| -|``C_{drag}`` | rain drop drag coefficient | - | ``0.55`` | ``C_{drag}`` is such that the mass averaged terminal velocity is close to [Grabowski1996](@cite) | |``\tau_{acnv\_rain}`` | cloud liquid to rain water autoconversion timescale | ``s`` | ``10^3`` | eq (5a) [Grabowski1996](@cite) | |``\tau_{acnv\_snow}`` | cloud ice to snow autoconversion timescale | ``s`` | ``10^2`` | | |``q_{liq\_threshold}`` | cloud liquid to rain water autoconversion threshold | - | ``5 \cdot 10^{-4}`` | eq (5a) [Grabowski1996](@cite) | @@ -356,7 +241,9 @@ The mass weighted terminal velocity ``v_t`` (following [Ogura1971](@cite)) is: \label{eq:vt} \end{equation} ``` -Integrating the default 1-moment ``m(r)`` and ``v_{term}(r)`` relationships +See [here](https://clima.github.io/CloudMicrophysics.jl/dev/TerminalVelocity.html) + for discussion of the different parameterizations of ``v_{term}``. +Integrating the 1-moment ``m(r)`` and power-law ``v_{term}(r)`` relationships over the assumed Marshall-Palmer distribution results in group terminal velocity: ```math \begin{equation} @@ -365,31 +252,21 @@ Integrating the default 1-moment ``m(r)`` and ``v_{term}(r)`` relationships {\Gamma(m_e + \Delta_m + 1)} \end{equation} ``` - -Integrating [Chen2022](@cite) formulae for rain and ice - over the assumed Marshall-Palmer size distribution, - results in group terminal velocity: -```math -\begin{equation} - 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} -``` -Finally, integrating [Chen2022](@cite) formulae for snow - over the assumed Marshall-Palmer distribution, - results in group terminal velocity: -```math -\begin{equation} - 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: +Integrating the Chen et al. [Chen2022](@cite) formulae over the assumed Marshall-Palmer size distribution + results in the group terminal velocity (eq. 20 in [Chen2022](@cite)): ```math \begin{equation} -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}} + v_t = \phi_{avg}^\kappa \sum_{i=1}^{j} \frac{a_i \lambda^{\delta} \Gamma(b_i + \delta)}{(\lambda + c_i)^{b_i + \delta} \; \Gamma(\delta)}, \end{equation} ``` -and $k = 3$. +where ``\delta = 4`` for the case of an exponential size distribution and the mass-weighted mean. +For snow, for simplicity, we first compute the + mass-weighted mean aspect ratio over the size distribution of particles ``\phi_{avg}`` + and then treat this as constant when computing the group terminal velocity. + +!!! note + For snow, we only use the B4 coefficients from [Chen2022](@cite). + We should switch to doing partial integrals and include also the B2 coefficients. ## Rain autoconversion @@ -759,7 +636,9 @@ If ``T > T_{freeze}``: ## Rain radar reflectivity -The rain radar reflectivity factor (``Z``) is used to measure the power returned by a radar signal when it encounters rain particles, and it is defined as the sixth moment of the rain particles distribution: +The rain radar reflectivity factor Z is used to measure the power + returned by a radar signal when it encounters rain particles. +It is defined as the 6th moment of the rain particle size distribution: ```math \begin{equation} Z = {\int_0^\infty r^{6} \, n(r) \, dr}. @@ -776,20 +655,19 @@ where: - ``n_{0}^{rai}`` - rain drop size distribution parameter, - ``\lambda`` - as defined in eq. 7 -By dividing ``Z`` with the equivalent return of a ``1 mm`` drop in a volume of a meter cube (``Z_0``) and applying the decimal logarithm to the result, we obtains the logarithmic rain radar reflectivity ``L_Z``, which is the variable that is commonly used to refer to the radar reflectivity values: +By dividing ``Z`` by the radar reflectivity factor ``Z_0`` for a drop of radius ``1 mm`` in a volume of ``1 m^3`` and applying the decimal logarithm to the result, we obtain the normalized logarithmic rain radar reflectivity ``L_Z``, which is the variable that is commonly referenced for radar reflectivity values: ```math \begin{equation} -L_Z = {10 \, \log_{10}(\frac{Z}{Z_0})}. +L_Z = {10 \, \log_{10} \left( \frac{Z}{Z_0} \right)}. \end{equation} ``` -The resulting logarithmic dimensionless unit is decibel relative to ``Z``, or ``dBZ``. +The resulting logarithmic dimensionless unit is decibel relative to ``Z_0``. ## Example figures ```@example include("plots/Microphysics1M_plots.jl") ``` -![](terminal_velocity.svg) ![](autoconversion_rate.svg) ![](accretion_rate.svg) ![](accretion_rain_sink_rate.svg) diff --git a/docs/src/Microphysics2M.md b/docs/src/Microphysics2M.md index 92cb72810..aefb12c6c 100644 --- a/docs/src/Microphysics2M.md +++ b/docs/src/Microphysics2M.md @@ -302,61 +302,6 @@ The default free parameter values are: !!! note In the paper for ``\overline{D}_{eq} < \overline{D}_r`` the equation ``\Phi_{br}(\Delta \overline{D}_r) = 2 exp(\kappa_{br} \Delta \overline{D}_r) -1`` is given. This equations seems to be missing parentheses as the equation must be continuous at ``\Delta \overline{D}_r = 0`` as shown in Fig. 2 of the paper. -### Terminal velocity - -For the two moment scheme which is based on number density and mass, it is straightforward to model sedimentation of particles by using number- and mass-weighted mean terminal velocities. For rain water these terminal velocities are obtained by calculating the following integral: -```math -\begin{equation} - \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, ``k=0`` for number density and ``k=1`` for mass. -The individual terminal velocity of particles is approximated by -```math -v(x) = \left(\rho_0/\rho\right)^{\frac{1}{2}} [a_R - b_R exp(-c_R D_r)] -``` -where ``a_R``, ``b_R`` and ``c_R`` are three free parameters and ``D_r`` is the particle diameter. -Evaluating the integral results in the following equation for terminal velocity: -```math -\begin{equation} - \overline{v}_{r,\, k} = \left(\frac{\rho_0}{\rho}\right)^{\frac{1}{2}}\left[a_R - b_R \left(1+\frac{c_R}{\lambda_r}\right)^{-(3k+1)}\right], -\label{eq:SBTerminalVelocity} -\end{equation} -``` -where ``\lambda_r`` is the raindrops size distribution parameter (based on diameter): ``\lambda_r = (\phi \rho_w/\overline{x}_r)^{1/3}``. To avoid numerical instabilities, especially when ``N_{rai} \rightarrow 0`` and ``q_{rai} \rightarrow 0``, ``\lambda_r`` is bounded. The limiting algorithm is as follows: -```math -\begin{align} -\widetilde{x}_r &= max \left(\overline{x}_{r,\, min} , min \left(\overline{x}_{r,\, max} , \frac{\rho q_{rai}}{N_{rai}}\right)\right),\\ -N_0 &= max \left(N_{0,\, min} , min \left(N_{0,\, max} , N_{rai}\left(\frac{\pi \rho_w}{\widetilde{x}_r}\right)^{\frac{1}{3}}\right)\right),\\ -\lambda_r &= max \left(\lambda_{min} , min \left(\lambda_{max} , \left(\frac{\pi \rho_w N_0}{\rho q_{rai}}\right)^{\frac{1}{4}}\right)\right). -\end{align} -``` - -When the limiting algorithm is not applied, the terminal velocity given by eq. (\ref{eq:SBTerminalVelocity}) can become negative for small mean radius values. This occurs because the equation for the individual particle's terminal velocity may predict negative values for small particles. To avoid negative terminal velocities (or a sudden transition to zero velocity if we return zero instead of negative values), we need to adjust the integration bounds. Specifically, the integrals should be evaluated from the radius at which the individual terminal velocity is zero to infinity. This adjustment leads to the following equation for the terminal velocity: - -```math -\begin{align} - \overline{v}_{r,\, k} &= \frac{1}{M_r^k} \int_{r_c}^\infty x^k f_r(x) v(x) dx \nonumber\\ - &= \left(\frac{\rho_0}{\rho}\right)^{\frac{1}{2}}\left[a_R \frac{\Gamma(3k+1, 2r_c\lambda_r)}{\Gamma(3k+1)} - b_R \frac{\Gamma(3k+1, 2r_c(\lambda_r + c_R))}{\Gamma(3k+1)} \left(1+\frac{c_R}{\lambda_r}\right)^{-(3k+1)}\right], -\label{eq:SBModifiedTerminalVelocity} -\end{align} -``` -where ``r_c = -ln(a_R / b_R)/(2c_R)``. - -The default free parameter values are: - -| symbol | default value | -|----------------------------|-------------------------------------| -|``a_R`` | ``9.65 \, m \cdot s^{-1}`` | -|``b_R`` | ``10.3 \, m \cdot s^{-1}`` | -|``c_R`` | ``600 \, m^{-1}`` | -|``\overline{x}_{r,\, min}`` | ``6.54 \times 10^{-11} \, m`` | -|``\overline{x}_{r,\, max}`` | ``5 \times 10^{-6} \, m`` | -|``N_{0,\, min}`` | ``3.5 \times 10^{5} \, m^{-4}`` | -|``N_{0,\, max}`` | ``2 \times 10^{10} \, m^{-4}`` | -|``\lambda_{min}`` | ``1 \times 10^{3} \, m^{-1}`` | -|``\lambda_{max}`` | ``4 \times 10^{4} \, m^{-1}`` | - ### Rain evaporation The parametrization of rain evaporation is obtained by considering the time scale of evaporation of individual raindrops: @@ -427,6 +372,87 @@ include("plots/RainEvapoartionSB2006.jl") ``` ![](SB2006_rain_evaporation.svg) +### Terminal velocity + +The number- and mass-weighted terminal velocities for rain water are obtained by calculating the following integral: +```math +\begin{equation} + \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. + +The individual terminal velocity of particles is approximated by +```math +v(D) = \left(\rho_0/\rho\right)^{\frac{1}{2}} [a_R - b_R exp(-c_R D)] +``` +where ``a_R``, ``b_R`` and ``c_R`` are free parameters and ``D`` is the particle diameter. +Evaluating the integral results in +```math +\begin{equation} + \overline{v}_{r,\, k} = \left(\frac{\rho_0}{\rho}\right)^{\frac{1}{2}}\left[a_R - b_R \left(1+\frac{c_R}{\lambda_r}\right)^{-(3k+1)}\right], +\label{eq:SBTerminalVelocity} +\end{equation} +``` +where ``\lambda_r`` is the raindrop size distribution parameter (based on diameter): ``\lambda_r = (\phi \rho_w/\overline{x}_r)^{1/3}``. To avoid numerical instabilities, especially when ``N_{rai} \rightarrow 0`` and ``q_{rai} \rightarrow 0``, ``\lambda_r`` is bounded. The limiting algorithm is as follows: +```math +\begin{align} +\widetilde{x}_r &= \max \left[ \overline{x}_{r,\, \min} , \min \left(\overline{x}_{r,\, \max} , \frac{\rho q_{rai}}{N_{rai}}\right)\right] ,\\ +N_0 &= \max \left[ N_{0,\, \min} , \min \left(N_{0,\, \max} , N_{rai}\left(\frac{\pi \rho_w}{\widetilde{x}_r}\right)^{\frac{1}{3}}\right)\right],\\ +\lambda_r &= \max \left[ \lambda_{\min} , \min \left(\lambda_{\max} , \left(\frac{\pi \rho_w N_0}{\rho q_{rai}}\right)^{\frac{1}{4}}\right)\right]. +\end{align} +``` + +When the limiting algorithm is not applied, the terminal velocity given by eq. (\ref{eq:SBTerminalVelocity}) can become negative for small mean radius values. This occurs because the equation for the individual particle's terminal velocity may predict negative values for small particles. To avoid negative terminal velocities (or a sudden transition to zero velocity if we return zero instead of negative values), we need to adjust the integration bounds. Specifically, the integrals should be evaluated from the radius at which the individual terminal velocity is zero to infinity. This adjustment leads to the following equation for the terminal velocity: +```math +\begin{align} + \overline{v}_{r,\, k} &= \frac{1}{M_r^k} \int_{r_c}^\infty x^k f_r(x) v(x) dx \nonumber\\ + &= \left(\frac{\rho_0}{\rho}\right)^{\frac{1}{2}}\left[a_R \frac{\Gamma(3k+1, 2r_c\lambda_r)}{\Gamma(3k+1)} - b_R \frac{\Gamma(3k+1, 2r_c(\lambda_r + c_R))}{\Gamma(3k+1)} \left(1+\frac{c_R}{\lambda_r}\right)^{-(3k+1)}\right], +\label{eq:SBModifiedTerminalVelocity} +\end{align} +``` +where ``r_c = -\ln(a_R / b_R)/(2c_R)``. + +The default parameter values are: + +| symbol | default value | +|----------------------------|-------------------------------------| +|``a_R`` | ``9.65 \, m \cdot s^{-1}`` | +|``b_R`` | ``10.3 \, m \cdot s^{-1}`` | +|``c_R`` | ``600 \, m^{-1}`` | +|``\overline{x}_{r,\, \min}`` | ``6.54 \times 10^{-11} \, m`` | +|``\overline{x}_{r,\, \max}`` | ``5 \times 10^{-6} \, m`` | +|``N_{0,\, \min}`` | ``3.5 \times 10^{5} \, m^{-4}`` | +|``N_{0,\, \max}`` | ``2 \times 10^{10} \, m^{-4}`` | +|``\lambda_{\min}`` | ``1 \times 10^{3} \, m^{-1}`` | +|``\lambda_{\max}`` | ``4 \times 10^{4} \, m^{-1}`` | + +Below we compare number-weighted (left) and mass-weighted (right) terminal velocities + for the original parameterization of SB2006 [SeifertBeheng2006](@cite) + and the modifed parameterzation given by eq. (\ref{eq:SBModifiedTerminalVelocity}), + without the limiting distribution shape factor $\lambda_r$. +```@example +include("plots/TerminalVelocity2M.jl") +``` +![](2M_terminal_velocity_comparisons.svg) + +For the Chen et al. [Chen2022](@cite) terminal velocity parameterization, + the number- and mass-weighted group terminal velocities are: +```math +\begin{equation} + \overline{v_k} = \frac{\int_0^\infty v(D) \, D^k \, n(D) \, dD} + {\int_0^\infty D^k \, n(D) \, dD} + = (\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 gamma function, $\delta = k + 1$, and +$\lambda$ is the size distribution parameter. +The velocity $\overline{v_k}$ is the number-weighted mean terminal velocity when k = 0, +and the mass-weighted mean terminal velocity when k = 3. + +## Diagnostics + ### Radar reflectivity The radar reflectivity factor (``Z``) is used to measure the power returned by a radar signal when it encounters atmospheric particles (cloud and rain droplets), and it is defined as the sixth moment of the particles distributions (``n(r)``) : @@ -509,61 +535,8 @@ where: ## Additional 2-moment microphysics options -### Terminal Velocity - -[Chen2022](@cite) provides a terminal velocity parameterisation -based on an empirical fit to a high accuracy model. -It consideres the deformation effects of large rain drops, -as well as size-specific air density dependence. -The fall speed of individual raindrops $v(D)$ is parameterized as: -```math -\begin{equation} - v(D) = (\phi)^{\kappa} \Sigma_{i = 1,3} \; a_i D^{b_i} e^{-c_i*D} -\end{equation} -``` -where D is the diameter of the particle, -$a_i$, $b_i$, and $c_i$ are free parameers that account for deformation at larger sizes and air density dependance, -$\phi$ is the aspect ratio (assumed to be 1 for spherical droplets), and -$\kappa$ is 0 (corresponding to a spherical raindrop). -$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$ refers to air density. -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$ | - -Assuming the same size distribution as in [SeifertBeheng2006](@cite), -the number- and mass-weighted mean terminal velocities are: -```math -\begin{equation} - \overline{v_k} = \frac{\int_0^\infty v(D) \, D^k \, n(D) \, dD} - {\int_0^\infty D^k \, n(D) \, dD} - = (\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 gamma function, $\delta = k + 1$, and -$\lambda$ is the size distribution parameter. -$\overline{v_k}$ corresponds to the number-weighted mean terminal velocity when k = 0, -and mass-weighted mean terminal velocity when k = 3. - -Below, in the top-left panel we compare the individual terminal velocity formulas for Chen2022 [Chen2022](@cite), SB2006 [SeifertBeheng2006](@cite) and data from [Gunn1949](@cite). -In the top-right panel, we compare bulk number weighted [ND] and mass weighted [M] -terminal velocities for both formulas integrated over the size distribution from SB2006 [SeifertBeheng2006](@cite). -We also show the mass weighted terminal velocity from the 1-moment scheme. In the bottom panels we compare number-weighted (left) and mass-weighted (right) terminal velocities for the original parameterization of SB2006 [SeifertBeheng2006](@cite) and the modifed parameterzation given by eq. (\ref{eq:SBModifiedTerminalVelocity}) without limiting distribution shape factor $\lambda_r$. -```@example -include("plots/TerminalVelocityComparisons.jl") -``` -![](2M_terminal_velocity_comparisons.svg) - -### Accretion and Autoconversion - -The other autoconversion and accretion rates in the `Microphysics2M.jl` module are implemented after Table 1 from [Wood2005](@cite) - and are based on the works of +The other autoconversion and accretion rates in the `Microphysics2M.jl` module + are implemented after Table 1 from [Wood2005](@cite) and are based on the works of [KhairoutdinovKogan2000](@cite), [Beheng1994](@cite), [TripoliCotton1980](@cite) and diff --git a/docs/src/MicrophysicsNonEq.md b/docs/src/MicrophysicsNonEq.md index 22603d973..acd9bf95d 100644 --- a/docs/src/MicrophysicsNonEq.md +++ b/docs/src/MicrophysicsNonEq.md @@ -92,26 +92,62 @@ where: - ``c_p`` is the specific heat of air at constant pressure, - ``R_v`` is the gas constant of water vapor, - ``L_v`` and ``L_s`` is the latent heat of vaporization and sublimation. - - Note that these forms of condensation/sublimation and deposition/sublimation are equivalent to those described in the adiabatic parcel model with some rearrangements and assumptions. It is just necessary to use the definitions of ```\tau```, ```q_{sl}```, and the thermal diffusivity ```D_v```: - - ```math - \begin{equation} - \tau = 4 \pi N_{tot} \bar{r}, \;\;\;\;\;\;\;\; - q_{sl} = \frac{e_{sl}}{\rho R_v T}, \;\;\;\;\;\;\;\; - D_v = \frac{K}{\rho c_p} - \end{equation} - ``` - and if we assume that the supersaturation S can be approximated by specific humidities (this is only exactly true for mass mixing ratios): - ```math - \begin{equation} - S_l = \frac{q_{vap}}{q_{sl}} - \end{equation} - ``` - then we can write: - ```math - \begin{equation} - q_{vap} - q_{sl} = q_{sl}(S_l - 1) - \end{equation} - ``` - and ```Gamma``` is equivalent to the ```G``` function used in our parcel and parameterizations. \ No newline at end of file + +Note that these forms of condensation/sublimation and deposition/sublimation + are equivalent to those described in the adiabatic parcel model with some rearrangements and assumptions. +To see this, it is necessary to use the definitions of ``\tau``, ``q_{sl}``, and the thermal diffusivity ``D_v``: + +```math +\begin{equation} + \tau = 4 \pi N_{tot} \bar{r}, \;\;\;\;\;\;\;\; + q_{sl} = \frac{e_{sl}}{\rho R_v T}, \;\;\;\;\;\;\;\; + D_v = \frac{K}{\rho c_p}. +\end{equation} +``` +If we then assume that the supersaturation ``S`` can be approximated by the specific humidities (this is only exactly true for mass mixing ratios): +```math +\begin{equation} + S_l = \frac{q_{vap}}{q_{sl}}, +\end{equation} +``` +we can write +```math +\begin{equation} + q_{vap} - q_{sl} = q_{sl}(S_l - 1). +\end{equation} +``` +``\Gamma_l`` and ``\Gamma_i`` then are equivalent to the ``G`` function used in our parcel model and parameterizations. + +## Cloud condensate sedimentation + +We use the Chen et al. [Chen2022](@cite) parameterization for cloud liquid and cloud ice sedimentation velocities. +In the 1-moment precipitation scheme, we assume that cloud condensate is a continuous field + and don't introduce an explicit particle size distribution. +For simplicity, we assume a monodisperse size distribution + and compute the group terminal velocity based on the volume radius + and prescribed number concentration: + +```math +\begin{equation} + D_{vol} = \frac{\rho_{air} q}{N \rho} +\end{equation} +``` +where: + - ``\rho_{air}`` is the air density, + - ``q`` is the cloud liquid or cloud ice specific humidity, + - ``N`` is the prescribed number concentration (``500/cm^3`` by default), + - ``\rho`` is the cloud water or cloud ice density. + +The sedimentation velocity then is +```math +\begin{equation} + v_t = v_{term}(D_{vol}). +\end{equation} +``` + +!!! note + We are using the B1 coefficients from Chen et al. [Chen2022](@cite) to compute + the cloud condensate velocities. They were fitted for larger particle sizes. + To mitigate the resulting errors, we multiply by a correction factor. + We should instead find a parameterization that was designed for the cloud droplet + size range. diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index 1f03e4b37..a760dff3f 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -182,24 +182,12 @@ include("plots/P3LambdaErrorPlots.jl") ## Terminal Velocity -We use the [Chen2022](@cite) velocity parameterization: -```math -V(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D} -``` -where -```math -\phi = (16 \rho_{ice}(D)^2 A(D)^3) / (9 \pi m(D)^2) -``` -is the aspect ratio, and ``\kappa``, ``a_i``, ``b_i`` and ``c_i`` are the free parameters. - -The aspect ratio of a spheroid is defined ``\phi = \frac{a}{c}``, where ``a`` is the equatorial radius and ``c`` is the distance from the pole to the center. - In terms of ``a`` and ``c``, a spheroid's volume can be represented as ``V(a, c) = \frac{4}{3} \pi a^2 c``, and its cross-sectional area can be assumed ``A(a, c) = \pi a c``. - We use ``m(D)`` and ``A(D)`` from P3, so by substituting ``m(a, c) = \rho_{ice} V(a, c)``, ``A(a, c)`` for ``m(D)``, ``A(D)`` into the formulation of aspect ratio above, - we demonstrate agreement with the definition ``\phi = \frac{a}{c}``. +We use the Chen et al. [Chen2022](@cite) velocity parameterization, + see [here](https://clima.github.io/CloudMicrophysics.jl/dev/TerminalVelocity.html#Chen-et.-al.-2022) for details. +The figure below shows the implied aspect ratio for different particle size regimes. Note that ``\phi = 1`` corresponds to spherical particles (small spherical ice (``D < D_{th}``) and graupel (``D_{gr} < D < D_{cr}``)). -The plot provided below helps to visualize the transitions between spherical and nonspherical regimes. ```@example include("plots/P3AspectRatioPlot.jl") ``` diff --git a/docs/src/TerminalVelocity.md b/docs/src/TerminalVelocity.md new file mode 100644 index 000000000..4d61becc8 --- /dev/null +++ b/docs/src/TerminalVelocity.md @@ -0,0 +1,145 @@ +# Terminal velocity + +`CloudMicrophysics.jl` offers three parameterizations of the relationship between +particle size and terminal velocity: + - A simple power-law used in the 1-moment microphysics scheme, + - 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). + +The above terminal velocities need to be averaged over the assumed + particle size distributions, to get the mass- or number-weighted + group terminal velocities used in the bulks schemes: + - In the non-equilibrium microphysics, we use the Chen et al. [Chen2022](@cite) parameterization and assume a + monodisperse size distribution to obtain the cloud liquid and cloud ice sedimentation velocities. + - In 1M microphysics, we assume an exponential size distribution for rain and snow + 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. + - 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 + for more details. + +## Power-law relationship + +We assume that rain drop and snow terminal velocities can be parameterized as +```math +v_{term}(r) = \chi_v \, v_0 \left(\frac{r}{r_0}\right)^{v_e + \Delta_v} +``` +where: + - ``r`` is the particle radius, + - ``r_0`` is the typical particle radius used to nondimensionalize, + - ``v_0, \, v_e \,`` are the default coefficients, + - ``\chi_v``, ``\Delta_v`` + are the coefficients that can be used during model calibration to expand + around the default values. + Without calibration all ``\chi`` parameters are set to 1 + and all ``\Delta`` parameters are set to 0. + +The rain drop terminal velocity coefficients are obtained from the + balance between the gravitational acceleration (taking into account + the density difference between water and air) and the drag force: +```math +\begin{equation} +v_0^{rai} = + \left( + \frac{8}{3 \, C_{drag}} \left( \frac{\rho_{w}}{\rho} -1 \right) + \right)^{1/2} (g r_0^{rai})^{1/2} +\label{eq:vdrop} +\end{equation} +``` +where: + - ``g`` is the gravitational acceleration, + - ``C_{drag}`` is the drag coefficient, + - ``\rho`` is the density of air, + - ``\rho_w`` is the density of water. + +Assuming a constant drag coefficient for rain drops is a big approximation, + as it should be size and flow + [dependent](https://www1.grc.nasa.gov/beginners-guide-to-aeronautics/drag-of-a-sphere/). +The coefficients for the snow terminal velocity are empirical. + +| symbol | definition | units | default value | reference | +|---------------|-----------------------------------------|-----------------|-----------------------|-----------| +|``C_{drag}`` | rain drop drag coefficient | - | ``0.55`` | ``C_{drag}`` is such that the mass averaged terminal velocity is close to [Grabowski1996](@cite) | +|``r_0^{rai}`` | typical rain drop radius | ``m`` | ``10^{-3} `` | | +|``v_e^{rai}`` | exponent in ``v_{term}(r)`` for rain | - | ``0.5`` | | +|``r_0^{sno}`` | typical snow crystal radius | ``m`` | ``10^{-3} `` | | +|``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) | + +## Seifert and Beheng 2006 + +Seifert and Beheng [SeifertBeheng2006](@cite) uses an empirical relationship between the rain drop size + and its terminal velocity +```math +v_{term}(D) = \left(\frac{\rho_0}{\rho}\right)^{\frac{1}{2}} [a_R - b_R e^{-c_R D}] +``` +where: + - ``a_R``, ``b_R`` and ``c_R`` are free parameters, + - ``D`` is the rain drop diameter. + +The default parameter values are + +| symbol | default value | units | +|----------------------------|---------------|-----------------| +|``a_R`` | ``9.65`` | ``\frac{m}{s}`` | +|``b_R`` | ``10.3`` | ``\frac{m}{s}`` | +|``c_R`` | ``600`` | ``\frac{1}{m}`` | + +## Chen et. al. 2022 + +Chen et al. [Chen2022](@cite) provide a terminal velocity parameterization + based on an empirical fit to a high-accuracy model. +The terminal velocity depends on particle shape, size, and density, + and it considers the deformation of large rain drops, + as well as a size-specific air density dependence. +The fall speed of individual particles 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, ``j=2``, and for rain, ``j=3``, to account for deformation at larger drop sizes. +The rain parameterization is valid for ``D > 100 \; \mu m``. +The free parameter values for solid particles are different for small and large sizes, + with a cutoff diameter ``D = 625 \; \mu m``. + +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. +In the non-equilibrium cloud microphysics scheme, we assume that both cloud + droplets and ice crystals are spherical and ``\phi=1``. +In all bulk schemes, we assume that rain drops are spherical and ``\phi=1``. +In the 1-moment bulk scheme for snow and in the P3 scheme, we assume that ``\kappa = -1/3``. +We compute the aspect ratio based on the assumed ``m(r)`` and ``a(r)`` relationships, + assuming that the particle volume can be represented as ``V_p = 4\pi/3 \; a^2 c``, + and particle cross section area as ``A_p = \pi a c``. + +## Example figures + +In the top row of the figure below, we show different particle terminal velocities + in cloud (left panel) and precipitation (middle panel) size range. +In the top right panel, we also show the snow particle aspect ratio. +In the bottom row, we show the group sedimentation velocities for cloud condensate (left panel) + and precipitation (right panel). + +```@example +include("plots/TerminalVelocity.jl") +``` +![](1M_individual_terminal_velocity_comparisons.svg) diff --git a/docs/src/plots/Microphysics1M_plots.jl b/docs/src/plots/Microphysics1M_plots.jl index 5aee2ecb2..73fcaf4cf 100644 --- a/docs/src/plots/Microphysics1M_plots.jl +++ b/docs/src/plots/Microphysics1M_plots.jl @@ -22,18 +22,6 @@ const Chen2022 = CMP.Chen2022VelType(FT) const Blk1MVel = CMP.Blk1MVelType(FT) const ce = CMP.CollisionEff(FT) -# eq. 5d in [Grabowski1996](@cite) -function terminal_velocity_empirical( - q_rai::DT, - q_tot::DT, - ρ::DT, - ρ_air_ground::DT, -) where {DT <: Real} - rr = q_rai / (DT(1) - q_tot) - vel = DT(14.34) * ρ_air_ground^DT(0.5) * ρ^-DT(0.3654) * rr^DT(0.1346) - return vel -end - # eq. 5b in [Grabowski1996](@cite) function accretion_empirical(q_rai::DT, q_liq::DT, q_tot::DT) where {DT <: Real} rr = q_rai / (DT(1) - q_tot) @@ -68,26 +56,6 @@ function rain_evap_empirical( return DT(1) / (DT(1) - q.tot) * S * F * G end -function terminal_velocity_chen_large( - (; pdf, mass)::CMP.CloudIce{FT}, - vel::CMP.Chen2022VelTypeSnowIce{FT}, - ρ::FT, - q::FT, -) where {FT} - fall_w = FT(0) - if q > FT(0) - # coefficients from Table B1 from Chen et. al. 2022 - aiu, bi, ciu = CO.Chen2022_vel_coeffs_large(vel, ρ) - # size distribution parameter - λ::FT = CM1.lambda(pdf, mass, q, ρ) - # eq 20 from Chen et al 2022 - fall_w = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 3)) - # It should be ϕ^κ * fall_w, but for rain drops ϕ = 1 and κ = 0 - fall_w = max(FT(0), fall_w) - end - return fall_w -end - # example values q_liq_range = range(1e-8, stop = 5e-3, length = 100) q_ice_range = range(1e-8, stop = 5e-3, length = 100) @@ -96,84 +64,6 @@ q_snow_range = range(1e-8, stop = 5e-3, length = 100) ρ_air, ρ_air_ground = 1.2, 1.22 q_liq, q_ice, q_tot = 5e-4, 5e-4, 20e-3 -PL.plot( - q_rain_range * 1e3, - [ - CM1.terminal_velocity(rain, Blk1MVel.rain, ρ_air, q_rai) for - q_rai in q_rain_range - ], - linewidth = 3, - label = "Rain-1M-default", - color = :blue, - xlabel = "q_rain/ice/snow [g/kg]", - ylabel = "terminal velocity [m/s]", -) -PL.plot!( - q_rain_range * 1e3, - [ - CM1.terminal_velocity(rain, Chen2022.rain, ρ_air, q_rai) for - q_rai in q_rain_range - ], - linewidth = 3, - label = "Rain-Chen", - color = :blue, - linestyle = :dot, -) -PL.plot!( - q_rain_range * 1e3, - [ - terminal_velocity_empirical(q_rai, q_tot, ρ_air, ρ_air_ground) for - q_rai in q_rain_range - ], - linewidth = 3, - label = "Rain-Empirical", - color = :green, -) -PL.plot!( - q_ice_range * 1e3, - [ - CM1.terminal_velocity(ice, Chen2022.snow_ice, ρ_air, q_ice) for - q_ice in q_ice_range - ], - linewidth = 3, - label = "Ice-Chen", - color = :pink, - linestyle = :dot, -) -PL.plot!( - q_ice_range * 1e3, - [ - terminal_velocity_chen_large(ice, Chen2022.snow_ice, ρ_air, q_ice) for - q_ice in q_ice_range - ], - linewidth = 3, - label = "Ice-Chen-Large", - color = :lightpink1, - linestyle = :dot, -) -PL.plot!( - q_snow_range * 1e3, - [ - CM1.terminal_velocity(snow, Blk1MVel.snow, ρ_air, q_sno) for - q_sno in q_snow_range - ], - linewidth = 3, - label = "Snow-1M-default", - color = :red, -) -PL.plot!( - q_snow_range * 1e3, - [ - CM1.terminal_velocity(snow, Chen2022.snow_ice, ρ_air, q_sno) for - q_sno in q_snow_range - ], - linewidth = 3, - label = "Snow-Chen", - color = :red, - linestyle = :dot, -) -PL.savefig("terminal_velocity.svg") # hide - T = 273.15 PL.plot( q_liq_range * 1e3, diff --git a/docs/src/plots/TerminalVelocity.jl b/docs/src/plots/TerminalVelocity.jl new file mode 100644 index 000000000..9e0f83bf6 --- /dev/null +++ b/docs/src/plots/TerminalVelocity.jl @@ -0,0 +1,199 @@ +import Plots as PL +using Measures + +import ClimaParams as CP + +FT = Float64 + +import CloudMicrophysics.MicrophysicsNonEq as CMNe +import CloudMicrophysics.Microphysics1M as CM1 +import CloudMicrophysics.Microphysics2M as CM2 +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.Common as CMO + +const rain = CMP.Rain(FT) +const liquid = CMP.CloudLiquid(FT) +const ice = CMP.CloudIce(FT) +const snow = CMP.Snow(FT) + +const SB2006 = CMP.SB2006(FT) +const SB2006_no_lim = CMP.SB2006(FT, false) + +const Chen2022 = CMP.Chen2022VelType(FT) +const SB2006Vel = CMP.SB2006VelType(FT) +const Blk1MVel = CMP.Blk1MVelType(FT) + +function aspect_ratio_snow_1M(snow::CMP.Snow, D::FT) where {FT <: Real} + (; r0, m0, me, χm, Δm) = snow.mass + (; a0, ae, Δa, χa) = snow.area + ρᵢ = snow.ρᵢ + + aᵢ = χa * a0 * (D / 2 / r0)^(ae + Δa) + mᵢ = χm * m0 * (D / 2 / r0)^(me + Δm) + + return 16 * ρᵢ^2 * aᵢ^3 / (9 * FT(π) * mᵢ^2) +end + +""" + rain_terminal_velocity_individual_Chen(velo_scheme ρ, D_r) + + - `velo_scheme` - structs with free parameters + - `ρ` - air density + - `D_r` - diameter of the raindrops + +Returns the fall velocity of a raindrop from Chen et al 2022 +""" +function rain_terminal_velocity_individual_Chen( + velo_scheme::CMP.Chen2022VelTypeRain, + ρ::FT, + D_r::FT, #in m +) where {FT <: Real} + ai, bi, ci = CMO.Chen2022_vel_coeffs_B1(velo_scheme, ρ) + + v = 0 + for i in 1:3 + v += (ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i])) + end + return v +end + +function ice_terminal_velocity_individual_Chen( + velo_scheme::CMP.Chen2022VelTypeSnowIce, + ρ::FT, + D_r::FT, #in m +) where {FT <: Real} + ai, bi, ci = CMO.Chen2022_vel_coeffs_B2(velo_scheme, ρ) + + v = 0 + for i in 1:2 + v += (ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i])) + end + return v +end + +function snow_terminal_velocity_individual_Chen( + snow::CMP.Snow, + velo_scheme::CMP.Chen2022VelTypeSnowIce, + ρ::FT, + D_r::FT, #in m +) where {FT <: Real} + ai, bi, ci = CMO.Chen2022_vel_coeffs_B4(velo_scheme, ρ) + + ϕ = aspect_ratio_snow_1M(snow, D_r) + + v = 0 + for i in 1:2 + v += ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i]) + end + v_term = ifelse(ϕ > 0, ϕ^(-1 / 3) * v, ϕ^(1 / 6) * v) + return max(FT(0), v_term) +end + +""" + rain_terminal_velocity_individual_SB(param_set, ρ, D_r) + + - `param_set` - set with free parameters + - `ρ` - air density + - `D_r` - diameter of the raindrops + +Returns the fall velocity of a raindrop from Seifert and Beheng 2006 +""" +function rain_terminal_velocity_individual_SB( + (; ρ0, aR, bR, cR)::CMP.SB2006VelType, + ρ::FT, + D_r::FT, +) where {FT <: Real} + return v = (ρ0 / ρ)^(1 / 2) * (aR - bR * exp(-cR * D_r)) +end + +""" + terminal_velocity_individual_1M(velo_scheme, ρ, D_r) + + - `velo_scheme` - set with free parameters + - `ρ` - air density + - `D_r` - particle diameter + +Returns the fall velocity of a raindrop or snow from 1-moment scheme +""" +function terminal_velocity_individual_1M( + velo_scheme::Union{CMP.Blk1MVelTypeRain, CMP.Blk1MVelTypeSnow}, + ρ::FT, + D_r::FT, +) where {FT <: Real} + (; χv, ve, Δv, r0) = velo_scheme + v0 = CM1.get_v0(velo_scheme, ρ) + vt = χv * v0 * (D_r / (2 * r0))^(Δv + ve) + return vt +end + +ρ_air = 1.2 +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) + +#! format: off +# velocity values for cloud particle sizes +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] +Ch_liq_small = [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_air, D_r) for D_r in D_r_range_small] +Ch_ice_small = [ice_terminal_velocity_individual_Chen(Chen2022.snow_ice, ρ_air, D_r) for D_r in D_r_range_small] +Ch_rain_small = [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_air, D_r) for D_r in D_r_range_small] +Ch_snow_small = [snow_terminal_velocity_individual_Chen(snow, Chen2022.snow_ice, ρ_air, D_r) for D_r in D_r_range_small] +# velocity values for precip particle sizes +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] +Ch_liq = [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_air, D_r) for D_r in D_r_range] +Ch_ice = [ice_terminal_velocity_individual_Chen(Chen2022.snow_ice, ρ_air, D_r) for D_r in D_r_range] +Ch_rain = [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_air, D_r) for D_r in D_r_range] +Ch_snow = [snow_terminal_velocity_individual_Chen(snow, Chen2022.snow_ice, ρ_air, D_r) for D_r in D_r_range] +# obs data +D_Gunn_Kinzer = [0.0, 0.078, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8] .* 1e-3 +u_Gunn_Kinzer = [0.0, 18.0, 27, 72, 117, 162, 206, 247, 287, 327, 367, 403, 464, 517, 565, 609, 649, 690, 727, 757, 782, 806, 826, 844, 860, 872, 883, 892, 898, 903, 907, 909, 912, 914, 916, 917] ./ 100 +D_Gunn_Kinzer_small = [0.0, 0.078, 0.1] .* 1e-3 +u_Gunn_Kinzer_small = [0.0, 18.0, 27] ./ 100 +# aspect ratio plot +Aspect_Ratio = [aspect_ratio_snow_1M(snow, D_r) for D_r in D_r_range] +# group velocity values +bM1_rain = [CM1.terminal_velocity(rain, Blk1MVel.rain, ρ_air, q) for q in q_range] +bM1_snow = [CM1.terminal_velocity(snow, Blk1MVel.snow, ρ_air, q) for q in q_range] +bCh_rain = [CM1.terminal_velocity(rain, Chen2022.rain, ρ_air, q) for q in q_range] +bCh_snow = [CM1.terminal_velocity(snow, Chen2022.snow_ice, ρ_air, q) 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.snow_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, 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) +p1 = PL.plot!(D_r_range_small * 1e6, Ch_ice_small * 1e2, linewidth = 3, xlabel = "D [um]", ylabel = "terminal velocity [cm/s]", label = "Ice Chen", color = :orange) +# individual particle comparison - precip size range +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, 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_ice, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", legend=false, color = :orange) +# snow aspect ratio +p3 = PL.plot(D_r_range * 1e3, Aspect_Ratio, linewidth = 3, xlabel = "D [mm]", ylabel = "aspect ratio", label = "Snow aspect ratio", color = :darkviolet) + +# 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") + +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") +p5 = PL.plot!(q_range * 1e3, bCh_rain, linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [m/s]", color = :blue, label="Rain Chen") +p5 = PL.plot!(q_range * 1e3, bCh_snow, linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [m/s]", color = :darkviolet, label="Snow Chen") +p5 = PL.plot!(q_range * 1e3, bCh_liq, linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [m/s]", color = :green, label="Liq Chen") +p5 = PL.plot!(q_range * 1e3, bCh_ice, linewidth = 3, xlabel = "q [g/kg]", ylabel = "terminal velocity [m/s]", color = :orange, label="Ice Chen") + +PL.plot(p1, p2, p3, p4, p5, layout = (2, 3), size = (1200, 750), dpi=400) +PL.savefig("1M_individual_terminal_velocity_comparisons.svg") + +#! format: on diff --git a/docs/src/plots/TerminalVelocity2M.jl b/docs/src/plots/TerminalVelocity2M.jl new file mode 100644 index 000000000..c2cc6b602 --- /dev/null +++ b/docs/src/plots/TerminalVelocity2M.jl @@ -0,0 +1,43 @@ +import Plots as PL +using Measures + +import ClimaParams as CP +import CloudMicrophysics.Microphysics2M as CM2 +import CloudMicrophysics.Parameters as CMP + +FT = Float64 + +const rain = CMP.Rain(FT) +const SB2006 = CMP.SB2006(FT) +const SB2006_no_lim = CMP.SB2006(FT, false) +const SB2006Vel = CMP.SB2006VelType(FT) + +ρ_air = 1.2 +# random q and N values for plotting vt vs. mean radius +q_rain_random = rand(10000) .* 5e-3 +N_rain_random = 10.0 .^ (1.0 .+ 6.0 .* rand(10000)) +r_mean = + ((ρ_air .* q_rain_random ./ N_rain_random) / 1000.0 * 3 / 4 / pi) .^ + (1.0 / 3) .* 1e6 + +#! format: off + +# SB2006 terminal velocities for random combinations of q and N +SB_rain_bN_rm = [CM2.rain_terminal_velocity(SB2006, SB2006Vel, q_rain_random[i], ρ_air, N_rain_random[i])[1] for i in 1:length(q_rain_random)] +SB_rain_bM_rm = [CM2.rain_terminal_velocity(SB2006, SB2006Vel, q_rain_random[i], ρ_air, N_rain_random[i])[2] for i in 1:length(q_rain_random)] +SB_rain_bN_rm_nolim = [CM2.rain_terminal_velocity(SB2006_no_lim, SB2006Vel, q_rain_random[i], ρ_air, N_rain_random[i])[1] for i in 1:length(q_rain_random)] +SB_rain_bM_rm_nolim = [CM2.rain_terminal_velocity(SB2006_no_lim, SB2006Vel, q_rain_random[i], ρ_air, N_rain_random[i])[2] for i in 1:length(q_rain_random)] + +# SB2006 group velocities vs. mean radius +p3 = PL.scatter( r_mean./1000, SB_rain_bN_rm, ms = 1.0, markerstrokewidth=0, xlabel = "r_mean [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [ND]", margin=10mm) +p3 = PL.scatter!(r_mean./1000, SB_rain_bN_rm_nolim, ms = 1.0, markerstrokewidth=0, xlabel = "r_mean [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [ND], no lim") +p3 = PL.plot!(xlim = [0, 2], ylim = [0, 5.5]) + +p4 = PL.scatter( r_mean./1000, SB_rain_bM_rm, ms = 1.0, markerstrokewidth=0, xlabel = "r_mean [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [M]") +p4 = PL.scatter!(r_mean./1000, SB_rain_bM_rm_nolim, ms = 1.0, markerstrokewidth=0, xlabel = "r_mean [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [M], no lim") +p4 = PL.plot!(xlim = [0, 2], ylim = [0, 10]) +# save plot +PL.plot(p3, p4, layout = (1, 2), size = (900, 450), dpi = 300) +PL.savefig("2M_terminal_velocity_comparisons.svg") + +#! format: on diff --git a/docs/src/plots/TerminalVelocityComparisons.jl b/docs/src/plots/TerminalVelocityComparisons.jl deleted file mode 100644 index ec061d156..000000000 --- a/docs/src/plots/TerminalVelocityComparisons.jl +++ /dev/null @@ -1,280 +0,0 @@ -import Plots as PL - -import ClimaParams as CP - -FT = Float64 - -import CloudMicrophysics.Microphysics1M as CM1 -import CloudMicrophysics.Microphysics2M as CM2 -import CloudMicrophysics.Parameters as CMP -import CloudMicrophysics.Common as CMO - -const rain = CMP.Rain(FT) -const liquid = CMP.CloudLiquid(FT) -const ice = CMP.CloudIce(FT) -const snow = CMP.Snow(FT) - -const SB2006 = CMP.SB2006(FT) -const SB2006_no_lim = CMP.SB2006(FT, false) - -const Chen2022 = CMP.Chen2022VelType(FT) -const SB2006Vel = CMP.SB2006VelType(FT) -const Blk1MVel = CMP.Blk1MVelType(FT) - -""" - rain_terminal_velocity_individual_Chen(velo_scheme ρ, D_r) - - - `velo_scheme` - structs with free parameters - - `ρ` - air density - - `D_r` - diameter of the raindrops - -Returns the fall velocity of a raindrop from Chen et al 2022 -""" -function rain_terminal_velocity_individual_Chen( - velo_scheme::CMP.Chen2022VelTypeRain, - ρ::FT, - D_r::FT, #in m -) where {FT <: Real} - ai, bi, ci = CMO.Chen2022_vel_coeffs_small(velo_scheme, ρ) - - v = 0 - for i in 1:3 - v += (ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i])) - end - return v -end - -""" - rain_terminal_velocity_individual_SB(param_set, ρ, D_r) - - - `param_set` - set with free parameters - - `ρ` - air density - - `D_r` - diameter of the raindrops - -Returns the fall velocity of a raindrop from Seifert and Beheng 2006 -""" -function rain_terminal_velocity_individual_SB( - (; ρ0, aR, bR, cR)::CMP.SB2006VelType, - ρ::FT, - D_r::FT, -) where {FT <: Real} - return v = (ρ0 / ρ)^(1 / 2) * (aR - bR * exp(-cR * D_r)) -end - -""" - terminal_velocity_individual_1M(velo_scheme, ρ, D_r) - - - `velo_scheme` - set with free parameters - - `ρ` - air density - - `D_r` - particle diameter - -Returns the fall velocity of a raindrop or snow from 1-moment scheme -""" -function terminal_velocity_individual_1M( - velo_scheme::Union{CMP.Blk1MVelTypeRain, CMP.Blk1MVelTypeSnow}, - ρ::FT, - D_r::FT, -) where {FT <: Real} - (; χv, ve, Δv, r0) = velo_scheme - v0 = CM1.get_v0(velo_scheme, ρ) - vt = χv * v0 * (D_r / (2 * r0))^(Δv + ve) - return vt -end - -""" - ice_terminal_velocity_individual_Chen(velo_scheme, ρ, D_r) - - - `velo_scheme` - set with free parameters - - `ρ` - air density - - `D_r` - diameter of the raindrops - -Returns the fall velocity of an ice particle from Chen et al 2022 -""" -function ice_terminal_velocity_individual_Chen( - velo_scheme::CMP.Chen2022VelTypeSnowIce, - ρ::FT, - D_r::FT, #in m -) where {FT <: Real} - - (; As, Bs, Cs, Es, Fs, Gs) = velo_scheme - - ai = [Es * ρ^As, Fs * ρ^As] - bi = [Bs + ρ * Cs, Bs + ρ * Cs] - ci = [0, Gs] - - D_r = D_r * 1000 #D_r is in mm in the paper --> multiply D_r by 1000 - v = 0 - for i in 1:2 - v += (ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i])) - end - return v -end - -""" - ice_terminal_velocity_large_individual_Chen(velo_scheme, ρ, D_r) - - - `velo_scheme` - set with free parameters - - `ρ` - air density - - `D_r` - diameter of the raindrops - -Returns the fall velocity of a large ice particle (D > 0.625mm) from Chen et al 2022 -""" -function ice_terminal_velocity_large_individual_Chen( - velo_scheme::CMP.Chen2022VelTypeSnowIce, - ρ::FT, - D_r::FT, #in m -) where {FT <: Real} - - (; Al, Bl, Cl, El, Fl, Gl, Hl) = velo_scheme - - ai = (Bl * ρ^Al, El * ρ^Al * exp(Hl * ρ)) - bi = (Cl, Fl) - ci = (FT(0), Gl) - - D_r = D_r * 1000 #D_r is in mm in the paper --> multiply D_r by 1000 - v = 0 - for i in 1:2 - v += (ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i])) - end - return v -end - -""" - snow_terminal_velocity_individual_Chen(precip, velo_scheme, ρ, D_r) - - - `precip`, `velo_scheme` - structs with free parameters - - `ρ` - air density - - `D_r` - particle diameter - -Returns the fall velocity of snow from Chen et al 2022 -""" -function snow_terminal_velocity_individual_Chen( - precip::CMP.Snow, - velo_scheme::CMP.Chen2022VelTypeSnowIce, - ρ::FT, - D_r::FT, -) where {FT <: Real} - (; r0, m0, me, Δm, χm) = precip.mass - (; a0, ae, Δa, χa) = precip.area - - D_r = D_r * 1000 - m0_comb = m0 * χm - a0_comb = a0 * χa - me_comb = me + Δm - ae_comb = ae + Δa - α = FT(-1 / 3) - - (; As, Bs, Cs, Es, Fs, Gs, ρᵢ) = velo_scheme - ai = [Es * ρ^As, Fs * ρ^As] - bi = [Bs + ρ * Cs, Bs + ρ * Cs] - ci = [0, Gs] - - aspect_ratio = - ((16 * (a0_comb)^3 * (ρᵢ)^2) / (9 * π * (m0_comb)^2)) * - (D_r / (2 * r0 * 1000))^(3 * ae_comb - 2 * me_comb) - aspect_ratio = aspect_ratio^(α) - vt = 0 - for i in 1:2 - vt = vt + aspect_ratio * ai[i] * ((D_r))^(bi[i]) * exp(-1 * ci[i] * D_r) - end - return vt -end - -function aspect_ratio_snow_1M(snow::CMP.Snow, D_r::FT) where {FT <: Real} - (; r0, m0, me, χm, Δm) = snow.mass - (; a0, ae, Δa, χa) = snow.area - ρᵢ = snow.ρᵢ - m0_comb = m0 * χm - a0_comb = a0 * χa - me_comb = me + Δm - ae_comb = ae + Δa - - aspect_ratio = - ((16 * (a0_comb)^3 * (ρᵢ)^2) / (9 * π * (m0_comb)^2)) * - (D_r / (2 * r0))^(3 * ae_comb - 2 * me_comb) - return aspect_ratio -end - -ρ_air, ρ_air_ground = 1.2, 1.22 -q_liq, q_ice, q_tot = 5e-4, 5e-4, 20e-3 -N_rai = 1e7 #this value matters for group terminal velocity -q_rain_range = range(1e-8, stop = 5e-3, length = 100) -D_r_range = range(1e-6, stop = 6e-3, length = 1000) -N_d_range = range(1e6, stop = 1e9, length = 1000) -D_r_ar_range = range(1e-6, stop = 6.25e-4, length = 1000) -# random q and N values for plotting vt vs. mean radius -q_rain_random = rand(10000) .* 5e-3 -N_rain_random = 10.0 .^ (1.0 .+ 6.0 .* rand(10000)) -r_mean = - ((ρ_air .* q_rain_random ./ N_rain_random) / 1000.0 * 3 / 4 / pi) .^ - (1.0 / 3) .* 1e6 - -#! format: off - -SB_rain_bN = [CM2.rain_terminal_velocity(SB2006, SB2006Vel, q_rai, ρ_air, N_rai)[1] for q_rai in q_rain_range] -Ch_rain_bN = [CM2.rain_terminal_velocity(SB2006, Chen2022.rain, q_rai, ρ_air, N_rai)[1] for q_rai in q_rain_range] -SB_rain_bM = [CM2.rain_terminal_velocity(SB2006, SB2006Vel, q_rai, ρ_air, N_rai)[2] for q_rai in q_rain_range] -Ch_rain_bM = [CM2.rain_terminal_velocity(SB2006, Chen2022.rain, q_rai, ρ_air, N_rai)[2] for q_rai in q_rain_range] -M1_rain_bM = [CM1.terminal_velocity(rain, Blk1MVel.rain, ρ_air, q_rai) for q_rai in q_rain_range] -# SB2006 terminal velocities for random combinations of q and N -SB_rain_bN_rm = [CM2.rain_terminal_velocity(SB2006, SB2006Vel, q_rain_random[i], ρ_air, N_rain_random[i])[1] for i in 1:length(q_rain_random)] -SB_rain_bM_rm = [CM2.rain_terminal_velocity(SB2006, SB2006Vel, q_rain_random[i], ρ_air, N_rain_random[i])[2] for i in 1:length(q_rain_random)] -SB_rain_bN_rm_nolim = [CM2.rain_terminal_velocity(SB2006_no_lim, SB2006Vel, q_rain_random[i], ρ_air, N_rain_random[i])[1] for i in 1:length(q_rain_random)] -SB_rain_bM_rm_nolim = [CM2.rain_terminal_velocity(SB2006_no_lim, SB2006Vel, q_rain_random[i], ρ_air, N_rain_random[i])[2] for i in 1:length(q_rain_random)] - -SB_rain = [rain_terminal_velocity_individual_SB(SB2006Vel, ρ_air, D_r) for D_r in D_r_range] -Ch_rain = [rain_terminal_velocity_individual_Chen(Chen2022.rain, ρ_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] -Ch_snow = [snow_terminal_velocity_individual_Chen(snow, Chen2022.snow_ice, ρ_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] -Ch_ice = [ice_terminal_velocity_individual_Chen(Chen2022.snow_ice, ρ_air, D_r) for D_r in D_r_range] -Ch_ice_large = [ice_terminal_velocity_large_individual_Chen(Chen2022.snow_ice, ρ_air, D_r) for D_r in D_r_range] -D_Gunn_Kinzer = [0.0, 0.078, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, - 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, - 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8] .* 1e-3 -u_Gunn_Kinzer = [0.0, 18.0, 27, 72, 117, 162, 206, 247, 287, 327, 367, 403, 464, 517, - 565, 609, 649, 690, 727, 757, 782, 806, 826, 844, 860, 872, 883, - 892, 898, 903, 907, 909, 912, 914, 916, 917] ./ 100 - -Aspect_Ratio = [aspect_ratio_snow_1M(snow, D_r) for D_r in D_r_ar_range] - -# 2 Moment Scheme Figures - -# single drop comparison -p1 = PL.plot(D_r_range *1e3, SB_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006", color = 1) -p1 = PL.plot!(D_r_range * 1e3, Ch_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen2022", color = 2) -p1 = PL.plot!(D_r_range * 1e3, M1_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-1M", color = 3) -p1 = PL.scatter!(D_Gunn_Kinzer * 1e3, u_Gunn_Kinzer, ms = 3, color = 4, label = "Gunn and Kinzer (1949)") -# group velocity comparison -p2 = PL.plot(q_rain_range * 1e3, SB_rain_bN, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [ND]", linestyle = :dot, color = 1) -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 = 2) -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 = 1) -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 = 2) -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 = 3) -# SB2006 group velocities vs. mean radius -p3 = PL.scatter(r_mean./1000, SB_rain_bN_rm, ms = 1.0, markerstrokewidth=0, xlabel = "r_mean [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [ND]") -p3 = PL.scatter!(r_mean./1000, SB_rain_bN_rm_nolim, ms = 1.0, markerstrokewidth=0, xlabel = "r_mean [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [ND], no lim") -p3 = PL.plot!(xlim = [0, 2], ylim = [0, 5.5]) -p4 = PL.scatter(r_mean./1000, SB_rain_bM_rm, ms = 1.0, markerstrokewidth=0, xlabel = "r_mean [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [M]") -p4 = PL.scatter!(r_mean./1000, SB_rain_bM_rm_nolim, ms = 1.0, markerstrokewidth=0, xlabel = "r_mean [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [M], no lim") -p4 = PL.plot!(xlim = [0, 2], ylim = [0, 10]) -# save plot -PL.plot(p1, p2, p3, p4, layout = (2, 2), size = (800, 800), dpi = 300) -PL.savefig("2M_terminal_velocity_comparisons.svg") - -# 1 Moment Scheme Figures - -# individual particle comparison -p1 = PL.plot(D_r_range * 1e3, Ch_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen", color = :blue, linestyle = :dot) -p1 = PL.plot!(D_r_range * 1e3, M1_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-default-1M", color = :blue) -p1 = 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) -p1 = PL.plot!(D_r_range * 1e3, M1_snow, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Snow-default-1M", color = :red) -p1 = PL.plot!(D_r_range * 1e3, Ch_ice, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Ice-Chen", color = :pink, linestyle = :dot) -p1 = PL.plot!(D_r_range * 1e3, Ch_ice_large, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Ice-large-Chen", color = :lightpink, linestyle = :dot) -# snow aspect ratio plot -p2 = PL.plot(D_r_ar_range * 1e3, Aspect_Ratio, linewidth = 3, xlabel = "D [mm]", ylabel = "aspect ratio", label = "Aspect Ratio", color = :red) -# save plot -PL.plot(p1, p2, layout = (1, 2), size = (800, 500), dpi = 300) -PL.savefig("1M_individual_terminal_velocity_comparisons.svg") - -#! format: on diff --git a/src/Common.jl b/src/Common.jl index 856d5acf5..ae45ae9e3 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -16,9 +16,11 @@ export H2SO4_soln_saturation_vapor_pressure export a_w_xT export a_w_eT export a_w_ice -export Chen2022_vel_add -export Chen2022_vel_coeffs_small -export Chen2022_vel_coeffs_large +export Chen2022_monodisperse_pdf +export Chen2022_exponential_pdf +export Chen2022_vel_coeffs_B1 +export Chen2022_vel_coeffs_B2 +export Chen2022_vel_coeffs_B4 """ G_func(air_props, tps, T, Liquid()) @@ -206,15 +208,15 @@ function a_w_ice(tps::TPS, T::FT) where {FT} end """ - Chen2022_vel_coeffs_small(precip_type, velo_scheme, ρ) + Chen2022_vel_coeffs_B1(ρ) - - velo_scheme - type for terminal velocity scheme (contains free parameters) - ρ - air density + - velo_scheme - a struct with terminal velocity free parameters -Returns the coefficients from Appendix B in Chen et al 2022 +Returns the coefficients from Table B1 Appendix B in Chen et al 2022 DOI: 10.1016/j.atmosres.2022.106171 """ -function Chen2022_vel_coeffs_small( +function Chen2022_vel_coeffs_B1( velo_scheme::CMP.Chen2022VelTypeRain, ρ::FT, ) where {FT} @@ -232,7 +234,17 @@ function Chen2022_vel_coeffs_small( return (aiu, bi, ciu) end -function Chen2022_vel_coeffs_small( + +""" + Chen2022_vel_coeffs_B2(ρ) + + - ρ - air density + - velo_scheme - a struct with terminal velocity free parameters + +Returns the coefficients from Table B2 Appendix B in Chen et al 2022 +DOI: 10.1016/j.atmosres.2022.106171 +""" +function Chen2022_vel_coeffs_B2( velo_scheme::CMP.Chen2022VelTypeSnowIce{FT}, ρ::FT, ) where {FT} @@ -249,15 +261,15 @@ function Chen2022_vel_coeffs_small( end """ - Chen2022_vel_coeffs_large(velo_scheme, ρ) + Chen2022_vel_coeffs_B4(ρ) - - velo_scheme - type for terminal velocity scheme (contains free parameters) - ρ - air density + - velo_scheme - a struct with terminal velocity free parameters -Returns the coefficients from Appendix B (table B4 for large particles) in Chen et al 2022 +Returns the coefficients from Table B4 Appendix B in Chen et al 2022 DOI: 10.1016/j.atmosres.2022.106171 """ -function Chen2022_vel_coeffs_large( +function Chen2022_vel_coeffs_B4( velo_scheme::CMP.Chen2022VelTypeSnowIce{FT}, ρ::FT, ) where {FT} @@ -273,8 +285,27 @@ function Chen2022_vel_coeffs_large( return (aiu, bi, ciu) end +# Wrapper to cast types from SF.gamma +# (which returns Float64 even when the input is Float32) +# TODO - replace with parameterization of our own +Γ(a::FT) where {FT <: Real} = FT(SF.gamma(a)) + """ - Chen2022_vel_add(a, b, c, λ, k) + Chen2022_monodisperse_pdf(a, b, c, D) + + - a, b, c, - free parameters defined in Chen etl al 2022 + - D - droplet diameter + +Returns the addends of the bulk fall speed of rain or ice particles +following Chen et al 2022 DOI: 10.1016/j.atmosres.2022.106171 in [m/s]. +Assuming monodisperse droplet distribution. +""" +function Chen2022_monodisperse_pdf(a, b, c, D) + return a * D^b * exp(-c * D) +end + +""" + Chen2022_exponential_pdf(a, b, c, λ, k) - a, b, c, - free parameters defined in Chen etl al 2022 - λ - size distribution parameter @@ -282,11 +313,12 @@ end Returns the addends of the bulk fall speed of rain or ice particles following Chen et al 2022 DOI: 10.1016/j.atmosres.2022.106171 in [m/s]. -We are assuming exponential size distribution and hence μ=0. +Assuming exponential size distribution and hence μ=0. """ -function Chen2022_vel_add(a::FT, b::FT, c::FT, λ::FT, k::Int) where {FT} +function Chen2022_exponential_pdf(a::FT, b::FT, c::FT, λ::FT, k::Int) where {FT} μ = 0 # Exponential instead of gamma distribution δ = FT(μ + k + 1) - return a * λ^δ * SF.gamma(b + δ) / (λ + c)^(b + δ) / SF.gamma(δ) -end + return a * λ^δ * Γ(b + δ) / (λ + c)^(b + δ) / Γ(δ) end + +end # module end diff --git a/src/Microphysics1M.jl b/src/Microphysics1M.jl index ae7f8040c..1a0eea1a4 100644 --- a/src/Microphysics1M.jl +++ b/src/Microphysics1M.jl @@ -93,7 +93,7 @@ function lambda( return q > FT(0) ? ( - χm * m0 * n0 * SF.gamma(me + Δm + FT(1)) / ρ / q / r0^(me + Δm) + χm * m0 * n0 * CO.Γ(me + Δm + FT(1)) / ρ / q / r0^(me + Δm) )^FT(1 / (me + Δm + 1)) : FT(0) end @@ -153,28 +153,26 @@ function terminal_velocity( # size distrbution λ = lambda(pdf, mass, q, ρ) - return χv * - v0 * - (λ * r0)^(-ve - Δv) * - SF.gamma(me + ve + Δm + Δv + FT(1)) / SF.gamma(me + Δm + FT(1)) + return χv * v0 * (λ * r0)^(-ve - Δv) * CO.Γ(me + ve + Δm + Δv + FT(1)) / + CO.Γ(me + Δm + FT(1)) else return FT(0) end end function terminal_velocity( - (; pdf, mass)::Union{CMP.Rain{FT}, CMP.CloudIce{FT}}, - vel::Union{CMP.Chen2022VelTypeRain{FT}, CMP.Chen2022VelTypeSnowIce{FT}}, + (; pdf, mass)::CMP.Rain{FT}, + vel::CMP.Chen2022VelTypeRain{FT}, ρ::FT, q::FT, ) where {FT} fall_w = FT(0) if q > FT(0) # coefficients from Table B1 from Chen et. al. 2022 - aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(vel, ρ) + aiu, bi, ciu = CO.Chen2022_vel_coeffs_B1(vel, ρ) # size distribution parameter λ::FT = lambda(pdf, mass, q, ρ) # eq 20 from Chen et al 2022 - fall_w = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 3)) + fall_w = sum(CO.Chen2022_exponential_pdf.(aiu, bi, ciu, λ, 3)) # It should be ϕ^κ * fall_w, but for rain drops ϕ = 1 and κ = 0 fall_w = max(FT(0), fall_w) end @@ -187,37 +185,31 @@ function terminal_velocity( q::FT, ) where {FT} fall_w = FT(0) + # For now, we assume the B4 table coefficients for snow + # and B2 table coefficients for cloud ice. + # TODO - we should do partial integrals + # from D=125um to D=625um using B2 and D=625um to inf using B4. if q > FT(0) - - (; r0, m0, me, Δm, χm) = mass - (; a0, ae, Δa, χa) = area + # coefficients from Table B4 from Chen et. al. 2022 + aiu, bi, ciu = CO.Chen2022_vel_coeffs_B4(vel, ρ) + #aiu, bi, ciu = CO.Chen2022_vel_coeffs_B2(vel, ρ) + # size distribution parameter λ::FT = lambda(pdf, mass, q, ρ) - m0c = m0 * χm - a0c = a0 * χa - mec = me + Δm - aec = ae + Δa + # As a next step, we could keep ϕ(r) under the integrals + # ϕ(r) = 16 * ρᵢ^2 * aᵢ(r)^3 / (9 * π * mᵢ(r)^2) - # coefficients from Appendix B from Chen et. al. 2022 - aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(vel, ρ) + # compute the mass weighted average aspect ratio + (; r0, m0, me, Δm, χm) = mass + (; a0, ae, Δa, χa) = area ρᵢ = vel.ρᵢ - κ = FT(-1 / 3) #oblate - k = 3 # mass weighted - - tmp = - λ^(k + 1) * - ((16 * a0c^3 * ρᵢ^2) / (9 * π * m0c^2 * r0^(3 * aec - 2 * mec)))^κ - ci_pow = - (2 .* ciu .+ λ) .^ - (.-(3 .* aec .* κ .- 2 .* mec .* κ .+ bi .+ k .+ 1)) - - ti = tmp .* aiu .* FT(2) .^ bi .* ci_pow - - Chen2022_vel_add_sno(t, b, aec, mec, κ, k) = - t * SF.gamma(3 * κ * aec - 2 * κ * mec + b + k + 1) / - SF.gamma(k + 1) + α = 3 * (ae + Δa) - 2 * (me + Δm) + ϕ₀ = 16 * ρᵢ^2 / 9 / FT(π) * (χa * a0)^3 / (χm * m0)^2 / r0^α + ϕ_avg = ϕ₀ / λ^α * CO.Γ(α + 3 + 1) / CO.Γ(3 + 1) - fall_w = sum(Chen2022_vel_add_sno.(ti, bi, aec, mec, κ, k)) + # eq 20 from Chen 2022 + κ = FT(-1 / 3) # assuming oblate + fall_w = ϕ_avg^κ * sum(CO.Chen2022_exponential_pdf.(aiu, bi, ciu, λ, 3)) fall_w = max(FT(0), fall_w) end return fall_w @@ -340,7 +332,7 @@ function accretion( accr_rate = q_clo * E * n0 * a0 * v0 * χa * χv / λ * - SF.gamma(ae + ve + Δa + Δv + FT(1)) / (λ * r0)^(ae + ve + Δa + Δv) + CO.Γ(ae + ve + Δa + Δv + FT(1)) / (λ * r0)^(ae + ve + Δa + Δv) end return accr_rate end @@ -386,7 +378,7 @@ function accretion_rain_sink( accr_rate = E / ρ * n0 * n0_ice * m0 * a0 * v0 * χm * χa * χv / λ_ice / λ * - SF.gamma(me + ae + ve + Δm + Δa + Δv + FT(1)) / + CO.Γ(me + ae + ve + Δm + Δa + Δv + FT(1)) / (r0 * λ)^FT(me + ae + ve + Δm + Δa + Δv) end return accr_rate @@ -443,11 +435,11 @@ function accretion_snow_rain( accr_rate = FT(π) / ρ * n0_i * n0_j * m0_j * χm_j * E_ij * abs(v_ti - v_tj) / r0_j^(me_j + Δm_j) * ( - FT(2) * SF.gamma(me_j + Δm_j + FT(1)) / λ_i^FT(3) / + FT(2) * CO.Γ(me_j + Δm_j + FT(1)) / λ_i^FT(3) / λ_j^(me_j + Δm_j + FT(1)) + - FT(2) * SF.gamma(me_j + Δm_j + FT(2)) / λ_i^FT(2) / + FT(2) * CO.Γ(me_j + Δm_j + FT(2)) / λ_i^FT(2) / λ_j^(me_j + Δm_j + FT(2)) + - SF.gamma(me_j + Δm_j + FT(3)) / λ_i / λ_j^(me_j + Δm_j + FT(3)) + CO.Γ(me_j + Δm_j + FT(3)) / λ_i / λ_j^(me_j + Δm_j + FT(3)) ) end return accr_rate @@ -502,7 +494,7 @@ function evaporation_sublimation( b_vent * (ν_air / D_vapor)^FT(1 / 3) / (r0 * λ)^((ve + Δv) / FT(2)) * (FT(2) * v0 * χv / ν_air / λ)^FT(1 / 2) * - SF.gamma((ve + Δv + FT(5)) / FT(2)) + CO.Γ((ve + Δv + FT(5)) / FT(2)) ) end # only evaporation is considered for rain @@ -541,7 +533,7 @@ function evaporation_sublimation( b_vent * (ν_air / D_vapor)^FT(1 / 3) / (r0 * λ)^((ve + Δv) / FT(2)) * (FT(2) * v0 * χv / ν_air / λ)^FT(1 / 2) * - SF.gamma((ve + Δv + FT(5)) / FT(2)) + CO.Γ((ve + Δv + FT(5)) / FT(2)) ) end return evap_subl_rate @@ -592,7 +584,7 @@ function snow_melt( b_vent * (ν_air / D_vapor)^FT(1 / 3) / (r0 * λ)^((ve + Δv) / FT(2)) * (FT(2) * v0 * χv / ν_air / λ)^FT(1 / 2) * - SF.gamma((ve + Δv + FT(5)) / FT(2)) + CO.Γ((ve + Δv + FT(5)) / FT(2)) ) end return snow_melt_rate diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl index d6ebc3652..0244cd9b5 100644 --- a/src/Microphysics2M.jl +++ b/src/Microphysics2M.jl @@ -524,7 +524,7 @@ function rain_particle_terminal_velocity( Chen2022::CMP.Chen2022VelTypeRain, ρₐ::FT, ) where {FT} - (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρₐ) + (ai, bi, ci) = CO.Chen2022_vel_coeffs_B1(Chen2022, ρₐ) v = sum(@. sum(ai * D^bi * exp(-ci * D))) return v @@ -574,13 +574,13 @@ function rain_terminal_velocity( N_rai::FT, ) where {FT} # coefficients from Table B1 from Chen et. al. 2022 - aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(vel, ρ) + aiu, bi, ciu = CO.Chen2022_vel_coeffs_B1(vel, ρ) # size distribution parameter λ = pdf_rain_parameters(pdf_r, q_rai, ρ, N_rai).λr # eq 20 from Chen et al 2022 - vt0 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 0)) - vt3 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 3)) + vt0 = sum(CO.Chen2022_exponential_pdf.(aiu, bi, ciu, λ, 0)) + vt3 = sum(CO.Chen2022_exponential_pdf.(aiu, bi, ciu, λ, 3)) vt0 = N_rai < eps(FT) ? FT(0) : max(FT(0), vt0) vt3 = q_rai < eps(FT) ? FT(0) : max(FT(0), vt3) diff --git a/src/MicrophysicsNonEq.jl b/src/MicrophysicsNonEq.jl index 2998ef8c2..c8ca6d53f 100644 --- a/src/MicrophysicsNonEq.jl +++ b/src/MicrophysicsNonEq.jl @@ -9,6 +9,8 @@ module MicrophysicsNonEq import Thermodynamics as TD import Thermodynamics.Parameters as TDP +import CloudMicrophysics.Common as CO + import ..Parameters as CMP export τ_relax @@ -114,4 +116,62 @@ function conv_q_vap_to_q_liq_ice_MM2015( return (qᵥ - qᵥ_sat_ice) / (τ_relax * Γᵢ) end +""" + terminal_velocity(sediment, vel, ρ, q) + + - `sediment` - a struct with sedimentation type (cloud liquid or ice) + - `vel` - a struct with terminal velocity parameters + - `ρ` - air density + - `q` - cloud liquid or ice specific humidity + +Returns the mass weighted average terminal velocity assuming a +monodisperse size distribution with prescribed number concentration. +The fall velocity of individual particles is parameterized following +Chen et. al 2022, DOI: 10.1016/j.atmosres.2022.106171 +""" +function terminal_velocity( + (; ρw)::CMP.CloudLiquid{FT}, + vel::CMP.Chen2022VelTypeRain{FT}, + ρ::FT, + q::FT, +) where {FT} + fall_w = FT(0) + if q > FT(0) + # TODO: Coefficients from Table B1 from Chen et. al. 2022 are only valid + # for D > 100mm. We should look for a different parameterization + # that is more suited for cloud droplets. For now I'm just multiplying + # by an arbitrary correction factor. + aiu, bi, ciu = CO.Chen2022_vel_coeffs_B1(vel, ρ) + # The 1M scheme does not assume any cloud droplet size distribution. + # TODO - For now I compute a mean volume radius assuming a fixed value + # for the total number concentration of droplets. + N = FT(500 * 1e6) + D = cbrt(ρ * q / N / ρw) + corr = FT(0.1) + # assuming ϕ = 1 (spherical) + fall_w = sum(CO.Chen2022_monodisperse_pdf.(aiu, bi, ciu, D)) + fall_w = max(FT(0), corr * fall_w) + end + return fall_w +end +function terminal_velocity( + (; pdf, mass, ρi)::CMP.CloudIce{FT}, + vel::CMP.Chen2022VelTypeSnowIce{FT}, + ρ::FT, + q::FT, +) where {FT} + fall_w = FT(0) + if q > FT(0) + # Coefficients from Table B2 from Chen et. al. 2022 + aiu, bi, ciu = CO.Chen2022_vel_coeffs_B2(vel, ρ) + # See the comment for liquid droplets above + N = FT(500 * 1e6) + D = cbrt(ρ * q / N / ρi) + # assuming ϕ = 1 (spherical) + fall_w = sum(CO.Chen2022_monodisperse_pdf.(aiu, bi, ciu, D)) + fall_w = max(FT(0), fall_w) + end + return fall_w +end + end #module MicrophysicsNonEq.jl diff --git a/src/P3_helpers.jl b/src/P3_helpers.jl index f5a560d38..af18ae512 100644 --- a/src/P3_helpers.jl +++ b/src/P3_helpers.jl @@ -1,8 +1,3 @@ -# Wrapper to cast types from SF.gamma -# (which returns Float64 even when the input is Float32) -# TODO - replace with parameterization of our own -Γ(a::FT) where {FT <: Real} = FT(SF.gamma(a)) - # helper functions for Γ_approx function c₁(a::FT) where {FT <: Real} p1 = FT(9.4368392235e-03) @@ -81,7 +76,7 @@ function Γ_lower(a::FT, z::FT) where {FT <: Real} c₁(a) * z / a / (a + 1) + (c₁(a) * z)^2 / a / (a + 1) / (a + 2) ) * - (1 - W(a, z)) + Γ(a) * W(a, z) * (1 - c₄(a)^(-z)) + (1 - W(a, z)) + CO.Γ(a) * W(a, z) * (1 - c₄(a)^(-z)) end end @@ -91,7 +86,7 @@ end An approximated upper incomplete Gamma function computed as Γ(a) - Γ_lower(a, z) """ function Γ_upper(a::FT, z::FT) where {FT <: Real} - return Γ(a) - Γ_lower(a, z) + return CO.Γ(a) - Γ_lower(a, z) end """ diff --git a/src/P3_size_distribution.jl b/src/P3_size_distribution.jl index 848a47c5a..8da23715d 100644 --- a/src/P3_size_distribution.jl +++ b/src/P3_size_distribution.jl @@ -35,7 +35,7 @@ Returns the shape parameter N₀ from Eq. 2 in Morrison and Milbrandt (2015). """ function DSD_N₀(p3::PSP3, N::FT, λ::FT) where {FT} μ = DSD_μ(p3, λ) - return N / Γ(1 + μ) * λ^(1 + μ) + return N / CO.Γ(1 + μ) * λ^(1 + μ) end """ @@ -159,7 +159,7 @@ function L_over_N_gamma( D_th = D_th_helper(p3) λ = exp(log_λ) - N = Γ(1 + μ) / (λ^(1 + μ)) + N = CO.Γ(1 + μ) / (λ^(1 + μ)) return ifelse( F_rim == FT(0), (p3_F_liq_average( diff --git a/src/P3_terminal_velocity.jl b/src/P3_terminal_velocity.jl index c20569a22..805b7d72a 100644 --- a/src/P3_terminal_velocity.jl +++ b/src/P3_terminal_velocity.jl @@ -23,9 +23,9 @@ function ice_particle_terminal_velocity( use_aspect_ratio = true, ) where {FT} if D <= Chen2022.cutoff - (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρₐ) + (ai, bi, ci) = CO.Chen2022_vel_coeffs_B2(Chen2022, ρₐ) else - (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρₐ) + (ai, bi, ci) = CO.Chen2022_vel_coeffs_B4(Chen2022, ρₐ) end v = sum(@. sum(ai * D^bi * exp(-ci * D))) diff --git a/src/parameters/Microphysics1M.jl b/src/parameters/Microphysics1M.jl index f0cd17279..7a9d7ad92 100644 --- a/src/parameters/Microphysics1M.jl +++ b/src/parameters/Microphysics1M.jl @@ -113,20 +113,25 @@ The parameters and type for cloud liquid water condensate # Fields $(DocStringExtensions.FIELDS) """ -struct CloudLiquid{FT} <: CloudCondensateType{FT} +Base.@kwdef struct CloudLiquid{FT} <: CloudCondensateType{FT} "condensation evaporation non_equil microphysics relaxation timescale [s]" τ_relax::FT + "water density [kg/m3]" + ρw::FT end CloudLiquid(::Type{FT}) where {FT <: AbstractFloat} = CloudLiquid(CP.create_toml_dict(FT)) function CloudLiquid(toml_dict::CP.AbstractTOMLDict) - name_map = (; :condensation_evaporation_timescale => :τ_relax) - (; τ_relax) = + name_map = (; + :condensation_evaporation_timescale => :τ_relax, + :density_liquid_water => :ρw, + ) + parameters = CP.get_parameter_values(toml_dict, name_map, "CloudMicrophysics") FT = CP.float_type(toml_dict) - return CloudLiquid{FT}(τ_relax) + return CloudLiquid{FT}(; parameters...) end """ @@ -148,6 +153,8 @@ struct CloudIce{FT, PD, MS} <: CloudCondensateType{FT} r_ice_snow::FT "deposition sublimation non_equil microphysics relaxation timescale [s]" τ_relax::FT + "ice density [kg/m3]" + ρi::FT end CloudIce(::Type{FT}) where {FT <: AbstractFloat} = @@ -168,7 +175,7 @@ function CloudIce(toml_dict::CP.AbstractTOMLDict = CP.create_toml_dict(FT)) FT = CP.float_type(toml_dict) P = typeof(pdf) M = typeof(mass) - return CloudIce{FT, P, M}(pdf, mass, p.r0, p.r_ice_snow, p.τ_relax) + return CloudIce{FT, P, M}(pdf, mass, p.r0, p.r_ice_snow, p.τ_relax, p.ρi) end function ParticleMass(::Type{CloudIce}, td::CP.AbstractTOMLDict) diff --git a/test/common_functions_tests.jl b/test/common_functions_tests.jl index 5ff805a4d..62f89ebda 100644 --- a/test/common_functions_tests.jl +++ b/test/common_functions_tests.jl @@ -153,7 +153,7 @@ function test_Chen_coefficients(FT) Ch2022 = CMP.Chen2022VelType(FT) TT.@testset "Chen terminal velocity rain (B1)" begin - aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(Ch2022.rain, ρ) + aiu, bi, ciu = CO.Chen2022_vel_coeffs_B1(Ch2022.rain, ρ) TT.@test all( isapprox.( @@ -179,7 +179,7 @@ function test_Chen_coefficients(FT) end TT.@testset "Chen terminal velocity small ice (B2)" begin - aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(Ch2022.snow_ice, ρ) + aiu, bi, ciu = CO.Chen2022_vel_coeffs_B2(Ch2022.snow_ice, ρ) TT.@test all( isapprox.( @@ -195,7 +195,7 @@ function test_Chen_coefficients(FT) end TT.@testset "Chen terminal velocity large ice (B4)" begin - aiu, bi, ciu = CO.Chen2022_vel_coeffs_large(Ch2022.snow_ice, ρ) + aiu, bi, ciu = CO.Chen2022_vel_coeffs_B4(Ch2022.snow_ice, ρ) TT.@test all( isapprox.( diff --git a/test/microphysics1M_tests.jl b/test/microphysics1M_tests.jl index 5d4beece6..1441c9699 100644 --- a/test/microphysics1M_tests.jl +++ b/test/microphysics1M_tests.jl @@ -84,23 +84,6 @@ function test_microphysics1M(FT) TT.@test v_bigger > vt_rai end - TT.@testset "1M_microphysics - Chen 2022 ice terminal velocity" begin - #setup - ρ = FT(1.2) - q_ice = FT(5e-4) - N_rai = FT(1e4) - - #action - vt_ice = CM1.terminal_velocity(ice, Ch2022.snow_ice, ρ, q_ice) - v_bigger = CM1.terminal_velocity(ice, Ch2022.snow_ice, ρ, q_ice * 2) - - #test - TT.@test vt_ice ≈ 2.2309476335899388 rtol = 2e-6 - TT.@test CM1.terminal_velocity(ice, Ch2022.snow_ice, ρ, FT(0)) ≈ 0 atol = - eps(FT) - TT.@test v_bigger > vt_ice - end - TT.@testset "1M_microphysics - Chen 2022 snow terminal velocity" begin #setup ρ = FT(1.1) @@ -112,7 +95,7 @@ function test_microphysics1M(FT) v_bigger = CM1.terminal_velocity(snow, Ch2022.snow_ice, ρ, q_sno * 2) #test - TT.@test vt_sno ≈ 1.4722373984934243 rtol = 2e-6 + TT.@test vt_sno ≈ 0.7529349430245454 rtol = 2e-6 TT.@test CM1.terminal_velocity(snow, Ch2022.snow_ice, ρ, FT(0)) ≈ 0 atol = eps(FT) TT.@test v_bigger > vt_sno diff --git a/test/microphysics_noneq_tests.jl b/test/microphysics_noneq_tests.jl index 094b5383a..ba395a8ec 100644 --- a/test/microphysics_noneq_tests.jl +++ b/test/microphysics_noneq_tests.jl @@ -12,6 +12,7 @@ function test_microphysics_noneq(FT) ice = CMP.CloudIce(FT) liquid = CMP.CloudLiquid(FT) tps = TD.Parameters.ThermodynamicsParameters(FT) + Ch2022 = CMP.Chen2022VelType(FT) TT.@testset "τ_relax" begin TT.@test CMNe.τ_relax(liquid) ≈ FT(10) @@ -87,6 +88,38 @@ function test_microphysics_noneq(FT) #! format: on end + + TT.@testset "Cloud condensate sedimentation - liquid" begin + #setup + ρ = FT(1.1) + q_liq = FT(1 * 1e-3) + + #action + vt_zero = CMNe.terminal_velocity(liquid, Ch2022.rain, ρ, FT(0)) + vt_liq = CMNe.terminal_velocity(liquid, Ch2022.rain, ρ, q_liq) + v_bigger = CMNe.terminal_velocity(liquid, Ch2022.rain, ρ, q_liq * 2) + + #test + TT.@test vt_zero == FT(0) + TT.@test vt_liq ≈ 0.004249204292098458 rtol = sqrt(eps(FT)) + TT.@test v_bigger > vt_liq + end + + TT.@testset "Cloud condensate sedimentation - ice" begin + #setup + ρ = FT(0.75) + q_ice = FT(0.5 * 1e-3) + + #action + vt_zero = CMNe.terminal_velocity(ice, Ch2022.snow_ice, ρ, FT(0)) + vt_ice = CMNe.terminal_velocity(ice, Ch2022.snow_ice, ρ, q_ice) + v_bigger = CMNe.terminal_velocity(ice, Ch2022.snow_ice, ρ, q_ice * 2) + + #test + TT.@test vt_zero == FT(0) + TT.@test vt_ice ≈ 0.005587793323944559 rtol = sqrt(eps(FT)) + TT.@test v_bigger > vt_ice + end end test_microphysics_noneq(Float32)