Abstract

Many two-moment bulk schemes use a three-parameter gamma distribution of the form N(D) = N0DαeλD to describe the size spectrum of a given hydrometeor category. These schemes predict changes to the mass content and the total number concentration thereby allowing N0 and λ to vary as prognostic parameters while fixing the shape parameter, α. As was shown in Part I of this study, the shape parameter, which represents the relative dispersion of the hydrometeor size spectrum, plays an important role in the computation of sedimentation and instantaneous growth rates in bulk microphysics schemes. Significant improvement was shown by allowing α to vary as a diagnostic function of the predicted moments rather than using a fixed-value approach. Ideally, however, α should be an independent prognostic parameter.

In this paper, a closure formulation is developed for calculating the source and sink terms of a third moment of the size distribution—the radar reflectivity. With predictive equations for the mass content, total number concentration, and radar reflectivity, α becomes a fully prognostic variable and a three-moment parameterization becomes feasible. A new bulk microphysics scheme is presented and described. The full version of the scheme predicts three moments for all precipitating hydrometeor categories.

Simulations of an idealized hailstorm in the context of a 1D kinematic cloud model employing the one-moment, two-moment, and three-moment versions of the scheme are compared. The vertical distribution of the hydrometeor mass contents using the two-moment version with diagnostic-α relations are much closer to the three-moment than the one-moment simulation. However, the evolution of the surface precipitation rate is notably different between the three-moment and two-moment schemes.

1. Introduction

Given the increasing importance of bulk microphysics parameterization schemes in atmospheric models, it is important to develop detailed yet computationally efficient schemes and to understand the strength and limitation of various approaches. Many two-moment bulk schemes use a three-parameter gamma size distribution of the form N(D) = N0DαeλD. The schemes predict changes to the mass content and total number concentration thereby allowing N0 and λ to vary as prognostic parameters while fixing the spectral shape parameter, α. In Milbrandt and Yau (2005a, hereafter Part I), we examined the impact of changes to the relative dispersion of the hydrometeor size spectrum, as represented by the shape parameter. It was shown that α plays an important role on the computation of sedimentation and instantaneous growth rates in bulk microphysics schemes. A method that diagnoses α as a monotonically increasing function of the mean-mass diameter was introduced and shown to yield notable improvement over the standard fixed-value approach. It was also demonstrated, however, that α should be a fully prognosed parameter by predicting three moments of the size distribution.

In this paper, a closure formulation is developed for calculating the source and sink terms of a third moment of the size distribution, the radar reflectivity. With predictive equations for the mass content, total number concentration, and radar reflectivity, α becomes a fully prognostic variable and a three-moment parameterization becomes feasible. In view of the analysis in Part I, a new multimoment bulk scheme, with options to diagnose or prognose the shape parameter, has been developed and is presented here. In the following section, the proposed three-moment closure is described. Section 3 presents an overview of the new microphysics scheme followed by a complete description of the source/sink terms in section 4. In section 5, simulations of an idealized hailstorm using the various versions of the scheme in a 1D kinematic cloud model are discussed. Concluding remarks are given in section 6.

2. Prognostic relation for α—A proposed closure for the three-moment approach

From the results presented in Part I, it is clear that there are advantages to extend the bulk approach to allow the spectral shape parameter αx, as well as the intercept parameter N0x and the slope parameter λx, of a size distribution for hydrometeor species x to be independent. To predict the evolution of αx, a prognostic equation for another moment such as the radar reflectivity Zx with the form

 
formula

must be added to the continuity equations for the mass and total number concentration [(8) and (9) in Part I], along with the closure equation [(6) in Part I]. The readers are referred to appendix C for a complete list of symbols.

The source term for Zx in (1) is computed as the sum of the individual tendencies of Zx for each microphysical process listed in appendix A. What remains, therefore, is to derive the source terms for Zx for each process A in each hydrometeor category x.

We classify the source terms for Zx into three types. For the first type, it is assumed that the change in αx due to the particular process A is negligible. By differentiating the closure Eq. (6) in Part I, we can relate the tendency of Zx to the tendencies of qx and NTx as

 
formula

Equation (2) is applied to collection, diffusional growth, and melting. The assumption that αx does not change because of a particular process is analogous to holding λx constant for a process in order to relate the fractional change in NTx to the fractional change in qx [e.g., Murakami (1990, hereafter M90); Ferrier (1994, hereafter F94); Harrington et al. (1995); Reisner et al. (1998, hereafter RRB)]. Since αx is a measure of the relative dispersion of the size spectrum, our assumption can be interpreted physically to mean that the change in spectral width is negligible due to that particular process. While strictly speaking this assumption may not be completely valid, the analysis in Part I indicates that sedimentation can change αx very quickly and hence if the effect of sedimentation is correctly modeled, the errors in neglecting the change in α by a process A will likely be small.

The second type of source terms is related to processes in which hydrometeors are being initiated in a category, such as during nucleation. In this case, the initial Zxinit = qxinit = NTxinit = 0 at the beginning of the time step, and qxfin and NTxfin at the end of the time step are known from solving (8) and (9) in Part I [hereafter equation (x) in Part I will be referred to as equation (I.x)]. To obtain the tendency terms for Zx, the value of αx at the end of the time step for process A (αxA) has to be specified. Values for different initiation processes are given in appendix C. Applying (I.6) at the end of a time step, one obtains

 
formula

By writing (dZx/dt)|A = (ZfinxZinitxt), substituting (3) for Zxfin, and taking the limit as Δt → 0, we obtain

 
formula

The third type of source terms occur when one category x is converted into another category y. In this case,

 
formula

An example is the probabilistic freezing of rain to hail, where the tendency for Zr is computed from (2) while the tendency of Zh is computed from (5). The total radar reflectivity due to such a process is thus conserved.

3. Overview of the new scheme

A new multimoment bulk scheme has been developed. The scheme consists of six hydrometeor categories. As is standard in bulk parameterizations, the liquid water spectrum is partitioned into cloud, consisting of small nonsedimenting droplets, and rain, consisting of sedimenting drops. It has been shown by McCumber et al. (1991) that to properly model the ice phase, at least four frozen hydrometeor categories should be included. In view of this (and following F94), the proposed scheme includes ice, snow, graupel, and hail. The ice category represents pristine ice crystals. The snow category includes large crystals (with radii greater than 100 μm) and aggregates. The graupel category includes moderate-density graupel, formed from heavily rimed ice or snow. The hail category includes high-density hail and frozen raindrops. The size distribution of cloud droplets is described by (I.1) with a fixed value of νc = 3, following Cohard and Pinty (2000a; hereafter CP00a). The size distributions of rain and of all frozen hydrometeors are described by (I.3)–(I.6). All hydrometeor categories x, with x ∈ [c, r, i, s, g, h] referring to cloud, rain, ice, snow, graupel, and hail, respectively, have the mass–diameter relationship mx(Dx) = cxDxdx (with dx = 3). Except for ice, all particles are assumed to be spherical with cx = (π/6), ρx where ρx is the bulk density of the particles summarized in Table 1. Ice crystals are assumed to be bullet rosettes, which are believed to be the dominant crystal habit in thunderstorms, and have ci = 440 kg m−3 (F94).

Table 1.

Bulk densities and minimum particle sizes of hydrometeor categories in the proposed scheme.

Bulk densities and minimum particle sizes of hydrometeor categories in the proposed scheme.
Bulk densities and minimum particle sizes of hydrometeor categories in the proposed scheme.

There are three main versions of the proposed scheme. The single-moment version predicts qx using (I.8) allowing λx defined by (I.5) to vary prognostically, and has fixed values for NTc and N0x for the precipitating categories and αx for all categories. The two-moment version includes predictive equations for qx and also for NTx, given by (I.9), where the corresponding prognostic parameters are N0x, defined by (I.4), and λx. The parameter αx can be either fixed or diagnosed from (I.13) and (I.14). The constants for these equations are listed in Table 1 of Part I. The three-moment version also has the predictive Eq. (1) for Zx for all categories except cloud water, and αx is obtained by solving (I.6).

Sedimentation is computed following the equations described in the appendix of Part I. The terminal fall velocity parameters for each category are listed in Table 2. The microphysical source/sink terms in the continuity equations listed in appendix A were taken and adapted from various existing schemes. The equations for the warm-rain processes follow closely those of Cohard et al. (1998, hereafter CPB98) and CP00a. The expressions for ice phase processes are mainly adapted from Cotton et al. (1986, hereafter C86), F94, Kong and Yau (1997, hereafter KY97), Lin et al. (1983, hereafter LFO), Meyers et al. (1997, hereafter M97), and M90. Various modifications and simplifications to many of the parameterizations were made, as detailed below.

4. Description of source/sink terms

The source/sink terms described in this section are listed in appendix A. The notation for the terms involving two interacting categories are denoted by VABy,x, where V is the prognostic variable under consideration (Q for mixing ratio, N for total concentration, or Z for radar reflectivity), AB represents the microphysical process (CL for collection, CN for conversion, FZ for freezing, IM for ice multiplication, ML for melting, NU for nucleation, SH for shedding, VD for diffusional growth), and the subscript y, x indicates that mass is being transferred from category y to x [for x, y ∈ (c, r, i, s, g, h, υ), where υ denotes water vapor]. If the tendency for a prognostic variable is not the same for the interacting categories, then a subscript representing the affected category is added to the prognostic variable. For example, VxABy,x is the rate of change of variable V in category x due to ABy,x while VyABy,x is the rate of change of variable V in category y for the same process. If the tendencies are the same for the two interacting categories, no subscript is added to the prognostic variable.

For source/sink terms involving three interacting categories like three-component freezing where the destination category may be different from the two interacting categories, a three-letter subscript follows the microphysical process. For example, NCLy,x,z is the tendency for total number concentration due to the collection growth involving categories x and y giving rise to category z.

a. Cloud nucleation (NNUυc)

CPB98 developed an expression relating the concentration of nucleated cloud droplets NCCN in ascending air to the maximum supersaturation. CP00a described an iterative procedure to solve for the maximum supersaturation as a function of updraft speed w, temperature T, and pressure p. Given the total number concentration of cloud droplets NTc at a certain time, the cloud nucleation rate is

 
formula

To reduce computational cost, we performed a least-square fit to NCCN(w, T, p) of CPB98 solved by the iteration procedure of CP00a.

b. Condensation and evaporation of cloud and rain (QVDυc, QVDυr, NVDυc, NVDυr, ZVDυr)

The rate of change in cloud mixing ratio due to nucleation and subsequent growth by condensation, QVDυc, is parameterized following the saturation adjustment technique of KY97. For the rate of evaporation of rain in subsaturated air, we write

 
formula

where

 
formula

is the bulk ventilation coefficient for rain and

 
formula

is the thermodynamic function (Byers 1965).

During evaporation, NTc and NTr decrease at the rates of NVDυx, which is given by an equation of the form of (B1) in appendix B for x ∈ [c, r] and Zr decreases at the rate of ZVDυr, given by an equation of the form of (B3).

c. Warm-rain collection (QCNcr, QCLcr, NcCNcr, NrCNcr, NCLrr, ZCNcr, ZCLcr, ZCLrr)

The source terms for the mixing ratios and number concentrations of cloud and rain due to the warm-rain collection process come from CP00a, which are based on analytic solutions to the stochastic collection equation (SCE) using Long’s (1974) polynomial approximation for the collection kernel. A minor modification was made to the way autoconversion is computed (J.-P. Pinty 2001, personal communication).

The qx tendency for autoconversion, QCNcr, is adopted from CP00a. To compute the change in concentration of rain particles due to autoconversion, NrCNcr, we note that in general, for a microphysical process generating new particles of mean diameter Dmr_new, its NTx tendency is related to its qx tendency by (B2) where Dx0 = Ddxmx_new. To determine Dmr_new of rain due to autoconversion, a threshold mean diameter Dmr_aut is first computed using the empirical formula of Berry and Reinhardt (1974),

 
formula

where σc is the standard deviation of the size of the cloud droplets (see CP00a). If the existing mean-drop diameter Dmr = [(ρqr/crNTr)]1/3 exceeds Dmr_aut, then Dmr_new = Dmr; otherwise Dmr_new = Dmr_aut. Thus, NrCNcr is given by (B2) with Dx0 = D3mr_new and the equation for ZCNcr is of the form of (B4) with αxAB = αrAUT.

For the qx and NTx tendencies due to accretion (QCLcr and NCLcr) and self-collection (NcCNcr and NCLrr), we exactly follow CP00a. Note that the NTc tendency due to autoconversion is parameterized as cloud droplet self-collection. Spontaneous breakup of rain is parameterized by including a mean-drop-size limiter in the form of a size-dependent reduction factor in the bulk collection efficiency in NCLrr (see Cohard and Pinty 2000b). The Zr tendencies for accretion (ZCLcr) and self-collection (ZCLrr) are identical to the form of (B3).

d. Collection involving ice particles

1) General collection equations

The rates of collection amongst frozen and liquid particles are calculated by integrating the following SCEs (Walko et al. 1995, hereafter W95; M97) for particle category x collecting category y:

 
formula
 
formula

Bulk collection efficiencies, Exy, are used. For the collection of cloud droplets by graupel and hail, an approximate empirical formula from Macklin and Bailey (1966) is used,

 
formula

where x ∈ [g, h]. For dry collection among ice particles, we adopted the following temperature-dependent collection efficiencies following Ferrier et al. (1995):

 
formula
 
formula

However, during wet growth of hail, we set Eih = Esh = 1.

The values for the other bulk collection efficiencies are Eci = Eri = 1, Eii = 0.1 Eis, Ecs = 1, Ers = 1, Ess = 0.1 Eis, Erg = 1, Egg = 0, Erh = 1, Egh = 0, and Ehh = 0.

(i) Collection amongst rain and frozen categories (QCLri, QCLir, QCLrs, QCLsr, QCLrg, QCLgr, QCLrh, QCLis, QCLig, QCLih, QCLsh, NCLri, NCLir, NCLrs, NCLsr, NCLrg, NCLgr, NCLrh, NCLis, NCLig, NCLih, NCLsh, ZCLri, ZCLir, ZCLrs, ZCLsr, ZiCLis, ZsCLis, ZiCLig, ZgCLig, ZiCLih, ZhCLih, ZsCLsh, ZhCLsh, ZCLrg, ZCLgr, ZrCLrh, ZhCLrh)

For collection amongst sedimenting categories, M90 proposed the following approximation:

 
formula

Using (16) and the bulk collection efficiencies, (11) and (12) can be integrated analytically to yield

 
formula

and

 
formula

and where ZCLyx is of the form of (B3).

Equations (17), (18), and the equation for ZCLyx apply to the collection of rain by all frozen categories (CLrx), collection of ice, snow, or graupel by rain (CLxr) and collection of one frozen category by another frozen category (CLxy), except for hail when it is undergoing wet growth (see section 4c). Equations of the form of (B3) are also used to compute ZCLcr and ZCLrr (with QCLrr = 0).

(ii) Collection of cloud water by frozen categories (QCLci, QCLcs, QCLcg, QCLch, NCLci, NCLcs, NCLcg, NCLch, ZCLci, ZCLcs, ZCLcg, ZCLch)

For the collection of cloud by any frozen category, the approximation for ΔV need not be made because cloud droplets are assumed to have negligible terminal fall velocity. Thus,

 
formula
 
formula

and ZCLcx is calculated using (B3) with NCLyx = 0, x ∈ [i, s, g, h], and y = c.

2) Aggregation for snow (NCLss, ZCLss)

For snow aggregation, the rate of decrease in NTs (NCLss,) follows F94. The corresponding rate of change in radar reflectivity, ZCLss, is of the form of (B3) with QCLss = 0.

3) Wet growth of hail (NSHhr, ZSHhr, QCLrh)

In nature, a hailstone grows mainly by accreting liquid water. However, its surface temperature can increase if the hailstone cannot dissipate all of the latent heat released due to freezing and deposition. When the surface of the hailstone warms to the point that all of the accreted water cannot be frozen, it is said to have reached the Shumann–Ludlam limit (SLL; Young 1993). Beyond the SLL, the hailstone enters the wet growth mode and begins to shed some of the accreted water.

For dry growth, the mass growth rate of hail is equal to QCLch + QCLrh + QCLih + QCLsg, with each of the terms calculated using (18). For wet growth, the approach of LFO is followed to determine the criteria for wet growth and the rate for mass increase, QhWET, based on the heat balance equation (Musil 1970). The actual growth rate is chosen to be the smaller of the two rates. If the wet growth rate is chosen, then QCLch, QCLih, and QCLsh are recomputed using collection efficiencies of 1 to obtain QCL′ch, QCL′ih, and QCL′sh. Likewise, the NTx tendencies NCLch, NCLih, and NCLsh are also recalculated using collection efficiencies of 1 and ZCLih and ZCLsh are computed using the appropriate rates in equations of the form of (B3). Finally, the collection rate for rain is recomputed as the difference between QhWET and the sum of the new (wet) growth rates for collecting cloud, ice and snow

 
formula

Also the difference between the dry collection rate QCLrh and the wet collection rate QCL′rh is assumed to be the rate that the collected water mass by the hailstone is shed to form rain (see LFO for discussion).

The change in NTr because of shedding is determined by assuming a mean size for the shed drops which has been observed to be between 0.5 and 2.0 mm in diameter with a modal size of 1 mm (Rasmussen and Heymsfield 1987; Lesins et al. 1980). Applying (11) with Dmr_new = Dshed = 1 mm, we obtain

 
formula

The Zr tendency due to shedding is given by

 
formula

where αrSH = 2.

e. Ice nucleation

Ice crystals are initiated via three modes: primary nucleation, rime splintering, and homogeneous freezing of cloud droplets.

1) Primary nucleation (QNUυi, NNUυi, ZNUυi)

Primary nucleation includes contact nucleation, deposition nucleation, and condensation-freezing nucleation. The NTi tendency due to contact nucleation, NuCONT, is parameterized following Young (1974) as described in C86 and W95. The combination of deposition and condensation-freezing nucleation is parameterized by the empirical formula of Meyers et al. (1992), which gives NTi as a function of the saturation ratio with respect to ice, Si:

 
formula

The NTi tendency due to deposition/condensation-freezing nucleation is thus

 
formula

A method of estimating the supersaturation with respect to ice, as a function of vertical velocity, temperature, and pressure only, and that is not influenced by the model time step is also used, thus reducing excessive ice nucleation. The nucleation rate for primary ice nucleation, NNUυi, is the sum of NuCONT and NuDEPSOR. To prevent overdepletion of water vapor in one time step, a maximum nucleation rate, analogous to the maximum depositional growth rate is imposed (following KY97),

 
formula

Also, primary nucleation can result in an increase in the number of ice crystals present, but will never result in a decrease. Thus,

 
formula

The mixing ratio and reflectivity factor tendencies are calculated using an assumed nucleated ice crystal mass mi0 = 10−12 kg (KY97; RRB), giving

 
formula

and, where ZNUυi is of the form of (B4), with αxAB = αiNU = 0.

2) Ice multiplication (QIMsi, QIMgi, NIMii, NIMsi, NIMgi, ZIMii, ZIMsi, ZIMgi)

Ice multiplication, or rime splintering, for riming of ice, snow, and graupel at temperatures between −3° and −8°C is based on Hallet and Mossop (1974) and is parameterized following the equations of RRB for the qx and NTx tendencies and with the Zx tendencies given by equations following the form of (B4), with αxAB = αiIM and y ∈ [i, s, g]. Note that for rime splintering due to riming of ice, the value of QIMii used to compute ZIMii is formulated the same as for QIMsi and QIMgi (see RRB) but there is no change in the ice mixing ratio [i.e., QIMii = 0 and does not appear in (A4)].

3) Homogenous freezing of cloud droplets (QFZci, NFZci, ZFZci)

Another source of ice crystals is the homogenous freezing of cloud droplets at temperatures below −30°C. DeMott et al. (1994) give the number of droplets ΔNfreeze that freeze in time Δt at a given temperature as

 
formula

where

 
formula

and V is the droplet volume. By substituting V with V = (π/6)Dmc3, the mean-droplet volume in (29), we obtain

 
formula

the fraction of the total cloud concentration that freezes in one time step. Hence

 
formula
 
formula

and the equation for ZFZci follows the form of (B3). For Tc > −30°C, ffr = 0, while for Tc < −50°C, ffr = 1. For continental cloud condensation nuclei (CCN), this approach produces appreciable amounts of freezing at temperatures several degrees warmer than −40°C, which is sometimes used as a threshold temperature for homogeneous freezing (e.g., KY97; RRB).

f. Deposition/sublimation (QVDυi, QVDυs, QVDυg, QVDυh, NVDυi, NVDυs, NVDυg, NVDυh, ZVDυi, ZVDυs, ZVDυg, ZVDυh)

All frozen categories undergo deposition (sublimation) in an environment supersaturated (subsaturated) with respect to ice. The diffusional growth rate for frozen category x ∈ [i, s, g, h] is calculated by

 
formula

where

 
formula

is the thermodynamic function and

 
formula

is the mass-weighted ventilation factor (F94). The second term on the right hand side of (34) represents the decease (increase) in the deposition (sublimation) rate due to the latent heating of the surface of the frozen particle due to accretion of liquid water. To prevent excessive supersaturation or subsaturation resulting from excessive diffusional growth in one time step (which is possible if the time step is large), a maximum diffusional growth rate VDmax is computed similar to that for condensation/evaporation (see KY97).

For particles undergoing sublimation, the number concentrations are reduced at the rate NVDυr given by (B1), following F94, but no change in the concentrations of frozen particles occurs during depositional growth. For both deposition and sublimation the reflectivity factor is changed by ZVDυr, which is of the form of (B3).

g. Freezing of rain

1) Probabilistic freezing (QFZrh, NFZrh, ZrFZrh, ZhFZrh)

When the ambient air temperature is below 0°C, rain can undergo spontaneous or probabilistic freezing to form small, frozen drops or hail embryos (e.g., Ziegler 1985, hereafter Z85). The rate of change of number concentration is given by Bigg (1953) as

 
formula

and the rate of change in mass by

 
formula

where A′ = 0.66 K−1 and B′ = 100 m−3 s−1. The reflectivity change rates are given by ZrFZrh, which is of the form of (B3), and ZhFZrh, which is of the form of (B5).

2) Collisional (three-component) freezing (NCLirg, NCLirh, NCLsrs, NCLsrg, NCLsrh, ZCLirg, ZCLirh, ZCLsrs, ZCLsrg, ZCLsrh)

When T < 0°C, rain drops also freeze when they come into contact with particles in a frozen category x. A simplification of F94’s three-component freezing method is used. A new frozen category z is produced and the density of which (ρz) is calculated as

 
formula

where x ∈ [i, s, g], the destination category z ∈ [s, g, h], and Dmz = max(Dmx, Dmr). We assume here that during contact, water is evenly distributed throughout the volume of the colliding frozen particle to increase its mass (and bulk density) but not its volume. The destination particle is then classified as belonging to the frozen category with the closest bulk density to ρz. It is classified as snow if ρz ≤ 0.5(ρs + ρg), as graupel if 0.5(ρs + ρg) < ρz ≤ 0.5(ρg + ρh), and as hail if ρz > (ρg + ρh).

The sink terms for the mixing ratios, number concentrations, and reflectivities of rain and category x are computed by the collection equations described above. The terms for rain collecting category x are QCLxr, NCLxr, and ZCLxr. The terms for category x collecting rain are QCLrx, NCLrx, and ZCLrx. The source term for the mixing ratio of the destination category z is the sum of the mixing ratio sink terms for the two colliding categories. Thus the delta function variable δxrz = 1 for the destination category z, but is 0 for the other frozen categories. For example, if rain and snow collide to form graupel, then δsrg = 1 while δsrs = δsrh = 0 so the rate of change of the graupel mixing ratio is (QCLrg + QCLgr). The NTz tendency for the destination category is calculated using the mean-mass diameter Dmz. Thus,

 
formula

where ρz is the actual density of category z, not the density calculated from (39) used to determine the destination category. The Zz tendency for the resulting category is given by

 
formula

where αzCL = 0 is the assumed αz for the newly formed particles.

h. Conversion processes

1) Ice to snow (QCNis, NiCNis, NsCNis, ZiCNis, ZsCNis)

Ice is converted to snow by growth of ice crystals from riming and deposition to the size of “embryo” snow particles and also by aggregation. Assuming ms0 = 4.4 × 10−10 kg as the mass of an embryo snow particle, the total conversion rate is the sum of the conversion rates by deposition, riming, and aggregation (M90; RRB),

 
formula

The change in NTs is calculated with the assumption that all snow particles converted from ice have an initial mass ms0, giving the equation for NsCNis of the form of (B2). The change in NTi due to conversion to snow by deposition and riming is computed as the negative rate of change in NTs due to conversion from ice deposition and riming. For conversion by ice aggregation, the NTi tendency is given by M90,

 
formula

where

 
formula

and Dmi is the mean ice crystal diameter, Vi is the fall velocity of ice crystals, Eii is the (dry, temperature dependent) collection efficiency amongst ice crystals and Xdisp is the dispersion of the fall velocity spectrum of the ice crystals (assumed for simplicity to be 0.25, following M90). The total NTi tendency due to conversion to snow is therefore

 
formula

The reflectivity changes are given by ZiCNis, which is of the form of (B3) with x = i, and ZsCNis, which is of the form of (B4) with αxAB = αsCN = 0 and with x = s.

2) Ice to graupel (QCNig, NiCNig, NgCNig, ZiCNig, IgCNig)

Ice is converted to graupel when the riming rate of ice crystals exceeds its depositional growth rate (following RRB). The tendency of ice mass converted to graupel is given by twice the difference between the tendency of growth by riming and that by deposition; that is

 
formula

The NTg tendency is calculated based on an assumed embryo graupel mass mg0 = 1.6 × 10−10 kg, giving the equation for NgCNig of the form of (B2), while NTi is reduced following M90 by NiCNig, which is of the form of (B1). The reflectivity tendencies are given by ZiCNig, which is of the form of (B3) with x = i, and ZgCNig, which is of the form of (B4) with αxAB = αgCN = 0 and with x = g.

3) Snow to graupel (QCNsg, NCNsg, ZsCNsg, ZgCNsg)

Snow is converted to graupel by riming. The snow-to-graupel conversion rate is given by the sum of the riming rate (collection of cloud water by snow) and the rate at which a portion of snow is converted to graupel in unit time following M90. Conversion occurs if the riming rate QCLcs exceeds the depositional growth rate QVDυs. It is assumed that the change in snow particle size due to riming is negligibly small. The conversion rate is

 
formula

The concentration tendency NCNsg is given by (B2) assuming that newly converted graupel particles are embryos with mass mg0 (RRB). This is a source for NTg and a sink for NTs. The reflectivity tendencies are given by ZsCNsg, which is of the form of (B3) with x = s, and ZgCNsg, which is of the form of (B5) with x = g.

4) Graupel to hail (QCNgh, NCNgh, ZgCNgh, ZhCNgh)

When a frozen particle growing by accreting liquid water first reaches the SLL, it is termed a hailstone (Young 1993). The conversion of graupel to hail is parameterized with the premise that a graupel particle growing by accretion becomes redefined as hail the moment it first reaches the SLL. The following exponential function approximates the size of a particle at the SLL, or the embryo diameter of a hailstone, as a function of the environmental temperature and water and ice contents (Z85)

 
formula

For a spectrum of graupel particles, all particles smaller than Dh0 will continue to undergo dry growth while all particles larger than Dh0 undergo wet growth and convert to hail. The conversion rate is thus parameterized by

 
formula

However, (49) needs to be modified as Dh0 can be much smaller than the mean graupel diameter when, for example, the ambient temperature is relatively warm and/or the LWC is relatively high. To prevent the situation that the total mass that is converted to hail becomes larger than the total graupel mass plus the water and ice mass that the graupel accretes, the following limit is placed on the conversion rate:

 
formula

When the ambient temperature is relatively cold and/or the LWC is low, Dh0 is very large and wet growth of graupel and conversion to hail are not expected. A lower limit of 0.1 is therefore placed on the ratio Dmg/Dh0 below which the conversion rate is set to zero. This value was found appropriate in preventing graupel from converting to hail when wet growth is not expected.

The rate of change in number concentrations due to conversion of graupel to hail is calculated from (B2) with Dx0 = D3h0. The NTg tendency is set to the negative of the NTh tendency. The Zx tendencies are given by ZgCNgh, which is of the form of (B3) with x = g, and ZhCNgh, which is of the form of (B5) with x = h.

(i) Melting of frozen particles (QMLir, QMLsr, QMLgr, QMLhr, NMLir, NMLsr, NMLgr, NMLhr, ZiMLir, ZrMLir, ZsMLsr, ZrMLsr, ZgMLgr, ZrMLgr, ZhMLhr, ZrMLhr)

Ice melts instantaneously to rain upon falling into warm (T > 0°C) air. Thus,

 
formula

and

 
formula

For snow, graupel, and hail, the melting rate is based on a heat balance with cooling by melting offset by heating from conduction and convection at the particle surface, latent heat of condensation/evaporation of water to/from the particle surface, and sensible heating from the collected cloud and rainwater (Wisner et al. 1972; LFO). The qx and NTx tendencies are

 
formula

where VENTx is the mass-weighted ventilation factor, and NMLxr, given by (B1). For all frozen categories, the Zx tendencies are given by ZxMLxr, which is of the form of (B3), and ZyMLxr, which is of the form of (B5) with the x and y indices reversed.

5. Results in a 1D kinematic model

For testing and demonstration, the three main versions of the scheme—the one-moment version (SM), the two-moment with the diagnosed-αx (DIAG), and the three-moment (TM) versions—have been interfaced with a 1D kinematic column model and used to simulated an idealized hailstorm. The model was initialized using a conditionally unstable sounding (not shown) to provide the initial temperature and water vapor mixing ratio at each level in a 12-km vertical domain. The 0°C isotherm is at approximately 3.5 km in the initial sounding. A time-dependent, kinematic updraft (shown in the inset of Fig. 1) is prescribed, growing sinusoidally from an initial peak value of 2 m s−1 to a maximum of 20 m s−1 over 15 min and then decreasing to zero. At each time step, all quantities are advected upward and the microphysics scheme is invoked, computing changes to the various hydrometeor moments, the temperature, the water vapor content (as described in section 4 and appendix A) and also sedimentation (as described in the appendix of Part I). To incorporate a constant low-level moisture supply into the domain, the water vapor mixing ratio in the lowest 1 km is not permitted to decrease below half of its initial value. The column model uses a two-time-level semi-Lagrangian advection scheme and the following simulations were performed with a vertical grid spacing of 240 m and a time step of 20 s.

Fig. 1.

Total surface precipitation rates vs time for idealized simulations with a 1D kinematic column model using the one-moment (heavy dashed), two-moment (solid), and three-moment (heavy dot–dashed) versions of the microphysics scheme. (The peak precipitation rate for the one-moment simulation is 340 mm h−1.) Inset shows the prescribed updraft profiles every 5 min.

Fig. 1.

Total surface precipitation rates vs time for idealized simulations with a 1D kinematic column model using the one-moment (heavy dashed), two-moment (solid), and three-moment (heavy dot–dashed) versions of the microphysics scheme. (The peak precipitation rate for the one-moment simulation is 340 mm h−1.) Inset shows the prescribed updraft profiles every 5 min.

The time series for the total surface precipitation rates for the three simulations during the first 50 min are shown in Fig. 1. Figure 2 depicts instantaneous profiles of the mass concentrations of each hydrometeor at 20 and 30 min for each run. The surface precipitation shown in all of these simulations comes from a combination of hail and rain (from melted hail) as can be seen in Figs. 2d,e,f (at times later than 50 min other particle types eventually precipitate out as well). The peak precipitation rate in SM is much larger than that of DIAG and TM, with values of 340 mm h−1 compared to 23 and 38 mm h−1, respectively. This is largely related to the fact that the bulk hail fall velocity, VQh, in SM is monotonically related to the hail mass content, Qh, and an unrealistically large accumulation zone results. At 20 min, the peak updraft velocity is approximately 16 m s−1 (Fig. 1 inset). With the constant supply of low-level water vapor, there is continuous condensation to cloud water that is then quickly accreted by hail. Figure 2a shows that the hail mass content is ∼5 to 7.5 g m−3 between 4 and 6 km in the vertical. The bulk fall velocity of hail turns out to be approximately 14–18 m s−1, which closely matches the updraft speed. As the updraft decreases in strength, the large quantity of hail sediments resulting very quickly in a sudden spike of very large precipitation rates at the surface. This is not the case in DIAG and TM, where NTh and VQh are not monotonically related to Qh and where the size-sorting mechanism allows hail with high number concentrations, and hence relatively lower bulk fall velocities, to be advected aloft and away from the region of high LWC (Figs. 2b,c). On the other hand, hail with lower number concentrations and higher bulk fall velocities is able to fall through the updraft. Overall, it can be observed that the general prediction of the hydrometeor mass distributions between the two-moment and three-moment versions of the scheme are quite similar, though differences in the timing and characteristics of the surface precipitation are apparent. In contrast, the results for the single-moment version of the scheme are considerably different.

Fig. 2.

Profiles of the mass contents of cloud (solid), rain (thick gray), ice (dashed), snow (thick dashed), graupel (dot–dashed), and hail (thick solid) for idealized 1D kinematic column model simulations using the (a), (d) one-moment, (b), (e) two-moment (with diagnosed αx), and (c), (f) three-moment versions of the proposed microphysics scheme at (top) 20 and (bottom) 30 min.

Fig. 2.

Profiles of the mass contents of cloud (solid), rain (thick gray), ice (dashed), snow (thick dashed), graupel (dot–dashed), and hail (thick solid) for idealized 1D kinematic column model simulations using the (a), (d) one-moment, (b), (e) two-moment (with diagnosed αx), and (c), (f) three-moment versions of the proposed microphysics scheme at (top) 20 and (bottom) 30 min.

6. Conclusions

A three-moment closure formulation for bulk microphysics parameterizations has been proposed. By introducing a tendency equation for a third moment, the radar reflectivity, as a function of the tendencies of the other two predicted moments for each microphysical process, the spectral shape parameter, α, of the gamma size distribution becomes an independent prognostic variable. In view of the advantages of the variable α approaches, discussed in Part I, a new multimoment bulk scheme has been developed. The two-moment version of the scheme predicts the zeroth and the third moments of the size spectrum of each hydrometeor category and uses the diagnostic equations for the shape parameter, based on the mean-particle size, introduced in Part I. The three-moment version also predicts α by an additional prognostic equation for radar reflectivity.

It was shown in Part I that the two-moment approach is superior to the one-moment approach, the diagnosed-α two-moment method is an improvement over the fixed-α method, and the fully prognosed-α (three-moment) technique is the best. The 1D kinematic model simulations presented in this paper show that while the two-moment (diagnosed α) version of the new scheme does a much better job at reproducing the results of the full three-moment version than the one-moment version, there are still notable differences in certain aspects such as the surface precipitation rate. Thus, there appears to be a gain in skill for the prognosed-α over the diagnosed-α approach.

Note, however, that additional computational costs are involved with both of these variable-α methods. For the diagnostic-α compared to the fixed-α two-moment approach, although no additional prognostic variable needs to be advected, many of the growth equations contain expressions involving the gamma function with αx as an argument (see section 4); these expressions must be computed at each time step and grid point if αx is a variable, whereas they may be precomputed and hardwired if αx is a constant. Thus the diagnostic-α two-moment approach involves additional computational cost compared to the fixed-α approach, even though all of the parametric equations are identical.

Since sedimentation and the source/sink terms feed back onto one another and can significantly affect the evolution of the microphysics and dynamics of a modeled storm system, detailed examination of the overall affects of multiple-moment bulk schemes should be carried out in the context of a full 3D dynamical model. Further examination of the validity of the proposed diagnostic relations for the shape parameter, and possible improvements to these equations, should also be done in that context since processes other than sedimentation will act to change the shape parameter in a three-moment scheme. The new microphysics scheme has been implemented into the Canadian Mesoscale Compressible Community (MC2) model (Benoit et al. 1997) and applied to successfully simulate a severe hailstorm using a horizontal grid spacing of 1 km. Future papers in this series (Milbrandt and Yau 2005b, c, manuscripts submitted to J. Atmos. Sci.) will present the results of the control simulation, using the full three-moment approach, along with sensitivity experiments using the various versions of the microphysics scheme.

Table 2.

Terminal fall velocity parameters for hydrometeor category x (see Ferrier 1994 for references).

Terminal fall velocity parameters for hydrometeor category x (see Ferrier 1994 for references).
Terminal fall velocity parameters for hydrometeor category x (see Ferrier 1994 for references).

List of Symbols

List of Symbols
List of Symbols

List of Symbols—Continued

List of Symbols—Continued
List of Symbols—Continued

List of Symbols—Continued

List of Symbols—Continued
List of Symbols—Continued

Acknowledgments

We thank the three anonymous reviewers for their suggestions for improving the original manuscript. In particular, we thank one reviewer for the helpful suggestions in improving the presentation of the scheme description. We would also like to express our appreciation to Dr. J.-P. Pinty for the valuable discussion and clarifications and also for providing us with the source code for their warm-rain scheme. This research was supported by the Canadian Foundation for Climate and Atmospheric Science.

REFERENCES

REFERENCES
Benoit
,
R.
,
J. M.
Desgagné
,
P.
Pellerin
,
S.
Pellerin
,
Y.
Chartier
, and
S.
Desjardins
,
1997
:
The Canadian MC2: A semi-Lagrangian, semi-implicit wideband atmospheric model suited for finescale process studies and simulation.
Mon. Wea. Rev.
,
125
,
2382
2415
.
Berry
,
E.
, and
R.
Reinhardt
,
1974
:
An analysis of cloud drop growth by collection. Part II: Single initial distributions.
J. Atmos. Sci.
,
31
,
1825
1831
.
Bigg
,
E. K.
,
1953
:
The supercooling of water.
Proc. Phys. Soc. London
,
B66
,
668
694
.
Byers
,
H. R.
,
1965
:
Elements of Cloud Physics.
The University of Chicago Press, 191 pp
.
Cohard
,
J-M.
, and
J-P.
Pinty
,
2000a
:
A comprehensive two-moment warm microphysical bulk scheme. I: Description and tests.
Quart. J. Roy. Meteor. Soc.
,
126
,
1815
1842
.
Cohard
,
J-M.
, and
J-P.
Pinty
,
2000b
:
A comprehensive two-moment warm microphysical bulk scheme. II: 2D experiments with a non-hydrostatic model.
Quart. J. Roy. Meteor. Soc.
,
126
,
1843
1859
.
Cohard
,
J-M.
,
J-P.
Pinty
, and
C.
Bedos
,
1998
:
Extending Twomey’s analytical estimate of nucleated cloud droplet concentrations from CCN spectra.
J. Atmos. Sci.
,
55
,
3348
3357
.
Cotton
,
W. R.
,
G. J.
Tripoli
,
R. M.
Rauber
, and
E. A.
Mulvihill
,
1986
:
Numerical simulation of the effects of varying ice crystal nucleation rates and aggregation processes on orographic snowfall.
J. Climate Appl. Meteor.
,
25
,
1658
1680
.
DeMott
,
P. J.
,
M. P.
Meyers
, and
W. R.
Cotton
,
1994
:
Parameterization and impact of ice initiation processes relevant to numerical model simulations of cirrus clouds.
J. Atmos. Sci.
,
51
,
77
90
.
Ferrier
,
B. S.
,
1994
:
A two-moment multiple-phase four-class bulk ice scheme. Part I: Description.
J. Atmos. Sci.
,
51
,
249
280
.
Ferrier
,
B. S.
,
W-K.
Tau
, and
J.
Simpson
,
1995
:
A two-moment multiple-phase four-class bulk ice scheme. Part II: Simulations of convective storms in different large-scale environments and comparisons with other bulk parameterizations.
J. Atmos. Sci.
,
52
,
1001
1033
.
Hallet
,
J.
, and
S. C.
Mossop
,
1974
:
Production of secondary ice particles during the riming process.
Nature
,
249
,
26
28
.
Harrington
,
J. Y.
,
M. P.
Meyers
,
R. L.
Walko
, and
W. R.
Cotton
,
1995
:
Parameterization of ice crystal conversion processes due to vapor deposition for mesoscale models using double-moment basis functions. Part I: Basic formulation and parcel model results.
J. Atmos. Sci.
,
52
,
4344
4366
.
Kong
,
F.
, and
M. K.
Yau
,
1997
:
An explicit approach to microphysics in MC2.
Atmos.
,
Ocean
,
33
,
257
291
.
Lesins
,
G.
,
R.
List
, and
P.
Joe
,
1980
:
Ice accretions. Part I: Testing of new atmospheric icing concepts.
J. Rech. Atmos.
,
14
,
347
356
.
Lin
,
Y-L.
,
R. D.
Farley
, and
H. D.
Orville
,
1983
:
Bulk parameterization of the snow field in a cloud model.
J. Climate Appl. Meteor.
,
22
,
1065
1092
.
Long
,
A. B.
,
1974
:
Solutions to the droplet collection equation for polynomial kernels.
J. Atmos. Sci.
,
31
,
1040
1052
.
Macklin
,
W. C.
, and
I. H.
Bailey
,
1966
:
On the critical liquid water concentrations of large hailstones.
Quart. J. Roy. Meteor. Soc.
,
92
,
297
300
.
McCumber
,
M.
,
W-K.
Tao
, and
J.
Simpson
,
1991
:
Comparison of ice-phase microphysical parameterization schemes using numerical simulations of tropical convection.
J. Appl. Meteor.
,
30
,
985
1004
.
Meyers
,
M. P.
,
P. J.
DeMott
, and
W. R.
Cotton
,
1992
:
New primary ice-nucleation parameterizations in an explicit cloud model.
J. Climate Appl. Meteor.
,
31
,
708
721
.
Meyers
,
M. P.
,
R. L.
Walko
,
J. Y.
Harrington
, and
W. R.
Cotton
,
1997
:
New RAMS cloud microphysics. Part II: The two-moment scheme.
Atmos. Res.
,
45
,
3
39
.
Milbrandt
,
J. A.
, and
M. K.
Yau
,
2005a
:
A multimoment bulk microphysics parameterization. Part I: Analysis of the role of the spectral shape parameter.
J. Atmos. Sci.
,
62
,
3051
3064
.
Murakami
,
M.
,
1990
:
Numerical modeling of dynamical and microphysical evolution of an isolated convective cloud—The 19 July 1981 CCOPE cloud.
J. Meteor. Soc. Japan
,
68
,
107
128
.
Musil
,
D. J.
,
1970
:
Computer modeling of hailstone growth in feeder clouds.
J. Atmos. Sci.
,
27
,
474
482
.
Rasmussen
,
R. M.
, and
A. J.
Heymsfield
,
1987
:
Melting and shedding of graupel and hail. Part II: Sensitivity study.
J. Atmos. Sci.
,
44
,
2764
2782
.
Reisner
,
J.
,
R. M.
Rasmussen
, and
T.
Bruintjes
,
1998
:
Explicit forecasting of supercooled liquid water in winter storms using the MM5 mesoscale model.
Quart. J. Roy. Meteor. Soc.
,
124
,
1071
1107
.
Walko
,
R. L.
,
W. R.
Cotton
,
M. P.
Meyers
, and
J. Y.
Harrington
,
1995
:
New RAMS cloud microphysics. Part I: The one-moment scheme.
Atmos. Res.
,
38
,
29
62
.
Wisner
,
C.
,
R. D.
Orville
, and
C.
Myers
,
1972
:
A numerical model of a hail-bearing cloud.
J. Atmos. Sci.
,
29
,
1160
1181
.
Young
,
K. C.
,
1974
:
The role of contact nucleation in ice phase initiation in clouds.
J. Atmos. Sci.
,
31
,
1735
1748
.
Young
,
K. C.
,
1993
:
Microphysical Processes in Clouds.
Oxford University Press, 427 pp
.
Ziegler
,
C. L.
,
1985
:
Retrieval of thermal and microphysical variables in observed convective storms. Part 1: Model development and preliminary testing.
J. Atmos. Sci.
,
42
,
1497
1509
.

APPENDIX A

Source and Sink Terms

The microphysical source/sink terms of the continuity equations of the proposed scheme are listed below. The following equations are for the full three-moment version of the scheme. For the two-moment version, (A14)–(A18) are not used, nor are any equations for tendencies of Zx. For the one-moment version, (A8)(A18) are not used, nor are any equations for the tendencies of NTx and Zx.

The tendencies for the mass mixing ratios are

 
formula
 
formula
 
formula
 
formula
 
formula
 
formula
 
formula

The tendencies for the total number concentrations are

 
formula
 
formula
 
formula
 
formula
 
formula
 
formula

The tendencies for the radar reflectivity factors are

 
formula
 
formula
 
formula
 
formula
 
formula

The temperature change equation is

 
formula

APPENDIX B

General Forms of the Tendency Equations for NTx and Zx

Many of the tendency equations for NTx can be described by one of the two following general equations, using the notation for the source/sink terms described in section 4 for hydrometeor category x interacting with category y through process AB. The first type of equation,

 
formula

is based on the assumption that the mean particle mass does not change due to process AB. The second type of equation,

 
formula

assumes that particles of mass mx0 are being initiated into category x. In some cases, the initial particle mass is specified with a prescribed initial size Dx0, where mx0 = cxDdxx0 in (B2).

The three types of tendency equations for Zx are given by (2), (4), and (5). Using the notation in section 4, the equations can be expressed as

 
formula
 
formula

respectively, and

 
formula

Note that the negative sign that appears in (5) has been dropped in (B5); the signs are applied appropriately in (A14)(A18).

APPENDIX C

List of Symbols

Footnotes

Corresponding author address: Dr. Jason A. Milbrandt, Meteorological Research Branch, Environment Canada, 2121 Trans-Canada Highway, 5th Floor, Dorval, QC H9P 1J3, Canada. Email: jason.milbrandt@mcgill.ca