## Abstract

Many two-moment bulk schemes use a three-parameter gamma distribution of the form *N*(*D*) = *N*_{0}*D ^{α}e*

^{−}

*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*

^{λD}*N*

_{0}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*) = *N*_{0}*D ^{α}e*

^{−λD}. The schemes predict changes to the mass content and total number concentration thereby allowing

*N*

_{0}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

*N*

_{0x}and the slope parameter

*λ*, of a size distribution for hydrometeor species

_{x}*x*to be independent. To predict the evolution of

*α*, a prognostic equation for another moment such as the radar reflectivity

_{x}*Z*with the form

_{x}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 *Z _{x}* in (1) is computed as the sum of the individual tendencies of

*Z*for each microphysical process listed in appendix A. What remains, therefore, is to derive the source terms for

_{x}*Z*for each process

_{x}*A*in each hydrometeor category

*x*.

We classify the source terms for *Z _{x}* into three types. For the first type, it is assumed that the change in

*α*due to the particular process

_{x}*A*is negligible. By differentiating the closure Eq. (6) in Part I, we can relate the tendency of

*Z*to the tendencies of

_{x}*q*and

_{x}*N*as

_{Tx}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

*λ*constant for a process in order to relate the fractional change in

_{x}*N*to the fractional change in

_{Tx}*q*[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

_{x}*α*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 *Z _{x}*

^{init}=

*q*

_{x}^{init}=

*N*

_{Tx}^{init}= 0 at the beginning of the time step, and

*q*

_{x}^{fin}and

*N*

_{Tx}^{fin}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

*Z*, the value of

_{x}*α*at the end of the time step for process

_{x}*A*(

*α*) 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

_{xA}By writing (*dZ _{x}*/

*dt*)|

_{A}= (

*Z*

^{fin}

_{x}−

*Z*

^{init}

_{x}/Δ

*t*), substituting (3) for

*Z*

_{x}^{fin}, and taking the limit as Δ

*t*→ 0, we obtain

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

## 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

*m*(

_{x}*D*) =

_{x}*c*(with

_{x}D_{x}^{dx}*d*= 3). Except for ice, all particles are assumed to be spherical with

_{x}*c*= (

_{x}*π*/6),

*ρ*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

_{x}*c*= 440 kg m

_{i}^{−3}(F94).

There are three main versions of the proposed scheme. The single-moment version predicts *q _{x}* using (I.8) allowing

*λ*defined by (I.5) to vary prognostically, and has fixed values for

_{x}*N*and

_{Tc}*N*

_{0x}for the precipitating categories and

*α*for all categories. The two-moment version includes predictive equations for

_{x}*q*and also for

_{x}*N*, given by (I.9), where the corresponding prognostic parameters are

_{Tx}*N*

_{0x}, defined by (I.4), and

*λ*. 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

_{x}*Z*for all categories except cloud water, and

_{x}*α*is obtained by solving (I.6).

_{x}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 *VAB _{y,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,

*V*AB

_{x}*is the rate of change of variable*

_{y,x}*V*in category

*x*due to AB

*while*

_{y,x}*V*AB

_{y}*is the rate of change of variable*

_{y,x}*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, *N*CL* _{y,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 (N*NU_{υc})

_{υc})

CPB98 developed an expression relating the concentration of nucleated cloud droplets *N*_{CCN} 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 *N _{Tc}* at a certain time, the cloud nucleation rate is

### b. *Condensation and evaporation of cloud and rain (Q*VD_{υc}, QVD_{υr}, NVD_{υc}, NVD_{υr}, ZVD_{υr})

_{υc}, Q

_{υr}, N

_{υc}, NVD

_{υr}, Z

_{υr})

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

where

is the bulk ventilation coefficient for rain and

is the thermodynamic function (Byers 1965).

During evaporation, *N _{Tc}* and

*N*decrease at the rates of

_{Tr}*N*VD

*, which is given by an equation of the form of (B1) in appendix B for*

_{υx}*x*∈ [

*c*,

*r*] and

*Z*

_{r}decreases at the rate of

*Z*VD

*, given by an equation of the form of (B3).*

_{υr}### c. *Warm-rain collection (Q*CN_{cr}, QCL_{cr}, N_{c}CN_{cr}, N_{r}CN_{cr}, NCL_{rr}, ZCN_{cr}, ZCL_{cr}, ZCL_{rr})

_{cr}, Q

_{cr}, N

_{c}

_{cr}, N

_{r}

_{cr}, N

_{rr}, Z

_{cr}, Z

_{cr}, Z

_{rr})

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 *q _{x}* tendency for autoconversion,

*Q*CN

*, is adopted from CP00a. To compute the change in concentration of rain particles due to autoconversion,*

_{cr}*N*CN

_{r}*, we note that in general, for a microphysical process generating new particles of mean diameter*

_{cr}*D*

_{mr_new}, its

*N*tendency is related to its

_{Tx}*q*tendency by (B2) where

_{x}*D*

_{x0}=

*D*

^{dx}_{mx_new}. To determine

*D*

_{mr_new}of rain due to autoconversion, a threshold mean diameter

*D*

_{mr}_{_aut}is first computed using the empirical formula of Berry and Reinhardt (1974),

where *σ _{c}* is the standard deviation of the size of the cloud droplets (see CP00a). If the existing mean-drop diameter

*D*= [(

_{mr}*ρq*/

_{r}*c*)]

_{r}N_{Tr}^{1/3}exceeds

*D*

_{mr_aut}, then

*D*

_{mr_new}=

*D*; otherwise

_{mr}*D*

_{mr_new}=

*D*

_{mr_aut}. Thus,

*N*CN

_{r}*is given by (B2) with*

_{cr}*D*

_{x0}=

*D*

^{3}

_{mr_new}and the equation for

*Z*CN

*is of the form of (B4) with*

_{cr}*α*

_{xAB}=

*α*

_{rAUT}.

For the *q _{x}* and

*N*tendencies due to accretion (

_{Tx}*Q*CL

*and*

_{cr}*N*CL

*) and self-collection (*

_{cr}*N*CN

_{c}*and*

_{cr}*N*CL

*), we exactly follow CP00a. Note that the*

_{rr}*N*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

_{Tc}*N*CL

*(see Cohard and Pinty 2000b). The*

_{rr}*Z*tendencies for accretion

_{r}*(Z*CL

*) and self-collection (*

_{cr}*Z*CL

*) are identical to the form of (B3).*

_{rr}### 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:*

Bulk collection efficiencies, *E _{xy}*, are used. For the collection of cloud droplets by graupel and hail, an approximate empirical formula from Macklin and Bailey (1966) is used,

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

However, during wet growth of hail, we set *E _{ih}* =

*E*= 1.

_{sh}The values for the other bulk collection efficiencies are *E _{ci}* =

*E*= 1,

_{ri}*E*= 0.1

_{ii}*E*,

_{is}*E*= 1,

_{cs}*E*= 1,

_{rs}*E*= 0.1

_{ss}*E*,

_{is}*E*= 1,

_{rg}*E*= 0,

_{gg}*E*= 1,

_{rh}*E*= 0, and

_{gh}*E*= 0.

_{hh}##### (i) *Collection amongst rain and frozen categories (Q*CL_{ri}, QCL_{ir}, QCL_{rs}, QCL_{sr}, QCL_{rg}, QCL_{gr}, QCL_{rh}, QCL_{is}, QCL_{ig}, QCL_{ih}, QCL_{sh}, NCL_{ri}, NCL_{ir}, NCL_{rs}, NCL_{sr}, NCL_{rg}, NCL_{gr}, NCL_{rh}, NCL_{is}, NCL_{ig}, NCL_{ih}, NCL_{sh}, ZCL_{ri}, ZCL_{ir}, ZCL_{rs}, ZCL_{sr}, Z_{i}CL_{is}, Z_{s}CL_{is}, Z_{i}CL_{ig}, Z_{g}CL_{ig}, Z_{i}CL_{ih}, Z_{h}CL_{ih}, Z_{s}CL_{sh}, Z_{h}CL_{sh}, ZCL_{rg}, ZCL_{gr}, Z_{r}CL_{rh}, Z_{h}CL_{rh})

_{ri}, Q

_{ir}, Q

_{rs}, Q

_{sr}, Q

_{rg}, Q

_{gr}, Q

_{rh}, Q

_{is}, Q

_{ig}, Q

_{ih}, Q

_{sh}, N

_{ri}, N

_{ir}, N

_{rs}, N

_{sr}, N

_{rg}, N

_{gr}, N

_{rh}, N

_{is}, N

_{ig}, N

_{ih}, N

_{sh}, Z

_{ri}, Z

_{ir}, Z

_{rs}, Z

_{sr}, Z

_{i}

_{is}, Z

_{s}

_{is}, Z

_{i}

_{ig}, Z

_{g}

_{ig}, Z

_{i}

_{ih}, Z

_{h}

_{ih}, Z

_{s}

_{sh}, Z

_{h}

_{sh}, Z

_{rg}, Z

_{gr}, Z

_{r}

_{rh}, Z

_{h}

_{rh})

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

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

and

and where *Z*CL* _{yx}* is of the form of (B3).

Equations (17), (18), and the equation for *Z*CL* _{yx}* apply to the collection of rain by all frozen categories (CL

*), collection of ice, snow, or graupel by rain (CL*

_{rx}*) and collection of one frozen category by another frozen category (CL*

_{xr}*), except for hail when it is undergoing wet growth (see section 4c). Equations of the form of (B3) are also used to compute*

_{xy}*Z*CL

*and*

_{cr}*Z*CL

*(with*

_{rr}*Q*CL

*= 0).*

_{rr}##### (ii) *Collection of cloud water by frozen categories (Q*CL_{ci}, QCL_{cs}, QCL_{cg}, QCL_{ch}, NCL_{ci}, NCL_{cs}, NCL_{cg}, NCL_{ch}, ZCL_{ci}, ZCL_{cs}, ZCL_{cg}, ZCL_{ch})

_{ci}, Q

_{cs}, Q

_{cg}, Q

_{ch}, N

_{ci}, N

_{cs}, N

_{cg}, N

_{ch}, Z

_{ci}, Z

_{cs}, Z

_{cg}, Z

_{ch})

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,

and *Z*CL* _{cx}* is calculated using (B3) with

*N*CL

*= 0,*

_{yx}*x*∈ [

*i*,

*s*,

*g*,

*h*], and

*y*=

*c*.

#### 2) Aggregation for snow *(N*CL_{ss}, ZCL_{ss})

_{ss}, Z

_{ss})

#### 3) Wet growth of hail (*N*SH_{hr}, ZSH_{hr}, QCL_{rh})

_{hr}, Z

_{hr}, Q

_{rh}

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 *Q*CL* _{ch}* +

*Q*CL

*+*

_{rh}*Q*CL

*+*

_{ih}*Q*CL

*, 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,*

_{sg}*Q*WET, 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

_{h}*Q*CL

*,*

_{ch}*Q*CL

*, and*

_{ih}*Q*CL

*are recomputed using collection efficiencies of 1 to obtain*

_{sh}*Q*CL′

*,*

_{ch}*Q*CL′

*, and*

_{ih}*Q*CL′

*. Likewise, the*

_{sh}*N*tendencies

_{Tx}*N*CL

*,*

_{ch}*N*CL

*, and*

_{ih}*N*CL

*are also recalculated using collection efficiencies of 1 and*

_{sh}*Z*CL

*and*

_{ih}*Z*CL

*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*

_{sh}*Q*WET and the sum of the new (wet) growth rates for collecting cloud, ice and snow

_{h}Also the difference between the dry collection rate *Q*CL* _{rh}* and the wet collection rate

*Q*CL′

*is assumed to be the rate that the collected water mass by the hailstone is shed to form rain (see LFO for discussion).*

_{rh}The change in *N _{Tr}* 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

*D*

_{mr_new}=

*D*

_{shed}= 1 mm, we obtain

The *Z _{r}* tendency due to shedding is given by

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 (*Q*NU_{υi}, NNU_{υi}, ZNU_{υi})

_{υi}, N

_{υi}, Z

_{υi})

Primary nucleation includes contact nucleation, deposition nucleation, and condensation-freezing nucleation. The *N _{Ti}* 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

*N*as a function of the saturation ratio with respect to ice,

_{Ti}*S*:

_{i}The *N _{Ti}* tendency due to deposition/condensation-freezing nucleation is thus

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, *N*NU* _{υ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),

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

The mixing ratio and reflectivity factor tendencies are calculated using an assumed nucleated ice crystal mass *m*_{i0} = 10^{−12} kg (KY97; RRB), giving

and, where *Z*NU* _{υi}* is of the form of (B4), with

*α*

_{xAB}=

*α*

_{iNU}= 0.

#### 2) Ice multiplication (*Q*IM_{si}, QIM_{gi}, NIM_{ii}, NIM_{si}, NIM_{gi}, ZIM_{ii}, ZIM_{si}, ZIM_{gi})

_{si}, Q

_{gi}, N

_{ii}, N

_{si}, N

_{gi}, Z

_{ii}, Z

_{si}, Z

_{gi})

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 *q _{x}* and

*N*tendencies and with the

_{Tx}*Z*tendencies given by equations following the form of (B4), with

_{x}*α*

_{xAB}=

*α*and

_{iIM}*y*∈ [

*i*,

*s*,

*g*]. Note that for rime splintering due to riming of ice, the value of

*Q*IM

*used to compute*

_{ii}*Z*IM

*is formulated the same as for*

_{ii}*Q*IM

*and*

_{si}*Q*IM

*(see RRB) but there is no change in the ice mixing ratio [i.e.,*

_{gi}*Q*IM

*= 0 and does not appear in (A4)].*

_{ii}#### 3) Homogenous freezing of cloud droplets (*Q*FZ_{ci}, NFZ_{ci}, ZFZ_{ci})

_{ci}, N

_{ci}, Z

_{ci})

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 Δ*N*_{freeze} that freeze in time Δ*t* at a given temperature as

where

and *V* is the droplet volume. By substituting *V* with *V* = (*π*/6)*D _{mc}*

^{3}, the mean-droplet volume in (29), we obtain

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

and the equation for *Z*FZ* _{ci}* follows the form of (B3). For

*T*> −30°C,

_{c}*f*= 0, while for

_{fr}*T*< −50°C,

_{c}*f*= 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).

_{fr}### f. *Deposition/sublimation (Q*VD_{υi}, QVD_{υs}, QVD_{υg}, QVD_{υh}, NVD_{υi}, NVD_{υs}, NVD_{υg}, NVD_{υh}, ZVD_{υi}, ZVD_{υs}, ZVD_{υg}, ZVD_{υh})

_{υi}, Q

_{υs}, Q

_{υg}, Q

_{υh}, N

_{υi}, N

_{υs}, N

_{υg}, N

_{υh}, Z

_{υi}, Z

_{υs}, Z

_{υg}, Z

_{υ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

where

is the thermodynamic function and

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 VD_{max} is computed similar to that for condensation/evaporation (see KY97).

For particles undergoing sublimation, the number concentrations are reduced at the rate *N*VD* _{υ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

*Z*VD

*, which is of the form of (B3).*

_{υr}### g. Freezing of rain

#### 1) Probabilistic freezing (*Q*FZ_{rh}, NFZ_{rh}, Z_{r}FZ_{rh}, Z_{h}FZ_{rh})

_{rh}, N

_{rh}, Z

_{r}

_{rh}, Z

_{h}

_{rh})

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

and the rate of change in mass by

#### 2) Collisional (three-component) freezing (*N*CL_{irg}, NCL_{irh}, NCL_{srs}, NCL_{srg}, NCL_{srh}, ZCL_{irg}, ZCL_{irh}, ZCL_{srs}, ZCL_{srg}, ZCL_{srh})

_{irg}, N

_{irh}, N

_{srs}, N

_{srg}, N

_{srh}, Z

_{irg}, Z

_{irh}, Z

_{srs}, Z

_{srg}, Z

_{srh})

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

where *x* ∈ [*i*, *s*, *g*], the destination category *z* ∈ [*s*, *g*, *h*], and *D _{mz}* = max(

*D*,

_{mx}*D*). 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

_{mr}*ρ*. It is classified as snow if

_{z}*ρ*≤ 0.5(

_{z}*ρ*+

_{s}*ρ*), as graupel if 0.5(

_{g}*ρ*+

_{s}*ρ*) <

_{g}*ρ*≤ 0.5(

_{z}*ρ*+

_{g}*ρ*), and as hail if

_{h}*ρ*> (

_{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 *Q*CL* _{xr}*,

*N*CL

*, and*

_{xr}*Z*CL

*. The terms for category*

_{xr}*x*collecting rain are

*Q*CL

*,*

_{rx}*N*CL

*, and*

_{rx}*Z*CL

*. The source term for the mixing ratio of the destination category*

_{rx}*z*is the sum of the mixing ratio sink terms for the two colliding categories. Thus the delta function variable

*δ*= 1 for the destination category

_{xrz}*z,*but is 0 for the other frozen categories. For example, if rain and snow collide to form graupel, then

*δ*= 1 while

_{srg}*δ*=

_{srs}*δ*= 0 so the rate of change of the graupel mixing ratio is (

_{srh}*Q*CL

*+*

_{rg}*Q*CL

*). The*

_{gr}*N*tendency for the destination category is calculated using the mean-mass diameter

_{Tz}*D*

_{m}_{z}. Thus,

where *ρ _{z}* is the actual density of category

*z*, not the density calculated from (39) used to determine the destination category. The

*Z*tendency for the resulting category is given by

_{z}where *α*_{zCL} = 0 is the assumed *α _{z}* for the newly formed particles.

### h. Conversion processes

#### 1) Ice to snow (*Q*CN_{is}, N_{i}CN_{is}, N_{s}CN_{is}, Z_{i}CN_{is}, Z_{s}CN_{is})

_{is}, N

_{i}

_{is}, N

_{s}

_{is}, Z

_{i}

_{is}, Z

_{s}

_{is}

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 *m*_{s0} = 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),

The change in *N _{Ts}* is calculated with the assumption that all snow particles converted from ice have an initial mass

*m*

_{s0}, giving the equation for

*N*CN

_{s}*of the form of (B2). The change in*

_{is}*N*due to conversion to snow by deposition and riming is computed as the negative rate of change in

_{Ti}*N*due to conversion from ice deposition and riming. For conversion by ice aggregation, the

_{Ts}*N*tendency is given by M90,

_{Ti}where

and *D _{mi}* is the mean ice crystal diameter,

*V*is the fall velocity of ice crystals,

_{i}*E*is the (dry, temperature dependent) collection efficiency amongst ice crystals and

_{ii}*X*

_{disp}is the dispersion of the fall velocity spectrum of the ice crystals (assumed for simplicity to be 0.25, following M90). The total

*N*tendency due to conversion to snow is therefore

_{Ti}#### 2) Ice to graupel (*Q*CN_{ig}, *N*_{i}CN_{ig}, *N*_{g}CN_{ig}, *Z*_{i}CN_{ig}, I_{g}CN_{ig})

_{ig}

_{i}

_{ig}

_{g}

_{ig}

_{i}

_{ig}

_{g}

_{ig}

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

The *N _{Tg}* tendency is calculated based on an assumed embryo graupel mass

*m*

_{g0}= 1.6 × 10

^{−10}kg, giving the equation for

*N*CN

_{g}*of the form of (B2), while*

_{ig}*N*is reduced following M90 by

_{Ti}*N*CN

_{i}*, which is of the form of (B1). The reflectivity tendencies are given by*

_{ig}*Z*CN

_{i}*, which is of the form of (B3) with*

_{ig}*x*=

*i*, and

*Z*CN

_{g}*, which is of the form of (B4) with*

_{ig}*α*

_{xAB}=

*α*

_{gCN}= 0 and with

*x*=

*g*.

#### 3) Snow to graupel (*Q*CN_{sg}, NCN_{sg}, Z_{s}CN_{sg}, Z_{g}CN_{sg})

_{sg}, N

_{sg}, Z

_{s}

_{sg}, Z

_{g}

_{sg}

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 *Q*CL* _{cs}* exceeds the depositional growth rate

*Q*VD

*. It is assumed that the change in snow particle size due to riming is negligibly small. The conversion rate is*

_{υs}The concentration tendency *N*CN* _{sg}* is given by (B2) assuming that newly converted graupel particles are embryos with mass

*m*

_{g0}(RRB). This is a source for

*N*and a sink for

_{Tg}*N*. The reflectivity tendencies are given by

_{Ts}*Z*CN

_{s}*, which is of the form of (B3) with*

_{sg}*x*=

*s*, and

*Z*CN

_{g}*, which is of the form of (B5) with*

_{sg}*x*=

*g*.

#### 4) Graupel to hail (*Q*CN_{gh}, NCN_{gh}, Z_{g}CN_{gh}, Z_{h}CN_{gh})

_{gh}, N

_{gh}, Z

_{g}

_{gh}, Z

_{h}

_{gh}

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)

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

However, (49) needs to be modified as *D*_{h0} 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:

When the ambient temperature is relatively cold and/or the LWC is low, *D*_{h0} 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 *D _{mg}*/

*D*

_{h0}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 *D*_{x0} = *D*^{3}_{h0}. The *N _{Tg}* tendency is set to the negative of the

*N*tendency. The

_{Th}*Z*tendencies are given by

_{x}*Z*CN

_{g}*, which is of the form of (B3) with*

_{gh}*x*=

*g*, and

*Z*, which is of the form of (B5) with

_{h}CN_{gh}*x*=

*h*.

##### (i) *Melting of frozen particles (Q*ML_{ir}, QML_{sr}, QML_{gr}, QML_{hr}, NML_{ir}, NML_{sr}, NML_{gr}, NML_{hr}, Z_{i}ML_{ir}, Z_{r}ML_{ir}, Z_{s}ML_{sr}, Z_{r}ML_{sr}, Z_{g}ML_{gr}, Z_{r}ML_{gr}, Z_{h}ML_{hr}, Z_{r}ML_{hr})

_{ir}, Q

_{sr}, Q

_{gr}, Q

_{hr}, N

_{ir}, N

_{sr}, N

_{gr}, N

_{hr}, Z

_{i}

_{ir}, Z

_{r}

_{ir}, Z

_{s}

_{sr}, Z

_{r}

_{sr}, Z

_{g}

_{gr}, Z

_{r}

_{gr}, Z

_{h}

_{hr}, Z

_{r}

_{hr})

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

and

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 *q _{x}* and

*N*tendencies are

_{Tx}## 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.

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, *V _{Qh}*, in SM is monotonically related to the hail mass content,

*Q*, and an unrealistically large accumulation zone results. At 20 min, the peak updraft velocity is approximately 16 m s

_{h}^{−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

*N*and

_{Th}*V*are not monotonically related to

_{Qh}*Q*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.

_{h}## 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

*α*is a variable, whereas they may be precomputed and hardwired if

_{x}*α*

_{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.

## 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

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### 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 *Z _{x}*. For the one-moment version, (A8)–(A18) are not used, nor are any equations for the tendencies of

*N*and

_{Tx}*Z*.

_{x}The tendencies for the mass mixing ratios are

The tendencies for the total number concentrations are

The tendencies for the radar reflectivity factors are

The temperature change equation is

### APPENDIX B

#### General Forms of the Tendency Equations for NTx and Zx

Many of the tendency equations for *N _{Tx}* 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,

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

assumes that particles of mass *m*_{x0} are being initiated into category *x*. In some cases, the initial particle mass is specified with a prescribed initial size *D*_{x0}, where *m*_{x0} = *c _{x}D^{dx}*

_{x0}in (B2).

The three types of tendency equations for *Z _{x}* are given by (2), (4), and (5). Using the notation in section 4, the equations can be expressed as

respectively, and

### APPENDIX C

#### List of Symbols

^{0101}^{0102}^{0103}

## 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