Bulk microphysics parameterizations that are used to represent clouds and precipitation usually allow only solid and liquid hydrometeors. Predicting the bulk liquid fraction on ice allows an explicit representation of mixed-phase particles and various precipitation types, such as wet snow and ice pellets. In this paper, an approach for the representation of the bulk liquid fraction into the predicted particle properties (P3) microphysics scheme is proposed and described. Solid-phase microphysical processes, such as melting and sublimation, have been modified to account for the liquid component. New processes, such as refreezing and condensation of the liquid portion of mixed-phase particles, have been added to the parameterization. Idealized simulations using a one-dimensional framework illustrate the overall behavior of the modified scheme. The proposed approach compares well to a Lagrangian benchmark model. Temperatures required for populations of ice crystals to melt completely also agree well with previous studies. The new processes of refreezing and condensation impact both the surface precipitation type and feedback between the temperature and the phase changes. Overall, prediction of the bulk liquid fraction allows an explicit description of new precipitation types, such as wet snow and ice pellets, and improves the representation of hydrometeor properties when the temperature is near 0°C.
The passage of a warm front can produce favorable environmental conditions for many winter precipitation types such as rain, freezing rain, ice pellets, wet snow, and snow (Stewart 1985). In such situations, a temperature (T) inversion, characterized by a warm layer aloft with T > 0°C (called the melting layer) and cold layer near the surface with T < 0°C (called the refreezing layer), occurs (e.g., Carmichael et al. 2011; Thériault et al. 2006; Gyakum and Roebber 2001; Hanesiak and Stewart 1995; Lin and Stewart 1986). The precipitation types formed when the temperature is near 0°C involve several microphysical processes including melting, refreezing, vapor deposition, collection, and wet growth. For example, when T > ~0°C, wet snow is formed from the partial melting of snow during which the melted water tends to accumulate on the ice core to form mixed-phase particles (Fujiyoshi 1986). Ice pellets are formed by refreezing in cold layers of the partially melted particles from aloft (Carmichael et al. 2011; Thériault and Stewart 2010; Gibson and Stewart 2007). Freezing rain corresponds to the surface icing of supercooled raindrops present in cold layers originating from the total melting of snow in warm layers aloft. Often, several types of precipitation can coexist, and, for example, the amount of surface freezing rain can be greatly inhibited because of the collection in the cold layers between supercooled raindrops and solid hydrometeors, such as ice pellets (Barszcz 2017; Carmichael et al. 2011; Hogan 1985). When T < 0°C, wet growth of graupel/hail and the shedding of accumulated liquid water can occur in conditions of high liquid water content, also involving the formation of mixed-phase particles (e.g., Phillips et al. 2014). In turn, the phase changes occurring with some of these processes impact the environmental temperature as, for example, the cooling induced by melting and evaporation (e.g., Lackmann et al. 2002) or the warming induced by refreezing and condensation (e.g., Thériault et al. 2006).
Most bulk microphysics schemes allow only for the representation of solid-phase hydrometeors (e.g., snow, graupel, hail) and liquid-phase hydrometeors (e.g., rain, cloud) (e.g., Morrison et al. 2009; Thompson et al. 2008; Milbrandt and Yau 2005). Mixed-phase hydrometeors, such as wet snow, are typically not represented. For example, instead of forming wet snow from partial melting, most schemes immediately transfer the melted water directly into the rain category. To represent the evolution of mixed-phase particles in microphysical parameterizations explicitly, it is necessary to track in space and time the liquid fraction of mixed-phase hydrometeors, defined as the ratio of liquid water mass to the total particle mass.
Some studies have experimented with an explicit parameterization of the liquid fraction. Based on the experimental and theoretical relationships developed by Mitra et al. (1990, hereafter M90), the study of Szyrmer and Zawadzki (1999) predicted the liquid fraction in a bulk microphysics module. This scheme was coupled with a two-dimensional nonhydrostatic, fully compressible dynamic framework to evaluate the effects of partial melting on the simulated “brightband” parameters (Fabry and Szyrmer 1999; Braun and Houze 1995; Fabry and Zawadzki 1995; Yokoyama and Tanaka 1984; Austin and Bemis 1950). They showed that the description of the particle habit prior to melting (unrimed or rimed and spherical or nonspherical) influences the melting rate and the liquid fraction evolution. They also pointed out the relevance of predicting the liquid fraction for radar applications. Based on the Szyrmer and Zawadzki (1999) parameterization, Thériault and Stewart (2010) modified the Milbrandt and Yau (2005) scheme to include several winter precipitation type categories. The new categories were wet snow, almost completely melted snow, refrozen wet snow, and two types of ice pellets. They tested the scheme using a one-dimensional (1D) cloud model and showed good agreement with observations for several winter storms (Thériault and Stewart 2010; Thériault et al. 2010, 2006). More recently, and also based on Szyrmer and Zawadzki (1999), Frick et al. (2013) implemented partial melting in a bulk microphysics scheme coupled with the Consortium for Small-Scale Modeling (COSMO; Doms and Schättler 2002) mesoscale model. They clearly showed the potential of their parameterization to simulate wet snow; however, their parameterization did not include the refreezing process, which is necessary to produce ice pellets. The implementation of an explicit liquid fraction has also been tested using bin microphysics parameterizations (Reeves et al. 2016; Geresdi et al. 2014; Phillips et al. 2007). Walko et al. (2000) developed an algorithm to diagnose the liquid fraction in the Regional Atmospheric Modeling System, which can be used with larger numerical time steps and regional modeling.
Recently, the predicted particle properties (P3) bulk microphysics scheme described in Morrison and Milbrandt (2015, hereafter MM15) and Milbrandt and Morrison (2016, hereafter MM16) was introduced. The P3 scheme completely abandons the use of predefined ice-phase categories and introduces the idea of “free” ice categories. With four prognostic variables per free category, important physical properties of ice can evolve realistically and smoothly in time and space, and thus, a wide range of ice particle types can be represented. Despite the conceptual improvement over traditional fixed-category schemes for representing ice, P3 is still limited in that it cannot represent ice particles that accumulate liquid water. To address this limitation, in this study, the original P3 scheme is modified to include the prediction of the bulk liquid fraction of ice categories and thus allow the representation of mixed-phase particles.
This paper is organized as follows. Section 2 briefly describes the original P3 scheme. Section 3 describes the implementation of the bulk liquid fraction into P3. Section 4 shows a validation of the proposed melting parameterization compared to a Lagrangian benchmark model. Section 5 shows results from two idealized experiments using a one-dimensional modeling framework to illustrate the overall behavior of the new parameterization and compare with the original one. Section 6 gives concluding remarks.
2. Overview of the original P3 microphysics scheme
This section gives a brief overview of the original P3 microphysics scheme (referred to as P3_ORIG); further details are in MM15 and MM16. As in other bulk schemes, P3 has two liquid-phase categories, cloud and rain, both of which are two moment, with prognostic mass and number mixing ratios for each. There is a user-specified number of free ice-phase categories, each with four prognostic variables. In this study, only the single-ice category configuration is discussed, but the modifications described below are general and can apply to multi-ice category configurations. The four prognostic ice variables are total ice mass qi,tot (kg kg−1), total ice number Ni,tot (kg−1), rime mass qi,rim (kg kg−1), and rime volume Bi,rim (m3 kg−1) mixing ratios. One can calculate the vapor deposition ice mass qi,dep from qi,tot minus qi,rim. Several relevant bulk properties, including the rime mass fraction (Fi,rim = qi,rim/qi,tot), the bulk ice density, the bulk rime density (ρi,rim = qi,rim/Bi,rim), the mean particle size, and the mean number- and mass-weighted fall speeds are derived from these four conserved prognostic variables.
The total number and ice mass mixing ratios are, respectively, given by
where D is the maximum dimension of the ice particles, md(D) is the mass–dimension relationship, and the particle size distribution (PSD) is assumed to follow a gamma function (all symbols for variables and parameters used in the paper are defined in Table 1). The gamma distribution is described by N0, λ, and μ, respectively, being the intercept, slope, and shape parameters. In P3, the shape parameter follows from the observational study of Heymsfield (2003): μ = 0.001 91λ0.8 − 2, where λ has units of per meter. For a given combination of prognostic variables (Ni,tot, qi,tot), to solve for the PSD intercept and slope parameters, the md(D) is needed.
where ρi = 917 kg m−3 is the density of bulk ice and ρg is the density of fully rimed ice (graupel/hail). The regimes are bounded by three threshold diameters,
that can evolve within the PSD as a function of the ice regime. The parameters αva and βva are empirical constants, and herein, we use default P3 values of and βva = 1.9 following Brown and Francis (1995), modified for the correction proposed by Hogan et al. (2012). Diameter Dth in (4a) smoothly divides small spherical ice particles from unrimed nonspherical particles and is found by equating (3a) and (3b). This threshold is needed because extrapolation of the md(D) relationship for unrimed nonspherical ice to sizes smaller than Dth would give particle densities larger than ρi, which is unphysical. Diameter Dgr, given by (4b), smoothly separates dense, unrimed, nonspherical ice particles from graupel/hail and is found by equating (3b) and (3c). This is the size at which the masses of a completely rime-filled particle (graupel/hail) and of an unrimed nonspherical ice particle are equal for a given diameter D. Finally, Dcr in (4c) smoothly divides graupel/hail from partially rimed nonspherical ice and is obtained by equating (3c) and (3d). Physically, this threshold originates directly from the assumption that the rime mass fraction does not vary with D, giving a critical size Dcr below which particles are completely filled in with rime.
The value of ρg is an Fi,rim-weighted average of the rime density ρi,rim and the density underlying the unrimed structure of the particle ρur (see MM15 for details). Because of their complicated interrelationship, the equations for Dgr, Dcr, ρg, and ρur are solved together by iteration (MM15; Dietlicher et al. 2018). The projected area–diameter relationship Ad(D) also follows specific power laws for various size-dependent regimes consistent with the md(D) relationship (MM15). The terminal velocity–diameter relationship Vt(D) is computed using md(D) and Ad(D) according to Mitchell and Heymsfield (2005).
Once the parameters for the PSD (N0, λ, and μ) and the md(D), Ad(D), and Vt(D) relationships are obtained, all the microphysical process rates are integrated. This is done offline, and the values are stored in lookup tables as a function of the normalized total ice mass qi,tot/Ni,tot, the bulk rime density ρi,rim, and the bulk rime mass fraction Fi,rim.
3. Parameterization description of the bulk liquid fraction
a. Overview of the liquid mass mixing ratio
To simulate mixed-phase particles and Fi,liq explicitly in P3, a new conserved prognostic variable, the liquid mass mixing ratio on ice particles qi,liq, has been added to P3_ORIG. The modified scheme will be referred to as P3_MOD. In P3_MOD,
where qi,ice = qi,rim + qi,dep. The bulk rime mass and liquid mass fractions are thus defined by
respectively. The conservation equation for qi,liq is
where t is time, ρa is the air density, u is the 3D wind vector, z is height, Vm is the mass-weighted fall speed (see appendix A), Δ*(qi,liq) is a subgrid-scale mixing operator, and is a source–sink term that includes various microphysical processes. The qi,liq microphysical tendency is
where Qi,mlt is the melting, Ql,wgrth the wet growth, Ql,coll,r (Ql,coll,c) the collection of rain (cloud droplets), Ql,shd the shedding, Ql,frz the refreezing, and Ql,dep is the vapor transfer. The latter can be a source term (condensation) or a sink term (evaporation). The source–sink terms are described in appendix A. Note, the modified definitions of qi,tot and Fi,rim in (5) and (6), respectively, are consistent with P3_ORIG since, for Fi,liq = 0, these variables simply revert back to the original formulations in MM15.
b. Main assumptions
A major question when predicting Fi,liq in bulk schemes is how the particle size evolves during melting. In traditional multimoment bulk parameterizations, the mean size of ice particles typically does not change during melting (e.g., Morrison et al. 2009; Thompson et al. 2008; Milbrandt and Yau 2005). This is a common closure assumption in order to calculate the decrease in number concentration during melting. However, observations show a decrease of size during melting for individual particles (Fujiyoshi 1986). This decrease can be captured by predicting Fi,liq. Past studies approximated the spherical drop equivalent diameter Deq as a function of Fi,liq and utilized a “melting critical diameter” D*, such as Thériault and Stewart (2010) and Szyrmer and Zawadzki (1999). The melting critical diameter determines the largest particles in the PSD that will melt completely within one time step, and it depends on Fi,liq. Thus, this method requires iteration to solve Deq because of the interdependence of D* and Fi,liq. This is one reason why this approach for the implementation of Fi,liq has not been widely used in traditional bulk schemes designed with several solid precipitation categories.
Past theoretical and experimental studies have described the melting behavior of ice particles depending on their size and type. Rasmussen et al. (1984a) and Rasmussen and Pruppacher (1982) found that small spherical ice particles with D < ~0.1 cm melt quickly, in less than 1 min, into liquid drops. Fujiyoshi (1986) showed that the liquid water produced by the melting of unrimed nonspherical ice particles accumulates to form wet snow. Rasmussen et al. (1984b) found that a portion of the liquid produced by the melting of large rimed ice particles (D > 0.1 cm) is shed, while the rest is accumulated around the ice core. Leinonen and von Lerber (2018) showed that the melting of lightly rimed crystals is similar to the melting of unrimed aggregates (no shedding), whereas the melting of moderately rimed crystals is similar to the melting of graupel/hail (with shedding).
Based on Rasmussen et al. (1984a), all liquid water produced by the melting of small spherical ice particles (D ≤ Dth) in P3_MOD is transferred to the rain category in one time step. Thus, it is assumed that Dth corresponds to the melting critical diameter in P3_MOD. For the parameter values of αva and βva used here, Dth is ~66 μm. Based on the studies discussed above, shedding during melting is considered when Fi,rim > 0, detailed in appendix A.
The simplest approach to describe the melting process is to assume that the liquid water is uniformly distributed around an ice core. For simplicity and because of a lack of detailed observations, ice cores are assumed to have the same properties (mass, projected area, capacitance, ventilation coefficient, and so on) as in P3_ORIG. Melting decreases ice-core mass qi,ice but not the total ice number Ni,tot, except for particles with D ≤ Dth that are transferred to the rain category as described above; thus, the mean size of ice cores decreases during melting. Because the bulk ice particle density increases with decreasing particle size in P3, this assumption is physically reasonable and implies a tendency toward small spherical ice cores with a density equal to ρi as particles melt.
According to the assumptions below, the most straightforward approach for the implementation of Fi,liq into P3 is to separate processes based on whether they apply to the ice core embedded within the particle or the whole mixed-phase particle (liquid and ice). This way, each process can be integrated over the appropriate size distribution based on the microphysical process computed. It is assumed that melting depends directly on properties of the ice core embedded in mixed-phase particles. Thus, the bulk melting rate (detailed in appendix A) is calculated by integration over the PSD computed using the ice-core diameter Dd. Ice-core PSD parameters N0,core, μcore, and λcore are computed using qi,ice = (1 − Fi,liq)qi,tot and md(D) following (3a)–(3d). Note that sublimation/deposition processes also depend directly on properties of the ice core embedded in the whole particle, as explained in appendix A.
Other microphysical processes, such as refreezing, condensation/evaporation of the liquid component, self-collection, shedding, and collection with other particle categories depend on properties of the whole particle (both the ice and liquid components). Bulk rates for these processes (also detailed in appendix A) are calculated by integrating over the PSD computed using the full mixed-phase particle diameter Dp. To compute the whole particle PSD parameters N0,p, μp, and λp, the mass–diameter relationship mt(Dp, Fi,liq) is needed. For simplicity and because of lack of observations, it is assumed that at a given Dp, mt(Dp, Fi,liq) is calculated by a linear interpolation based on Fi,liq:
Here, md(Dp) is given by (3a) to (3d) and (ρw = 1000 kg m−3), evaluated over the full mixed-phase PSD (from 0 to ∞). The relationship mt(Dp, Fi,liq) is used to solve the PSD parameters N0,p, μp, and λp over qi,tot defined by
The linear interpolation of Fi,liq for mt(Dp, Fi,liq) is consistent with a density of the mixed-phase particle being equal to a Fi,liq-weighted average of the ice and liquid parts, which was also used by Szyrmer and Zawadzki (1999). Also, this approach is consistent with 1) assuming a constant Fi,liq with diameter over the size distribution, 2) physically consistent behaviors in the limits of Fi,liq = 0 and Fi,liq = 1, 3) a mixed-phase particle density that must always be less than or equal to the density of a liquid drop with the same diameter, and 4) a straightforward computation of the size distribution parameters as a function of Fi,liq.
An approach similar to (10) is made for the projected area At(Dp, Fi,liq) and for the terminal velocity Vt(Dp, Fi,liq) relationships. These are, respectively,
where Ad(Dp) is the P3_ORIG projected-area relationships, , and
where Vt(Dp, Fi,liq = 0) and Vt(Dp, Fi,liq = 1) are given in appendix B [(B4) and (B6), respectively, but with Dd replaced by Dp]. These expressions follow Mitchell and Heymsfield (2005) and Simmel et al. (2002), respectively.
The Vt(Dp, Fi,liq) of unrimed ice (Fi,rim = 0) and fully rimed ice (Fi,rim = 1) at four different values of Fi,liq are shown in Figs. 1a and 1b, respectively. The terminal velocity increases with both the particle size Dp and Fi,liq in agreement with the observations of M90 for unrimed ice (Fig. 1a). Direct comparison with the observed and experimental terminal velocities of M90 is, however, difficult because the M90 relationship can only be applied to unrimed ice particles. Also, the M90 fall speed relationship used a different ice-core mass–dimension relationship compared to P3. There are no observations of terminal velocities of partially melted rimed (partially or fully) ice particles. Based on Rasmussen et al. (1984b), it is assumed that for large particles, accumulated liquid water is shed. However, for unrimed and partially rimed ice, the liquid water is retained and forms a shell, which flattens as the particles fall, thereby limiting the increase of fall speed with Dp similar to the behavior of large purely liquid drops. This is seen in Fig. 1b for particles larger than 0.4 cm.
4. Comparison with a benchmark melting model
a. Experimental design
A benchmark comparison of P3_MOD with a detailed Lagrangian model using the M90 microphysics parameterization is presented. It aims at illustrating differences between P3_MOD and M90 in the melting behavior of a population of ice particles. M90, Yokoyama and Tanaka (1984), and Matsuo and Sayso (1981) have evaluated the distance below the 0°C level for a single falling ice crystal to melt completely using a Lagrangian model with their respective melting equations. Here, this distance is evaluated for a population of unrimed nonspherical ice particles to calculate the evolution of melting particles as they fall, with comparison between the P3_MOD scheme and a Lagrangian model we developed with the M90 parameterization. Two different populations of ice particles, shown in Fig. 2a, provide input boundary conditions. The first, PSD 1, includes a moderate number of large ice particles characterized by a total ice mass mixing ratio of qi,tot = 0.6 g kg−1 and a total number mixing ratio of ice of Ni,tot = 10 550 kg−1. The second, PSD 2, includes a larger number of small ice particles characterized by qi,tot = 0.65 g kg−1 and Ni,tot = 41 550 kg−1. Rimed particles are not included in this experiment because they are not accounted for in the M90 parameterization.
Both the Lagrangian benchmark and P3_MOD bulk model use the melting equation developed in appendix A [see (A3)]. To solve (A3), both models use a temperature variation with height of 0.6°C (100 m)−1 as shown in Fig. 2b. For both models, melting is the only microphysical process included, and feedback between temperature and latent cooling during melting is neglected; the temperature profile is constant in time. The parameters Le, Lf, ka, Dυ, qs,0, qυ, and ρa are computed with the air temperature and pressure, assuming hydrostatic balance. The atmosphere is saturated with respect to water.
The Lagrangian model solves (A3) using the M90 parameterization. It is numerically integrated as an initial-value problem using a simple Euler forward scheme with a time step of 1 s. The model is initialized by dividing the particle population into 10 000 particle sizes with varying Dd, covering a size range from 0.0001 to 2 cm. This provides the evolution of particle properties along vertical trajectories for each of the 10 000 initial sizes. Bulk properties of the model, such as Fi,liq and qi,ice, are calculated from these trajectories every meter below the 0°C level. The particle fall speed follows from M90 [Geresdi et al. 2014; see their (3)].
P3_MOD solves bulk equations for the 1D (vertical) model prognostic variables as they evolve in time. Sedimentation is calculated using a simple first-order upwind approach using the mean number- and mass-weighted bulk fall speeds described in appendix A. The equations are solved using a time step of 1 s and a vertical grid spacing of 1 m. Results for comparison with the Lagrangian model are taken after 10 min of simulation time.
Plots of Fi,liq and the qi,ice as a function of the distance below the 0°C level are shown in Figs. 3a and 3b, respectively. Overall, the Lagrangian benchmark model and P3_MOD give similar results. For a given T, differences in Fi,liq between P3_MOD and the Lagrangian model are always smaller than 0.15 and relatively larger for small Fi,liq. The melting process in P3_MOD is somewhat slower than the Lagrangian model for Fi,liq < 0.6–0.7, leading to a slightly larger qi,ice in P3_MOD. This is compensated by slightly faster melting in P3_MOD for larger Fi,liq so that the total distance for melting is nearly the same as in the Lagrangian model.
The small differences between P3_MOD and the Lagrangian model seen in Fig. 3 can be explained mainly by differences in the parameterization used for the melting equation associated with the capacitance C(Dd, Fi,liq) and ventilation coefficient F(Dd, Fi,liq) relationships. For the Lagrangian model, the M90 relationships are based on Pruppacher and Klett (1997). For P3_MOD, the relationships are given in appendix B [(B1) and (B2), respectively]. Figure 4 shows comparisons between M90 and P3_MOD of the capacitance (Fig. 4a) and the ventilation coefficient (Fig. 4b) for unrimed ice (Fi,rim = 0). At a given Dd, the capacitance of P3_MOD is slightly smaller than M90, especially at low Fi,liq. Larger differences occur for the ventilation coefficient owing to its dependence on the terminal velocity parameterization in (A3). For a given Dd, the ventilation coefficient of P3_MOD is greater (smaller) than M90 at low (high) Fi,liq. As detailed in appendix B, the ventilation coefficient for P3_MOD is computed by linear interpolation as a function of Fi,liq between the ventilation parameters for dry ice in P3_ORIG when Fi,liq = 0 and those corresponding to a liquid drop when Fi,liq= 1. On the other hand, in M90, the ventilation coefficient is computed using their experimental terminal velocity, which is small at low Fi,liq and increases quickly at Fi,liq > 0.7. Also seen in Fig. 3, the differences between P3_MOD and the Lagrangian model using M90 are larger for PSD 1 than for PSD 2, especially when Fi,liq is small. This is because PSD 1 has a higher concentration of large particles, and differences in the ventilation coefficient parameterization between P3_MOD and M90 are greater for larger particles than smaller ones (Fig. 4b).
A comparison between M90 Lagrangian model and P3_MOD for the evolution of the PSDs as the two populations of ice particles fall below the 0°C is shown in Fig. 5. Note that the PSDs in Fig. 5 are shown as a function of the ice-core diameter Dd. For both initial PSDs 1 and 2 from P3_MOD, the shape parameter of the gamma PSD remains unchanged and is equal to 0 at small liquid fractions. However, when the liquid fraction increases, the slope parameter increases, which leads to an increase of the shape parameter. This behavior is seen with PSD 2 at Fi,liq = 0.9. The intercept parameter also increases with Fi,liq for both PSDs. In general, the three parameters of the PSD increase with the increasing Fi,liq, and their changes depend on the rime mass fraction and the values of qi,tot and Ni,tot that characterize the PSD. The evolution of the PSDs for the P3_MOD simulations compares well to that using the M90 Lagrangian model, especially at low Fi,liq, and reflects the decrease of particle size and shift of the PSD toward smaller sizes during melting (Thériault and Stewart 2010; Thériault et al. 2006). Overall, these simulations suggest that P3_MOD behaves realistically compared to M90.
Instead of comparing with the Lagrangian M90 model, the P3_MOD melting process for Fi,rim = 0 could also be compared to the Thériault and Stewart (2010) bulk microphysics scheme. It is believed that results will be comparable, since Thériault and Stewart (2010) used the same parameterizations as M90. However, as the parameterizations are very different between Thériault and Stewart (2010) and P3_MOD, and the Thériault and Stewart (2010) scheme does not include the prediction of rime fraction, P3_MOD will only be compared to P3_ORIG in the following.
5. Idealized simulations using P3_MOD
a. Experimental design
Characteristics of the precipitation type are investigated using P3_MOD coupled with a simple 1D kinematic model. P3_MOD is compared to P3_ORIG to illustrate differences in the precipitation properties with the prediction of the Fi,liq and the overall behavior of the modified scheme.
All 1D simulations are initialized with vertical profiles of temperature and water vapor mass mixing ratio for 50 vertical levels evenly spaced. The grid spacing is 60 m. The simulated period is 2 h, and the time step is 1 s. Hydrostatic balance is assumed. The precipitation type characteristics are analyzed after 90 min of simulation time. A total ice mass mixing ratio of qi,tot = 0.265 g kg−1 and total ice number mixing ratio of Ni,tot = 5000 kg−1 provide boundary conditions at the domain top. This corresponds to a snowfall rate of 1 mm h−1 when Fi,rim = 0. Sensitivity of the precipitation properties to variation of the Fi,rim value specified at the domain top is also investigated. The initial qi,rim, with a fixed bulk rime density of 900 kg m−3, is systematically modified to increase the snowfall rate at the top of the column such that it reaches 2.7 mm h−1 when Fi,rim = 1; the increased snowfall rate for a given qi,tot occurs because the mass-weighted mean fall speed increases with Fi,rim. The vertical air motion is zero throughout the column.
Two observed vertical profiles of temperature and dewpoint temperature were used to initialize the 1D simulations, shown in Fig. 6. Case 1 (Fig. 6a) has a vertical profile that was measured by a gondola during the Science of Nowcasting Olympic Weather for Vancouver 2010 (SNOW-V10; Isaac et al. 2014) field campaign. It is characterized by melting near the surface. A rain–snow boundary was observed along Whistler Mountain at around 2300 UTC 7 March 2010 (Thériault et al. 2014), with a mixture of wet snow and rain at the base of the mountain. In this profile, the near-surface melting layer is 500 m deep, and the environmental conditions are subsaturated with respect to liquid water in cold layers and at the melting-layer top.
The case 2 (Fig. 6b) thermodynamic profile is based on observations obtained on 2315 UTC 1 February 1992 at St. John’s, Newfoundland, Canada (Hanesiak and Stewart 1995). It is associated with a melting layer atop a refreezing layer below. A mixture of ice pellets and some needles were reported at the surface around 2317 UTC, and a mixture of ice pellets, freezing rain, and needles at around 2345 UTC. In the refreezing layer of this profile, the environmental conditions are near saturation with respect to liquid water, meaning that it is supersaturated with respect to ice. The top of the melting layer is subsaturated with respect to liquid water.
b. Case 1: Melting layer near the surface
Vertical profiles of the temperature, mass and number mixing ratios, liquid fraction, and mass-weighted density and fall speed after 90 min are shown in Fig. 7 for both P3_ORIG and P3_MOD. These results are for unrimed ice; that is, Fi,rim = 0 at the domain top. The main microphysical processes (not shown) for both simulations are sublimation in the cold layers (T < 0°C) and melting in the warm layers. Collection of rain and evaporation occur in the warm layer, with evaporation of melting ice neglected in P3_ORIG. In P3_MOD, these processes are sources and sinks for the liquid component, qi,liq. Collected rain in P3_ORIG is shed back to rain at temperatures above freezing, while in P3_MOD, the collected rain is added to qi,liq. The evaporation of qi,liq in P3_MOD cools the air near the top of the melting layer compared to P3_ORIG (Fig. 7a). The surface precipitation type formed in P3_ORIG is rain (Fig. 7b), while in P3_MOD, there is a mixture of rain and almost completely melting ice (Fig. 7b), the latter characterized by Fi,liq = 0.98 (Fig. 7d). The total number mixing ratio of rain (Fig. 7c) is higher in P3_MOD compared to P3_ORIG because the source term for rain number concentration is parameterized differently. In P3_ORIG, the rain number source from melting is proportional to the respective changes in qi,tot, while it is proportional to qi,ice in P3_MOD ( appendix A). Moreover, P3_ORIG calculates the number of raindrops formed by melting using a scaling factor of 0.2 to account for the rapid evaporation of small melting particles; that is, it assumes a proportion of one raindrop formed per five melted ice particles. No such scaling is applied in P3_MOD. Both mass-weighted density (Fig. 7e) and fall speed (Fig. 7f) are larger in P3_MOD compared to P3_ORIG in the melting layer, which has significant impacts on the melting time of particles; faster fall speeds mean less time to melt over a given distance.
Results from varying the rime mass fraction Fi,rim specified at the domain top are shown in Figs. 8 and 9. Vertical profiles of Fi,liq after 90 min (Fig. 8a) illustrate that at low Fi,rim, particles melt almost completely before reaching the surface. This is reflected by the percentages of surface precipitation types reaching the surface (Fig. 8b). The surface precipitation type is only rain using P3_ORIG for an initial Fi,rim ≤ 0.2, while P3_MOD produces a mixture of rain and very wet ice of Fi,liq > 0.95. For initial values of Fi,rim > 0.5, both P3_ORIG and P3_MOD show a mixture of surface precipitation types. Differences between P3_ORIG and P3_MOD are mainly due to collection of rain, shedding, and the production of faster-falling particles by P3_MOD. For example, the collection of rain by partially melted ice increases with the specified Fi,rim at the domain top because there is more shedding, which reduces the amount of rain reaching the surface in P3_MOD. Also, as seen in Fig. 9, an increase of the mean mass-weighted density (Fig. 9a) and the mean mass-weighted fall speeds (Fig. 9c) occurs in P3_MOD compared to P3_ORIG (Figs. 9b and 9d, respectively), which, as mentioned before, impacts the time spent by particles in the warm layer. The main added value of P3_MOD for this case is the ability to produce wet snow at the surface.
c. Case 2: Melting layer aloft
Results for case 2, specifying Fi,rim = 1 at the domain top, are shown in Fig. 10. As for case 1, only melting cools the environmental air in P3_ORIG, while there is also cooling from evaporation of qi,liq in P3_MOD at the top of the melting layer (Fig. 10a). In both P3_ORIG and P3_MOD, a mixture of supercooled raindrops and rimed ice is produced in the refreezing layer (Fig. 10b). Note that supercooled raindrops represent freezing rain as the surface temperature is <0°C in both P3_MOD and P3_ORIG. Note also that since Fi,rim = 1, the lines for qi,tot and qi,rim are superimposed for P3_ORIG. In P3_MOD, rimed ice is produced from the complete refreezing of partially melted ice particles aloft, which is the process forming ice pellets. The proportion of ice to the total precipitation is higher in P3_MOD compared to P3_ORIG because of the partial melting and the refreezing of partially melted ice, which is neglected in P3_ORIG. Values of Fi,liq for partially melted ice entering the cold layer are ~0.8 and decrease below 1.2 km because of refreezing (Fig. 10d). The rime mass fraction (Fig. 10d) remains close to 1 and the rime density (Fig. 10e) to 900 kg m−3 at all heights, which follows from the parameterization of melting and refreezing as described in appendix A. The mass-weighted mean density (Fig. 10e) and fall speed (Fig. 10f) increase in the warm layer and slightly decrease in the cold layer, the latter associated with ice pellets.
The temporal evolutions of P3_MOD and P3_ORIG air temperature are shown in Figs. 11a and 11b, respectively. Cooling due to melting and warming due to refreezing are the main processes affecting temperature, giving profiles that tend toward a 0°C isothermal layer. P3_MOD is generally colder in the melting layer and warmer in the refreezing layer than P3_ORIG (Fig. 11c). The formation of ice pellets by refreezing in P3_MOD warms the air compared to P3_ORIG. Differences in the melting layer are associated with faster-falling particles when Fi,liq > 0 after the onset of melting in P3_MOD, which vertically extends the region where particles melt compared to P3_ORIG. In the cold layer, temperature differences are larger because of refreezing in P3_MOD, which does not occur in the P3_ORIG simulation.
Percentages of surface precipitation types from varying Fi,rim at the domain top for case 2 are shown in Fig. 12. For Fi,rim ≤ 0.6 at the domain top, both P3_MOD and P3_ORIG produce only freezing rain at the surface over the entire period. For Fi,rim ≥ 0.6 at the model top, both P3_MOD and P3_ORIG show a mixture of ice and freezing rain. However, ice pellets formed from refreezing are the dominant precipitation type in P3_MOD when Fi,rim > 0.65 at the model top, whereas freezing rain is the dominant type in P3_ORIG. The refreezing process in P3_MOD has a major impact on the precipitation types reaching the surface because the generation of ice pellets from refreezing leads to an increase in the collection of supercooled raindrops, in turn reducing freezing rain at the surface consistent with Barszcz et al. (2018).
6. Summary and conclusions
A new parameterization approach is proposed to predict the bulk liquid mass fraction of mixed-phase particles, Fi,liq, in the predicted particle properties (P3) bulk microphysics scheme. The modified scheme, P3_MOD, can explicitly simulate the evolution of bulk mixed-phase particle properties, thus improving the representation of key microphysical processes such as melting, evaporation/condensation, and refreezing. It also allows the explicit prediction of several winter precipitation types, such as freezing rain, ice pellets, and wet snow because of the addition of a new prognostic variable, the bulk liquid mass mixing ratio accumulated on ice qi,liq.
P3_MOD produced comparable melting rates to a Lagrangian benchmark model based on the M90 melting parameterization, which was developed from observations. This supports the viability of P3_MOD and its bulk representation of the melting behavior of ice particles. Prediction of the liquid fraction affects the mean fall speed and density of hydrometeors falling into a melting layer, which in turn impacts distributions of latent cooling as well as other microphysical processes such as collection and condensation/evaporation, compared to P3_ORIG. Predicting the liquid fraction also allows the refreezing of partially melted ice particles to be explicitly represented, which increases the formation of ice pellets and the ratio of ice pellets to supercooled rain compared to P3_ORIG. Finally, the sensitivity of surface precipitation characteristics to riming differs between P3_MOD and P3_ORIG. For a melting layer 500 m deep and a surface temperature of 2°C, an increase in rime mass fraction leads to an increase in the ratio of snow to liquid precipitation in P3_MOD compared to P3_ORIG.
Overall, implementation of the bulk liquid fraction in P3 is a step forward toward better prediction of precipitation type and distribution when the temperature is near 0°C. Important precipitation types, such as ice pellets and hail, involve tracking mixed-phase particles. Forecasting these precipitation types will therefore benefit from the prediction of the bulk liquid fraction. Although the focus of this work is on improving the representation of wintertime mixed-phase precipitation, the parameterized microphysics in P3_MOD is general and thus the inclusion the predicted liquid fraction should also improve the simulation of hail through better representation of shedding during wet growth and melting. The effects of P3_MOD on the simulation of hail will be examined in a future study. Since P3_MOD only involves one additional prognostic variable per ice category, the additional computational cost is relatively small. It could be used in numerical weather prediction (NWP) as well as in convection-permitting climate models (CPCM).
Although the main objective of the paper was to describe the new approach for implementing Fi,liq into P3, and we show that the approach is comparable to the detailed Lagrangian model, comparison with observations is needed to validate the new parameterization. Therefore, in future work, P3_MOD will be tested by simulating a freezing rain and ice pellets storm using a three-dimensional atmospheric model, with comparison to detailed observations.
The Microphysical Process Formulations
The melting source/sink is divided into two terms. The first term is the meltwater transferred to the rain category Qr,mlt, and the second is the meltwater accumulated on ice particles Qi,mlt:
where Dd is the maximum dimension of the ice core and is the melting rate of a single particle. During melting, qi,rim is reduced to maintain a constant rime mass fraction. The ice (rain) number sink (source) terms are proportional to the respective changes in qi,ice: Nr,mlt = Qr,mltNi,tot/qi,ice. All remaining ice is transferred to the rain category when Fi,liq > 0.99.
M90 developed an equation for the melting behavior of an individual nonspherical unrimed ice crystal based on observations (Fujiyoshi 1986), using the thermodynamic model of Mason (1956). Two processes control the melting rate. The first is latent heating (cooling) associated with condensation (evaporation) of water vapor at the particle’s surface. The second is heat transfer from the environment to the particle surface and from the surface to the embedded ice core. During the melting process, it is assumed that the particle surface temperature is T0 = 273.15 K. The change of ice-core mass from melting is
where Dυ is the diffusivity of water vapor in air, Lf (Le) is the latent heat of fusion (evaporation), ka is the thermal conductivity of air, qυ is the water vapor mass mixing ratio, qs,0 is the saturated water vapor mass mixing ratio at the surface of the particle, C(Dd, Fi,liq) is the capacitance [ appendix B; (B1)], and F(Dd, Fi,liq) is the ventilation coefficient [ appendix B; (B2)]. The melting rate is treated in a simple way following most parameterizations (e.g., M90) that assumes the same ventilation coefficient and capacitance for both the latent term and the heat diffusion term in (A3).
The refreezing process describes how the liquid water surrounding the ice core freezes when mixed-phase particles are transported to regions where T < 0°C. The refreezing source/sink term Ql,frz is
where Dp is the full mixed-phase particle diameter and is the refreezing rate of a particle. The term Ql,frz represents a mass transfer from qi,liq to qi,rim, and the rime density of refreezing ice is assumed to be near solid bulk ice (900 kg m−3), following other drop freezing process in P3_ORIG.
From Pruppacher and Klett (1997), assuming a quasi-steady state during the refreezing process and a surface particle temperature of T0 = 273.15 K, the freezing rate equation is a balance between the latent heat of freezing and the conduction through the ice shell formed around the whole particle. The conduction through the ice shell is given by the conduction of heat between the surrounding air and the particle and the heat exchange by evaporation/condensation at the particle surface:
c. Collection with liquid-phase categories
It is assumed that rain Ql,coll,r and cloud water Ql,coll,c mass collected when T ≥ 0°C is transferred into qi,liq. At T < 0°C and in the dry growth regime, the collected rain and cloud masses are assumed to freeze instantaneously and are transferred to qi,rim, as in P3_ORIG (MM15). At T < 0°C and in the wet growth regime, calculated following Musil (1970), the total collected rain and cloud mass Ql,wgrth is assumed to be a source for qi,liq. Wet growth of hail is not the focus of the paper; however, prediction of Fi,liq in P3_MOD presents an interesting possibility for improving the prediction of wet growth and shedding from hail.
The rain and cloud collection rates are parameterized using a collection kernel derived from the projected area and the terminal velocity relationships for rain and ice or mixed-phase particles, numerically integrated over the respective PSDs. The fall speed of cloud droplets is neglected in the kernel equation.
d. Vapor transfer
For simplicity, deposition/sublimation of ice in P3_MOD is allowed only when Fi,liq = 0 because liquid water is assumed to be distributed around the ice core when Fi,liq > 0. The calculation of deposition/sublimation follows from P3_ORIG (see MM15). Note that the Fi,liq threshold for sublimation/deposition versus condensation/evaporation processes can be modified easily in P3_MOD. For example, it could be increased to 0.2 as in Thériault and Stewart (2010). If so, then the sublimation/deposition process would be calculated using the ice-core properties as for the melting process.
Vapor transfer of qi,liq (denoted Ql,dep), which represents the liquid mass of mixed-phase particles lost by evaporation or gained by condensation, is also computed with the quasi-analytic formulation for supersaturation of MM15, using liquid-phase thermodynamic parameters. The Ql,dep source/sink term when Fi,liq > 0 is
where δt=0 is the initial supersaturation, Δt is the time step, Γl the psychrometric correction associated with the latent heating/cooling, and Ai,wet is the change in δ due to vertical motion, turbulent mixing, and radiation. The overall supersaturation relaxation time scale in conditions when Fi,liq > 0 (τ) is . The supersaturation relaxation time scale associated with mixed-phase particles τi,wet is
where τc and τr are the relaxation time scales for cloud and rain, respectively (Morrison and Grabowski 2008), and C(Dp, Fi,liq) and F(Dp, Fi,liq) are the capacitance [ appendix B; (B7)] and ventilation coefficient [ appendix B; (B8)] relationships. During the evaporation process (when Ql,dep < 0), which is computed gradually as a function of the liquid fraction, it is assumed that Ni,tot decreases proportionally to the change in qi,tot: Nl,evp = Ql,depNi,tot/qi,tot.
Based on Rasmussen et al. (1984b), it is assumed that only ice particles with diameter larger than 9 mm within the PSD shed. The mass of qi,liq due to shedding Ql,shd is given by the total integrated liquid mass of particles with Dp > 9 mm within the PSD and interpolating as a function of Fi,rim, with no shedding when Fi,rim = 0 and all liquid mass (for Dp > 9 mm) shed when Fi,rim = 1. The increase in the rain number mixing ratio from shedding assumes a mean raindrop diameter of 1 mm.
The variables qi,liq, qi,tot, qi,rim, and Bi,rim use the total mass-weighted fall speed Vm for their sedimentation, while the number mixing ratio Ni,tot uses the number-weighted fall speed VN, given by, respectively,
All integrations are done offline, and the values are stored in lookup tables as a function of 50 values of normalized total ice mass qi,tot/Ni,tot, five bulk rime densities ρi,rim (50, 250, 450, 650, and 900 kg m−3), four bulk rime mass fractions Fi,rim (0, 0.333, 0.667, and 1), and four bulk liquid mass fractions Fi,liq (0, 0.333, 0.667, and 1).
Capacitance and Ventilation Coefficient
a. Processes depending on the ice-core properties
The capacitance C(Dd, Fi,liq) accounts for nonspherical shape of ice particles undergoing melting or deposition/sublimation. The capacitance is calculated simply by a linear interpolation based on Fi,liq between the capacitance of an ice particle Cd(Dd) when Fi,liq = 0 and that of a spherical drop Cliq(Dd) = 0.5Dd when Fi,liq = 1, for a given Dd:
where Cd(Dd) is a function of the ice particle properties for each regime of the distribution as detailed in MM15.
The ventilation coefficient F(Dd, Fi,liq) for the melting calculation is also given by linear interpolation over Fi,liq between the value for an ice particle Fd(Dd) when Fi,liq = 0 from P3_ORIG and that of a liquid drop Fliq(Dd) when Fi,liq = 1:
where Fd(Dd) is (Thorpe and Mason 1966)
and , Red = Vt(Dd, Fi,liq = 0)Ddρa/υ is the Reynolds number, υ is the dynamic viscosity of air, and Sc is the Schmidt number Sc = υ/(ρa Dυ). The terminal velocity Vt(Dd, Fi,liq = 0) follows Mitchell and Heymsfield (2005) with Fi,liq = 0, as in P3_ORIG:
The term Fliq(Dd) is (Pruppacher and Klett 1997)
where and ρw,g = 1 g m−3 is the density of water.
Equations (B1) and (B2) are the formulations to express the capacitance and ventilation coefficient as a function of the ice-core diameter, in contrast to those used in M90 that also depend on the real particle diameter. Although the capacitance and ventilation coefficient used for the latent heat term in the melting (A3) should also depend on the liquid part of the particle, the melting equation is treated following most parameterizations (e.g., Milbrandt and Yau 2005; M90) that assume the same ventilation coefficient and capacitance for both the latent heat term and the heat diffusion term in (A3). Note also that the relationships for capacitance and ventilation coefficient differ for Fi,liq = 0 between M90 and P3 and the M90 coefficients are only for unrimed particles, which is why we use the P3 coefficients rather than M90.
b. Processes depending on properties of the entire particle
For condensation/evaporation and refreezing of mixed-phase particles with Fi,liq > 0, the capacitance C(Dp, Fi,liq) is similar to (B1) but uses the full particle diameter Dp:
The ventilation coefficient F(Dp, Fi,liq) is (Thorpe and Mason 1966)
The authors would like to thank the Canadian Network for Regional Climate and Weather Processes (CNRCWP) and the NSERC Discovery grant for funding this research. One of the authors (Mélissa Cholette) would like to thank the Fonds de Recherche du Québec—Nature et Technologies (FQRNT) and the UQAM Faculty of Science for graduate fellowships. We would like to thank the three anonymous reviewers for their constructive comments.
Denotes content that is immediately available upon publication as open access.