If one introduces $N \geq 1$ tracers (i.e. parts of the air that do not belong to dry air), of which $1 \leq N_C \leq N - 1$ are condensed, one needs $N$ additional continuity equations. If one makes the assumption of local thermodynamic equilibrium (all condensates have the same temperature as the gas phase), no additional first laws need to be introduced. The basic ansatz for the dynamic and thermodynamic coupling of the gaseous and condensed tracers to the system of equations for dry air was established in Sect. 10.3.
Let $\rho$ be a mass density. In preparation, the Crank-Nicolson method is written down with $n$ as the current time step:
\[ \begin{align} \rho^{(n + 1)} = \rho^{(n)} + \Delta t\frac{1}{2}\left(\frac{\partial\rho^{(n)}}{\partial t} + \frac{\partial\rho^{(n + 1)}}{\partial t}\right) \end{align} \]
At this point, one is only concerned with vertical convergence; the other terms are calculated explicitly:
\[ \begin{align} \rho^{(n + 1)} = \rho^{(n)} + \Delta t\newdot{\rho}^{(n)}_{\text{expl}} - \Delta t\frac{1}{2}\left(\frac{\partial\left(\rho^{(n)}w^{(n + 1)}\right)}{\partial z} + \frac{\partial\left(\rho^{(n + 1)}w^{(n + 1)}\right)}{\partial z}\right) \end{align} \]
Only the generalized density is treated in the form of the Crank-Nicolson method; the velocity is taken from the new time step for reasons of stability. The part
\[ \begin{align} -\frac{1}{2}\frac{\partial\left(\rho^{(n)}w^{(n + 1)}\right)}{\partial z} \end{align} \]
is absorbed into the explicit tendency; that is, one redefines
\[ \begin{align} \newdot{\rho}^{(n)}_{\text{expl}} \to \newdot{\rho}^{(n)}_{\text{expl}} - \frac{1}{2}\frac{\partial\left(\rho^{(n)}w^{(n + 1)}\right)}{\partial z} \end{align} \]
Vertical discretization requires a layer index $1 \leq i \leq N_L$:
\[ \begin{align} \rho_i^{(n + 1)} = \rho_i^{(n)} + \Delta t\newdot{\rho_i}^{(n)}_{\text{expl}} - \frac{\Delta t}{2}\frac{A_{i - 1/2}\newoverline{\rho_k^{(n + 1)}}^{(i - 1/2)}w_{i - 1/2}^{(n + 1)} - A_{i + 1/2}\newoverline{\rho_k^{(n + 1)}}^{(i + 1/2)}w_{i + 1/2}^{(n + 1)}}{V_i}\tag{36.5}\label{eq:adv_gen_density_deriv_0} \end{align} \]
For the vertical averages of the density, one writes in general
\[ \begin{align} \newoverline{\rho_k}^{(i + 1/2)} = \alpha_{i + \frac{1}{2}, i}\rho_{i} + \alpha_{i + \frac{1}{2}, i + 1}\rho_{i + 1} \end{align} \]
with
\[ \begin{align} \alpha_{i + \frac{1}{2}, i} + \alpha_{i + \frac{1}{2}, i + 1} = 1. \end{align} \]
Substituting this into Eq. (36.5), one obtains
\[ \begin{align} &\rho_i^{(n + 1)} = \rho_i^{(n)} + \Delta t\newdot{\rho_i}^{(n)}_{\text{expl}}\nonumber\\ &- \frac{\Delta t}{2}\left(\alpha_{i - \frac{1}{2}, i - 1}\rho_{i - 1}^{(n + 1)} + \alpha_{i - \frac{1}{2}, i}\rho_i^{(n + 1)}\right)\frac{A_{i - 1/2}w_{i - 1/2}^{(n + 1)}}{V_i} + \frac{\Delta t}{2}\left(\alpha_{i + \frac{1}{2}, i}\rho_i^{(n + 1)} + \alpha_{i + \frac{1}{2}, i + 1}\rho_{i + 1}^{(n + 1)}\right)\frac{A_{i + 1/2}w_{i + 1/2}^{(n + 1)}}{V_i}. \end{align} \]
One now defines
\[ \begin{align} W_{i + 1/2} \coloneqq A_{i + 1/2}w_{i + 1/2}^{(n + 1)}, \end{align} \]
thus one obtains
\[ \begin{align} &\rho_i^{(n + 1)} = \rho_i^{(n)} + \Delta t\newdot{\rho_i}^{(n)}_{\text{expl}}\nonumber\\ &- \frac{\Delta t}{2}\left(\alpha_{i - \frac{1}{2}, i - 1}\rho_{i - 1}^{(n + 1)} + \alpha_{i - \frac{1}{2}, i}\rho_i^{(n + 1)}\right)\frac{W_{i - 1/2}^{(n + 1)}}{V_i} + \frac{\Delta t}{2}\left(\alpha_{i + \frac{1}{2}, i}\rho_i^{(n + 1)} + \alpha_{i + \frac{1}{2}, i + 1}\rho_{i + 1}^{(n + 1)}\right)\frac{W_{i + 1/2}^{(n + 1)}}{V_i}. \end{align} \]
At $i = 1$, $w_{i - 1/2} = 0$:
\[ \begin{align} &\rho_1^{(n + 1)} = \rho_1^{(n)} + \Delta t\newdot{\rho}^{(n)}_{1, \text{expl}} + \frac{\Delta t}{2}\left(\alpha_{\frac{3}{2}, 1}\rho_1^{(n + 1)} + \alpha_{\frac{3}{2}, 2}\rho_{2}^{(n + 1)}\right)\frac{W_{3/2}^{(n + 1)}}{V_1} \end{align} \]
At $i = N_L$, $w_{i + 1/2} = 0$:
\[ \begin{align} \rho_{N_L}^{(n + 1)} = \rho_{N_L}^{(n)} + \Delta t\newdot{\rho}^{(n)}_{N_L, \text{expl}} - \frac{\Delta t}{2}\left(\alpha_{N_L - \frac{1}{2}, N_L - 1}\rho_{N_L - 1}^{(n + 1)} + \alpha_{N_L - \frac{1}{2}, N_L}\rho_{N_L}^{(n + 1)}\right)\frac{W_{N_L - 1/2}^{(n + 1)}}{V_{N_L}}. \end{align} \]
The vector $\mathbf{x}$ of the unknowns is defined by
\[ \begin{align} \mathbf{x} = \left(\begin{array}{c} \rho^{(n + 1)}_{1}\\ \rho^{(n + 1)}_{2}\\ \vdots\\ \rho^{(n + 1)}_{N_L} \end{array}\right), \end{align} \]
for which a linear system of equations holds
\[ \begin{align} A\cdot\mathbf{x} = \mathbf{r} \end{align} \]
with a matrix $A$ and a right-hand side $\mathbf{r}$. For this right-hand side, one has
\[ \begin{align} \mathbf{r} = \left(\begin{array}{c} \rho^{(n)}_{1} + \Delta t\newdot{\rho}^{(n)}_{1, \text{expl}}\\ \rho^{(n)}_{2} + \Delta t\newdot{\rho}^{(n)}_{2, \text{expl}}\\ \rho^{(n)}_{3} + \Delta t\newdot{\rho}^{(n)}_{3, \text{expl}}\\ \vdots\\ \rho^{(n)}_{N_L - 1} + \Delta t\newdot{\rho}^{(n)}_{N_L - 1, \text{expl}}\\ \rho^{(n)}_{N_L} + \Delta t\newdot{\rho}^{(n)}_{N_L, \text{expl}} \end{array}\right). \end{align} \]
For the matrix $A$ one obtains
\[ \begin{align} A &= \left(\begin{array}{cccc} d_1 & e_1 & \dots & 0 \\ b_1 & d_2 & e_2 \hspace{2 cm}\dots & 0 \\ \vdots & \hspace{2 cm}\ddots & \ddots & \vdots \\ 0 & \dots & b_{N_L - 1} & d_{N_L} \end{array}\right) \end{align} \]
with vectors $\mathbf{c}, \mathbf{e} \in \mathbb{R}^{N_L - 1}$, $\mathbf{e} \in \mathbb{R}^{N_L}$. For these, one obtains
\[ \begin{align} c_i &= \frac{\Delta t}{2}\frac{\alpha_{i + \frac{1}{2}, i}W_{i + 1/2}^{(n + 1)}}{V_{i + 1}},\\ d_1 &= 1 - \frac{\Delta t}{2}\alpha_{\frac{3}{2}, 1}\frac{W_{3/2}^{(n + 1)}}{V_1},\\ d_i &= 1 + \frac{\Delta t}{2}\frac{\alpha_{i - \frac{1}{2}, i}W_{i - 1/2}^{(n + 1)} - \alpha_{i + \frac{1}{2}, i}W_{i + 1/2}^{(n + 1)}}{V_i},\\ d_{N_L} &= 1 + \frac{\Delta t}{2}\frac{\alpha_{N_L - \frac{1}{2}, N_L}W_{N_L - 1/2}^{(n + 1)}}{V_{N_L}},\\ e_i &= -\frac{\Delta t}{2}\alpha_{i + \frac{1}{2}, i + 1}\frac{W_{i + 1/2}^{(n + 1)}}{V_i}. \end{align} \]
Limiters are used to enforce bounds on fluxes. In practice, these bounds usually reflect the requirement that mass densities resulting from the fluxes must not become negative, since negative densities are unphysical. For tracers with strongly varying densities, this can occur through advection alone. Because limiters are numerical interventions that cannot be derived purely from physical laws, the numerics should be designed so that limiters are invoked as rarely as possible and, when invoked, apply only minimal corrections.
After advection, diffusion, and phase transitions, negative densities may still occur despite the limiter. These must be reset to zero at the end of the time step. Although this violates mass conservation, there is no practical alternative.
If the air is supersaturated at time $t$, a phase transition from gaseous to liquid or from gaseous to solid occurs until the air is saturated. If the air is undersaturated and cloud water or cloud ice is still present, the phase transition takes place in the other direction. For modeling it is assumed that after a time step of length $\Delta t$ the air is saturated, i.e
\[ \begin{align} p_v\left(t + \Delta t\right) = p_S\left(t + \Delta t\right),\tag{36.22}\label{eq:model_phase_trans_deriv_1} \end{align} \]
Here $p_v$ is the vapor pressure and $p_S$ is, as usual, the saturation vapor pressure. From this the density $\Delta\rho_v$ of water vapor should be determined, which must condense or resublimate to reach this state. This is not simply $\rho_S\left(t\right) - \rho_v\left(t\right),$ where $\rho_S$ is the saturation density, as this would leave the phase transition heat unaccounted for. From Eq. (36.22) follows
\[ \begin{align} \rho_v\left(t + \Delta t\right)R_vT\left(t + \Delta t\right) &= p_S\left[T\left(t + \Delta t\right)\right]\nonumber\\ \Leftrightarrow\left(\rho_v + \Delta\rho_v\right)R_vT\left(t + \Delta t\right) &= p_S\left[T\left(t + \Delta t\right)\right], \end{align} \]
Here, quantities without explicitly stated time dependence denote quantities at time $t$, i.e. $\rho_v\coloneqq\rho_v\left(t\right).$ From this it further follows, with $\rho_t$ as the total density,
\[ \begin{align} \left(\rho_v + \Delta\rho_v\right)R_v\left(T - \Delta\rho_v\frac{L_v}{\rho_tc_p}\right) &= p_S\left(T - \Delta\rho_v\frac{L_v}{\rho_tc_p}\right). \end{align} \]
If one multiplies out the left-hand side and performs a first-order Taylor expansion on the right-hand side, one obtains
\[ \begin{align} p_v + \Delta\rho_vR_vT - R_v\Delta\rho_v\rho_v\frac{L_v}{\rho_tc_p} - \left(\Delta\rho_v\right)^2R_v\frac{L_v}{\rho_tc_p} &= p_S - p_S'\Delta\rho_v\frac{L_v}{\rho_tc_p}\nonumber\\ \Leftrightarrow -\left(\Delta\rho_v\right)^2R_v\frac{L_v}{\rho_tc_p} + \Delta\rho_v\left(R_vT - R_v\rho_v\frac{L_v}{\rho_tc_p} + p_S'\frac{L_v}{\rho_tc_p}\right) + p_v - p_S &= 0\nonumber\\ \Leftrightarrow a\left(\Delta\rho_v\right)^2 + b\Delta\rho_v + c &= 0\nonumber\\ \Leftrightarrow \left(\Delta\rho_v\right)^2 + p\Delta\rho_v + q &= 0 \end{align} \]
with
\[ \begin{align} a &\coloneqq -R_v\frac{L_v}{\rho_tc_p} < 0,\\ b &\coloneqq R_vT - R_v\rho_v\frac{L_v}{\rho_tc_p} + p_S'\frac{L_v}{\rho_tc_p},\\ c &\coloneqq p_v - p_S,\\ p &\coloneqq \frac{b}{a},\\ q &\coloneqq \frac{c}{a}. \end{align} \]
This implies two solutions
\[ \begin{align} \left(\Delta\rho_v\right)_{1,2} = -\frac{p}{2} \pm \sqrt{\frac{p^2}{4} - q}. \end{align} \]
Now it remains to be determined which of the two signs is the physically correct one. In the case of undersaturation, $p_v<p_S,$ which implies $c<0$, so $q>0.$ Furthermore, $\left(\Delta\rho_v\right)^2>0$ and $\Delta\rho_v>0,$ which, because $p\Delta\rho_v = -q - \left(\Delta\rho_v\right)^2 < 0$, implies the relation $p<0$. For $\Delta\rho_v>0$ to hold, the negative sign must be chosen in this case. In the case of supersaturation, $p<0$ also holds, since the saturation vapor pressure does not appear in the expression for $p$, but $q<0$, which also implies the negative sign, since otherwise $\Delta\rho_v\to 0$ does not follow for $q\to 0$. Thus, in any case,
\[ \begin{align} \Delta\rho_v = -\frac{p}{2} - \sqrt{\frac{p^2}{4} - q}. \end{align} \]
If phase transitions are calculated using the averaged thermodynamic variables, the subscale variability considered in Sect. 17.9.3 is neglected. Here the fundamental considerations made in Sect. 17.9.3 are to be made usable for a model. Let $T = T\left(\mathbf{r}\right)$ be the temperature field valid at a fixed time $t$ in a grid box $V\subseteq\mathbb{R}^3$. Let $T_\mathrm{min}$ and $T_\mathrm{max}$ be the minimum and maximum temperatures within this set. Then, due to the strict monotonicity of $p_S\left(T\right)$, the minimum and maximum saturation vapor pressure is derived:
\[ \begin{align} p_{S,\text{min}} &= p_S\left(T_\text{min}\right),\\ p_{S,\text{max}} &= p_S\left(T_\text{max}\right). \end{align} \]
Let $p_v$ be the average vapor pressure within the grid box under consideration. Now one makes the assumption that the vapor pressure $p_v$ is homogeneous within $V$. The basic idea is that condensation begins at $p_v>p_{S,\mathrm{min}}$, but initially in a very small volume fraction, the cloud-filled volume fraction of the grid box $c_f$. For $c_f$, one makes the following ansatz
\[ \begin{align} c_f\left(p_v\leq p_{S,\text{min}}\right) &= 0,\\ c_f\left(p_v\geq p_{S,\text{max}}\right) &= 1. \end{align} \]
This is correct in any case. The simplest way to fulfill this for $p_{S,\mathrm{min}}\leq p_v\leq p_{S,\mathrm{max}}$ is a linear interpolation:
\[ \begin{align} c_f\left(p_v\right) &= \frac{p_v - p_{S,\text{min}}}{p_{S,\text{max}} - p_{S,\text{min}}}. \end{align} \]
This would be the case, for example, if the temperature were uniformly distributed between $T_\mathrm{min}$ and $T_\mathrm{max}$ and the saturation vapor pressure in this interval depended linearly on the temperature, which is satisfied to first order. The idea is now to express the vapor pressure with which the phase-transition rates are calculated, in the case $p_v>p_{S,\mathrm{min}}$ (the only case in which cloud formation can take place at all), via the ansatz
\[ \begin{align} p_v^{(\text{eff})} = p_{S,\text{min}} + \xi\left(p_v - p_{S,\text{min}}\right)\tag{36.38}\label{eq:p_v_eff} \end{align} \]
What remains to be determined is $\xi$. It is calculated via
\[ \begin{align} \xi = \frac{p_v^{(\text{eff})} - p_{S,\text{min}}}{p_v - p_{S,\text{min}}} = \frac{\int_V\text{max}\left(p_v - p_S\left(\mathbf{r}\right),0\right)d^3r}{V\left(p_v - p_{S,\text{min}}\right)}. \end{align} \]
Under the assumption already made above that $T$ and $p_S$ are uniformly distributed within $V$, in the case $p_{S,\mathrm{min}}\leq p_v\leq p_{S,\mathrm{max}}$ it follows
\[ \begin{align} \xi &= \frac{\int_{p_{S,\text{min}}}^{p_v}p_v - p_S'dp_S'}{\left(p_{S,\text{max}} - p_{S,\text{min}}\right)\left(p_v - p_{S,\text{min}}\right)} = \frac{p_v\left(p_v - p_{S,\text{min}}\right) - \frac{1}{2}\left(p_v^2 - p_{S,\text{min}}^2\right)}{\left(p_{S,\text{max}} - p_{S,\text{min}}\right)\left(p_v - p_{S,\text{min}}\right)}\nonumber\\ &= \frac{p_v\left(p_v - p_{S,\text{min}}\right) - \frac{1}{2}\left(p_v + p_{S,\text{min}}\right)\left(p_v - p_{S,\text{min}}\right)}{\left(p_{S,\text{max}} - p_{S,\text{min}}\right)\left(p_v - p_{S,\text{min}}\right)} = \frac{p_v - \frac{1}{2}\left(p_v + p_{S,\text{min}}\right)}{p_{S,\text{max}} - p_{S,\text{min}}} = \frac{p_v - p_{S,\text{min}}}{2\left(p_{S,\text{max}} - p_{S,\text{min}}\right)} = \frac{c_f}{2}. \end{align} \]
In the case $p_v>p_{S,\text{max}}$, one has
\[ \begin{align} c_f &= 1, \end{align} \]
i.e. the grid box is 100 percent covered with clouds. In this case it follows for $\xi$
\[ \begin{align} \xi &= \frac{\int_{p_{S,\text{min}}}^{p_{S,\text{max}}}p_v - p_S'dp_S'}{\left(p_{S,\text{max}} - p_{S,\text{min}}\right)\left(p_v - p_{S,\text{min}}\right)} = \frac{p_v\left(p_{S,\text{max}} - p_{S,\text{min}}\right) - \frac{1}{2}\left(p_{S,\text{max}}^2 - p_{S,\text{min}}^2\right)}{\left(p_{S,\text{max}} - p_{S,\text{min}}\right)\left(p_v - p_{S,\text{min}}\right)}\nonumber\\ &= \frac{p_v\left(p_{S,\text{max}} - p_{S,\text{min}}\right) - \frac{1}{2}\left(p_{S,\text{max}} + p_{S,\text{min}}\right)\left(p_{S,\text{max}} - p_{S,\text{min}}\right)}{\left(p_{S,\text{max}} - p_{S,\text{min}}\right)\left(p_v - p_{S,\text{min}}\right)} = \frac{p_v - \frac{1}{2}\left(p_{S,\text{max}} + p_{S,\text{min}}\right)}{p_v - p_{S,\text{min}}}. \end{align} \]
The previous derivation in this section refers to the formation of clouds (dealing with supersaturation, also known as saturation adjustment). The fact that only a certain proportion $c_f$ of the grid boxes is filled with clouds also has consequences for the interactions and evaporation of condensates. The approach is as follows:
Divide all condensate densities by $c_f$. If $c_f<0.1$, use 0.1 for stability. This is done because cloud physics only takes place within the cloud-filled region of the grid box and the densities there are higher by a factor of $1/c_f$.
Calculate the cloud physics with the rescaled densities.
Multiply the resulting source densities by $\text{max}\left(c_f;0,1\right)$ to obtain the tendencies for the averaged densities (which are the prognostic variables).
For the fraction $c_p$ of the grid box filled with precipitation, one makes the ansatz:
\[ \begin{align} c_p = \begin{cases} \frac{1}{3}, c_f = 0 \land \rho_s + \rho_r + \rho_g>0,\\ c_f, \:\text{otherwise}. \end{cases} \end{align} \]
It is therefore assumed that the fraction filled with precipitation is equal to the fraction filled with clouds, unless there is only precipitation but no clouds in the grid box. This is typically the case under clouds from which precipitation falls.
This scheme takes the following four classes of condensate into account:
liquid cloud water (index $l$)
cloud ice (index $i$)
liquid precipitation (index $r$)
solid precipitation (index $s$)
This scheme is suitable for global models. It is assumed that at and below $-35^\circ$ C the cloud water is completely frozen and at and above $-10^\circ$ C the cloud water is completely liquid and there is a linear transition in between. This can be justified by observations [35].
The size distributions of the condensates are relevant for calculating the fall velocities, the interactions with the radiation field and their mutual interactions. Therefore, a formula for the effective radius must be specified for each condensate class.
An effective radius of 12 micrometers is assumed for cloud droplets.
For the effective radius of cloud ice, the formula developed in [36] is one of the most common approximations. It reads:
\[ \begin{align} d_{\text{eff},i} &= \left(1,2351 + 0,0105\cdot\left(t - t_0\right)\right)\cdot\left(45,8966\cdot w^{0,2214} + 0,7957\cdot w^{0,2535}\cdot\left(t - 83,15\right)\right).\tag{36.44}\label{eq:sun_rikus} \end{align} \]
Here, $t\coloneqq T/\text{K}$ is the measure of the temperature in Kelvin, $t_0\coloneqq 273.15$, $w \coloneqq 1000\left(\rho_i + \rho_s\right)/\left(\text{kg}/\text{m}^3\right)$ is the measure of the cloud ice density (including snow) in $\text{g}/\text{m}^3$ and $d_\text{eff}$ the measure of the effective diameter in micrometers. The result of Eq. (36.44) is restricted in the following way:
\[ \begin{align} d_\text{eff} \to \text{max}\left(\text{min}\left(d_\text{eff},155\right),20 + 40\cdot\cos\left(\phi\right)\right) \end{align} \]
with $\phi$ as the geographical latitude. For the effective radius $r_\text{eff}$, the relation
\[ \begin{align} r_{\text{eff},i} = \frac{3\sqrt{3}}{8}d_{\text{eff},i}\cdot 10^{-6}\text{m} \end{align} \]
is derived in (36.44). The fact that the factor is not simply 1/2 is because ice crystals are generally not spherical.
It is assumed that the size distribution of the raindrops corresponds to a Marshall-Palmer distribution, Eq. (23.95).
It is assumed that the size distribution of the snowflakes corresponds to a Marshall-Palmer distribution, Eq. (23.95), according to Sect. 23.3.1.1.
The assumption of Marshall-Palmer distributions implies that the process of breaking up condensates, which limits the spectrum upwards, does not have to be explicitly simulated, but is already parameterized via the distribution.
To determine the fall velocities of the condensates, one follows Sect. 23.2.3, but with some simplifications:
A fall velocity of zero is assumed for cloud droplets.
A fall velocity is assumed for cloud ice.
For snow, the empirical formula
\[ \begin{align} v_S = a\cdot r^b = 9,356\frac{\text{m}^{0,7}}{\text{s}}\cdot\left(2\cdot r_\text{eff}\right)^{0,3} = 11,519\frac{\text{m}^{0,7}}{\text{s}}\cdot\left(r_\text{eff}\right)^{0,3}\tag{36.47}\label{eq:fall_snow} \end{align} \]
[16] is assumed.
There is turbulent flow around raindrops; for them, Eq. (23.90) with $c_w = 0.8$ is used.
If the fall velocity of a precipitation class is $v_{D,P}\left(D\right)$, it has a volume mixing ratio $r_P$ and a size spectrum $n\left(D\right)$, then a mass-weighted fall velocity can be calculated via
\[ \begin{align} v_P\left(D\right) = \frac{1}{r_P}\int_0^\infty\frac{\pi D^3}{6}v_{D,P}\left(D\right)n\left(D\right)dD \end{align} \]
Assuming that $v_{D,P}$ has the form
\[ \begin{align} v_{D,P}\left(D\right) = aD^b \end{align} \]
and that the Marshall-Palmer distribution, Eq. (23.95), holds, it follows
\[ \begin{align} v_P\left(D\right) &= \frac{1}{r_P}\int_0^\infty\frac{\pi D^3}{6}aD^bn_0e^{-\lambda D}dD = \frac{an_0\pi}{6r_P}\int_0^\infty D^{b + 3}e^{-\lambda D}dD = \frac{an_0\pi}{6r_P\lambda^{b + 4}}\int_0^\infty D^{b + 3}e^{-D}dD\nonumber\\ &= \frac{an_0\pi\Gamma\left(b + 4\right)}{6r_P\lambda^{b + 4}}. \end{align} \]
From Eq. (23.99), one obtains
\[ \begin{align} \frac{1}{\lambda^4} = \frac{r_P}{\pi n_0}. \end{align} \]
From this it follows
\[ \begin{align} v_P\left(D\right) &= \frac{a\Gamma\left(b + 4\right)}{6\lambda^b} \stackrel{\href{ch-22-cloud-microphysics.html#eq:marshall-palmer_d_lambda}{\text{Eq. (23.100)}}}{=} \frac{\Gamma\left(b + 4\right)}{6}aD^b. \end{align} \]
The simplest understanding of how rain occurs is to assume that it occurs when cloud water $\rho_w$ exceeds a critical mixing ratio $q_{w,\mathrm{krit}}$. A simple mathematical ansatz based on this for the source density of rainwater reads:
\[ \begin{align} Q_r = \frac{\text{max}\left(\rho_w - q_{w,\text{krit}}\cdot\rho_h,0\right)}{\tau}. \end{align} \]
This is the so-called Kessler approach [25]. Typical values are $q_{w,\mathrm{krit}} = 0.2\cdot 10^{-3}\mathrm{ kg}/\mathrm{m}^3$, $\tau = 1000\mathrm{ s}$.
The Kessler approach presented in Sect. 36.5.3.1 can also be applied to the formation of snow, typical values are $q_{i,\text{krit}} = 1.0\cdot 10^{-3}\text{ kg}/\text{m}^3$, $\tau = 1000\text{ s}$.
The Kessler approach presented in Sect. 36.5.3.1 can also be applied to the formation of graupel, typical values are $q_{s,\text{krit}} = 0.6\cdot 10^{-3}\text{ kg}/\text{m}^3$, $\tau = 1000\text{ s}$.
In this scheme, a fifth condensate class is added to the four-condensate scheme, namely graupel. The five condensate classes used here are:
liquid cloud water (index $l$),
cloud ice (index $i$),
liquid precipitation (index $r$),
snow (index $s$),
graupel (index $g$).
This scheme is suitable for regional models. It corresponds to the scheme presented in Sect. 36.5 with the following modifications:
An effective radius of 0.5 millimeters is assumed for graupel particles.
There is turbulent flow around graupel particles; for them, Eq. (23.90) with $c_w = 0.8$ and the replacement $\rho_l' \to \rho_i'$ is used.
In the schemes considered so far, each condensate class is described by a prognostic variable, namely the mass density or the specific density. This would only be a complete description if all condensates of the respective class had the same radius and temperature as the gas phase. However, neither is generally the case. Two extensions are available:
Each condensate class can be assigned its own temperature.
The size distribution within the condensate classes can be taken into account. This happens with the so-called multi-moment schemes, in which the higher moments of the size distributions (or comparable quantities describing the shape of the distributions) are introduced as additional prognostic variables.