Skip to content

Commit

Permalink
cleanup in docs
Browse files Browse the repository at this point in the history
  • Loading branch information
@trontrytel
trontrytel committed Aug 11, 2023
1 parent 1ff2095 commit 474ce0b
Show file tree
Hide file tree
Showing 2 changed files with 107 additions and 150 deletions.
245 changes: 101 additions & 144 deletions docs/src/Microphysics1M.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -259,17 +339,16 @@ 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}
{\int_0^\infty n(r) \, m(r) \, dr}
\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}
Expand All @@ -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

Expand Down
12 changes: 6 additions & 6 deletions docs/src/TerminalVelocityComparisons.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")

Expand Down

0 comments on commit 474ce0b

Please sign in to comment.