Parameterization and Explicit Modeling of Cloud Microphysics: Approaches, Challenges, and Future Directions

Clouds, precipitation, radiation, and their interactions are some of the most critical yet uncertain aspects of weather and climate prediction (Yu et al., 2019; Meehl et al., 2020), because related processes and their interactions are complex, nonlinear, and multiscale in nature. Microphysical processes occur at the smallest end of the multiscale hierarchy (Liu, 2019) and thus must be parameterized not only in large-scale global circulation models (GCMs) used in weather and climate prediction but also in various higher-resolution limited-area models such as cloud-resolving models (CRMs) and large-eddy simulation (LES) models. Accurate microphysical parameterizations have recently proven to also be essential to improving solar energy forecasting (Jimenez et al., 2016; Jiménez et al., 2016), a relatively new research area with growing potential in contributing to a sustainable economy and society that harnesses renewable sources of energy (Liu et al., 2021, 2022). Compared to weather and climate prediction, solar forecasting poses further challenges because it requires not only the total solar irradiance, as in weather and climate models, but also the detailed partition of total irradiance into direct and diffuse irradiances, and at much shorter time horizons (e.g., minutes). The total horizontal irradiance (GHI) is key to solar panels, while the direct normal irradiance (DNI) is essential to concentrating solar technologies (Weinstein et al., 2015; Xie et al., 2020). The need for improving microphysical parameterizations is especially acute for investigating mesoscale convective systems and high-impact extreme events (Cao et al., 2022).

This review is devoted to bulk microphysics parameterizations (BMPs) and cloud modeling with explicit cloud microphysics, with an emphasis on the general themes, understudied topics, outstanding challenges, emerging areas, and future directions. More general discussions on related topics can be found in other reviews (e.g., Khain et al., 2015; Morrison et al., 2020), two excellent books (Stensrud, 2007; Straka, 2009), and the special issue of J. Geophys. Res. on fast physics in large-scale models (https://agupubs.onlinelibrary.wiley.com/doi/toc/10.1002/(ISSN)21698996.FASTPHYS).

Specially, after briefly describing microphysical processes and the variety of hydrometeor types involved in formation, evolution, dissipation, and interactions in section 2, section 3 discusses BMPs, which have been the main way to represent various cloud microphysical processes, with an emphasis on general characterization, main approaches to obtaining BMPs, theoretical foundations, outstanding challenges, and future directions. Section 4 introduces and compares the two types of explicit microphysical models that resolve the details of hydrometeor sizes and hydrometeor types and have been often used as bench-marking tools to facilitate development and evaluation of BMPs: spectral bin microphysics (bin microphysics hereafter) and Lagrangian particle-based microphysics (particle-based microphysics hereafter). Traditionally, bin microphysics has been mainly used to represent explicit cloud microphysics in research with CRMs or LES models due to its high computational cost (Khain and Pinsky, 2018). Particle-based models have been recently emerging as a promising alternative to the traditional bin microphysics for representing explicit cloud microphysics (Grabowski et al., 2019). Section 5 concentrates on particle-resolved direct numerical simulation (DNS) models, which not only resolve the smallest turbulent eddies but also track individual cloud particles, with applications to understand turbulence–microphysics interactions. Section 6 provides a summary and outlook on future research.

Clouds in the atmosphere can be grossly classified into liquid clouds and ice clouds according to the phase of water in clouds, with mixed-phase clouds being included in the category of ice clouds for convenience. Depending on the phase and stage of clouds, multiple microphysical processes can occur and interact, generating hydrometeors of different sizes and of different types. As briefly illustrated in Fig. 1, liquid phase processes include activation of cloud condensation nuclei (CCN) into cloud droplets; diffusional growth/shrinkage of cloud droplets and raindrops by condensation/evaporation; further drop growth by collision-coalescence (collection) process; spontaneous and collisional breakups of raindrops; and differential sedimentation of drops of different sizes. Ice-phase processes are more complicated, including various mechanisms of formation of ice crystals (e.g., homogeneous freezing of cloud drops into ice crystals, contact freezing, sorption and immersion freezing, and deposition nucleation); diffusional growth/shrinkage of ice crystals by deposition/sublimation; further aggregation growth of ice crystals into snow aggregates; riming and density change of ice particles; graupel initiation by riming process; graupel and frozen drops becoming hail embryos by riming; wet and dry growth of hailstones; melting of solid phase hydrometeors; Bergeron–Findeisen process in mixed-phase clouds; and secondary production of ice crystals (ice multiplication) by rime-splintering, ice fracturing, and shedding during hail growth. These individual processes do not occur in isolation; instead, they interact closely with one another to determine the properties of a hydrometeor population. Note that the processes listed in the yellow box influence all the other processes (e.g., turbulence and sedimentation). More detailed discussions on these microphysical processes can be found in textbooks on cloud physics (e.g., Roger and Yau, 1989; Pruppacher and Klett, 1997).

Figure 1.  Illustration of main microphysical processes. The processes in the yellow box likely occur or affect all the other processes.

Without loss of generality, a hydrometeor population is characterized by two fundamental properties: number concentration distribution with respect to hydrometeor sizes (hydrometeor size distributions, HSD hereafter) and hydrometeor types. For liquid-phase clouds, hydrometeors are largely homogeneous water spheres and can be classified into small cloud droplets and large raindrops according to the particle diameter or radius. Solid hydrometeors are more complicated, with multiple ice types characterized by different particle morphologies, different densities, and varying riming degrees. A cloud microphysics parameterization essentially involves developing adequate representation of HSD and hydrometeor type in association with the microphysical processes.

Bulk microphysics parameterizations (BMPs) seek to simplify the representation of microphysical processes and their associated HSDs and hydrometer types with a minimal number of prognostic quantities/parameters. Kessler (1969) pioneered the development of BMPs for warm rain processes by partitioning the whole liquid drop population into small cloud droplets and large raindrops. He assumed that the initial formation of embryonic raindrops was represented through the so-called autoconversion process that is further parameterized as a linear function of cloud liquid water content. Raindrops follow the Marshall–Palmer negative exponential size distribution and grow by collecting cloud droplets through continuous accretion. The Kessler BMP is a simple one-moment scheme for warm rain processes that only predicts liquid water mixing ratio. Subsequent studies have applied similar ideas to develop BMPs for ice processes (Wisner et al., 1972;; Cotton et al., 1982; Lin et al., 1983; Rutledge and Hobbs, 1984). Early ice BMPs assume that all hydrometeors are of spherical shapes, and these BMPs simply partition the complex solid hydrometeor types into ice and snow (e.g., Lin et al., 1983).

The early one-moment BMPs have been extended to predict two moments (zeroth moment for number concentration in addition to water mixing ratio) since the late 1980’s and early 1990’s (e.g., Hu and He, 1987; Murakami, 1990; Ferrier, 1994; Meyers et al., 1997; Reisner et al., 1998; Khairoutdinov and Kogan, 2000; Seifert and Beheng, 2001; Lou et al., 2003; Chen and Liu, 2004; Milbrandt and Yau, 2005a; Morrison et al., 2005; Morrison and Grabowski, 2007; Phillips et al., 2007; Lim and Hong, 2010; Thompson and Eidhammer, 2014 ). Motivation for the development of two-moment BMPs stemmed from scientific interests in aerosol indirect effects on climate and climate change (Twomey, 1977; Albrecht, 1989), together with the recognition of one-moment BMPs suffering from obvious shortcomings, and growing computational power. The prediction of both number concentration and water mixing ratio permits effective radius and other cloud properties being represented in a more realistic manner, allowing for better consideration of aerosol–cloud interactions compared to one-moment BMPs that only predict water mixing ratio. Such two-moment BMPs have been increasingly used in CRMs and LES models (e.g., Ferrier, 1994; Meyers et al., 1997; Seifert and Beheng, 2001; Morrison et al, 2005) and GCMs (e.g., Ghan et al., 1997; Lohmann et al., 1999; Ming et al., 2007; Morrison and Gettelman, 2008). Although two-moment BMPs still dominate as of now, development and use of three-moment BMPs have started since 2005 and are the focus of this paper (see section 3.2 for detailed discussion). Considering more moments in BMPs improves the process representation in principle on one hand, and on the other hand, it poses the necessity to consider more processes/quantities and thus increases the computational cost at the same time. Take liquid-phase clouds as an example. According to the idea of partitioning the total drop population into cloud droplets and raindrops, four different conversion mechanisms are involved with the collection equation: (1) autoconversion process as the formation of embryonic raindrops by collection of cloud droplets among themselves; (2) accretion as the process of raindrops collecting cloud droplets; (3) self-collection of cloud droplets whereby the resultant droplets remain in the category of cloud droplets; and (4) self-collection of raindrops whereby the resultant drops remain in the category of raindrops. Only the processes of autoconversion and accretion need to be considered in predicting changes of water mixing ratios in one-moment BMPs since self-collections of cloud droplets and raindrops do not affect the corresponding water contents. However, self-collections change the number concentration and spectral shape of the corresponding hydrometeor populations and thus need to be seriously considered in multi-moment BMPs (Beheng, 1994).

Another extension of BMP development is concerned with the treatment of hydrometeor types. To represent the multiple hydrometeor types within the constraint of computational resources, most existing BMPs use a limited number of hydrometeor categories defined by a set of prescribed physical hydrometeor characteristics (e.g., morphology, bulk density, and terminal velocity) that broadly describe a given “typical” hydrometeor type. Most liquid phase BMPs have followed Kessler (1969) to partition the whole liquid phase into two categories of cloud droplets and raindrops, probably with the exception reported in Saleeby et al. (2015) and related studies whereby a second mode of large cloud droplets with diameters from 40 μm to 80 μm is added. The main complexity of representing hydrometeor types lies with solid-phase hydrometeors, with the trend of an ever-increasing number of ice categories being considered in BMPs. For example, the simplest BMPs assume two solid categories of ice and snow analogous to liquid clouds. A rimed ice category (“graupel” or “hail”) was added in Lin et al. (1983) and Rutledge and Hobbs (1984). Other BMPs considered graupel and hail as separate categories to allow for both slower- and faster-falling rimed ice particles (Ferrier, 1994; Milbrandt and Yau, 2005b; Mansell et al., 2010). Walko et al. (1995) and Meyers et al. (1997) further added separate snow and aggregates categories. Straka and Mansell (2005) developed a BMP with 10 different categories of solid particles: column ice crystal, plate ice crystal, rimed cloud ice, snow (ice crystal aggregates), three categories of graupel with different densities and intercepts, frozen drops, small hail, and large hail.

Regardless of their specific differences, BMPs have mostly followed a paradigm of two pillars. The first pillar assumes some analytical HSD function such as the three-parameter Gamma function to close the set of equations for predicting HSD moments. The second pillar assumes that the diverse variety of hydrometeors can be grouped into several hydrometeor types characterized by a set of relationships (e.g., power-laws) between hydrometeor properties (e.g., mass, area, mass density, and terminal velocity) and hydrometer maximum diameter. Note that maximum dimension has been commonly used to describe non-spherical hydrometeors in literature; but maximum radius is used throughout this paper to avoid confusion with fractal dimension discussed in section 3.5. Figure 2 presents a BMP complexity diagram based on this two-pillar paradigm to illustrate the different levels of BMP complexity. In this diagram, the x-axis and y-axis denote the number of HSD moments and the number of considered hydrometeor categories, respectively. Accordingly, we recommend a BMP notation system “M_iC_j” to classify a specific BMP wherein the subscripts “i” and “j” denote the number of prognostic moments and the number of hydrometeor categories, respectively. Use of this classification system facilitates a first-order differentiation among the ever-increasing number of BMPs in literature. Table 1 provides some examples of BMPs available in the WRF4.1.2 model that are coupled with the radiative scheme to provide key microphysical properties (e.g., liquid water path and effective radius) for calculating cloud radiative properties.

Figure 2.  Classification diagram to illustrate the complexities of bulk microphysics schemes. A bulk scheme is characterized by the number of distribution moments (x-axis) and number of hydrometeor categories (y-axis). The notation M_xC_y is used in representing the microphysics schemes. The scheme complexity increases with either or both numbers, with the extremes approaching the explicit microphysics and particle-resolved DNS discussed in sections 4 and 5, respectively. Note that classifications of the exemplary schemes here should be viewed as approximate because not all the schemes consider the same hydrometer categories with the same number of moments.

Three extreme “BMP” scenarios are worth mentioning. First, when the number of prognostic HSD moments is sufficiently high, the moment-based BMPs become essentially equivalent to explicit bin microphysics according to the equivalence between the HSD moment and bin discretization:

where n(r) denotes the HSD; r is the hydrometeor radius; and $ {M}_{p} $ is the p-th HSD moment. This extreme BMP is denoted by i = ∞ in the BMP classification system. Note that there exist mathematical situations where the size distributions are not uniquely determined by their moments, i.e., the so-called dissimilar iso-momental distributions (White, 1990); however, they are less of a concern in reality (McGraw et al., 1998). Second, there has been a recently emerging fundamental change in representing hydrometeor types that moves from the discrete classification of different hydrometeor types to representation of continuous hydrometeor properties (hereafter “continuous” representation for convenience; see section 3.4 for more discussions). Likewise, this continuous representation of hydrometeor types is denoted by j = ∞ in the BMP classification system. Evidently, the ultimate representation of cloud microphysics is a combination of explicit microphysics with “continuous” hydrometeor type, or explicit “continuous” microphysics for short.

These extreme scenarios of explicit microphysics are embodied in the schemes of bin microphysics and particle-based microphysics detailed in section 4, and further in the particle-resolved direct numerical simulation (PR-DNS) models detailed section 5. In addition, discussions on plausible theoretical foundations buttressing BMPs are deferred to section 3.5.

Similar to Milbrandt and Yau (2005a, b), the prognostic equation for the pth HSD moment of the jth hydrometeor category can be generically written as

where u is the 3D air velocity vector, and the terms on the right-hand side of Eq. (2) represent, respectively, advection/divergence, turbulent mixing representing the subgrid turbulent transport (note its difference from the turbulent entrainment-mixing processes discussed in section 3.3), sedimentation, and microphysical sources (or sinks when it is negative). The source term is computed as a sum of the tendencies resulting from individual microphysical processes. The p-th-moment weighted terminal velocity is given by:

For three-moment BMPs, the three common prognostic moments are number concentration (p = 0), water mixing ratio (p = 3), and radar reflectivity (p = 6), and the corresponding equations are often closed by assuming that the HSD follows a three-parameter analytical HSD (e.g., Gamma distribution n(r) = N₀r^μexp(−λr), Weibull distribution, or lognormal distribution) and relating the prognostic moments to the HSD parameters (see Appendix for the different forms of analytical HSDs and the corresponding relationships). For two-moment BMPs, the spectral shape parameter such as μ in the Gamma HSD is assumed as a prespecified constant. The addition of M₆ allows for explicit consideration of the effect of HSD spectral shape in principle, which is neglected in both one-moment and two-moment BMPs. In practice, to parameterize the effects of HSD spectral shape on individual microphysical processes in BMPs is not trivial, and existing studies are limited in scope, some of which are highlighted next.

Raindrop sedimentation is arguably one of the first topics on which the effects of spectral shape have been extensively investigated in efforts to address the issue of size sorting that results from large drops falling faster than small ones. It is obvious that one-moment BMPs are not able to catch the size sorting phenomenon due to using a single mass-weighted sedimentation. Wacker and Seifert (2001) mathematically analyzed the equations of raindrop sedimentation resulting from treatments of a one-moment BMP, a two-moment BMP, and a reference bin microphysics model, and investigated the resultant differences in vertical profiles of rainwater and surface precipitation. They showed that the two-moment BMP, while producing a solution closer to the reference bin microphysics than the one-moment BMP, overestimated sedimentation of large raindrops due to the assumed inverse-exponential HSD. They also found that the treatment of moments transformed the original linear equation for bin microphysics to quasi-linear moment equations, resulting in spurious shockwaves. Using a similar approach, Milbrandt and Yau (2005a) showed that the two-moment BMP could mimic the effects of size sorting when distinct number-weighted and mass-weighted mean fall speeds are used to sediment M₀ and M₃, respectively. However, when the HSD is constrained to an inverse-exponential function, the differential sedimentations led to excessive size sorting in the two-moment BMP. They also examined the limitation of an inverse exponential HSD by using a three-parameter Gamma HSD function, by comparing different fixed values of μ, and by diagnosing μ through an empirical relationship derived from a bin microphysics model. Milbrandt and Yau (2005b) further extended the study to include a prognostic equation for M₆ (radar reflectivity) (thus μ) and developed the first three-moment BMP. Their results demonstrated that by allowing μ to vary, either diagnostically or prognostically, the problems of excessive size sorting associated with two-moment BMPs became manageable, because the difference between the number- and mass-weighted fall speeds increased with broadening HSDs. It is evident that this important feedback between the spectral shape and different moment-weighted fall speeds cannot be adequately represented with one- or two-moment BMPs.

Although most BMPs assume that cloud droplets do not fall, some recent modeling studies revealed that cloud droplet sedimentation plays a crucial role in determining cloud properties as well, through affecting cloud-top entrainment and droplet evaporation. For example, Ackerman et al. (2004) found, by simulating marine stratocumulus clouds with an LES with bin-microphysics, that decreases in sedimentation and precipitation due to increased droplet concentration (thus reduced droplet sedimentation) enhanced cloud-top entrainment and suppressed surface precipitation, and that the response of cloud water to increased droplet concentrations was determined by competition between moistening from suppressed surface precipitation and drying from increased entrainment of overlying air. When the overlying air was dry, the response of cloud water to aerosol-induced suppression of precipitation was dominated by enhanced entrainment of overlying dry air, reducing cloud water as droplet concentration increases. A complete suppression of sedimentation accelerated entrainment, drying the boundary layer, and thinning the cloud layer as the cloud base rises faster than the cloud top. They also pointed out that the traditional concept of “non-precipitating clouds” could be misleading, since it was the change in the precipitation flux from droplet sedimentation within clouds that modulates the drying of the boundary layer by cloud-top entrainment. Ackerman et al. (2009) even found that liquid water path responded more strongly to cloud water sedimentation than to drizzle in the models. Bretherton et al. (2007) examined the effect of cloud droplet sedimentation on entrainment rate and liquid water path of a nocturnal non-drizzling stratocumulus layer using an LES with bulk microphysics. In agreement with Ackerman et al. (2004), they found that consideration of droplet sedimentation decreased entrainment rate and increased liquid water path. However, instead of attributing the sedimentation–entrainment link to boundary-layer turbulence as Ackerman et al., they suggested that droplet sedimentation reduced entrainment by removing liquid water from the entrainment zone and thus inhibiting two mechanisms that promote the sinking of entrained air into the cloud layer: entrainment-induced evaporative cooling and longwave radiative cooling. A sensitivity study shows that the radiative effect is less important than the reduced evaporation. Caldwell and Bretherton (2009) further found that the impact of droplet concentration on entrainment is primarily through droplet sedimentation feedback rather than through drizzle processes. These results suggest that the impacts of droplet sedimentation on cloud-top entrainment should be considered in climate model simulations of aerosol indirect effects.

Note that these studies on droplet sedimentation have assumed the Stokes terminal velocity and lognormal cloud droplet size distribution with a fixed geometric standard deviation (or relative dispersion). The effect of varying spectral shape of the cloud droplet size distribution through droplet sedimentation has not yet been explored, although some qualitative effect can be gleaned from the fact that σ_g has been used as a tunable parameter to prevent the model from over-entraining (too small droplet sedimentation). Together with the studies on raindrop sedimentation, these studies have demonstrated that adequate consideration of HSD spectral shape effect is important for treating sedimentation of not only large raindrops but also small cloud droplets. Further extension to solid hydrometers such as ice crystals merits investigation as well.

Autoconversion is another important process in which spectral shape has been demonstrated to play a critical role. Briefly, all the autoconversion parameterizations that have been developed so far can be generically written as a two-part function given by:

where R is the autoconversion rate; R₀ represents the conversion rate after the onset of the autoconversion process and is called rate function; and T represents the threshold function describing the threshold behavior of the autoconversion process. The rate function R₀ has been the primary focus of most studies of autoconversion over the last few decades (e.g., Kessler, 1969; Manton and Cotton, 1977; Liou and Ou, 1989; Baker et al., 1993; Liu and Daum, 2004).

The threshold function T, however, has received less attention and is thus focused on here. Liu et al. (2006a) pointed out that according to how the threshold function is specified, existing autoconversion parameterizations can be grouped into three major types: Kessler-type, Sundquist-type, and Berry-type, named after the key original researchers. Briefly, Kessler (1969) assumed that the autoconversion process exhibited an abrupt threshold behavior described by:

where the Heaviside step function H(LWC − LWC_c) indicates that no autoconversion occurs when the liquid water content (LWC) is less than the threshold value LWC_c. Later Kesser-type parameterizations (Manton and Cotton, 1977; Liou and Ou, 1989; Baker, 1993; Liu and Daum, 2004) replace L with some measure of droplet size such that:

where r_m and r_c denote the driving and critical radius, respectively. Although it has been agreed that the threshold process in the Kessler-type parameterizations is driven by some mean diameter, the exact r_m differs with different Kessler-type parameterizations. For example, r_m respectively represents the mean volume radius (r₃) in the parameterization of Manton and Cotton (1977), mean radius of the fourth moment (r₄) in the parameterizations of Liou and Ou (1989), Baker (1993), and Boucher et al. (1995), and mean radius of the sixth moment (r₆) in Liu and Daum (2004). The essential discrepancy among different driving radii lies in the role of representing the effect of HSD spectral shape; for monodispersed HSDs, all the different mean diameters are equal. Furthermore, the critical radius r_c had been treated as an empirical parameter tuned in modeling studies until Liu et al. (2004), in which a theoretical expression was derived to relate r_c to LWC and cloud droplet concentration. To overcome the deficiency of all-or-nothing autoconversion with the Kessler-type parameterization, Sundqvist (1978) proposed another ad hoc threshold function given by

Del Genio et al. (1996) introduced a slightly different threshold function

Equation (5d) exhibits a cloud-to-rain transition sharper than Eq. (5c), but still smoother than the Heaviside function. Liu et al. (2006b) presented a more general Sundquist-type equation that covers more threshold behaviors and explicitly considers droplet concentration as well:

It can be readily shown that Eq. (5e) covers the three types of autoconversion threshold functions by varying the empirical parameter ν ≥ 0: it approaches the Kessler-type and Berry-type when ν approaches $\infty $ and 0, respectively. Berry-type refers to the BMPs that do not consider the threshold behavior or implicitly assume T = 1 (Berry, 1967; Beheng, 1994; Khairoutdinov and Kogan, 2000). Nevertheless, these discrepant expressions lack clear physics, thus precluding a sound choice between them. Furthermore, they are only for the conversion rate of liquid water content (third moment). These deficiencies were eliminated in Liu et al. (2006a, 2007), in which a general autoconversion threshold function is theoretically derived such that:

where Γ and $\gamma $ represent the complete and incomplete Gamma functions, respectively, and x_c is the ratio of the critical to mean masses of the droplet population. The parameter q is the spectral shape parameter in the Weibull droplet size distribution and has a unique relationship with the relative dispersion ε_c of cloud droplet size distribution given by:

The parameter p denotes the order of HSD moment, and Eq. (6) describes the threshold behavior of droplet concentration, liquid water content, and radar reflectivity, respectively, when p = 0, 3, and 6. Note that the use of ε_c, instead of the HSD shape parameter such as q in the Weibull HSD, has the advantage of encompassing the representation of HSD spectral shape embodied in the commonly used different analytical functions (e.g., Gamma, Weibull, and lognormal distributions; see Appendix for details) and different processes. Following these ideas, Xie and Liu (2011) derived the formulation by use of the four-parameter modified Gamma function that covers the three-parameter Gamma and Weibull HSDs as special cases.

Inspection of Eq. (6) reveals that the autoconversion threshold behavior is in fact determined by two dimensionless numbers: critical mass ratio x_c and ε_c (Fig. 3a). The critical mass ratio x_c determines where the “abrupt” conversion occurs, and ε_c determines the conversion steepness or threshold types. Most importantly, Eq. (6) reveals that the different types of autoconversion threshold functions essentially reflect the variation of ε_c; as ε_c increases from 0 to $ \infty $, the autoconversion threshold function changes from the Kessler-type through the Sundqvist-type to the Berry-type. The theoretical autoconversion parameterization eliminates the need for tunable empirical parameters (e.g., r_c and ν) at microphysical scales. Different versions of the theoretical threshold formulation have been applied to climate modeling studies (Rotstayn and Liu, 2009; Xie et al., 2018), CRM studies (Guo et al., 2008), threshold radar reflectivity separating precipitating clouds from non-precipitating clouds (Liu et al., 2008a), the dependence of rain initiation height on aerosol properties and updraft velocity (Chen et al., 2018a), and fog investigation (Niu et al., 2010; Lu et al., 2013a). It is noteworthy that the current expression for critical radius (thus x_c) does not account explicitly for turbulence effect and ε_c; further research is needed to address these issues. Also worth mentioning is that the autoconversion process is conventionally thought to be something as a practical need/convenience to represent the droplet collection process in models; however, the theoretical formulation reveals that it is a real physical threshold process that links droplet condensation, evaporation, and collection processes.

Figure 3.  Illustration of effects of relative dispersion (ε). The left panel, adapted from Liu et al. (2008b), shows that the dispersion-induced error in cloud albedo (left y-axis) and cloud radiative forcing (right y-axis). The right panel, adapted from Liu et al. (2006a), shows that different ad hoc types of autoconversion parameterization reflect the variation of relative dispersion.

The impact of HSD spectral shape on cloud radiative properties has been investigated primarily through its effect on cloud droplet effective radius (r_e) defined as the ratio of the third to the second moment of cloud droplet size distribution (Hansen and Travis, 1974; Slingo, 1989; Liu and Daum, 2000). Like the autoconversion parameterization, early parameterizations expressed effective radius as either a linear or a cubic root function of LWC without considering the dependence on droplet concentration and ε_c (Stephen, 1978; Fouquart et al., 1989). A serious deficiency of such “one-moment” r_e parameterization is its inability to study the first aerosol indirect effect. Later, Bower and Choularton (1992) and Bower et al. (1994) proposed an expression that relates r_e to the ratio of LWC to droplet number concentration for monodisperse cloud droplet size distribution or ε_c = 0. However, monodisperse droplet size distributions seldom occur in real clouds; Martin et al. (1994) empirically estimated ε_c for continental and marine clouds; and a similar idea has been used to distinguish the difference in ε_c between continental and marine clouds in climate models (Ghan et al., 1997; Lohmann et al., 1999; Rotstayn, 1999). However, such two-moment r_e parameterization does not fully consider the influence of spectral shape. Subsequent studies (Liu and Hallett, 1997; Liu and Daum, 2000; Lu et al., 2007) have shown that r_e can be generally parameterized through the “1/3” power-law:

where N denotes the droplet number concentration, and the dimensionless parameter β_e is an increasing function of ε_c that has a unique correspondence with the commonly used Gamma, Weibull, and lognormal size distributions (see Appendix for a list of the expressions). A larger ε_c leads to a larger β_e, which in turn leads to a larger r_e, smaller cloud optical depth, smaller cloud albedo, and smaller cloud radiative forcing, other factors being the same. This bias induced by incorrect representation of ε_c is called dispersion bias by Liu et al. (2008b), which showed that an incorrect ε_c could lead to a bias in cloud radiative forcing comparable in magnitude to the climate forcing caused by doubled CO₂ (Fig. 3b).

The significance of considering spectral shape is further reinforced by the fact that a perturbation in aerosol properties can concurrently alter cloud droplet concentration and ε_c, and the aerosol-induced dispersion effect acts to buffer the better-studied aerosol indirect effect through aerosol-induced change in droplet concentration (number effect hereafter) such that it negates the strong number effect in the aerosol-limited regime (Liu and Daum, 2002; Wood et al., 2002; Peng and Lohmann, 2003; Rotstayn and Liu, 2003, 2005, 2009; Yum and Hudson, 2005; Lu et al., 2007; Peng et al., 2007; Chen et al., 2012; Ching et al., 2012; Pandithurai et al., 2012; Xie et al., 2018; Wang et al., 2020), but enhances the weak number effect in the updraft-limited regime (Martins and Dias, 2009; Ma et al., 2010; Berg et al., 2011; Hudson et al., 2012; Chen et al., 2016a, 2018a; Xie et al., 2017). Igel et al. (2017) showed that ε_c significantly affects cloud fraction as well.

Although the importance of considering HSD spectral shape in BMPs has been increasingly recognized, representing it adequately in BMPs remain extremely challenging. A few studies have examined the dispersion effect by use of empirical expressions derived from measurements in climate and weather models, including the Australian CSIRO Mk3.0 (Rotstayn and Liu, 2003, 2005, 2009), the ECHAM (Peng and Lohmann, 2003), the Community Atmosphere Model (CAM) (Xie et al., 2017; Wang et al., 2020), the Institute of Atmospheric Physics’s atmospheric GCM (IAP AGCM) (Xie et al., 2018), and the Weather Research and Forecasting (WRF) model (Xie and Liu, 2011, 2015; Wang et al., 2013, 2018; Xie et al., 2013). Among many findings from these modeling studies, two points are worth emphasizing. First, the discrepancy among the different empirical expressions is large enough to cause significant discrepancies in numerical simulations in clouds, precipitation, and radiative forcing. Second, these empirical expressions are based on limited sets of observational data collected from certain types of clouds. Their generalization is thus questionable.

To overcome such deficiencies, there have been two rare attempts at formulating theoretical expressions that can be used to diagnose relative dispersion in two-moment BMPs (Liu et al., 2006c; Liu and Li, 2015). Briefly, both formulations start with the approximate equation for regular condensational growth given by (Rogers and Yau, 1989; Pruppacher and Klett, 1997):

where r is the droplet radius, S is the fractional water vapor supersaturation, and G is a function of the air temperature and pressure. The Liu et al. (2006c) formulation is based on the relationship

along with a steady-state equation for diagnosing the mean radius r₁. The Liu-Li formulation first diagnoses the mean radius from the condensation-induced change of LWC and supersaturation S based on

Relative dispersion ${\varepsilon }_{c}$ is then diagnosed from the relationship between r₁ and mean volume radius r₃ for the Gamma droplet size distribution

Wang et al. (2020) compared the impacts of the two diagnostic formulations and other empirical representations on climate model simulations and found significant differences.

It is worth mentioning that several studies have tried to obtain empirical relationships between the parameters of rain drop size distributions by analyzing measurements (Ulbrich, 1983; Zhang et al., 2001; Brandes et al., 2003; Cao and Zhang 2009) or model simulations with bin microphysics (Milbrandt and Yau, 2005a; Seifert, 2008).

Although the above-mentioned diagnostic expressions are valuable to improve two-moment BMPs with fixed spectral shape, they are not universal and may not be generalized beyond the conditions or assumptions on which these diagnostic expressions are based. To some extent, this deficiency manifests itself in the diversity of the empirical expressions for $ {\varepsilon }_{c} $. To allow for a more general consideration of HSD spectral shape, increasing efforts have been devoted to developing prognostic three-moment BMPs, in which an additional moment (usually M₆, which is proportional to radar reflectivity for spherical particles of constant density) is added to predict the spectral shape parameter (Milbrandt and Yau, 2005b; Szyrmer et al., 2005; Shipway and Hill, 2012; Dawson et al., 2014; Loftus et al., 2014; Naumann and Seifert, 2016). MY2005b is a three-moment scheme with the six hydrometeor categories being cloud, rain, ice, snow, graupel, and hail. All hydrometeors but ice crystals are assumed to have the power-law mass–radius relationship with a power exponent of 3; ice crystals are assumed to be bullet rosettes following Ferrier (1994). The primary physical process with prognostic shape parameter is the effect of hydrometeor sedimentation to mitigate the problem of excessive size sorting of the two-moment scheme with prescribed constant value of μ. Treatment of other microphysical processes largely follows previous studies (e.g., Ferrier, 1994) to balance M₆ budgets. Naumann and Seifert (2016) developed a three-moment BMP for warm rain processes, including evaporation and collection processes. Paukert et al. (2019) designed a three-moment BMP for rain microphysics by adding prognostic M₆ to replace the empirical shape–slope relationship based on Cao et al. (2008) in the so-called P3 scheme (see section 3.4 for more discussion on the P3 scheme). A unique feature of this work, compared to existing three-moment schemes, is the explicit consideration of influences of raindrop self-collection and collisional breakup on HSDs. Chen and Tsai (2016) developed a three-moment BMP for depositional growth of cloud ice crystals. The crystal size distribution was described with a three-parameter Gamma function; the geometric shape of ice crystals was represented using the volume-weighted aspect ratio. It was found that failure to consider the ice crystal shape led to a 45% underestimation of mass growth. Using only two moments to describe the crystal size distribution led to a 37% underestimation in mass and 28% underestimation in the bulk aspect ratio of the ice crystals. The proposed scheme was able to capture the shape memory effect and the gradual adaptation of ice crystal aspect ratios to a new growth habit regime.

Note that the treatment of cloud droplets in most existing so-called three-moment BMPs is actually two-moment, for example, assuming a Gamma size distribution with a fixed value of μ = 3 (Cohard and Pinty, 2000). Clark (1974) might be the first attempt at developing a three-moment BMP for cloud droplet population under the assumption of Gamma or lognormal droplet size distribution; but he only tested the BMP in highly simplified one- and two-dimensional cloud models. Deng et al. (2018) recently took up a similar task and tested their three-moment BMP with the WRF model. These rare three-moment BMPs for cloud droplets likely face challenges in representing both droplet concentration and ε_c in the updraft-limited regime similar to the empirical or analytical expressions (Ghan et al., 2011; Chen et al., 2016a).

It is noteworthy that similar multi-moment modeling efforts with assumed mathematical size distribution functions can be found in studies of aerosol dynamics (Williams, 1986). It should be valuable to compare the works in the two fields, especially in seeking to improve representation of aerosol–cloud interactions (Wang et al., 2013).

The above-mentioned BMPs and investigations into spectral shape effects are based primarily on the adiabatic assumption without (adequately) considering the effects of turbulence–microphysics interactions, including turbulent entrainment-mixing processes and stochastic condensation to be discussed in this subsection.

In seeking explanations for the discrepancies between observed and predicted droplet size distributions and rain initiation, a few studies as far back as the 1970s proposed the concepts of various turbulent entrainment-mixing mechanisms and investigated their impacts on cloud droplet size distributions (Warner, 1973; Latham and Reed, 1977; Baker and Latham, 1979; Baker et al., 1980). Figure 4 schematically illustrates the general characteristics of the various processes in terms of the commonly used microphysical mixing diagram (Burnet and Brenguier, 2007). Briefly, in homogeneous entrainment-mixing, entrained dry air is assumed to quickly mix with cloudy air such that all droplets experience the same condition for evaporation, leading to decreased droplet sizes and unchanged droplet number in the cloudy volume but decreased droplet number concentration because entrained air increases cloud volume through dilution (positive correlation between the mean droplet size and droplet number concentration). In the extreme inhomogeneous entrainment-mixing process, turbulent mixing is so slow that some droplets can completely evaporate to saturate the entrained dry air before mixing occurs, leading to decreased droplet concentration but unchanged mean droplet sizes. Fewer droplets means less competition for water vapor; therefore, the remaining droplets can grow bigger if the cloud parcel is lifted upward again after entrainment-mixing. This process is described as inhomogeneous mixing with subsequent ascent and leads to negative correlation between the mean droplet size and droplet concentration. In real clouds, any scenario between these extremes can occur, depending on the coupled dynamical and microphysical conditions characterized by the dimensionless Damkoler number defined as:

Figure 4.  Microphysical mixing diagram of mean volume radius vs. droplet concentration to illustrate major turbulent entrainment-mixing mechanisms and the corresponding microphysical relationships. See text for detailed explanation.

where τ_mix, τ_r, ε, and L are the turbulent mixing time, microphysical time, turbulent dissipation rate, and entrained dry eddy size, respectively.

Several lines of progress have been recently made in unifying the quantification of all different entrainment-mixing processes and developing parameterizations to represent their microphysical effects. The first line of progress lies in finding some microphysical measures to unify the various entrainment-mixing processes. In a modeling study, Morrison and Grabowski (2008) proposed an equation to estimate the effects of turbulent entrainment-mixing processes:

where the subscript “0” denotes the variable after entrainment but before evaporation, and 0 ≤ α ≤ 1 is an empirical parameter determining the type of turbulent entrainment-mixing processes (α = 0 for homogeneous entrainment-mixing and α = 1 for extreme inhomogeneous entrainment-mixing). In a series of studies (Lu et al., 2013a, b; Lu et al., 2014a), we further defined various microphysical measures called homogeneous mixing degrees such that a larger homogeneous mixing degree indicates a stronger degree of homogeneous mixing process.

Recognizing that “Da” requires knowledge of entrained eddy sizes that is not unique in turbulent clouds and difficult to know, Lehmann et al. (2009) suggested using the transition length (L*) in calculation of turbulent mixing time, which is the turbulent eddy size that corresponds to Da = 1 and is given by,

They showed that smaller L* leads to stronger homogeneous mixing. Note that a similar transition length was actually introduced earlier (Kabanov et al., 1971; Baker et al., 1980); but Lehmann et al. (2009) were the first to conduct a systematic examination of it. Lu et al. (2011) further introduced the dimensionless transition scale number (N_L) given by:

where η is the Kolmogorov turbulence scale and ν is the kinematic viscosity. The simplification of the Damkoler number to transition scale length and then to transition scale number is the second line of progress toward developing a parameterization to represent all the turbulent entrainment-mixing mechanisms in atmospheric models.

The final line of progress lies in building empirical relationships between the homogenous mixing degree and transition scale number that can be used to parameterize microphysical influences of entrainment-mixing processes on cloud microphysical properties (Lu et al., 2013b, 2014a; Gao et al., 2018; Luo et al., 2020; Desai et al., 2021). A parameterization was recently implemented in the LES version of WRF-Solar, a WRF version tailored for solar irradiance forecast (Endo et al., 2015; Jimenez et al., 2016), and examined for its influences on simulated cloud and radiative properties (Xu et al., 2022). There have been a few modeling studies with different and more complicated treatments of entrainment-mixing processes (Jarecka et al., 2009, 2013; Hoffmann and Feingold, 2019).

Beside the framework of homogeneous/inhomogeneous mechanisms, other ideas based on vertical mixing have been proposed to explain the microphysical variations during entrainment-mixing processes, including entity-type entrainment-mixing (Telford, 1996) and vertical circulation mixing (Wang et al., 2009). However, it is not clear as to how to represent these vertical mixing ideas in atmospheric models and how to reconcile the vertical mixing ideas with the homogeneous/inhomogeneous mixing framework, despite some efforts (Yum et al., 2015; Yeom et al., 2017).

Another related topic is stochastic condensation, which considers the growth of a droplet population as a stochastic process and relates the cloud droplet size distribution to multiscale supersaturation fluctuations associated with turbulence. Stochastic condensation was pursued in the 1960s, especially by Chinese and former Soviet Union scientists (Gu, 1962; Zhou, 1963; Levin and Sedunov, 1966; Sedunov, 1974). These early studies often replaced the full growth equations with simplified versions of kinetic equations amenable to analytical analysis, assumed Gaussian fluctuations, and claimed that turbulence fluctuations lead to spectral broadening of droplet size distributions. On the contrary, by numerically solving the growth equations under Gaussian fluctuating environments generated by Monte-Carlo simulations, Warner (1969), and Bartlett and Jonas (1972) predicted that turbulent fluctuations only slightly broadened droplet size distributions. They argued that supersaturation and updraft were so closely coupled to one another that a droplet experiencing a higher supersaturation and growing faster was likely in a stronger updraft and had a shorter time to grow. Manton (1979) demonstrated that turbulent mixing could break the link between supersaturation and updraft, leading to spectral broadening. Khvorostyanov and Curry (1999a, b) extended these early studies by presenting a more general set of stochastic kinetic equations and pointed out that the early low-frequency theories of stochastic condensation generally yielded droplet size distributions of the Gaussian type while observations tended to follow positively skewed distributions. They derived a more general mean-field equation and showed that the Gamma distribution was the solution to their equation under certain assumptions in the low-frequency regime. Khvorostyanov and Curry (2008a, b, c) further extended their formulation to include ice particles. By formulating stochastic condensation in terms of the stochastic Langevin equation and Fokker-Planck equation with known supersaturation fluctuation, McGraw and Liu (2006) derived a special type of Weibull droplet size distribution.

Several points are noteworthy. First, although addressing the issue of spectral shape is a primary motivation to invoke different turbulent entrainment-mixing processes, previous studies have been focused most on LWC and droplet concentration and thus only modify two-moment BMPs. Incorporation of ε_c into the framework remains elusive (Korolev et al., 2016; Pinsky et al., 2016a, b; Luo et al., 2021). Second, it is not trivial to represent supersaturation fluctuation/variation in atmospheric models that commonly assume saturation adjustment. Finally, turbulent entrainment-mixing and stochastic condensation have been largely investigated in separation, although entrainment-mixing and fluctuations are actually closely connected and act together on droplets over a range of scales (Su et al., 1998). Research is needed to consider entrainment-mixing processes and stochastic condensation together. Fully addressing these issues calls for particle-resolved DNS models with sufficiently large model domain size (see section 5 for more discussion on particle-resolved DNS).

In addition to choosing the optimal number of HSD moments, the other challenge in designing BMPs confronts the other BMP pillar—diversity of above-mentioned hydrometeor categories. The traditional way to handle the hydrometeor category diversity is through introducing more hydrometeor categories into BMPs, as illustrated in Fig. 2. Although it seems natural, the ever increase of hydrometeor categories is clearly limited by computational resources and the additional load of parameterizing more physical processes involved in BMP developments. Moreover, this traditional paradigm likely suffers from some fundamental difficulties. Take simpler BMPs in GCMs as examples (e.g., Rasch and Kristjánsson, 1998; Morrison and Gettelman, 2008; Gettelman and Morrison, 2015). GCM BMPs often assume two ice categories: small “cloud ice” and large precipitating “snow”, with the effects of rimed particles (graupel and hail) being neglected. Conversion of mass between the cloud ice and snow categories is parameterized by analogy to warm bulk microphysics schemes. Ice-to-snow autoconversion is treated in an ad hoc way that varies among schemes, while accretion is formulated by assuming continuous collection with a gravitational collection kernel, neglecting the fall speed of cloud ice. One-moment BMPS (e.g., Rasch and Kristjánsson, 1998) represent ice autoconversion by converting mass in the cloud ice category to the snow category when it exceeds a threshold mass mixing ratio. Two-moment BMPs (e.g., Morrison and Gettelman, 2008; Gettelman and Morrison, 2015) represent ice autoconversion by converting both mass and number of cloud ice larger than a critical size. Like liquid clouds, model simulations are highly sensitive to the value of critical ice size (Zhang et al., 2013; Eidhammer et al., 2014). However, unlike warm clouds, the ice-to-snow autoconversion may not represent any real individual physical processes because small ice particles can grow to snow through several pathways (e.g., deposition, aggregation, and riming). Furthermore, the physical properties of cloud ice and snow (density, shape, and terminal fall speed) vary widely. Thus, the separation of solid hydrometeors into cloud ice and snow as distinct categories and introduction of autoconversion may not have as strong of a physical basis as the liquid counterparts, and the abrupt (discrete) transitions between ice properties may be physically unrealistic. The situation could change for three-moment BMPs if a theoretical ice autoconversion threshold function can be derived similar to the autoconversion process of warm rain discussed in section 3.2. However, no work has been reported in the literature, and it may not be possible in view of the fundamental problems just discussed. Another deficiency in the conventional treatment of hydrometeor category is that empirical hydrometeor mass–size, projected area–size, and fall speed–size relationships are often used, which may not be self-consistent. This inconsistency likely causes unphysical results. Obviously, self-consistency of these relationships and quantities is desired because they are physically coupled in nature (Mitchell et al., 2011).

An entirely different paradigm has been emerging since the early 2000s that continuously tracks hydrometeor properties in time and space, instead of separating solid hydrometeors into different predefined categories. This emerging paradigm is hereafter called as the “continuous paradigm” in contrast to the traditional “discrete paradigm”, and key studies are briefly discussed here. The idea of the continuous paradigm originated with the bin ice microphysics of Hashino and Tripoli (2007). Harrington et al. (2013a, b) and Sulia et al. (2014) developed a bulk scheme that predicts particle habit evolution by including the crystal a- and c-axis mixing ratios as prognostic variables, thereby allowing for prediction of crystal axis ratio from vapor depositional growth. Other schemes use a diagnostic approach to include variability in ice particle properties. Lin and Colle (2011) included separate categories for cloud and precipitating ice and diagnosed the degree of riming and ice particle properties (mass–size and fall speed–size relationships) for precipitating ice as a function of the ratio of the riming rate to the sum of riming and deposition growth rates. Such a diagnostic approach is computationally efficient because it does not require additional prognostic variables, but particle properties are calculated locally and are not tracked in time and space. Morrison and Milbrandt (2015) generalized the approach and presented a method to predict several bulk physical properties of ice particles, eliminating the need for artificial conversion between ice categories and forming the basis for a conceptually new bulk microphysics scheme (P3 scheme for Predicted Particle Properties; see also Milbrandt and Morrison (2016). To represent the evolution of various physical properties in space and time, the P3 scheme includes a single solid-phase category with four prognostic mixing ratio variables: the total ice mass q_i, ice number N_i, the ice mass from rime growth q_rim, and the bulk rime volume B_rim. They also made a distinction between prognostic variables that are tracked through all the important mechanisms of change and growth, including vapor dynamical tendencies from advection and subgrid-scale mixing and microphysical tendencies (growth/decay processes and sedimentation), and predicted variables that are derived from the prognostic variables, including the rime mass fraction, bulk density, and mean particle size. Eidhammer et al. (2017) implemented the P3 scheme in CAM5 by modifying the two-moment bulk scheme (Morrison and Gettelman, 2008; Gettelman and Morrison, 2015). The mass–size and projected area–size relationships vary with particle sizes following Erfani and Mitchell (2016) and Morrison and Milbrandt (2015). Different lookup tables were used for the integration over HSD due to the size-dependent mass– and area–size relationships. Riming effect on physical properties of ice particle was neglected in Eidhammer et al. (2017). In a similar effort, Zhao et al. (2017) merged cloud ice and snow into a single prognostic variable, i.e., total ice, following Lin and Colle (2011). They incorporated the shape and riming impacts on ice particle properties through an environment-dependent riming intensity. In this BMP, eight microphysical processes (autoconversion of cloud ice to snow, accretion of cloud ice by snow, accretion of cloud water by snow, Bergeron process between cloud water and snow, melting of snow, fast freezing of rain, deposition, and sublimation of snow) are no longer needed or consolidated.

For liquid-phase clouds, it is generally thought that category classification is much less of a problem since the partition of cloud droplets and rain drops has clear correspondence with distinct growth modes of vapor diffusion for cloud droplets and collision–coalescence for rain, and cloud droplets and small raindrops are approximated by liquid spheres reasonably well. Nevertheless, some studies have shown that additional drop modes may be needed to accurately represent warm-rain processes (Saleeby et al., 2015). Kogan and Belochitski (2012) applied a similar idea of one hydrometeor category to develop a BMP for liquid-phase processes by predicting five moments of the full drop size distribution that consists of both small cloud droplets and large raindrops. The five prognostic moments are the zero, second, third, fourth, and sixth, corresponding respectively to number concentration, surface area, water content, drizzle flux, and radar reflectivity. The process rates and key microphysical properties are then diagnosed from the empirical multiple power-law fits to the LES simulations with a detailed bin microphysics of some stratocumulus clouds. Thus, this effort can be characterized as M₅C₁ according to the parameterization classification suggested in section 2. It is interesting to note that to some extent, this scheme trades hydrometeor categories for HSD moments.

Although the continuous paradigm of hydrometeor category is conceptually attractive, and perhaps has practical benefits according to the above-mentioned studies, many questions remain to be answered. For example, as discussed above, the computational cost associated with any BMP is likely proportional to the product of the number of hydrometeor size distribution moments and the number of hydrometeor categories. Will the benefit of reducing the latter be at the cost of adding the former, and to what extent? More studies are recommended on this promising new paradigm, with the high hope of a transformative change like the celebrated Scientific Revolution in astronomy from the Ptolemaic geocentric system to the Copernican heliocentric system.

Regardless of their detailed differences, virtually all BMPs are built on two pillars: (1) HSD is described by some known function, and perhaps the most used function is the Gamma distribution; (2) the relationships between hydrometeor properties and hydrometer radius (e.g., mass–radius relationship) can be described by power-laws generically given by:

where a_z and b_z are the prefactor and power exponent for the corresponding variable z, respectively. However, the ansatz of the two pillars is largely empirical or heuristic. In searching for a theoretical foundation for the common form of empirical HSDs and why observed droplet size distributions are generally broader than those predicted by the classical theory, a theoretical formalism has been developed by integrating into cloud physics the ideas of statistical physics and information theory (Zhang and Zeng, 1994; Liu et al., 1995; Liu, 1995, 1997; Liu and Hallett, 1997; Yano et al., 2016; Wu and McFarquhar, 2018). This theoretical formalism considers atmospheric particles as a system instead of following individual particles and thus is referred to as systems theory. We argue that the systems theory can be used to provide a theoretical foundation for the HSD function assumed in BMPs, because a BMP essentially seeks to represent the collective behavior of stochastic subgrid microphysical processes over the model grid, which is many orders of magnitude larger than the scale at which cloud microphysics acts. Furthermore, the assumed power-laws are closely related to the self-affine fractal geometry of individual hydrometeors and scale invariance properties of cloud turbulence fields (Mandelbrot, 1983; Liu, 1997). Below is a summary of the key equations involved in the theoretical formulations presented in those studies. In the framework developed by Liu and his coauthors throughout the years (Liu formalism hereafter), the particle system is constrained by two equations:

where z is called the restriction variable and related to the physical conservation laws controlling the particle system; Z is the total amount of z per unit volume; N is the total particle concentration; n(z) is the particle concentration per unit volume per unit z interval; ρ(z) = n(z)/N can be considered as the probability that the particle of z occurs. As the Boltzmann energy distribution describes the most probable energy distribution of a molecular system and the Maxwell velocity distribution characterizes the most probable velocity distribution of a molecular system, there exists a characteristic HSD that occurs most probably among all the possible HSDs. The most probable HSD is obtained by maximizing the spectral entropy defined after the Shannon entropy for complex systems:

where k is a proportional constant that has no effect on the derivation of the most probable droplet size distribution. Maximization of the spectral entropy subject to the constraints given by Eqs. (18a, b) yields the most probable distribution with respect to z:

where α = Z/N represents the mean Z per particle (compared to the physical meaning of “K_BT” in the Boltzmann energy distribution representing the mean energy per molecule). A combination of Eqs. (17) and (20) leads to the general most probable HSD of Weibull form:

where the parameters N₀ = a_zb_z/α and λ = a_z/α.

It is worth emphasizing that the restriction variable z, and thus parameters a_z and b_z, is related to the conservation laws controlling the particle system such as turbulent entrainment-mixing processes and particle geometrical shapes. Wu and McFarquhar (2018) presented a similar formulation based on relative entropy (hereafter WM formalism), instead of the Shannon information entropy. They derived the four-parameter generalized gamma distribution as the most probably size distribution:

The spectral shape parameter μ is related to the prior distribution satisfying the power law:

Note that $ {\rho }_{0}\left(r\right) $ was called the invariant measure in Wu and McFarquhar (2018). The power-law form of the prior distribution in the WM formalism is derived from the scale invariance principle that two different cloud volumes have the same shape of prior distribution.

It can be readily shown that the generalized Gamma function covers the Weibull distribution and Gamma distribution as special cases when μ = b−1 and b = 1, respectively. The primary distinction between the Liu formalism and WM formalism lies in their use of the prior distribution. The Liu formalism implicitly assumes a uniform prior distribution whereas the WM formalism assumes that the prior distribution obeys a power-law relationship to hydrometer radius independent of the power law relationship associated with the conservation law. Thus, the Liu formalism is equivalent to the WM formalism under the assumption that the exponent in the power-law prior is related to the power-law relationship associated with the conservation law such that μ = b−1. It remains an open question whether this equation holds as a result of some deeper scale-invariance principles, in addition to the common laws of conservation used to build the kinetic equations.

The other BMP pillar is that various hydrometeor properties (e.g., mass, volume, surface area, density, and terminal velocity) assume power-law relationships with hydrometer radii. Most power-law relationships are readily understood for homogeneous spherical particles such as cloud droplets. After analyzing the empirical results scattered in literature since 1935 on atmospheric aerosol particles and hydrometeors, Liu (1995, 1997) showed that the variety of power-law relationships for atmospheric particles could be unified into self-affine fractals (Mandelbrot, 1983), which covers as special cases irregularly shaped self-similar fractals and regular Euclidean geometric shapes. In practice, unlike a self-similar fractal that requires only one fractal dimension to quantify its morphology, characterization of self-affine fractals requires at least two fractal dimensions and thus simultaneous quantification (e.g., measurements) of mass, surface area, and size of particles. More studies on the power-law relationships for fractal-like ice crystal aggregates (both measured or simulated) can be found in Schmitt and Heymsfield (2010), Ishimoto et al. (2012), Letu et al. (2016), Lawson et al. (2019), Schmitt et al. (2019), Li et al. (2022), and Przybylo et al. (2019, 2022a, b). However, the self-consistency of the different power-law relationships has not been rigorously investigated in context of self-affine fractals. A systematic investigation of the various power-laws and their relationships under different conditions also holds promise in shifting the paradigm from discrete BMPs to continuous BMPs (Lin and Colle, 2011; Lin et al., 2011; Zhao et al., 2021).

Heuristic analysis suggests that the power-law relationships between the restriction variable and hydrometer radius are likely related to the physical laws controlling the processes and hydrometeor morphology. All else being equal, more irregular (smaller fractal dimension) hydrometers (Liu, 1995; Schmitt and Heymsfield, 2010) and stronger turbulent entrainment-mixing (Liu et al., 2002) tend to have broader HSDs. For solid hydrometeors, hydrometer morphology or power-laws may be related to turbulence and entrainment-mixing processes as well; however, such a connection has not been investigated adequately. The link between the scale-invariance for individual fractal hydrometers and the scale-invariance for hydrometeor population merits further investigation.

It should be noted that the above discussion is for a system that is either constrained by a single physical law or the physical laws can be characterized by a single power-law between the restriction variable z and hydrometeor radius r, i.e., the single mode HSD case. The general formulation of the systems theory that covers multiple physical constraints and more complex relationships between the restriction variable z and particle radius can be found in Liu et al. (1995). Some preliminary attempts have also been made to extend the systems theory to address the least probable distribution and scale-dependence of cloud droplet size distributions in turbulent clouds (Liu and Hallett, 1998; Liu et al., 2002), and formulate a kinetic potential theory on rain initiation (McGraw and Liu, 2003, 2004). More research along these lines is warranted.

A literature survey indicates that six different approaches have been used by different researchers to develop BMPs. The first is based on either educated guess or empirical observations. The Kessler BMP is a typical example of this approach. Such empirical parameterizations suffer from a major deficiency of lacking clear physics. The second approach is theoretical formulation. The advantage of the theoretical approach lies in its clear underlying physics, and thus, it should be used whenever possible. For example, Manton and Cotton (1977) used simple dimensional analysis arguments to extend the Kessler BMP to include the influence of varying cloud droplet concentrations. Verlinde et al. (1990) showed that an analytical solution to the full stochastic collection equation could be obtained if the collection efficiencies are held constant. In a series of publications, Liu and his coauthors (Liu et al., 2004, 2005, 2006a, b, 2007) derived a sequence of theoretical autoconversion parameterizations. Their theoretical analysis revealed the implicit assumptions underlying some commonly used autoconversion schemes and eliminated some empirical parameters (e.g., critical radius). Despite the high desirability of clear physics and reduced number of tunable parameters, it remains challenging to derive analytical expressions describing all the microphysical processes and hydrometeor types involved. This is particularly true with more complex BMPs (e.g., increasing number of moments and/or hydrometeor categories). The third approach relies on curve-fitting of numerical simulations by detailed numerical models with explicit microphysics. Most BMPs are built on parcel or one-dimensional models with bin microphysics (Berry, 1967, 1968; Berry and Reinhardt, 1974a, b, c, and d; Beheng, 1994; Chen and Liu, 2004). Khairoutdinov and Kogan (2000) derived their BMP from the statistical analysis of simulated stratiform clouds by an LES with bin microphysics. Noh et al. (2018) explored using simulations with a particle-based cloud model to develop BMPs. The fourth approach is based on numerical simulations with bin microphysics as well. However, instead of establishing simple empirical expressions describing individual processes, this approach relies on look-up tables (LUT) generated directly from bin microphysics models to circumvent the difficulties of fitting complex non-linear relationships (Feingold et al., 1988, 1999; Saleeby et al., 2015). Feingold et al. (1988) replaced the autoconversion formulation in the Regional Atmospheric Modeling System (RAMS) with full stochastic collection solutions considering droplet self-collection and rain (drizzle) drop collection of cloud droplets with more realistic collection kernels from Long (1974) and Hall (1980) rather than assuming constant collection efficiencies used in the earlier versions of RAMS. Feingold et al. (1999) applied a similar LUT approach to emulate representation of collection and drop sedimentation. Walko et al. (1995) described the implementation of this approach in RAMS for prediction of hydrometeor mixing ratios. Meyers et al. (1997) extended this approach to mixing ratio and number concentration. Saleeby et al. (2015) added explicit nucleation of cloud droplets as well as a second cloud mode of large cloud droplets with diameters from 40 μm to 80 μm. The addition of a new droplet category changes the BMP in RAMS to a M₂C₈ BMP that considers small cloud droplets, large cloud droplets, rain, pristine ice, snow, aggregates, graupel, and hail. A similar LUT approach was used in Paukert et al. (2019) in their three-moment BMP for rain processes. The fifth approach—the machine learning approach—is emerging as a viable candidate as a result of the advances in computational sciences, machine learning, and big data. The machine learning approach can be viewed as an extension of the third approach based on statistical regression models, or using more advanced machine learning models to fit more complex nonlinear relationships, or as a mathematical formalization of the LUT approach. The data used to develop the machine learning models may come from numerical models (virtual data) and/or measurements. For example, Rothenberg et al. (2018) built a polynomial chaos expansion emulator of an adiabatic parcel model with bin microphysics that can be used to diagnose droplet number at the maximum supersaturation from a wide range of dynamic, thermodynamic, and aerosol conditions. They found that the emulator requires much less computational time to build, store, and evaluate than a high-dimensional lookup table. Liu et al. (2018) emulated their parcel model simulations with a multiple layer perceptron (MLP) model to estimate liquid water content, droplet concentration, and ε_c at levels of maximum supersaturation. Seifert and Rasp (2020) explored the use of neural networks for parameterizing the processes of autoconversion, accretion, and self-collection in a two-moment BMP based on particle-based simulations of the collision-coalescence equation. Chiu et al. (2021) based their neural network modeling on the numerical simulations of the stochastic collection equation with bin microphysics constrained by measured droplet size distributions. Among other things, both studies suggested that the autoconversion parameterization may be improved by incorporating information about rain, due perhaps to the fact that although rain has no direct impact on autoconversion by definition, it contains information on the spectral shape of cloud droplet size distribution known to affect autoconversion rate (e.g., Liu et al., 2006a). Gettelman et al. (2021) replaced the bin microphysical model with several neural networks designed to emulate the autoconversion and accretion rates produced by the bin microphysical model. Despite its advantages, the machine learning approach aggravates the shortcomings of unclear physics and interpretability fraught with the traditional empirical and LUT approaches. It is interesting to note that the activation functions commonly used in machine-learning models resemble the functions describing such microphysical threshold processes as activation of aerosol particles into cloud droplets and autoconversion of cloud droplets to raindrops. Such functional similarities may be valuable to develop interpretable machine learning BMPs. In practice, hybrid approaches that combine two or more above approaches have often been, and will likely continue to be, used as the sixth approach. For example, LUT is partly used for some complex processes in some popular BMP packages (Morrison and Gettelman, 2008; Thompson et al., 2008).

Despite many efforts and achievements, much remains to be done in the development of multi-moment schemes. The following points are particularly worth mentioning. First, existing studies have been concerned mainly with simple hydrometeor categories and a limited number of microphysical processes, and thus, more comprehensive and systematic investigation is needed. Second, as the number of predicted moments increases, the optimal choice of moments to use is worth exploring (Wacker and Lupkes, 2009; Milbrandt and McTaggart-Cowan, 2010; Shan, 2020). Third, so far, we have discussed BMPs without specific differentiation between the types and spatiotemporal resolutions of numerical models wherein the BMPs are used (e.g., GCM vs. CRM or LES). However, BMPs are expected to be scale-dependent from the fundamental nature of multiscale fluctuations associated with turbulent clouds (Liu et al., 2002). Although one often deals with such issues by tuning empirical parameters embedded in BMPs, the need for scale-aware BMPs will become more acute for models with an adaptive or unstructured mesh. Finally, as discussed above, detailed models with explicit microphysics, especially bin microphysics, have been instrumental in developing BMPs. However, many aspects of explicit microphysics modeling are in dire need of improvement. Furthermore, although it is well known that small-scale turbulence (e.g., sub-LES scales) and related processes play a critical role in determining cloud microphysics, the effects of turbulence-related processes have not been considered in most BMP_S, and serious knowledge gaps remain in understanding turbulence–microphysics interactions. Modeling explicit microphysics and the effects of turbulence are the foci of section 4 and section 5, respectively.

As alluded in the preceding sections, cloud models with explicit microphysics are instrumental in understanding various microphysical processes and developing BMPs by providing virtual data and benchmarking simulations. In explicit microphysics models, the dynamics of individual aerosol/cloud/precipitation particles and their stochastic collisions are “explicitly” expressed, and thus, there are fewer assumptions than in bulk models. Up to now, two different types of numerical approaches have existed for dealing with explicit microphysics: bin scheme and particle-based scheme. Bin schemes have about 70 years of tradition and have been the predominant representation of explicit microphysics since their inception. Excellent reviews can be found in the literature (e.g., Khain et al., 2015; Khain and Pinsky, 2018; Grabowski et al., 2019; Morrison et al., 2020). The less popular particle-based scheme has emerged as a promising methodology for treating explicit microphysics during the last decade (e.g., Shima et al., 2009). Figure 5 illustrates the conceptual difference between bin schemes and particle-based schemes. In the rest of this section, the basic characteristics of the two explicit microphysics approaches are elucidated, with slightly more emphasis on particle-based schemes.

Figure 5.  Conceptual difference between bin schemes and particle-based schemes [adopted from Fig. 4 of Grabowski et al. (2019), ©American Meteorological Society, used with permission]. (a) Scatter plot of real droplets on a plane of droplet size versus CCN size. (b) This droplet population can be expressed by a two-dimensional histogram. (c) Most of the bin schemes approximate the intrinsically multidimensional histogram by a one-dimensional histogram. (d) In particle-based schemes, the real droplet population is approximated by a much smaller number of computational particles, often called super-droplets or super-particles. Each super-droplet represents a multiple number of real droplets, which is denoted by the multiplicity variable $ \mathit{\xi } $. The color in the figure changes from dark blue to red as the multiplicity $ \mathit{\xi } $ increases.

This subsection formally provides the equations that describe the detailed dynamics and evolution of all the atmospheric particles. Briefly, the state of a particle is characterized by two types of vector variables: position $ \boldsymbol{x} $ and attribute $ \boldsymbol{a} $. The attributes consist of a number of variables, denoted by d, that represent the internal state of the particle (e.g., particle radius, mass of liquid water, mass of ice water, masses of other chemical substances, shape of ice particle, and amount of electric charge). Let us assign a unique index $ i $ to each particle such that the state of the cloud microphysical system is uniquely specified by the positions and attributes of all the particles $ \{\left({\boldsymbol{x}}_{i},{\boldsymbol{a}}_{i}\right),i=1,\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }2,\mathrm{ }\dots ,{N}_{\mathrm{r}}\} $, where $ {N}_{\mathrm{r}} $ is the total number of particles. Without loss of generality, various cloud microphysical processes can be categorized into two general types: individual particle dynamics without interactions with other particles and collision with particle–particle interactions. Individual dynamics can be formally expressed by a set of ordinary differential equations:

where $ {\boldsymbol{v}}_{i}\in {\boldsymbol{a}}_{i} $ is the velocity of particle $ i $ and $ {\boldsymbol{G}}_{i} $ is the state of ambient air. Equation (24a) represents the transport of particles by the flow and gravity; Equation (24b) describes various growth processes such as diffusional growth by condensation/deposition and chemical reactions. Collision of particles can be regarded as a memoryless Markovian stochastic process. Here, for simplicity, we only consider collision-coalescence, and ignore collision–breakup, though collision–breakup can be formulated similarly. Following Gillespie (1972), consider a region with volume $ \mathrm{\Delta }V $. If $ \mathrm{\Delta }V $ is sufficiently small, we can consider that the particles inside this region are “well-mixed”, e.g., by atmospheric turbulence. Then, all the pairs of particles in the volume have a possibility to collide and coalesce, and the probability that a particle pair $ j $ and $ k $ inside the small volume $ \mathrm{\Delta }V $ collide and coalesce within an infinitesimal time interval $ (t,t+{\rm{d}}t) $ is given by:

where the function $ K\left({\boldsymbol{a}}_{j},{\boldsymbol{a}}_{k};\boldsymbol{G}\right) $ is called the collision-coalescence kernel, and $ \boldsymbol{G} $ denotes the state of the moist air in $ \mathrm{\Delta }V $. The two sets of equations together describe the temporal evolution of the particle population $ \left\{\left({\boldsymbol{x}}_{i},{\boldsymbol{a}}_{i}\right),i=1,\mathrm{ }2,\mathrm{ }\dots ,{N}_{\mathrm{r}}\right\} $. However, because of the high number of particles in a typical cloud [on the order of ${10}^{8}\;{\mathrm{m}}^{-3}\times (10\;\mathrm{k}{\mathrm{m})}^{3}={10}^{20}]$ and the wide range of scales spanned by turbulent clouds and associated physical processes (e.g., from 1 mm to 10 km), it is practically impossible to conduct a full direct numerical simulation (DNS) of this fundamental model that resolves individual particles and the smallest turbulent eddies and covers a sufficiently large model domain. Approaches to simplifying the model are still necessary in practice, although it is referred to as explicit microphysics both in literature and in this paper. This section is focused on the bin scheme (section 4.2) and particle-based scheme (section 4.3) and their pros and cons (section 4.4). Studies on and with particle-resolved DNS models are deferred to section 5.

Instead of following all the particles, bin schemes approximate the mean distribution function of particles by a histogram (see Figs. 5b and c) and predict its temporal evolution. Let $ n(\boldsymbol{a},\boldsymbol{x},t) $ be the mean distribution function of particles, which represents the number density of particles with attribute $ \boldsymbol{a} $ at position $ \boldsymbol{x} $ at time $ t $. If we assume that the mean number densities of any two different attributes are not correlated, we can derive the following temporal evolution equation of $ n(\boldsymbol{a},\boldsymbol{x},t) $ (Gillespie, 1972):

where ${a}{{'}}{{'}}$ is such an attribute that it becomes $ {a} $ if it coalesces with ${a}\mathrm{{'}}$; the superscript d denotes the number of variables describing the internal state of the particle (e.g., droplet radius, masses of aerosol substances, mass and shape of ice crystals, and electric charge). This equation is often called the Smoluchowski equation (Smoluchowski, 1916), which is also referred to as the “stochastic collection equation” in the cloud microphysics community. The basic idea of bin schemes is to solve the Smoluchowski equation by approximating the mean distribution function $ n({{\boldsymbol{a}}},{{\boldsymbol{x}}},t) $ by a histogram. Due to its computational cost, most of the bin schemes approximate the multidimensional histogram by a one-dimensional histogram in terms of particle radius (see Figs. 5b and c). The development of bin schemes has been ongoing since the late 1940s (Howell, 1949; Mason and Ramanadham, 1954; Berry, 1967; Bleck, 1970; Clark, 1973; Berry and Reinhard, 1974a, b, c, d). They are now commonly used in three-dimensional cloud models for simulating various types of clouds or as a type of microphysical parameterizations; readers can find reviews on bin microphysics and its applications in literature (e.g., Khain et al., 2015; Khain and Pinsky, 2018; Grabowski et al., 2019; Morrison et al., 2020).

Bin microphysics (especially simplified versions) has been used for microphysical parameterizations as well as serving as a tool for improving physical understanding and generating synthetic data to build BMPs. However, the practical advantages of bin microphysics over BMPs are not conclusive based on some comparative studies. For example, Seifert and Beheng (2006) found that the two-moment bulk scheme can be adjusted to produce results consistent with a spectral bin scheme for an isolated convection cell and a squall line. Morrison and Grabowski (2007) suggested that warm cloud simulations with bulk schemes are comparable to the results from using bin microphysics, but the one-moment bulk scheme without predicting number concentrations of hydrometeors yields significant errors in the simulations of cloud microphysical properties. Other deficiencies in bulk schemes were reported by Li et al. (2009a, b), suggesting that an over-estimated rain evaporation rate in the bulk scheme produces a much stronger near-surface cold pool than the bin scheme simulation, which reproduces the distinct life cycles of convective clouds. Khain et al. (2009) and Ekman et al. (2011) pointed out that the bulk scheme with prescribed aerosol concentration fails to show the different precipitation sensitivities to aerosol concentrations in dry/humid environments. More recently, a modeling study by Fan et al. (2012) revealed striking differences in cloud microphysical properties and opposite aerosol effects on convection and heavy precipitation between the bin scheme and BMPs implemented in WRF, mainly because of pre-described aerosols in the bulk scheme. Gettelman et al. (2021) explored replacing the BMP with bin microphysics for warm rain formation process in a GCM and found improvements in the rain onset and frequency of light rain, compared to high resolution process models and observations and at the expense of a 400% computational slowdown. They then tested replacing the computationally demanding bin microphysical model with a machine learning emulator consisting of several neural networks.

In particle-based schemes, the population of all the aerosol/cloud/precipitation particles $ \left\{\left({\boldsymbol{x}}_{i},{\boldsymbol{a}}_{i}\right),i=1,\mathrm{ }2,\mathrm{ }\dots ,{N}_{r}\right\} $ is approximated by a much smaller number of computational particles, often called super-droplets or super-particles $ \{\left({\xi }_{i},{\mathbf{x}}_{i},{\boldsymbol{a}}_{i}\right),i=1,\mathrm{ }2,\mathrm{ }\dots ,{N}_{\mathrm{s}}\} $, where $ {N}_{s}\ll {N}_{r} $ is the number of super-particles (see Figs. 5a and d). A super-particle is characterized by multiplicity $ {\xi }_{i} $, position $ {\boldsymbol{x}}_{i} $, and attributes $ {\boldsymbol{a}}_{i} $. We consider that the $ i $-th super-particle represents $ {\xi }_{i} $ number of real particles $ \left({\boldsymbol{x}}_{i},{\boldsymbol{a}}_{i}\right) $. Particle-based schemes for cloud microphysics have been gaining popularity over the last decade (Jensen and Pfister, 2004; Paoli et al., 2004; Shirgaonkar and Lele, 2006; Andrejczuk et al., 2008, 2010; Shima et al., 2009, 2020; Sölch and Kärcher, 2010; Riechelmann et al., 2012; Grabowski and Abade, 2017; Abade et al., 2018; Brdar and Seifert, 2018; Grabowski et al., 2018; Jaruga and Pawlowska, 2018; Hoffmann et al., 2019; Seifert et al., 2019).

The individual dynamics of real droplets can be applied straightforwardly to super-particles without much complication. However, the treatment of the collision-coalescence was the major challenge in the early stage of development, because a naive implementation of collision-coalescence could lead to huge computational cost or errors. A direct calculation of binary collision costs $ O\left({N}_{s}^{2}\right) $. Another concern is the size and change of super-particle number $ {N}_{\mathrm{s}} $. To get enough statistics, too small numbers of super-particles have to be avoided, even though coalescence always decreases the number of real particles. On the other hand, a too large number of super-particles is not acceptable either in view of the computational cost. Shima et al. (2009) developed a stochastic algorithm referred to as the “super-droplet method” (SDM) that eliminates those problems. The SDM scales linearly with the number of super-particles $ O\left({N}_{\mathrm{s}}\right) $, keeping the number of super-particles $ {N}_{\mathrm{s}} $ unchanged during collision-coalescence. Unterstrasser et al. (2017, 2020), and Dziekan and Pawlowska (2017) showed that the SDM provides an efficient and accurate numerical algorithm for collision-coalescence. There are several other collision-coalescence algorithms being used in other areas, such as O’Rourke’s method (1981) and the no-time counter method modified for spray combustion (Schmidt and Rutland, 2000); the weighted flow algorithm for aerosol dynamics (DeVille et al., 2011); and Ormel and Spaans’ method (2008) and Johansen et al.'s method (2012) for astrophysics. Li et al. (2017) confirmed that the performance of the SDM is better than Johansen et al.'s method (2012), but performance comparison with other algorithms remains to be conducted.

Note also that there is some arbitrariness as to how to initialize super-particles, and the performance of the SDM is sensitive to it (Unterstrasser et al., 2017). Unterstrasser et al. (2017), and Dziekan and Pawlowska (2017) proposed a procedure using a grid to distribute the super-particles more uniformly along the particle size axis. However, it is known that the “discrepancy” of an axis-aligned grid decreases slowly in a high dimension, though it is optimal for one dimensional space (Niederreiter, 1978). Hence, their procedures would not be suitable for 3D cloud simulations with various particle attributes. More discussion on this point was provided in section 5.3 of Shima et al. (2020).

As well summarized in Grabowski et al. (2019), there are three known issues in bin schemes, but the SDM and other particle-based schemes can mitigate these problems. The first issue is the numerical diffusion. Because bin schemes discretize the physical and attributional space into a grid, numerical diffusion is unavoidable, causing unphysical broadening of particle size distributions (Clark, 1973; Morrison et al., 2018). On the other hand, most of the particle-based schemes do not need discretization, and hence, there is no numerical diffusion.

The second issue is the “curse of dimensionality”. Most existing bin schemes are one-dimensional and predict only the size distribution of particles. Although there are various motivations to consider multiple attributes, modeling with multiple attribute bins is computationally prohibitive. Briefly, if we need $ d $ number of attributes, then the bin space becomes $ d $ dimensional. Let us assume that we need 100 bins for each dimension. Then, the total number of bins needed yields $ {100}^{d} $. Further, to calculate the collision-coalescence process, all the combinations of the bins have to be investigated. Then, the computational cost scales to $ {100}^{2d} $, which increases rapidly as $ d $ increases; each additional dimension requires a 10⁴ times larger computational cost. This unfavorable feature is referred to as the curse of dimensionality of bin schemes. The computational cost of particle-based schemes also increases with increasing $ d $ because we need to use more super-particles to sample the higher-dimensional attribute space. Nevertheless, particle-based schemes can relax this issue. Shima et al. (2009) deduced from a theoretical argument that the SDM is more computationally efficient than bin schemes when the number of attributes is larger than a critical number of 2–4. The derivation is based on the general features of the SDM being that the computational cost is linear $ O\left({N}_{\mathrm{s}}\right) $ and that the number of super-particles $ {N}_{\mathrm{s}} $ is almost always conserved. They also assumed that the discretization error of bin schemes is in the range of 1st to 2nd order. For more details about the estimation, see Appendix B of Shima et al. (2009) and Shima (2008). Recently, Li et al. (2017) elucidated that the performance of the SDM is superior to bin schemes, at least under the specific condition they investigated.

The third issue is related to the limitation of the Smoluchowski equation itself, which is the equation for bin schemes to solve numerically. Derivation of the Smoluchowski equation assumes no correlation between the number densities of two different particle sizes. However, the correlation cannot be ignored in general. For example, when the well-mixed volume is small, the particles in the volume are so few that their discreteness is evident (see Dziekan and Pawlowska, 2017 and references therein). Furthermore, the size of a well-mixed volume varies with the atmospheric condition; when turbulence is weak or the time scale of collision-coalescence is short, the well-mixed volume becomes smaller. An additional limitation of the Smoluchowski equation is that it cannot account for fluctuation around the mean adequately, making it incapable of describing the rapid growth of lucky droplets that could be crucial for rapid precipitation onset (Telford, 1955; Kostinski and Shaw, 2005). In contrast, the SDM is independent of the no-correlation assumption. Therefore, as Dziekan and Pawlowska (2017) confirmed, the SDM provides the correct solution even when the well-mixed volume is small, and the Smoluchowski equation is not valid. Furthermore, the SDM can be used for simulating lucky droplets. However, particle-based schemes are expected to produce a higher variance than reality because the collision-coalescence process is resolved by a smaller number of samples than in nature. More work is needed to better understand the impact of statistical aspects on rain formation through collision-coalescence.

Despite its popularity, bin microphysics faces some numerical and conceptual problems that are difficult to overcome. Particle-based schemes can resolve or relax those issues, thereby paving the way toward a more accurate representation of explicit microphysics. Considering the growing interest in the methodology, great progress can be expected in the years to come. That said, much remains to be done for particle-based schemes and related applications. First, a quantitative understanding of elementary cloud microphysical processes is becoming more important. Particle-based schemes are already being used, not only for warm clouds, but also for in-cloud aqueous-phase chemistry and ice/mixed phase clouds. However, our process-level understanding of cloud microphysics, in particular for ice/mixed phase, is not sufficient, which is critical for particle-based models. Establishing a reliable particle-based algorithm for spontaneous/collisional breakup is also a high priority. Second, the computational cost of the current particle-based cloud models is at least one or two orders of magnitude larger than bulk models. Incorporating the impact of SGS turbulence fluctuation to particle-based microphysics is one way to reduce the computational cost. Similar to cloud models with bin microphysics, cloud models with particle-based microphysics typically have grid sizes of a few tens of meters or less with typical CRM and LES domain sizes. Therefore, introducing SGS turbulence models not only to atmospheric fluid dynamics but also to cloud microphysics should improve the grid convergence characteristics, enabling the use of larger grid spacing. For example, SGS models for collision-coalescence were proposed by Wang et al. (2008), Onishi and Seifertet al. (2016), and Chen et al. (2018b). SGS models for supersaturation fluctuation were proposed by Grabowski and Abade (2017), Abade et al. (2018), and Hoffmann et al. (2019). Another idea is to use the Twomey activation approach if cloud-processed CCNs are of secondary importance (Grabowski et al., 2018). A huge reduction of the computational cost can be expected, because the super-particles are used only inside the clouds, but adaptive load balancing of parallel computation is required. Adaptive adjustment and resampling of super-particles are also worth considering to optimize the distribution of super-particles in the physical and attributional space (Unterstrasser and Sölch, 2014; Schwenkel et al., 2018). Even with those improvements, it would not be practical to directly use particle-based cloud models in large-scale NWPs and GCMs; attempts to construct high-fidelity BMPs from particle-based cloud models need to be continued (Noh et al., 2018), as well as comparisons with bin microphysics (Chandrakar et al., 2022).

It is known that turbulence affects the evolution of a cloud on many scales; however, our knowledge on turbulence–microphysics interactions is extremely limited, and most microphysics parameterizations do not adequately consider such turbulence influences, including high resolution LES models with either bin microphysics (Khairoutdinov and Kogan, 2000) or particle-based microphysics (Noh et al., 2018, Hoffmann et al., 2019). Today, cloud–turbulence interaction represents a challenging area of research with limited knowledge, and this section focuses on the interaction of turbulence with cloud droplets at the finest spatial scales in a warm cloud, mostly from the viewpoint of very recent progress using particle-resolved direct numerical simulation (PR-DNS) models. Unlike LES, which relies on subgrid-scale (SGS) parameterizations to represent the processes at sub-LES scales, DNS resolves the smallest turbulent eddies at the Kolmogorov microscale (~1 mm in clouds) without any SGS parameterizations. PR-DNS models further track individual cloud particles affected by turbulence, instead of using the approximations of bin schemes or SDM schemes. Also noted is the limited accuracy of the so-called kinetic equations solved in most models with bin microphysics, which remain questionable themselves, especially in the context of stochastic condensation and multiscale fluctuations (Khvorostyanov and Curry, 2014). Thus, except for their limited model domain size due to computational power constraints, PR-DNS models constitute the microphysical models at the most fundamental level. The material presented here covers turbulent effects on condensation/evaporation, collision-coalescence, and the simultaneous condensation-collision-coalescence processes. It supplements some earlier reviews of the subject, including Vaillancourt and Yau (2000), Shaw (2003), Devenish et al. (2012), Grabowski and Wang (2013), and Gao et al. (2018).

For a rising air parcel containing cloud droplets, the equations governing the turbulent flow velocity, the velocity of the droplets, the disturbance flow exerted by neighboring droplets, the rate of condensational growth of droplets, and the changes in thermodynamic and water vapor variables are given in Vaillancourt et al. (2001), Franklin et al. (2005), and Chen et al. (2018b, c). Specifically, the turbulent velocity field U is governed by the incompressible Navier-Stokes equations:

where P is the perturbation pressure deviation from the hydrostatic parcel mean pressure $ {P}_{M} $, g is gravity, B is buoyancy in a form given in Vaillancourt et al. (2001), and F is the external forcing to maintain the turbulence following Chen et al. (2016b). The droplet velocity V(t) is governed by the drag force of the fluid and gravity of the form:

where $ {\tau }_{p} $ denotes the droplet response time and $\tilde{\boldsymbol{U}}$ is the flow velocity at the droplet center location, $ \boldsymbol{X}\left(t\right), $ represented by the sum of the turbulent flow velocity $ \boldsymbol{U} $ and the disturbance flow velocity $ {\boldsymbol{U}}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}} $ contributed by the neighboring droplets. The superposition method by Wang et al. (2005) is used to calculate the disturbance flow $ {\boldsymbol{U}}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}} $. The details are given in Chen et al. (2018b).

For the condensational growth of droplets, it is convenient to solve a set of microscopic equations and a set of macroscopic equations (Vaillancourt et al., 2001). The microscopic set includes the condensational growth equation of an individual droplet i with radius $ {r}_{i}, $ with the form:

where $ {{e}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is the saturation water vapor pressure; $ {f}_{v} $ refers to the droplet ventilation coefficient determined by the empirical formulas as functions of the Reynolds number and the Schmidt number obtained from the laboratory experiment of Beard and Pruppacher (1971); and S is the supersaturation in the grid cell where droplet i is located, defined as:

where $ {q}_{v} $ and $ {q}_{vs}(T,P) $ are the water vapor mixing ratio and the saturation water vapor mixing ratio, respectively. All droplets in the same grid cell are assumed to be exposed to the same supersaturation environment. The scaler fields of $ {{q}}_{\mathrm{v}} $ and temperature T can be decomposed into the parcel mean state and the perturbation state. The perturbations of temperature ($ {T}^{{{'}}} $) and water vapor mixing ratio ($ {q}_{v}^{{{'}}} $) are calculated as follows:

where ${W}^{\mathrm{{'}}}$ is the vertical perturbation velocity, ${{\varGamma }}_{d}$ is the dry adiabatic lapse rate, L is the latent heat of vaporization, ${C}_{p}=1004.0\;\mathrm{J}\;\mathrm{k}{\mathrm{g}}^{-1}\;{\mathrm{K}}^{-1}$ is the specific heat for dry air at constant pressure, ${D}_{t}=2.22\times {10}^{-5}\;{\mathrm{m}}^{2}\;{\mathrm{s}}^{-1}$ is the coefficient of thermal diffusivity, $ {D}_{v}=2.55\times {10}^{-5}{\mathrm{m}}^{2}{\mathrm{s}}^{-1} $ is the coefficient of water vapor diffusivity in the air, and ${C}_{d}^{{{'}}}={C}_{d}-{C}_{dM}$ is the differential condensation rate between the grid cell and the whole parcel. Given Eq. (30), the condensation rate inside the grid cell of mass $ {m}_{a} $ can be simplified as:

The macroscopic equations govern the time evolution of the parcel-mean temperature $ {T}_{M} $, water vapor mixing ratio $ {q}_{vM} $, pressure $ {P}_{M} $, and density $ {\rho }_{aM} $. All variables of parcel mean are denoted with a subscript M.

The gas constant for air is denoted by $ {R}_{a} $. The total fields of T and $ {q}_{v} $ are calculated by adding the macroscopic variables and the perturbation variables.

Analysis of cloud droplet size distributions in cumulus and stratocumulus clouds indicates that the measured spectra are broader than those obtained from adiabatic parcel calculations, even in quasi-adiabatic cloud cores (Brenguier and Chaumat, 2001; Pawlowska et al., 2006; Prabha et al., 2012). Srivastava (1989) first suggested that the random spatial distribution of droplets and the variations in the vertical air velocity could cause nonnegligible deviation of the supersaturation experienced by individual droplets. By assuming a Poisson spatial distribution of the droplets, he derived analytical estimates of the variance in the microscopic supersaturation and concluded that the variance could contribute significantly to the width of the droplet size distribution. However, the effects of small-scale turbulence and droplet sedimentation were neglected in Srivastava’s estimates. This shortcoming was alleviated in the work of Vaillancourt et al. (2002), which reported the first DNS study following the growth of a large number of cloud droplets by the diffusion of water vapor in homogeneous isotropic turbulence. They numerically solved the above set of equations within a cubic domain of 10 cm³ loaded with 50 000 droplets. When droplets were on a crossing trajectory, they were treated as ghost particles and passed through one other. Disturbance flows from neighboring drops were also neglected. Several DNS simulations were performed covering a range of turbulent eddy dissipation rates ε (1.9 cm² s⁻³, 14 cm² s⁻³, and 161 cm² s⁻³). Their results, depicted in Fig. 6, clearly show the spatial clustering of droplets in turbulence as a result of the finite inertia of the droplets, causing them to diverge from regions of high vorticity and preferentially converge in regions of low vorticity. However, the net effect of turbulence on droplet spectral broadening is small when sedimentation is taken into account. Figure 7a shows the standard deviation of the droplet size distribution σ as a function of time for an experiment with an eddy dissipation rate of 161 cm² s⁻³. The value of σ is largest in the simulation without sedimentation (No sed) and smallest in the case when droplet inertia is not allowed (No inertia). The value for the control experiment (full) lies somewhat in between. In general, the standard deviation of the size distribution in the control increases very slowly, at a rate of about 0.01 μm per minute. The underlying reason is that droplets do not stay long in regions of low concentration and high supersaturation. Faster growth at one instant is offset by slower growth at another instant when droplets sediment into regions of high concentration and low supersaturation such that the integrated supersaturation perturbations along the trajectory of the droplet become quite small. Similar results were also reported by Lanotte et al. (2009) and Sardina et al. (2015). However, another broadening mechanism – eddy hopping (Grabowski and Abade, 2017), which describes the effect of large turbulent eddies on broadening the droplet size spectrum through condensation – has gained increasing attention in recent years. The theory was originally proposed by Cooper (1989), suggesting that cloud droplets arriving at the same location inside a cloud may experience different Lagrangian supersaturation histories and form a broad size distribution. Grabowski and Abade (2017), using a turbulent adiabatic parcel model with a varying length scale, found that for a length scale of less than 10 m, the effect of turbulence is negligible, which is consistent with the DNS findings. However, as the length scale increases to 10–100 m, which is typical of current LES grid sizes and exceeds the DNS domain length, the turbulent eddies start to have an impact on the mean droplet radius and spectral width. The impact continues to amplify as the domain size and eddy dissipation rate increase (Fig. 7b), indicating the importance of large-scale turbulence on the condensational broadening.

Figure 6.  Cross section of field of vorticity magnitude (left panel) and positions of 20-μm droplets in the same cross section (right panel). Note clustering of droplets in low vorticity regions. Adapted from Vaillancourt et al. (2002), ©American Meteorological Society, used with permission.

Figure 7.  (a) Time evolution of the standard deviation of the droplet size distribution (from Vaillancourt et al., 2002). (b) Mean radius and spectral width at t = 1000 s as a function of the scale L for different eddy dissipation rates $\epsilon ({\rm{c}\rm{m}}^{2}\;{\rm{s}}^{-3} )$ (adapted from Grabowski and Abade, 2017, ©American Meteorological Society, used with permission).

It has long been postulated that small-scale turbulence enhances droplet collisional growth. There are three major effects affecting collision that are represented by three statistical parameters: a) the turbulent transport effect which increases the droplet radial relative velocity W_r through the local shear and air acceleration; b) the clustering effect measured by the radial distribution function g($ {r}_{1}+{r}_{2} $) which redistributes the droplets to cluster in regions of low vorticity and high shear because of their inertia, or in regions of low Lagrangian acceleration; c) the droplet–droplet interaction effect which affects the collision efficiency E($ {r}_{1},{r}_{2} $) or the probability that a droplet residing in a collector droplet’s effective swept-out volume will collide with the collector droplet. Specifically, W_r is the difference between the velocities of the colliding droplets in the direction of collision, and clustering occurs when g($ {r}_{1}+{r}_{2} $) >1. To calculate the collisional growth of cloud droplets using the continuous or stochastic collection equation, it is necessary to have knowledge of the collision kernel (see Rogers and Yau, 1989), which includes the product of the three parameters. Specifically, the effective collision kernel is defined as $ {{\varGamma }}^{\mathrm{E}\mathrm{F}\mathrm{F}} =$ $ {\varGamma }^{\mathrm{G}\mathrm{E}\mathrm{O}}E\left({r}_{1},{r}_{2}\right), $ where the geometric collision kernel is given by ${\varGamma }^{\mathrm{G}\mathrm{E}\mathrm{O}}=2\pi {\left({r}_{1}+{r}_{2}\right)}^{2}⟨\left|{W}_{r}\right|⟩g\left({r}_{1}+{r}_{2}\right)$ as formulated by Wang et al. (1998). Details for calculating W_r, g($ {r}_{1}+{r}_{2} $), and ${\varGamma }^{{\rm{GEO}}}$ in DNS can be found in Chen et al. (2016b).

Some recent DNS studies on droplet collision statistics with turbulence include those of Franklin et al. (2005), Wang et al. (2008), Ayala et al. (2008b), and Rosa et al. (2013); see Onishi and Seifert (2016) for a limited range of computational Taylor Reynolds numbers R_λ (related to number of grid points N) and eddy dissipation rates ε. More complete calculations covering R_λ from 63 (N =64³) to 589 (N=1024³) and ε values of 50 cm² s⁻³ , 200 cm² s⁻³, 500 cm² s⁻³, 1000 cm² s⁻³, and 1500 cm² s⁻³ are reported in Chen et al. (2016b). Figure 8 depicts the geometric collision kernel between 10-μm and 20-μm droplets as a function of R_λ for different values of ε. Even though the standard deviation of the statistics increases with R_λ which can be induced by the intensified intermittency of the turbulence as the domain size increases, the behavior of the curves and the calculation of the vorticity kurtosis, useful to quantify the level of intermittency, indicate that the effect of intermittency is secondary, at least over the simulated range of R_λ. The fluctuations in the curves were caused by slight variations in the dissipation rate in the simulations.

Figure 8.  Geometric collision kernel of 10–20-μm collisions at different turbulent dissipation rates ε (adapted from Chen et al., 2016b, ©American Meteorological Society, used with permission).

Figure 9 shows the collision statistics W_r, g($ {r}_{1}+{r}_{2} $), and $ {\varGamma }^{\mathrm{G}\mathrm{E}\mathrm{O}} $ between 10-μm and 20-μm droplets for different eddy dissipation rates. As indicated, all collision and pair statistics increase monotonically with ε. The fluctuation of statistics among different resolutions also grew slightly with ε, which can be related to the increasing intermittency of turbulence. The enhancements of the geometric collision kernel (relative to the non-turbulent gravitation collision kernel) for the 10-μm and 20-μm droplet pairs from different sources are compared in Fig. 10. The solid lines denote three DNS results from Franklin et al. (2005, F05), Ayala et al. (2008b, A08a), and Chen et al. (2016b, DNS-MPI). The dashed lines are the values from various parameterizations presented in Saffman and Turner (1956, ST56), Wang et al. (1998, W98), W98 with the Zhou et al. (2001) clustering parameterization (W98&Z01), Franklin et al. (2007, F07), Ayala et al. (2008a, A08b), and Chen et al. [2016, Formula (8)]. Formula (8) in Chen et al. (2016b) has the form:

Figure 9.  Statistics of (a) geometric collision kernel, (b) radial distribution function g, and (c) radial relative velocity between 10- and 20-μm droplets. The error bars show the standard deviation of the statistics with different values of R_λ. Adapted from Chen et al. (2016b), ©American Meteorological Society, used with permission.

Figure 10.  Enhancement of turbulent collision kernel for droplet pairs of 10–20 μm. Solid lines represent DNS simulation results. Dashed lines represent values given by different parameterizations. Adapted from Chen et al. (2016b), ©American Meteorological Society, used with permission.

where $ {\tau }_{pi} $ and $ \mathrm{S}{\mathrm{t}}_{i} $ are, respectively, the response time of the colliding droplets and the Stokes number defined as the ratio of the particle response time to the Kolmogorov time scale $ {\tau }_{k} $. The first term on the right-hand side is the gravitational collision kernel involving two droplets of radius r₁ and r₂; the term is $ \pi R^2\left|{V}_{t2}-{V}_{t1}\right| $, with R = r₁ + r₂ and $ {V}_{ti}=1.223\times {10}^{8}{r}_{i}^{2} $ being the terminal velocity of a droplet with radius r_i. The second term is the turbulence-enhanced component, which is the coupling of the clustering effect and the turbulent transport effect. This parameterized collision kernel is a function only of the droplet size and the dissipation rate ε (as the Kolmogorov length scale $ \eta $ and time scale $ {\tau }_{k} $ can be derived from ε). Unlike other parameterizations, the Taylor microscale Reynolds number $ {R}_{\lambda } $ and the turbulent root-mean-square velocity $u{'}$ are absent from the above formula. As demonstrated by Chen et al. (2016b), the quantities $ {R}_{\lambda } $ and $ u{'} $ in DNS depend on the computational domain size, with a larger domain corresponding to a larger $ {R}_{\lambda } $ and $u{'}$. For droplet collision, the relevant eddy sizes are those of the mean separation distance between droplets, which is of the same order of magnitude as the Kolmogorov length scale in clouds. As long as the Kolmogorov scale can be adequately resolved in DNS, the effects of larger eddies (or larger domain size) are not significant to the droplet collision problem. Therefore, the parameterization should be free of those domain-size dependent parameters.

In general, turbulence enhances the collision kernel. In terms of the DNS results, the enhancement in DNS-MPI is almost the same as A08a up to ε = 400 cm² s⁻³. The collision kernel from F05 is overestimated. For the parameterizations, Eq. (40) is the best, followed closely by A08b and W98&Z01. F07 is larger than all other results when ε > 200 cm² s⁻³. W98 and ST56 did not consider the clustering effect and underestimated the collision kernel over the range of ε examined (from 150 cm² s⁻³ to 1500 cm² s⁻³).

Knowledge of the collision efficiency is required in the effective collision kernel. In the past, the non-turbulent collision efficiencies were mainly obtained by numerical studies involving a pair of droplets interacting in still air. For turbulent collision efficiency, there have been various statistical modeling studies (e.g., De Almeida, 1979; Koziol and Leighton, 1996; Pinsky et al., 2007, 2008), but their results are not in general agreement because of different assumptions about the statistical description of the turbulent flow and the relative velocity of the droplets. Until very recently, DNS studies of turbulent collision efficiency have been few and have had their limitations. The results of Onishi et al. (2013) and Ayala et al. (2014) lack quantitative evaluation. The work of Wang et al. (2008) is limited to two values of ε, at 100 cm² s⁻³ and 400 cm² s⁻³. It also did not consider collector droplet size of less than 20 μm, and the cloud liquid water content specified was much higher than typical adiabatic values. Recently, Chen et al. (2018b) exploited the relative insensitivity of the collision statistics to the computational Reynolds number found in Chen et al. (2016b) and computed the complete collision efficiency covering the whole droplet size range (radius less than 30 μm) important for warm rain initiation. Their work investigates the disturbance flow of neighboring drops and covers the case with no turbulence and cases with turbulence at ε of 20 cm² s⁻³, 50 cm² s⁻³, 200 cm² s⁻³, and 500 cm² s⁻³. The cloud liquid water content specified was also in agreement with observations.

Figure 11 shows an example of the collision efficiency for collector droplets of size R colliding with collected droplets of size r in turbulence with different ε. The results from Pinsky et al. (2008) and Wang et al. (2008) are also included for comparison. For R = 10 μm and 15 μm, the collision efficiency increases monotonically with an increase in radii ratio r/R. For R = 20 μm and 25 μm, the collision efficiency increases with an increase in the radii ratio until r/R~0.6. Thereafter, the collision efficiency first decreases and then increases again, especially for r/R getting closer to unity. Quantitatively, the Chen et al. (2018b) efficiencies agree well with the results of Pinsky et al. (2008) in most of the cases but are slightly larger for the R = 20 μm case with (r/R) < 0.7. On the other hand, the results from Wang et al. (2008) are somewhat larger than the Chen et al. (2018b) results for R = 20 μm and equivalent eddy dissipation rates.

Figure 11.  Turbulent collision efficiency for different collector droplet radii R as a function of the radii ratio r/R with r being the size of the collected droplets. The solid curves with colors denote the results for different eddy dissipation rates, with the label Gravity denoting the no-turbulence case. The dashed curves are the results from Pinsky et al. (2008, P10 for ε=10 cm² s⁻³ and P200 for ε=200 cm² s⁻³) and Wang et al. (2008, W100 for ε=100 cm² s⁻³ and W400 for ε=400 cm² s⁻³) for comparison. Adapted from Chen et al. (2018b), ©American Meteorological Society, used with permission.

Figure 12 depicts the turbulence enhancement factor, defined as the collision efficiency normalized by its gravitational (non-turbulent) value for ε = 200 cm² s⁻³. Enhancement is greatest for similar-sized collisions, indicating that turbulence has its strongest influence in modifying the hydrodynamic interactions between droplets of similar sizes. In other words, turbulence effects exert their largest contribution in enhancing the collision rates between similar-sized droplets. Another finding is that the enhancement is highly sensitive to the radii ratio but only depends weakly on the size of the collector size R for a fixed radii ratio.

Figure 12.  Turbulence enhancement factor of collision efficiency for eddy dissipation rate of 200 cm² s⁻³. The results from Pinsky et al. (2008) and Wang et al. (2008) are included for comparison. Adapted from Chen et al. (2018b), ©American Meteorological Society, used with permission.

The initiation of warm rain in an observed time scale of 15–20 minutes has remained a classical cloud physics problem for decades. Classical condensation theory yields too long a time scale to form raindrops as the time rate of change of radius is inversely proportional to the radius in condensational growth of cloud droplets. As mentioned previously, (Vaillancourt et al., 2002) found that the inhomogeneous distribution of water vapor supersaturation, caused by the interaction of turbulence and the inhomogeneous spatial distribution of droplets in turbulence, produces too little broadening to explain the observed cloud droplet size distribution. However, they showed that the interaction between cloud droplets and turbulence can result in preferential concentration or clustering of cloud droplets that can enhance the growth of the particles by collision. Using DNS methodology, Chen et al. (2016b) showed that turbulence can significantly enhance the radial distribution function and the radial relative velocity of droplets in turbulent flows. Chen et al. (2018b) produced the most complete quantitative results on the turbulent enhancement of collision efficiency and showed that the enhancement is strongest among comparable-sized droplet collisions. Recognizing that condensational growth produces a narrow drop size distribution while turbulence enhances similar-sized droplet collisional growth, Chen et al. (2018c) performed the first sets of DNS experiments including both condensation and collision-coalescence growth in turbulence and examined their effects on the evolution of the cloud droplet size distributions. Three sets of experiments are performed: 1) condensation-only experiments, 2) collision-only experiments, and 3) condensation-collision experiments. The initial cloud droplet size spectrum is adopted from aircraft measurements in non-precipitating cumulus clouds (Raga et al., 1990) with an initial droplet concentration of 80 cm⁻³. The air parcel is assumed to rise with a constant mean upward velocity of 2.5 m s⁻¹. Figure 13 illustrates an example of the DNS droplet size spectra at 6.5 min of simulation time for an eddy dissipation rate of 500 cm³ s⁻². As stated in Chen et al. (2018c), “Consistent with past findings, the turbulence effect on droplet condensational growth is small. The condensation-only process produces the narrowest size distribution (red curve) among the three experiments and droplets grow no larger than 20 μm.” On the other hand, in both the collision-only experiment (blue curve) and the condensation–collision experiment (yellow curve), a substantial number of large droplets are found. Furthermore, compared to the collision-only simulation, the condensation–collision experiment generates substantially larger droplets both in number and in size. The largest radius reaches 100 μm compared to less than 65 μm in the collision-only case. Meanwhile, the number concentration of r >30 μm droplets in the condensation–collision case increases by a factor of 2.3 (0.35 cm⁻³ compared to 0.15 cm⁻³ in the collision-only case). A detailed analysis of the results indicated that the condensational process produces similar sized droplets, while turbulence enhances comparable-sized collisions thereby accelerating the broadening of the drop size distribution in a turbulent flow. This positive impact of condensation on turbulent collisions is most significant as it represents a plausible mechanism to explain the fast appearance of precipitation in observed warm clouds.

Figure 13.  Droplet size spectra at 6.5 min of DNS simulation for the condensation-only experiment (red), collision-only experiment (blue), and the condensation-collision experiment (yellow). The initial droplet size distribution is shaded. Adapted from Chen et al. (2018c), ©Author(s) 2018.

Another area of DNS applications is studying the turbulent entrainment-mixing processes and the between-clouds and surrounding environmental air and their subsequent effects on microphysical properties. In a series of publications (Andrejczuk et al., 2004, 2006, 2009), Andrejczuk and his coauthors conducted DNS studies of turbulent entrainment-mixing in decaying turbulence. They examined the effects of initial turbulence kinetic energy (TKE), cloud fraction, and droplet sizes, and explored the relationship of the mixing mechanism to the Damköhler number. De Lozar and Mellado (2013) added more features such as sedimentation and particle inertia in the bulk formulation. In Lanotte et al. (2009) and Celani et al. (2005), a model combining the Eulerian description of the turbulent velocity and supersaturation fields with a Lagrangian population of cloud droplets was used to study condensation and evaporation of cloud droplets in turbulent flows. Similar to Vaillancourt et al. (2002) and Lanotte et al. (2009), Kumar et al. (2012) developed a droplet-resolved DNS to study turbulent entrainment-mixing processes. In their work, a slab-like vapor field was adopted to mimic the supersaturated cloudy area and subsaturated environment, with an artificial isotropic volume forcing being introduced to maintain the turbulence (forced turbulence). The effects of temperature and buoyancy were ignored. Kumar et al. (2014) extended their previous work to both forced and decaying turbulence and found that the buoyancy due to droplet evaporation played a minor role in the mixing process in their simulations. Kumar et al. (2017) further used their DNS constrained by the thermodynamic conditions observed in monsoon convective clouds over the Indian subcontinent to investigate cloud-edge mixing and its impact on the cloud droplet size distribution. Kumar et al. (2018) explored the scale-dependence of entrainment-mixing processes by running DNS simulations with domain sizes from 12.8 cm to 2 m. They found that the mixing mechanisms tended to be more inhomogeneous as the domain size increased, similar to the observational study by Lu et al. (2014b). The shape of the droplet size distributions was also found to vary strongly with spatial scales. Gao et al. (2018) also reported a new droplet-resolved DNS and used it to investigate turbulent entrainment-mixing processes and their parameterization for larger scale models. Unlike most DNSs that have been based on pseudo-spectral methods with known limitations on dealing with realistic boundary conditions and sharp changes, this DNS uses the finite difference method coupled with the WENO (weighted essentially non-oscillatory) scheme (Jiang and Shu, 1996). Among other results, this study reinforced, via droplet-resolved DNS simulations under different conditions, the potential to unify the representation of various entrainment-mixing mechanisms in terms of the transition scale number originally studied by use of aircraft measurements (Lu et al., 2013b).

With ever-increasing computational power, DNS study of cloud physics problems is becoming increasingly feasible and realistic. As discussed above, considerable progress has been made in the understanding of warm rain processes. Future advances are expected in such areas as 1) adding ice particles and investigating the collision processes involving ice or mixed-phase particles in DNS to evaluate the model uncertainty caused by bin microphysics and BMP in LES sub-grid scales; 2) using DNS results to improve microphysical parameterizations, for example, applying the turbulent collision kernel and collision efficiency formula derived from Chen et al. (2016b, 2018b) in the bin scheme, and formulate the turbulent autoconversion rate into the BMP based on the study from Chen et al. (2020b) to account for turbulence effects; and 3) enlarging the DNS domain to bridge with LES simulations (Mellado et al., 2018). In particular, application of the SDM idea holds great potential for considering more physical processes and increasing DNS domain sizes (Li et al., 2017). More research is warranted in applying DNS results to advancing BMPs to consider turbulence–microphysics interactions and turbulent entrainment-mixing processes.

Instead of having a comprehensive review of existing microphysical parameterizations that have been developed over the years, this paper purposely concentrates on several key and emerging topics that, we believe, hold great potential for future advances of parameterization and modeling of cloud microphysics. Microphysical processes occurring in clouds and the types of hydrometeors involved are briefly discussed as well. The key topics are schematically summarized in Fig. 14.

Figure 14.  Schematic illustration of the main cloud microphysical processes in clouds (left panel) and the four ways of representation in numerical models ranging from large-eddy simulation (LES) to numerical weather prediction (NWP) to earth system models (right panel). The acronyms of BMP, SBM, SDM and DNS stand for bulk microphysics parameterization, spectral bin microphysics, super-droplet method, and particle-resolved direct numerical simulation, respectively. Use of DNS aims to understand and develop representation of turbulence-microphysics interactions.

On BMPs, we describe the general features of BMPs, overview the development of BMPs from one-moment to two-moment and three-moment schemes, elucidate the two general pillars in designing BMPs, and present a BMP classification system based on the number of HSD moments and the number of hydrometeor types assumed in BMPs. It is pointed out that although one-moment or two-moment BMPs still dominate microphysical representation in existing models ranging from LES to CRM to GCMs, three-moment BMPs capable of representing the effects of HSD spectral shape are expected to move to the central stage in near future. Studies in several exemplary areas are further highlighted from the perspectives of further improving two-moment BMPs and developing three-moment BMPs, including consideration of turbulent entrainment-mixing processes, and stochastic condensation, and the influences of HSD spectral shape on sedimentation of raindrops and cloud droplets, autoconversion threshold, and cloud radiative properties. Diagnostic expressions for cloud droplet relative dispersion and prognostic three-moment schemes are discussed. Also discussed is the representation of hydrometeor types with an emphasis on the emerging “continuous” representations that aim to replace the conventional separation of discrete types with a single hydrometeor type with continuously varying hydrometeor properties. Theoretical foundations for the two pillars of developing BMPs are presented, including the systems theory based on maximum entropy principle for the Gamma, Weibull, and modified Gamma HSDs, and the self-affine fractals to represent hydrometeor morphology. Finally, the approaches that have been used in developing BMPs are summarized and grouped into six general approaches: 1) educated guess and empirical analysis of observations; 2) theoretical derivation; 3) statistical curve-fitting of detailed model simulations with explicit cloud microphysics; 4) look-up tables from detailed model simulations with explicit cloud microphysics; 5) machine learning BMPs derived from detailed model simulations with explicit cloud microphysics and/or observations; and 6) hybrid approaches that combine two or more of the above five approaches.

After highlighting the essential roles of cloud models with explicit cloud microphysics to improve our understanding and work as references and synthetic data for developing BMPs, we review the traditional bin schemes and emerging particle-based schemes, and discuss their pros and cons. The particle-based cloud models with Lagrangian super-particles are advocated in view of their demonstrated advantages. Nevertheless, like models with bin schemes, most particle-based cloud models are CRMs or LES models wherein a significant portion of turbulence eddies that affects cloud microphysics remains unresolved and thus needs to be parameterized. To explicitly explore turbulence–microphysics interactions, we further introduce a type of studies that only began around 2000 on developing and using particle-resolved DNS models. Some interesting results that likely fill in our knowledge gaps are elucidated, and some efforts at using the DNS results to bring turbulence influences into BMPs are put forth.

Outstanding challenges and future research directions in the above-mentioned areas are also explored. Summarized below are several points worthy of special attention.

● Two-moment BMPs will continue to serve as the workhorse in the near future for representing cloud microphysics in models of many types (e.g., LES, CRM, and GCM), especially for operational applications. A few areas will likely receive increasing attention in future development, for example, improving representations of turbulent entrainment-mixing processes, stochastic condensation, and solid hydrometeor types.

● Representation of turbulent entrainment-mixing processes and stochastic condensation calls for improving representations of 1) subgrid supersaturation variability by going beyond the conventional treatment of saturation adjustment and 2) turbulent dissipation rate and entrained air properties.

● Our understanding and representation of hydrometeor morphology have been largely based on observations and semi-theoretical modeling of aggregates. It is obviously desirable to link the formation of hydrometeor morphology with the atmospheric fields of the particles embedded. Further development and application of particle-resolved DNS models hold great potential in this crucial area.

● Development of BMPs that treats three or more moments and thus the influences of HSD spectral shapes will become more important and find increasing use thanks to expected advances in computational power and physical understanding.

● To some extent, studies on BMPs, microphysical modeling, and numerical models (LES, CRM, and GCM) all started around the 1960s for different purposes following the advent of electronic computers. It is still anticipated that bin microphysics, particle-based microphysics, and particle-resolved DNS models will advance together and in parallel with multi-moment BMPs. Note that there were also studies that attempted to track a very limited number of realistic Lagrangian hydrometeors in 3D Eulerian cloud models (Xu and Duan, 2002), which can be viewed as particle-based microphysics with compromised representation of the number of hydrometeors in clouds. Comparative studies of BMPs, bin schemes, and particle-based schemes, and the discerning reasons underlying their discrepancies are valuable to further advancements of all three types of microphysical parameterizations and thus numerical simulations of clouds, weather, and climate.

● Machine learning techniques are expected to have an increasing role in developing multi-moment BMPs from detailed cloud model simulations with explicit bin and/or super-particle microphysics, and their coupling with cloud dynamics. Advances in machine learning and high-performance computing may transform the development of microphysical parameterization, its coupling with other processes such as convection, and detailed modeling itself (Schneider et al., 2017).

● BMPs have been largely used in various models of different types and with different spatiotemporal resolutions, without much differentiation except for concerns over computational resources and perhaps tuning the empirical parameters embedded in the BMP involved. However, the resolution dependence of BMP parameters has not been systematically investigated. Furthermore, there exist a potential conceptual mismatch in using explicit microphysics in large-scale models that cannot resolve cloud-scale motions, and scale mismatch in turbulent clouds whose properties are deemed scale-dependent (Liu et al., 2002). Thus, scale-aware BMPs will become increasingly vital as different types of models are unifying with adaptive or unstructured meshes (Chosson et al., 2014).

● In addition to applying existing particle-resolved DNS models to improve our physical understanding and parameterization of turbulence–microphysics interactions, further DNS development is necessary to consider more physical processes (e.g., ice processes) and enlarge the model domain sizes, e.g., from the current < 1 m³ to 10 m³. The increased domain size will provide unique opportunities to investigate the effects of the Reynolds number and bridge with LES models. Optimal integration between DNS and LES is itself an interesting research topic of increasing significance. Ideas of super-particles, treating realistic boundary conditions, and physics-informed machine learning likely find applications in this endeavor.

● Common to BMPs, explicit microphysics and particle-resolved DNS modeling are the critical issues of aerosol–cloud interactions. Addressing them demands adequate representations of not only cloud microphysics but also aerosol microphysics, aerosol activation into cloud droplets or ice crystals, and deactivation of cloud particles back into aerosol particles. This will further involve the other part of the loop – wet removal of aerosol particles. Also moving to the forefront are the issues of the aerosol–cloud continuum and the transitional zone between clouds and their surrounding environments (Charlson et al., 2007; Marshak et al., 2021), which may hold a key to reducing the stubborn uncertainties in aerosol indirect effects and in remote sensing of clouds and aerosols. Addressing these issues further calls for unified treatments of aerosol processes, cloud microphysics, and those processes that affect cloud edges, including turbulence–microphysics interactions (discussed in sections 3.3 and 5), and the potential radiative effects on cloud microphysics due to the conspicuous differences of radiative properties between a cloud and its surrounding atmosphere (Zeng et al., 2022).

● Although observational analysis and various modeling studies will continue to play essential roles in developing microphysical understanding and improve microphysical parameterizations, theoretical formulations are strongly recommended whenever possible to provide clear physics and to reduce the number of tunable parameters.

● Representation of cloud microphysics has been rigorously considered, mainly for large-scale stratiform clouds in most GCMs. The coupling of cloud microphysics with deep convection schemes has just recently started (Zhang and Song, 2016), and there has been little work on coupling cloud microphysics with shallow cumulus clouds. In view of the importance of cloud microphysics on aspects of parameterized convective clouds such as detrainment and role of entrained aerosols (Chen et al., 2020a), improvement and more research are definitely needed. This will become especially true for shallow cumulus convection as model resolutions increase.

● This review is focused on the topics of our own research interests, and thus, are likely biased. Other important yet not discussed topics include turbulence effects on ice processes, microphysics–convection coupling, cloud chemistry, cloud electricity, and radiative effects on cloud microphysics (Harrington et al., 2000; Marquis and Harrington, 2005; Zeng, 2022. Also not discussed are the advances in measurement technologies (both in situ and remote sensing), laboratory experiments, and model evaluation, which are equally essential to identifying parameterization deficiencies, tuning model parameters, obtaining physical understanding, and eventually advancing numerical models. Developing digital twins in parallel presents growing opportunity to transform smart integration of numerical models and measurements across different scales (Bauer et al., 2021).

Acknowledgements. Yangang LIU is supported by the US Department of Energy (DOE)’s Office of Science Atmospheric Systems Research (ASR) Program, the Office of Energy Efficiency and Renewable Energy (EERE) Solar Energy Technologies Office (SETO) award (33504), and the Brookhaven National Laboratory (BNL)’s Laboratory Directed Research & Development Program (LDRD) (22-065). The Brookhaven National Laboratory is operated by the Brookhaven Science Associates, LLC (BSA), for the US Department of Energy under Contract No. DESC0012704. The views expressed in this article do not necessarily represent the views of the Department of Energy or the US Government. Shin-ichiro SHIMA is supported by JSPS KAKENHI Grant No. 26286089 and MEXT KAKENHI Grant No. 18H04448.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

This Appendix summarizes and compares the key relationships between the HSD parameters and HSD moments for the commonly used Gamma size distribution, Weibull size distribution, and Lognormal size distribution. Also summarized are some empirical expressions.

For the modified Gamma size distribution given by n(r)=N₀r ^μexp(−λr^q), we have

where M_p (μ_p) denotes the (normalized) pth moment, and N₀, $\mu$, λ, and q are the distribution parameters. The Gamma and Weibull size distributions correspond to q = 1 and μ = q−1, respectively.

For the lognormal distribution given by $n\left(r\right)= $ $ \dfrac{N}{{\ln}(\chi )\sqrt{2\pi }} \dfrac{1}{r} \text{exp}\left(-\dfrac{1}{2}{\left(\dfrac{{\ln}\left(r\right)-{\ln}\left({r}_{0}\right)}{{\ln}(\chi )}\right)}^{2}\right)$, we have:

where the parameters r₀ and x are the geometric mean radius and standard deviation, respectively.

Application of the equations for the normalized moments and definitions of the commonly used statistics such as relative dispersion readily yields the most commonly used expressions in terms of the different distribution parameters for the Gamma, Weibull, and Lognormal distributions, which can be rewritten in terms of relative dispersion, a parameter of clear physical meaning. This variable transformation not only permits straightforward comparison between the parameterizations derived from these distinct assumed distribution forms, but also unifies the treatment of the spectral effects on microphysical processes or properties such as effective radius and autoconversion rate. As an example, Table A1 summarizes the key expressions used in parameterization of effective radius for the commonly used HSD functions. The monodisperse and Gaussian HSDs are also provided as references.