Growth of ice particle mass and projected area during riming

There is a long-standing challenge in cloud and climate models to simulate the process of ice particle riming realistically, partly due to the unrealistic parameterization of the growth of ice particle mass (m) and projected area (A) during riming. This study addresses this problem, utilizing ground-based measurements of m and ice particle maximum dimension (D) as well as theory to formulate simple expressions describing the dependence of m and A on riming. It was observed that β in the m−D power law m = αD appears independent of riming during the phase 1 (before the formation of graupel), with α accounting for the ice particle mass increase due to riming. This semi-empirical approach accounts for the degree of riming and renders a gradual and smooth ice particle growth process from unrimed ice particles to graupel, and thus avoids discontinuities in m and A during accretional growth. Once the graupel with quasispherical shape forms,D increases with an increase inm and A (phase 2 of riming). The treatment for riming is explicit, and includes the parameterization of the ice crystal–cloud droplet collision efficiency (Ec) for hexagonal columns and plates using hydrodynamic theory. In particular, Ec for cloud droplet diameters less than 10 μm are estimated, and under some conditions observed in mixed-phase clouds, these droplets can account for roughly half of the mass growth rate from riming. These physically meaningful yet simple methods can be used in models to improve the riming process.


Introduction
Milbrandt, 2015). 23 Note that riming occurs only when ice particles have a D greater than the riming threshold size 24 (D thres : the smallest ice crystal D for which riming can occur). Early observations (Harimaya,25 1975) and numerical studies (Pitter and Pruppacher, 1974;Pitter, 1977) determined a D thres being 26 around 300 µm. However, it was later shown by both observational (Bruntjes et al., 1987) and 27 numerical (WJ00) studies that such D thres is around 35 µm, 110 µm, and 200 µm for hexagonal 28 columns, hexagonal plates, and broad-branched crystals, respectively (note that all these 29 dimensions are along a-axis of crystals). Many models calculated ice particle mass by assuming that ice particles are spherical (e.g. Gettelman, 2008). However, this assumption is not realistic, and produces errors in the evolution 3 of snow-size spectra (Mitchell, 1988). Based on observations, several studies developed ice 4 particle mass-dimension (m-D) power law parameterizations to reduce the dimensionality of 5 complex ice particle shapes. For a specific ice particle shape or an environmental condition, this 6 relationship has the form: where both α and β are constants over a specific size range. They are determined via direct 8 measurements of ice particle mass and dimension (Locatelli and Hobbs, 1974; M90), or are 9 constrained through aircraft measurements of the ice particle size distribution (PSD) and IWC 10 (Heymsfield et al., 2010; Cotton et al., 2012). The prefactor α was considered to contain 11 information on particle density and thickness, whereas β was believed to have information on 12 particle shape. We will discuss the latter in Sect. 4.1 for the riming process. Similar power laws 13 have been developed for projected area-dimension (A-D) relationships: where γ and δ are constants over a specific size range derived by direct measurements of ice 15 particle projected area and dimension (M96). When comparing rimed particles with the same size, 16 lump graupel has the largest mass and area relative to cone-like graupel or hexagonal graupel, and expressions over the PSD because the size range for each microphysical step (ice crystal, graupel, 26 and transition) was calculated in a way to provide only continuous mass, and thus produced 27 discontinuities in projected area. JH15 developed a detailed ice growth model that simulates ice 28 particle habit and mass via the processes of vapor deposition and riming. This model is also a 29 single-category scheme, but it does not employ m-D and A-D power laws; instead, it computes the growth of ice particles along the major and minor axes of oblate or prolate spheroids (representing 1 hexagonal plates or columns). Therefore, the model is able to simulate simple ice particle shapes, 2 and also captures the temperature-dependency of vapor deposition and the riming processes (since 3 particle shape is a function of temperature and relative humidity; Magono and Lee 1966;4 Pruppacher and Klett 1997). The simulated ice particle shape, mass, and fallspeed are in good 5 agreement with observational data from wind tunnel experiments on ice crystal growth. 6 The computation of rime mass (an increase in ice particle mass by riming) in models is performed 7 by calculating the accretional mass growth rate (Heymsfield, 1982;Mitchell, 1995;JH15). When 8 an ice particle falls in a cloud of supercooled cloud droplets, the increase in its mass due to 9 accretion depends on ice particle features (e.g. fallspeed and projected area), droplet characteristics 10 (e.g. mass and number concentration of droplets), and the collision efficiency (E c ) between an ice 11 particle and droplet. More details on mass growth rates are provided in Sect. 6, and E c is described 12 in the next section.  15 One important factor in the modeling of riming is the calculation of the E c between ice particles 16 and cloud droplets (Pruppacher and Klett, 1997). E c was calculated as a function of ice particle D 17 and cloud droplet diameter (d) via both experimental measurements (Sasyo and Tokuue, 1973, 18 hereafter ST73; Kajikawa, 1974, hereafter K74; Murakami et al., 1985) and theoretical/numerical 19 calculations (Beard and Grover, 1974;Pitter and Pruppacher, 1974;Schlamp, 1975;Pitter, 1977;20 Wang and Ji, 2000, hereafter WJ00). The difference in E c between various studies is due to the 21 strong sensitivity of E c to the ice particle shape as well as the assumptions and limitations in 22 different studies. Experimental measurements of E c have been conducted in vertical wind tunnels.

23
Such studies are rare due to the difficulty and limitations of experiments, and were limited to only 24 planar ice crystals or circular disks with D > 1 mm (Reynolds number (Re) > 40). Murakami et al.

25
(1985) studied the E c between polystyrene latex spheres (d < 6 µm) and freely-falling planar ice 26 crystals (1.5 mm < D < 5 mm, and 70 < Re <300). ST73 investigated fixed hexagonal plates (5 27 mm < D < 20 mm) that were exposed to water droplets contained in airflow in a vertical wind for columnar crystals (Schlamp, 1975 density and surface roughness). To solve this issue, Rasmussen and Heymsfield (1985) suggested 6 that E c between graupel and cloud droplets can be calculated by modification of the results of 7 Beard and Grover (1974) for E c between raindrops and water droplets. Similar to Beard and 8 Grover (1974), they employed the superposition method for collision between particles, but they 9 assumed that the small cloud droplets do not change the graupel fallspeed, therefore they used 10 Stokes number instead of mixed Froude number in the non-dimensionalized momentum equation 11 (see Eqs. 1-6 in Rasmussen and Heymsfield, 1985 investigated experimental E c between freely falling spherical ice particles (initially 580 µm < D < 13 760 µm) and water droplets (20 µm < d < 40 µm) in a vertical wind tunnel with laminar flow.
14 They showed that collection kernels of ice particles are higher than that of raindrops, and therefore 15 calculated a correction factor to account for the error in E c , when assuming raindrops instead of 16 graupel.

17
The objective of this study is to develop various empirical and theoretical approaches to represent 18 the continuous and gradual growth of ice particle mass and projected area during riming in a 19 realistic and yet simple way, suitable for models. field campaign on cloud seeding from 1986 to 1988, and for one part of that project, natural ice particles were collected during snow storms in a polystyrene petri dish and then the particles were 1 photographed using a microscope equipped with a camera. Then a heat-lamp was used to melt 2 these ice particles, and immediately after melting another photograph was taken of the hemispheric 3 water drops (contact angle on polystyrene = 87.4 degrees). The images were used later in the lab to 4 measure the maximum dimension (D) of individual ice particles (defined as diameter of a 5 circumscribed circle around the particle). In addition, the diameter of the water hemispheres was 6 measured, and from this the volume and mass of individual ice particles were computed. Also 7 recorded were individual ice particle shapes, which were classified using the Magono and Lee 8 (1966) nomenclature scheme. The level of riming (i.e. light, moderate, heavy riming, or graupel) 9 was indicated based on this scheme, and the temperature range over which the observed ice 10 particle shape originated was recorded (e.g. for long columns, -8°C < T < -6°C). These riming

13
Photographic examples of these rimed particle types are shown in Fig. 2 of Locatelli and Hobbs 14 (1974). Software was developed to extract all combinations of particle shapes (for a detailed 15 explanation of sampling and measurements, see M90). SCPP is a unique dataset that measures 16 both ice particle size and mass and also determines the degree of riming. As explained in EM16, 17 the important problem with airborne measurements is that they are unable to measure single ice Particles In Cirrus (SPartICus) field campaign for D < 100 μm and a subset of SCPP data for D > 23 100 μm. Since CPI does not measure ice particle mass, EM16 developed a method that calculates 24 mass from the measurements of ice particle projected area, D, and aspect ratio by assuming that having graupel-like centers but with rimed branches extending outwards revealing the dendritic 12 origin. Also classified were total of 67 lump graupel (R4b), cone-like graupel (R4c), and 13 hexagonal graupel (R4a); R4b and R4c are graupel with non-discernable original habit, whereas 14 R4a forms just prior to R4b or R4c, with its hexagonal origin still recognizable. 15 In order to represent the natural variability of ice particle mass, all identifiable particles are 16 initially shown with their actual mass and maximum dimension. Thereafter, to quantify the 17 variability and to further investigate m-D power laws and the rimed-to-unrimed mass ratio, the ice 18 PSDs were divided into size bins with intervals of 100 μm between 100 and 1000 μm, and with 19 subsequent intervals of 200, 200, 400, 600, 600 and 1000 μm (up to 4000μm) at larger sizes to 20 supply sufficient sampling numbers in each size bin (for more details on bin intervals, see Table   21 1). In order to investigate the riming effect, all identifiable particles were sorted into either one of 22 three rimed categories or an unrimed category. Both unrimed and lightly-rimed ice particles are 23 included in the unrimed category, whereas the three rimed categories consist of densely rimed, 24 heavily rimed, and graupel particles. The purpose of this section is to demonstrate how the cold habit SCPP curve fit from EM16 (based 28 on unrimed ice crystals) compares with all the SCPP data, since this shows how the EM16 curve 29 fit appears representative for all ice particles sampled during SCPP and thus may be representative for Sierra Nevada snowfall. This comparison is shown in Fig. 1a for all ice particles that could be 1 classified (3781 ice particles). The curve fit appears to bisect the data well. Moreover, it is seen 2 that rimed ice particles tend to have larger mass on average, compared to unrimed ice particle of μm. In order to be even more quantitative, the percent difference between the total SCPP mean ice 19 particle mass in each size-bin of Fig A long-standing problem in cloud modeling is the treatment of α, β, γ and δ as a function of ice 3 particle riming. Since riming leads to graupel formation and graupel tends to be quasi-spherical, it 4 is intuitive to assume that β and δ will approach limiting values of 3 and 2, respectively 5 (corresponding to ice spheres), as more and more supercooled liquid water is accreted by an ice 6 particle to produce graupel. One common approach in many cloud models (that use an m-D 7 relationship) is to assume that β is equal to ~ 2 for unrimed crystals and is equal to ~ 3 for graupel. 8 This implies that riming enhances β. This assumption is tested in this section by using SCPP data 9 with the objective of developing observational-based guidelines for modeling the process of 10 riming. To test this assumption for β, the size-resolved masses of rimed and unrimed ice particles 11 from the same basic shape category are needed. In this section, we used heavily rimed dendrites 12 (R3a, R3b and R3c) and unrimed dendrites (P1e, P1d and P1f). In addition, this data was 13 partitioned into the same size-intervals described earlier to calculate the mean m and D in each 14 size-interval for unrimed and heavily rimed dendrite crystals, along with their ζ. All these results 15 are shown in Fig. 2. Size-intervals having less than 3 measurements are not represented. Most of 16 the data for unrimed crystals is associated with D > 600μm. One can see quantitatively how the 17 mean masses for rimed dendrites are substantially greater than those for unrimed dendrites on 18 average for the same size-interval, in agreement with the hypothesis of Heymsfield (1982). 19 Using only the size-intervals containing at least 3 measurements, the m-D power law for the 20 unrimed dendrites is: and for heavily rimed dendrites is: where all variables have cgs units. If the size-interval corresponding to the largest unrimed 23 dendrites is not used in the least-square fit calculation, the m-D expression for unrimed dendrites 24 becomes: having an exponent nearly identical to that in Eq. (4) for heavily rimed dendrites. This is contrary 1 to most cloud models that assume different ice categories (snowfall with β ~ 2 and graupel with β 2 ~ 3) and an abrupt increase in β upon a change in ice category (autoconversion). Based on SCPP 3 observations, it is apparent that the traditional hypothesis that β increases with riming is not 4 correct, at least not before the graupel onset. This can be understood by noting that β does not 5 necessarily indicate the morphology of an ice particle within a given size-interval, but rather 6 indicates the mass rate-of-change with respect to size (since β is the slope of the m-D line in log-7 log space). This can also be seen qualitatively in Fig. 2, where the rimed and unrimed data points 0.0078D m  that represents a slight increase in β for graupel which is 14 significantly less than spherical β (which is equal to 3). By assuming that initial graupel mass can where ρ i is a reduced density, such ρ i is equal to 0.18 g cm -3 for D 16 = 500 µm, which is lower than the ρ i for heavily rimed graupel in the dry growth regime (ρ i = 0.4 g 17 cm -3 ; Rutledge and Hobbs, 1984; Ferrier, 1994). This assumption would produce a fit parallel to 18 the ice spheres fit in Fig. 2, and is poorly fitted to the SCPP R4b and R4c data (R 2 = 0.67), 19 compared to the power law fit (R 2 = 0.94). 20 All these observations are in agreement with the experiment of Rogers (1974)  for both planar and columnar ice crystalsA similar assumption is also valid for hexagonal columns.

28
The impact of moderate to heavy riming on β for hexagonal columns was demonstrated in M90 29 (see their Table 1 and Sect. 4d). For these columnar crystals, riming had no effect on β (i.e., β was 1.8 for both rimed and unrimed columns), indicating that riming can be modeled by only 1 increasing α for these crystals. Thus, it appears justified to treat β as constant during the riming 2 process (until spherical shape) for both dendritic and columnar ice crystals: where subscript u denotes unrimed conditions. The IWC is defined as: where n(D) is number distribution. We explained that β and D do not change during riming. Also 5 unchanged is n(D), because the number of ice particles in each size bin is not affected by riming. 6 Therefore, the dependence of α on riming can be calculated by knowing the contribution of riming 7 to the IWC: Since β is essentially the same in Eqs. (4) and (5), their prefactor ratio (α in Eq. 4 divided by α in 9 Eq. 5, which is equal to 2.12) indicates that riming contributed slightly more than half the mass of deviation from the mean for 300 µm < D < 400 µm may be due to only a single unrimed ice crystal 15 of anomalous mass in this size bin. 16 Equations (4) and (5) also suggest a means of adapting the m-D curve fit in Fig. 1 for modeling the 17 riming process in mixed phase clouds. Since this curve fit is representative of ice particle 18 populations in frontal clouds (containing a mixture of unrimed and rimed particles), it can be 19 adapted for modeling the riming process in frontal clouds. Since β should be essentially the same 20 for both unrimed and the mixture of unrimed plus rimed SCPP ice particles, the ratio of their 21 corresponding prefactors (i.e. α u /α mix ) can be multiplied by the mass predicted by the curve fit 22 equation to yield the masses appropriate for unrimed particles. For the ice particles plotted in Fig.   23 1a, m u /m mix is equal to 0.650 (where m mix includes all these particles and m u /m mix was calculated by the same method that calculated m r /m u in Fig. 3). This implies that multiplying the mass predicted 1 by the curve fit in Fig. 1 by a factor of 0.65 will yield masses proper for unrimed ice particles. To 2 model the riming process in frontal clouds, these unrimed particles can be subjected to the riming 3 growth equations described below as well as Eq. (8). MG08 used different A-D power laws in each riming step, but such method led to discontinuities 7 in projected area during the transition from one ice category to another one. It seems that the A-D 8 and m-D that they used were not self-consistent (e.g. they were from different studies based on 9 different datasets). Here, we suggest an approach to avoid the discontinuity in projected area.
The reason for this can be explained by noting that the riming process often affects A but does not 18 change D (by filling the space between ice particle branches) significantly prior to graupel 19 formation. This is also evident from observations, as shown in Table 1 where A max is the maximum projected area due to riming in the phase 1 (which is the graupel A), 24 and R is the riming factor defined as:  found that their aspect ratio does not exceed 0.8. Using this value, JH15 showed good agreement 6 between their model and observational data from a wind tunnel. Based on such analysis, k is equal 7 to 0.8. Further observational data are needed to determine the value of A max more accurately. 8 So far, we discussed the phase 1 of riming growth (before the formation of graupel), where m and 9 A increases while D and therefore β and δ are conserved. Once the graupel stage is attained, phase 10 2 of riming starts and the graupel continues to grow through riming, and a different methodology 11 is required to describe riming growth at this growth stage, because graupel D increases by riming. 12 Once m = m max , then a graupel bulk density is defined as: where m is calculated as described in Sect. 6. As before, for a given D, , and in this way 17 riming growth is treated for all conditions.

Columnar ice crystals
1 Figure 4b represents m r /m u between graupel (R4b and R4c) and unrimed columnar crystals (N1e 2 and N2c) in order to determine m max for columnar crystals (initial graupel at the end of phase 1).

3
Relatively small variability of m r /m u (between 1.6 and 3) is found for D < 1400 µm, with larger 4 variability (from 1.4 to 9.4) found for larger ice particles, with the weighted average of m r /m u 5 equal to 2.4, and therefore m max ≈ 2.4×m u . The higher variability for D > 1400 μm is likely due to a 6 single graupel particle per size-bin. Based on SCPP dataset, we showed that D and β are constants 7 during the phase 1 of riming, and since initial graupel mass is 2.4 times larger than unrimed 8 column mass, these mean that graupel density is ~ 2.4 times larger than unrimed column density.  11 Some of the data shown in Fig. 2 describes an experiment investigating the ability of the Baker and 12 Lawson (2006) (hereafter BL06) m-A power law to reproduce the masses of unrimed dendrites that 13 presumably have relatively low area ratios (the ratio of the actual ice particle projected area to the 14 area of a circle having a diameter equal to D). A study by Avramov et al. (2011) found that this 15 power law overestimated the masses of low-density dendrites (P1b), high-density dendrites (P1c), 16 and low density dendrite aggregates, but that the BL06 power law yielded masses consistent with 17 high density dendrite aggregates at commonly observed sizes. It is important to understand the 18 potential limitations of this power law for dendrites due to their abundance in Arctic mixed phase 19 clouds and for the modeling of these clouds. BL06 used a subset of SCPP data (e.g. 865 ice 20 particles), of which 550 were identifiable, and 36% of such identifiable particles were moderately 21 or heavily rimed. They then developed a software to calculate ice particle projected area from their 22 magnified images. Thereafter, they calculated a m-A power law expression. Since BL06 used only  Although A is more strongly correlated with ice particle m than is D (based on BL06), inferring m 5 or volume from a 2-D measurement is still ambiguous since different crystal habits exhibit

Hexagonal plates 25
The numerical study of WJ00 is valid for unsteady flow, hexagonal ice plates with 1 < Re < 120 26 and 160 µm < D < 1700 µm, and water droplets having diameter d between 1 µm and 100 µm. Re 27 for hexagonal plates is calculated based on D (e.g., , where ε is kinematic viscosity). Since there is not sufficient agreement between the historical H80 relationship and the 1 data of WJ00, we provided best fits to the data of WJ00 that has the form of:  (15) where K is mixed Froude number of the system of water drop-ice particle, and is calculated as: where v is water drop fallspeed, and g is gravitational acceleration. Since cloud water drops are in 4 Stokes regime, v is calculated as the Stokes fallspeed (e.g., water density, a  is air density, and μ is dynamic viscosity), and K is the same as the Stokes 6 number in this flow regime. K crit is the critical value of K (where E c equals 0 in the third line in Eq. 7 15) and is expressed as a function of ice particle Re: .
Based on Eq. (15), E c in the third line is physically meaningful only when K ≥ K crit . When K < 9 K crit , E c in the third is imaginary and must be set to zero in order to avoid errors. K thres is the 10 threshold of K between small and large cloud droplets, and is calculated as: The resulting curve fits for E c (Fig. 5a) show that the provided equations can represent the data of 5 WJ00 very well in various ranges of K and Re. The percent error in E c between curve fits and 6 WJ00 data has a mean value of 6.65% with standard deviation of 3.67% for all Re and K. ≤ 640 are slightly lower than curve fit for Re = 120. This does not seem to be a discrepancy, 16 because it is observed from the curve fits (based on WJ00) that E c is not sensitive to Re when Re ≥ 17 60. This is also observed in K74 for large Re (their Fig. 14). E c from ST73 for Re = 97 is in good 18 agreement with the curve fit for K ~ 1.5, but is larger than the curve fit for K ~ 0.3. It is noteworthy 19 to explain the shortcomings of these experiments, as mentioned by Pruppacher and Klett (1997). 20 For the experiment of K74, when Re > 100, the flow is unsteady and leads to the eddy shedding 21 and formation of wakes at the top of the particle, which increases the uncertainty in fallspeed. For surrogates suitable for spherical raindrops.
Note that Eqs. (15)- (17) are derived for the range over which the data of WJ00 is valid (e.g., 1 < 1 Re < 120), and they should not be used for extrapolation to Re values larger or smaller than this 2 range. Since Re < 1 corresponds to ice particle smaller than D thres , it is justified to assume that E c = 3 0 in this Re range. When considering the range Re > 120, values of E c for Re = 120 should be 4 used; this is reasonable based on the experiments of K74 for 200 < Re <640, and the theoretical 5 study of WJ00 for 60 ≤ Re ≤ 120.  resolving models, except for the study of JH15 that calculated E c for prolate spheroids based on 11 their aspect ratio for use in models. In addition to hexagonal plates, WJ00 studied E c between 12 hexagonal columns (with width w between 47 and 292.8 µm, length l between 67.1 and 2440 µm 13 and 0.2 < Re < 20) and water drops of 1 µm < d < 100 µm. Note that WJ00 calculated Re for 14 columns in a different way than was done for plates. Re for columns was calculated from their 15 width, whereas Re for plates was computed from D (e.g.,  / Re wV columns  ). If the values of Re 16 were calculated from the column maximum dimension, they would have values comparable to 17 those for plates. In formulating E c for columns, we have followed the Re convention of WJ00. 18 Similar to hexagonal plates, we provide the best fits to the data of WJ00 for hexagonal columns:  (19) and r is a parameter related to the major radius of the ellipse fit and is determined as:  (20) and K thres is calculated as:  (21) The results are shown in Fig. 5b. Similar to hexagonal plates, the curve fits are able to represent unknown, but we suggest using E c for Re = 20 as a conservative underestimate of E c . 12 13 6 Mass growth rate by riming 14 In Sect. 4, the dependence of α on IWC was explained. Unrimed IWC can be derived from α and β 15 pertaining to unrimed ice crystals (see EM16). The riming rate for a single ice particle of size D 16 can be calculated by using the definition of riming mass growth rate, similar to Heymsfield (1982), 17 M95 and JH15: where t is time, d is diameter of a cloud droplet, A g (D,d) is the geometrical cross-section area of distribution, and d max is diameter of the largest cloud droplet. Note that the cloud droplet . For this equation, the riming rate is not sensitive to the droplet 5 distribution. . However, such a 10 relationship cannot represent the evolution of ice particle size and shape, and is often inconsistent 11 with the realistic dependence of V on the ice particle m/A ratio. This increases uncertainty in the developed an alternative method to improve M96 method, and calculated X as a function of m/A r 21 ratio, where A r is area ratio (defined as the ratio of ice particle projected area to the projected area 22 of a circumscribed circle around the particle; see Eq. 15 in Erfani and Mitchell, 2016).

23
Since the contribution of the cloud droplet projected area to A g (D,d) is negligible, A g (D,d) can be 24 approximated as the maximum ice particle cross-section area projected normal to the air flow. Ice 25 particles fall with their major axis perpendicular to the fall direction, therefore A g (D,d) is 26 approximated as the ice particle A, which is calculated in Sect. 4.2. The m(d) is calculated from 27 spherical geometry as: where E c was discussed in Sect. 5, and E s is the sticking efficiency (fraction of the water droplets that stick to the ice particle after collision), and is presumed to be unity since supercooled cloud droplets freeze and bond to an 1 ice particle upon collision. Conditions under which E s may be less than unity are addressed in 2 Pruppacher and Klett (1997). It is noteworthy that by using the above calculations, riming growth 3 will be represented in a self-consistent, gradual, and continuous way. Based on the explanations in 4 this section, Eq. (22) can be reduced to: second term on the RHS should be relatively small (riming has little impact on D prior to graupel 7 formation). Therefore, to a first approximation: and together with Eq. (23), a change in α due to riming can be determined. Since D and β do not 9 change by riming, dα/dt is linearly proportional to dm/dt. divides the droplet PSD mass into equal parts). E c is calculated from Eqs. (15) and (18) for 13 hexagonal plates and hexagonal columns, respectively, and a sub-exponential PSD is assumed for 14 cloud droplets that has the form: where λ is the PSD slope parameter, ν is the PSD dispersion parameter and N o is intercept 16 parameter. M95 used observational droplet spectra from Storm Peak lab (Steamboat Springs, 17 Colorado, USA), and calculated various PSD parameters:  Fig. 6a that 20 dm/dt for riming increases with increasing ice particle D. The dm/dt is linearly proportional to 21 LWC when MMD and D are constant. In addition, when LWC is constant, doubling MMD (from 8 22 to 16 µm) leads to a quadrupling of dm/dt. One important feature is the contribution of small droplets (d < 10 µm) to dm/dt, when K < 0.7 and E c < 0.3. It is seen in this figure that when MMD 1 is relatively small (= 8 µm), ignoring such small droplets results in values of (dm/dt) riming at the 2 largest crystal sizes that are ~ 40% (for plates) and ~ 70% (for columns) of those obtained when all 3 droplets are included. That is, small droplets contribute about 60% and 30%, respectively, to the 4 (dm/dt) riming values at the largest sizes. This surprising contribution from small droplets is partly 5 due to half of the LWC being associated with d < 8 μm. However, when MMD is larger (= 16 µm), 6 the contribution from small droplets is only ~ 5%. The size-dependence of dm/dt for hexagonal 7 columns (Fig. 6b) shows that dm/dt for columns is larger than that for hexagonal plates for a 8 specific crystal size when droplet MMD is 8 µm, partly because columns fall faster than plates (see  The collection kernel (K c ) can be calculated as A(D)V(D)E(D,d), which is alternatively equal to 13 dm/dt divided by the ice particle mass due to riming (see Eq. 23). MG08 approximated this 14 variable by using simple assumptions, and found that it is proportional to D 2 . Here, we showed by 15 more accurate analysis that K c has a form of a second-order polynomial fit, and is represented for multiplied by the riming fraction m r /m u or alternatively IWC/IWC u . A similar strategy could be 5 adopted for other ice particle shapes or shape mixtures in frontal clouds, as is done for columnar 6 particles in this study. By using this method, there is no discontinuity in the growth of m and A; 7 rather, the particles grow gradually during riming process. Phase 2 of riming starts when graupel 8 with quasi-spherical shape forms. In this phase, the increase in m and A causes an increase in D. 9 It is straightforward for models with multiple ice categories to utilize our new method. This can be 10 done by describing riming growth as two phases and removing the autoconversion process. Phase 11 1 simulates rimingthe growth from an unrimed ice crystal to the onset of graupel formation. In this 12 phase, ice particle mass and projected area gradually increase, but size is unchanged (Eqs. 6-11). 13 Phase 2 represents graupel growth. In this phase, the shape is unchanged, but mass, projected area, 14 and size gradually increase (Eqs. [12][13][14]. By using the method introduced in this study, the models 15 may still use multiple categories (e.g. ice crystal, rimed particle, graupel), but within each category 16 the rimed mass and projected area fraction can there is a gradually changeincrease, thus preventing 17 an abrupt change in ice particle attributes between from one categoriesy to another one. 18 Prior to this work, there was only none rigorous practical method ftor calculatinge the droplet size- 19 dependence of E c for use in models explicitly for columnar crystals. cloud droplet d and ice particle D in non-steady flow. In the future, this treatment of the riming 28 process will be employed in a new snow growth model that predicts the vertical evolution of ice 29 particle size spectra, mass, projected area, fallspeed, and snowfall rate in terms of the growth processes of vapor diffusion, aggregation and riming. These results will be compared with airborne 1 measurements from two spiral descents.        Additional curves (dashed red and dashed black curves) are produced by assuming that E c conforms to the ellipse 5 curves and is zero for smaller droplets.