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.
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
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
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
By writing (dZx/dt)|A = (Zfinx − Zinitx/Δt), substituting (3) for Zxfin, 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 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).
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
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
is the bulk ventilation coefficient for rain and
is the thermodynamic function (Byers 1965).
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),
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
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,
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 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:
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,
and ZCLcx is calculated using (B3) with NCLyx = 0, x ∈ [i, s, g, h], and y = c.
2) Aggregation for snow (NCLss, ZCLss)
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
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
The Zr 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 (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:
The NTi 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, 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),
Also, primary nucleation can result in an increase in the number of ice crystals present, but will never result in a decrease. Thus,
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
and V is the droplet volume. By substituting V with V = (π/6)Dmc3, 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 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
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 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
and the rate of change in mass by
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
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,
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
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),
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,
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
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
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
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)
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
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:
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,
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
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, 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.
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.
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.
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
The tendencies for the total number concentrations are
The tendencies for the radar reflectivity factors are
The temperature change equation is
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,
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 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).
List of Symbols
Corresponding author address: Dr. Jason A. Milbrandt, Meteorological Research Branch, Environment Canada, 2121 Trans-Canada Highway, 5th Floor, Dorval, QC H9P 1J3, Canada. Email: email@example.com