1. Introduction
Bulk microphysical schemes are widely used to model cloud microphysical processes. For bulk parameterizations, the droplet size distribution (DSD) of each hydrometeor category is approximated by a continuous function for which there are one or more free parameters. Usually, from the DSD one or more prognostic moments can be calculated, for which predictive equations for microphysical processes are computed.
Initially, bulk schemes (Kessler, 1969) incorporated only a single moment (usually the third moment regarding the diameter of the DSD, proportional to the liquid water content). As the complexity of the models increased, the zeroth moment with respect to the drop diameter (which is equal to the total number concentration) was incorporated in two-moment bulk schemes (e.g., Murakami, 1990; Ferrier 1994; Reisner et al., 1998; Seifert and Beheng, 2001; Morrison et al., 2005; Cohard and Pinty, 2000).
For parameterized microphysics (e.g., Khain et al., 2015; Milbrandt and Yau, 2005a; Seifert and Beheng, 2006; Morrison et al., 2009; Lim and Hong, 2010), a gamma distribution is usually assumed for the drop size distribution:
where D is drop diameter (in cm), n 0 is the intercept of the distribution (in cm -4 ), µ is the shape parameter (non-dimensional) and λ is the slope parameter (in cm -1 ). The moments of order 0 and order 3 with respect to radius or diameter (which are equal and proportional to the concentration and the liquid water content respectively) are usually predicted.
By fixing the shape parameter µ, the computational burden of the bulk approach is drastically reduced, but other drawbacks arise. For example, while calculating the sedimentation some authors reported an excess sorting (e.g., Wacker and Seifert, 2001; Milbrandt and McTaggart-Cowan, 2010; Shipway and Hill, 2012).
Size sorting can be observed in a polydisperse population of droplets, due to the fact that larger drops have larger terminal velocities and settle much more quickly, resulting in a spatial separation of droplets. It can be modeled in two-moment schemes, as M 3 sediments much faster than M 0 . However, when fixing µ, there is an over-prediction of M 3 , and an under-prediction of M 0 (when comparing with an explicit reference model) generating the excess size sorting reported by some authors.
For a scheme with a fixed value of µ, the excess sorting occurs due to the fact that the ratio of the moment-weighted sedimentation velocities (V k /V j , with k >j), which is a function of the shape parameter, is always positive and larger than 1 (Milbrandt and McTaggart-Cowan, 2010). Therefore, M k has always larger sedimentation rates than M j . Then, as the ratio M 3 /M 0 increases, the mean radius will also increase at the lower edge of the sedimentation profile. This behavior is a consequence of assuming a constant value of µ and, implicitly, a prescribed constant ratio of the bulk fall velocities.
For two-moment schemes, this problem can be mitigated by using a fixed but large value of µ, consequently obtaining moment-weighted fall velocities for M 3 and M 0 closer in value. An alternative solution (without changing µ) is to make the ratio V k /V j closer to 1 by increasing V j . However, these are only palliatives as the mean radius will always increase at the leading edge of the sedimentation profile. Milbrandt and Yau (2005a) proposed a diagnosing relationship for the shape parameter as a function of the mean droplet diameter µ = f(D m ) to control the size sorting. Also, in a more general approach, the shape parameter was parameterized as a function of M 3 and M 0 . Additionally, diagnosing relationships were extracted from bin microphysical models by fitting the explicit droplet size distributions to gamma distributions in order to calculate the shape parameter.
Three-moment parameterization schemes were introduced in order to obtain a more physically based representation of the droplet size distribution evolution. For these schemes, an additional moment (commonly the sixth moment, M 6 , which is proportional to the radar reflectivity) is added in order to calculate the shape parameter µ. As a result, they are more able to approximate the bin reference models for pure sedimentation (Milbrandt and Yau, 2005a). For the sedimentation process, three-moment parameterization schemes performed better than two-moment schemes, because they predict a larger reflectivity-weighted fall velocity than the mass-weighted fall speed, resulting in a larger sedimentation rate for M 6 . As a result, there is an increase in the shape parameter µ during size sorting and, consequently, the droplet size distribution narrows, thus further limiting the size sorting as the weighted-fall speeds for the different moments are closer. This feedback is the reason for three-moment schemes that provide a more realistic description of the sedimentation process (Dawson et al., 2014).
There are few three-moment schemes that consider other processes besides sedimentation (Szyrmer et al., 2005; Shipway and Hill, 2012; Dawson et al., 2014; Loftus et al., 2014; Naumann and Seifert, 2016). In Naumann and Seifert (2016), a three-moment rain scheme that includes the processes of sedimentation, evaporation, and collision-coalescence, was compared with a Lagrangian model obtaining a good correspondence between the bulk and the explicit models. Paukert et al. (2019) developed a three-moment scheme that included various microphysical processes, such as sedimentation, evaporation, self-collection, and collisional breakup.
The aim of this study is to advance further in this direction through the development and assessment of a three-moment microphysical scheme that includes the processes of sedimentation and collision-coalescence. To quantify the impact of these processes on drop size distribution and shape parameter evolution, simulations were performed varying the shape parameter of the initial distribution. The results obtained with the three-moment microphysical scheme are in good agreement with the explicit Eulerian reference model. It was found that the collision-coalescence process counteracts the sedimentation tendency to create narrow droplet size distributions.
The paper is organized as follows: section 2 presents the three-moment parameterization scheme. The numerical implementation of the sedimentation model for both the bulk and the explicit schemes is presented in section 3. Section 4 is devoted to the analysis of simulation results. Concluding remarks are given in section 5.
2. The bulk three-moment parameterization scheme
2.1 Obtaining the droplet size distribution shape as a function of moments
The three-parameter gamma distribution function (Eq. 1) has been widely used in bulk parameterization schemes (e.g., Seifert and Beheng, 2001; Milbrandt and Yau, 2005b; Morrison et al., 2005; Milbrandt and McTaggart-Cowan, 2010; Ziemer and Wacker, 2014). The parameters n 0 and λ must be positive, while µ can also be negative. From distribution parameters, the value of any moment M (i) can be computed analytically from the expression:
The main objective of a three-moment parameterization scheme is to obtain the three parameters of the gamma distribution (Eq. [1]) from the prognostic moments. The intersection n 0 and the scale parameter λ of the distribution in Eq. (1) are calculated from the moments M 0 and M 3 by using Eqs (3) and (4) (Milbrandt and MacTaggart-Cowan, 2010):
and the shape parameter µ is obtained from the solution of the cubic equation (Paukert et al., 2019):
Following Paukert et al. (2019), we set µ = 0 for K > 20, and µ = 20 for K < 1.46; then µ Є [0, 20]. For K Є [1.46, 20]. the equation can be solved either analytically by using Cardano’s formula (Press et al., 1992) or numerically. At each time step, in order to calculate the three parameters, the prognostic moments need first to be updated due to microphysical processes (sedimentation and collision-coalescence). Then, the new parameters of the distribution are calculated from Eqs. (2)-(5).
2.2 Updating the moments due to collision-coalescence process
Within the bulk approach, the DSD is decomposed in two parts. Drops smaller that a threshold radius, that is typically within a range from 20 µm (Khairoutdinov and Kogan, 2000; Wood and Blossey, 2005) to 41 µm (Cohard and Pinty, 2000; Beheng, 2010), are called cloud droplets, and larger droplets are called raindrops. In this paper, a threshold diameter of D = 82 µm was adopted following Cohard and Pinty (2000).
The following interactions between cloud droplets and raindrops, due to collision-coalescence, are assumed: cloud droplet-cloud droplet collisions, cloud droplet-raindrop collisions, and raindrop-raindrop collisions (Lee and Baik, 2017). For cloud droplet-cloud droplet interactions, two processes can be identified: autoconversion (if the resulting drop is a raindrop) and self-collection (if the formed drop is a cloud droplet). For cloud droplet-raindrop interactions and raindrop-raindrop collisions, the accretion and the self-collection processes can be identified. For both these latter processes, the result is a raindrop.
The tendencies for the number concentration, mass mixing ratio and radar reflectivity for both cloud droplets (N c , Q c , Z c ) and raindrops(N r , Q r , Z r ) due to autoconversion, accretion and self-collection are:
The moment tendencies are calculated from the former equations by noticing that N iT = M i0 , LWC = (6/µρ L )M 3 and Z i = M i6 . The equations for each process are described in detail in section S1 of the supplementary material. The changes of the reflectivity due to microphysical processes (autoconversion, accretion and self-collection) are calculated following Milbrandt and Yau (2005b), and only Type 1 tendency equations are considered. Type 2 and Type 3 tendencies are not considered in our model, as Type 2 tendencies represent changes in radar reflectivity when a new hydrometeor is initiated, and Type 3 tendencies represent conversion from one hydrometeor category to another (Milbrandt and Yau, 2005b). Type 1 tendency equations for the radar reflectivity can be obtained by taking the derivative of Eq. (12), considering that the shape parameter µ i is constant.
In Eqs. (12) and (13), i stands again for cloud or rain, c i = ρ i (µ/6), and ρ a is the density of air. Then, the reflectivity rates due to autoconversion, accretion and self-collection are calculated from the equations (Milbrandt and Yau, 2005b):
2.3 Implementation of the bulk model for sedimentation
For the sedimentation process, the parameterized microphysics consists of a system of three budget equations for each moment of the DSD:
where M k is the moment of the DSD, k is the order of the moment (k = 0.3 and 6 with respect to the drop diameter) and V k is the moment-weighted sedimentation velocity that is calculated from the equation:
In Eq. (21) Γ is the gamma function, and λ and µ are the slope and the shape parameter of the gamma distribution, respectively. For the numerical solution of Eq. (20), 80 vertical layers with Δz = 100m were defined. All equations are integrated in time by using a first-order Euler forward-scheme with a time step of Δt = 1s.
3. The bin microphysics reference model and methodology of comparison with the bulk microphysics scheme
3.1 The bin microphysics reference model
The spectral bin model, which is used as a reference, solves the partial differential equation:
where f(z, t, D) is the DSD (as a function of height z, time and droplet diameter D), V(D) is the terminal velocity, and (∂f(z, D)/∂t)ǀ coagulation is the source term due to the coagulation process defined by the kinetic collection equation (KCE). This equation, in its formulation for a size distribution f(x, t) with drop mass x (Pruppacher and Klett, 1997) has the form:
where K(x, y) is the hydrodynamic kernel, which is a symmetric function of the mass of the colliding droplets
where r x and r y are the radii of droplets with masses x and y, respectively, and E(r x ,r y ) are the Hall (1980) collection efficiencies. The KCE was solved using the flux method developed by Bott (1998), and a drop size range from 1 to 2500 µm was used during the simulations. The size distribution f(x) was defined the same as in Bott (1998) and Berry (1967), and is represented by 33 mass doubling categories, then mass M k in the category k is determined as M k = 2M k-1 .
The terminal velocity V (D) in Eq. (24) is given as a power-law relationship (Straka, 2009):
where a = 1300 cm 0.5 s -1 , b=0.5 (Gunn and Kinzer, 1949) and (ρ 0 /ρ) is the air density correction factor (with ρ and ρ 0 denoting the air density aloft and at the surface, respectively) that, for simplicity, is assumed to be 1 throughout this study. For the numerical solution of Eq. (22), a forward-in-time and upstream-in-space method was implemented with a time step of Δt = 1 s. The value of b is the same for Eqs. (21) and (25).
3.2 Methodology of comparison with the bulk microphysics scheme
As the size distribution f(x) for the bin microphysical model was defined the same as in Bott (1998), it has to be initialized from the condition (Berry, 1967)
Then, given a three-parameter gamma distribution N(D) with parameters n 0 , µ and λ, the corresponding droplet mass distribution f(x) has the form (with droplet mass in grams, and considering the water density ρ w =1):
Then, the mass distribution for the bin microphysical model must be initialized from Eq. (27).
4. Rainshaft model setup and simulation results
The performance of the three-moment parameterization was tested within a horizontally homogeneous environment (rainshaft model). Our computational domain is an 8 km height vertical column, discretized with 80 grid points with a spacing of Δz = 100 m, and the time step was set equal to Δt = 1 s. This simplified setup is very useful to evaluate microphysical processes and has been used by other authors (Seifert and Beheng, 2001; Ziemer and Wacker, 2014).
The only microphysical processes considered in our study are sedimentation and collision-coalescence, with the configurations SBM-S, SBM-SC, M3-S, M3-SC and M2-S (see Table I for the definitions).
Configuration | Description |
SBM-S | Spectral-bin microphysics with sedimentation. |
SBM-SC | Spectral-bin microphysics with sedimentation and collision-coalescence. |
M3-S | Three-moment (M3) parameterization scheme with sedimentation determining µ. |
M3-SC | Three-moment (M3) parameterization scheme with sedimentation and collision-coalescence determining µ. |
M2-S | Two-moment (M2) parameterization scheme with sedimentation and a fixed value of µ. |
4.1. Pure sedimentation case
The performance of three-moment schemes for pure sedimentation was analyzed in previous papers by Milbrandt and Yau (2005a), Milbrandt and McTaggart-Cowan (2010), and Ziemer and Wacker (2014). In this section, we also perform a comparison between M3-S, M2-S, and SBM-S before addressing the combined influence of collision-coalescence and sedimentation processes. For all the simulations, an initial 1.5 km thickness maritime cloud (which lies between 6000 and 7500 m) was assumed. Three different initial configurations were considered (with initial parameters taken from Ziemer and Wacker [2014]), with a liquid water content of L wc = 5 × 10-7 g cm-3, and a cloud droplet number concentration of n 0 = 3 × 10-3 cm-3. Only the reflectivity was varied (see values in Table II) in order to obtain different initial distributions with different shape parameters for the three schemes (two-moment, three-moment and the bin-scheme).
Z (cm3) | Shape parameter (µ) | n 0 | λ |
6.0793 × 10-9 | 0 | 7.9840 × 10-2 | 26.6134 |
3.7257 × 10-9 | 0.5 | 6.8793 × 10-1 | 34.5475 |
9.1052 × 10-10 | 4.8773 | 1.7501 × 107 | 100.0092 |
The shape parameters µ calculated from these initial configurations were found equal to 0, 0.5, and 4.8773, respectively (see Table II), and serve us to initialize both the bin and three-moment parameterization schemes, and as the shape parameter (that has been fixed to these values during the entire simulations) for the two-moment parameterization scheme.
The results obtained for the three experiments (with µ = 0, 0.5, and 4.8773) are displayed in Figures 1, 2, and 3 for three different times (200, 400, and 600 s). For the three experiments, the prognostic moments droplet concentration N, liquid water content L wc , and reflectivity Z for the M3-S model, perform much better than the M2-S model, which has a fixed shape parameter.
For µ = 0 there is a good coincidence between the three-moment and bin schemes for the prognostic moments liquid water content and reflectivity. However, the three-moment scheme slightly overestimates the maximum value for droplet number concentration at all times. For that case (µ = 0), the two-moment scheme overestimates the droplet concentration and underestimates the prognostic moments liquid water content and the reflectivity.
For the second experiment (with shape parameter µ = 0.5, Fig. 2), the M3-S model seems to slightly overestimate the maximum value for the number concentration at t = 600 s, but there is a good match between the three-moment and the spectral bin microphysical schemes for the prognostic moments liquid water content and the reflectivity. For the M2-S model, on the other hand, there is a marked underestimation of the prognostic moments liquid water content (L wc ) and reflectivity (Z), and an overestimation of the prognostic moment number concentration.
For the case with an initial distribution with µ = 4.8773, we saw that all the prognostic moments were well captured by the three-moment scheme. The M2-S model slightly underestimated the radar reflectivity at t = 600 s. We can conclude that the M3-S model, for all the prognostic moments (N, L, and Z) performs very well for narrower distributions (with µ = 0.5 and 4.8773) and outperforms the two-moment parameterization scheme.
For the two-moment scheme, the liquid water content L wc and reflectivity Z sediment faster than the M 0 . This is confirmed by the medium radius vertical profile, with much larger values than the reference spectral model. For the three-moment scheme, on the other hand, a good coincidence between the mean radius profile for the spectral reference and the parameterization was obtained for different values of the initial shape parameter of the DSD
As can be observed in Figures 4, 5, and 6, both the two-moment and three-moment schemes produce a size sorting effect; however, the two-moment scheme is not very accurate at reproducing the medium radius profile, while the three-moment scheme performs very well on predicting the medium radius at all heights for the three cases (with µ = 0, 0.5 and 4.8773), for all simulation times.
The cause of this discrepancy for two-moment schemes was already discussed in previous studies and outlined in the introduction (e.g., Milbrandt and MacTaggart-Cowan, 2010), and has its origin in the fact that for a scheme with a fixed value µ, the ratio of the moment-weighted sedimentation velocities (which is actually a function of the shape parameter) is always positive and larger than 1. The fact that a good correspondence was obtained between the three-moment and the spectral reference solution confirms that three-moment schemes give a physically based and more complete representation of the sedimentation processes.
Vertical profiles of shape parameter µ for three times (200, 400, and 600 s), and for three different initial values (µ = 0, 0.5, and 4.8773) obtained from the M3-S model, are compared with those from SBM-S, which serves as a benchmark (Figs. 7, 8, and 9). The figures were obtained for cloud droplet concentrations larger than 10-6 cm-3.
For this comparison, it is assumed that the DSD for the bin model (which evolves without restrictions for the bin microphysical model) follows a gamma distribution (Paukert et al., 2019). As can be observed, the shape parameter profiles from the parameterized model follow very closely the shape parameter obtained from the bin model under the assumption that the DSD is gamma. This leads us to the conclusion that the gamma distribution works well as an approximation of the DSD for the sedimentation case and that the bulk M3-S model is able to capture the evolution of the shape parameter.
To assess the impact of the sedimentation process on the DSD shape parameter, contour plots for two different initial values of the shape parameter (µ = 0.5 and 4.8773) were calculated (Figs. 10 and 11). As can be checked in these figures, there is an increase in the shape parameter as cloud height decreases, and consequently (as the shape parameter µ is related to the relative dispersion ε of the DSD from the relation µ = ε -2 -1) a decrease of the relative dispersion. The increase of µ is a result of size sorting, which tends to make DSD much narrower. For the two-moment scheme (with a constant value of the shape parameter), there is an excess size sorting and consequently an overestimation of the mean radius (Fig. 4, left panel). For the three-moment scheme, the DSD tends to be narrower, and the model reproduces quite accurately the mean radii obtained with the explicit model. We can conclude that the M3-S model (with variable shape parameter) captures well the narrowing size distribution resulting from size sorting (Fig. 4, right panel).
4.2 Combined effect of sedimentation and collision coalescence
Two simulations were performed with the configurations SBM-SC (spectral-bin microphysics with sedimentation and collision-coalescence) and M3-SC (three-moment parameterization scheme with sedimentation and collision-coalescence determining µ), the former serving as a benchmark. As for the pure sedimentation case, a 1.5 km thickness maritime cloud (which lies between 6000 and 7500 m) was assumed, with a cloud liquid water content of L wc = 2 × 10-6 g cm-3, a cloud drop number concentration of N 0c = 100 cm-3, and a reflectivity of Z 0c = 3.8213 × 10-13 cm3 (see Table III). For this combination of distribution moments, the initial value of the shape parameter was found equal to µ = 5.99. For the two experiments, the simulation time was set equal to t = 1800 sec.
Symbol | Description | Values |
L wc | Initial cloud liquid water content | 2 × 10-6 g cm-3 |
N 0c | Initial cloud droplet number concentration | 102 cm-3 |
Z 0c | Initial cloud reflectivity | 3.8213 × 10-13 cm3 |
The comparison between the SBM-SC and M3-SC cases was performed by plotting the evolution in time of the prognostic moments (concentration, liquid water content, and reflectivity). As can be observed in Figures 12 and 13, in general, there is a good agreement between the parameterized and the bin microphysics models.
The time-height distributions for cloud water drop number concentrations, cloud liquid water content, and radar reflectivity (Fig. 12) are very similar. For all cases, the slope in all the figures indicates that the sedimentation process is activated, and the hydrometeors (cloud and rain drops) are falling to the ground. At t = 0, we have the same distribution, and constant values between 6000 and 7500 m. As time evolves, the drop number concentration decreases due to the combined effect of the collision-coalescence and sedimentation processes.
However, a more critical analysis of Figs. 12 and 13 reveals some differences between the SBM-SC and the M3-SC models. Compared to SBM-SC, the three-moment scheme exhibits a slightly larger coalescence rate, a fact that becomes clear when observing the values of rain number concentration and liquid water content at the same time. For example, at t = 600 s, the parameterized scheme exhibits values of rain number concentration larger than 0.20 cm-3, while for the explicit model, the concentration values are in the order of 0.14 cm-3. Accordingly, the same behavior is observed for the rain liquid water content, with values of 5.5 × 10-9 g cm-3 and 5×10-9 g cm-3 for the M3-SC and the SBM-SC models, respectively.
These results for rainwater are consistent with a faster rate of decrease in cloud droplet concentration for the parameterized model. For example, at t = 1000 s for the explicit model, we can find cloud droplet concentrations as large as 90 cm-3 at 3400 m, while for the parameterized model cloud droplet concentrations are smaller than 90 cm-3 at all cloud levels. For cloud liquid water content and cloud reflectivity, in general, the three-moment parameterized model emulates very well the results obtained with the bin model, especially for cloud water reflectivity, with very similar time-height profiles for the two models.
Figures 12 and 13 show that the results obtained with the three-moment scheme for this case (with the combined effect of sedimentation and collision coalescence) are found to agree well with those obtained with the explicit model. Some differences are unavoidable, due to the complexity of the schemes and the incorporation of other processes. Overall, the model reproduces quite accurately the results obtained with the explicit model.
5. Discussion and conclusions
In this paper, a bulk three-moment parameterization scheme incorporating sedimentation and collision-coalescence was developed and evaluated. The performance of the parameterized scheme was assessed through a detailed comparison with a bin microphysical scheme within a one-dimensional kinematic setting. For the comparison, the predicted moments were number density N, liquid water content L wc , and radar reflectivity Z.
In our simulations, for the pure sedimentation case, results from the three-moment parameterized scheme are found to agree well with those from simulation using the bin microphysical scheme (for the three predicted moments) for initial values of the shape parameter µ = 0.5 and 4.8773. For µ = 0, the three-moment scheme slightly overestimates the maximum value for droplet concentration at all times. Better results are obtained for narrower initial distributions (with larger values of the shape parameter). The results obtained justified the adoption of the new procedure for calculating the shape parameter µ outlined in Paukert et al. (2019).
The mean radius vertical profiles calculated from the three-moment parameterized scheme match very well those obtained from simulations using the bin microphysical scheme, confirming (as discussed previously by Milbrandt and Yau [2005a], Milbrandt and McTaggart-Cowan [2010], and Ziemer and Wacker [2014]), that three-moment parameterization schemes outperform the two-moment schemes, and that the inclusion of a third moment in order to predict µ gives a more realistic description of the DSD evolution due to sedimentation.
When the collision-coalescence process is incorporated, overall, results from the three-moment scheme are found to agree qualitatively well with those obtained from simulation using the explicit scheme. However, the onset of precipitation occurs earlier in the M3-SC model, with a clear overestimation of the raindrop number concentration and the rain liquid water content.
The small differences between the prognostic moments profiles for the SBM-SC and the M3-SC models are clearly a result of the incorporation of the collision-coalescence process. Overall, results from the three-moment parameterized scheme are found to agree well with those from simulations using the bin microphysical scheme. However, we were unable to find the level of proximity to the reference solution that was obtained for the pure sedimentation case.
This can be explained by the complexity of the Kessler-type parameterizations that involved six prognostic moments (N c , Q c , Z c , and N r , Q r , Z r ) and three microphysical processes (autoconversion, accretion, and self-collection). Also, when developing parameterizations for the collection process (Cohard and Pinty, 2000), the polynomial form of the collection kernel developed by Long (1974) is considered to obtain analytical expressions for the integrals of the KCE (in the form of gamma functions). Then, the results will differ from those obtained from the integration of the KCE with the hydrodynamic kernel Eq. (24).
A further effort to improve the analytical approach (that uses an approximated form of the collection kernel) could lie in the adoption of the machine learning (ML) approach, which was used for the autoconversion process with good results by Alfonso and Zamora (2021). For the Kessler-type parameterizations, it could be extended in order to include the processes of accretion and self-collection, and embedded into a dynamic framework in order to check the performance through a direct comparison with a spectral model.