Ice Particle Properties Inferred From Aggregation Modelling

Abstract

We generated a large number 105,000 of aggregates composed of various monomer types and sizes using an aggregation model. Combined with hydrodynamic theory, we derived ice particle properties such as mass, projected area, and terminal velocity as a function of monomer number and size. This particle ensemble allows us to study the relation of particle properties with a high level of detail which is often not provided by in situ measurements. The ice particle properties change rather smoothly with monomer number. We find very little differences in all particle properties between monomers and aggregates at sizes below 1 mm which is in contrast to many microphysics schemes. The impact of the monomer type on the particle properties decreases with increasing monomer number. Whether, for example, the terminal velocity of an aggregate is larger or smaller than an equal-size monomer depends mostly on the monomer type. We fitted commonly used power laws as well as Atlas-type relations, which represent the saturation of the terminal velocity at large sizes (terminal velocity asymptotically approaching a limiting value) to the data set and tested the impact of incorporating different levels of complexity with idealized simulations using a 1D Lagrangian super particle model. These simulations indicate that it is sufficient to represent the monomer number dependency of ice particle properties with only two categories (monomers and aggregates). The incorporation of the saturation velocity at larger sizes is found to be important to avoid an overestimation of self-aggregation of larger snowflakes.

Introduction

The terminal velocity v_term of ice monomers and aggregated ice particles and its relation to size has manifold impacts on precipitation and radiative effects of ice containing clouds. For example, Morales et al. (2019) show that parameters describing v_term of aggregates have the largest impact on the precipitation of simulated orographic clouds. Experiments with global climate simulations revealed that also radiative fluxes are very sensitive to changes in v_term (Jakob, 2002). Sanderson et al. (2008) found that v_term of ice is the second most influential parameter for the climate sensitivity in their multimember perturbed physics General circulation model ensemble. Constraining v_term of cloud ice and aggregated ice particles can reduce the degrees of freedom in model tuning (e.g., to improve top of atmosphere radiative fluxes Schmidt et al., 2017) and improve the physical consistency in atmospheric models.

The importance of v_term of ice particle has been early recognized and has motivated first observational studies in the first third of the 20th century. Using initially manual observations and microphotography, pioneering studies such as (Brown, 1970; Kajikawa, 1972; Langleben, 1954; Locatelli & Hobbs, 1974; Nakaya & Terada, 1935; Zikmunda & Vali, 1972) investigated the relation of v_term to the particle's size for various ice particle habits and aggregates. In addition to the direct measurements of velocity, several studies started to investigate the principle relation between particle properties such as mass, size, and projected area to v_term which allows deriving v_term from these quantities (Cornford, 1965; Heymsfield, 1972). Due to the large efforts in performing these often manual measurements, the sample size of the derived relations is rather small. For example, some of the relations of the widely used relations by Locatelli and Hobbs (1974) are only based on 10–50 particles. One can assume that particles with ideal monomer types might have been subjectively chosen in order to easier associate the derived relationships to certain well-defined shapes. Nevertheless, the relations of size, mass, area, and v_term derived in these early studies are still used in microphysics parameterizations (e.g., the v_term-size relation of the snow category in Morrison and Milbrandt (2015) is taken from Locatelli and Hobbs (1974) mixed aggregates; see Figure 1). In Figure 1a, a selection of the aforementioned v_term relations is shown for their defined size range. The spread of velocities for different ice particle monomers is relatively high (e.g., Kajiwa, 1972) reported v_term to be about 0.2 m s⁻¹ for a dendrite but about 0.5 m s⁻¹ for a plate monomer. In contrast, v_term of aggregates of different monomer types appears to be relatively similar and always close to 1 m s⁻¹ in the reported size range.

Figure 1. (a) In situ measurements of vterm of monomers (separated by monomer type; blue Kajikawa, 1972) and aggregates composed of different monomers (green: LH74 Locatelli & Hobbs, 1974) and particle ensembles from the PIP-CARE data set (see section A0.1). (b) vterm of unrimed ice particles in two-moment microphysics schemes. The blue line represents the implementation of cloud ice (monomers), the green line the implementation for the snow (aggregates) category in (solid lines, SB Seifert & Beheng, 2006) and (dashed lines, Morr Morrison et al., 2005). The Predicted Particle Property (P3) scheme (Morrison & Milbrandt, 2015) assumes identical properties for all unrimed particles (yellow line).

Evolving computer technology allowed the realization of automated particle measurement systems such as the 2D Video Disdrometer (2DVD Kruger & Krajewski, 2002), the Snow Video Imager (SVI Newman et al., 2009), its successor the Particle Imaging Package (PIP Tiira et al., 2016), the Hydrometeor Velocity and Shape Detector (HVSD Barthazy et al., 2004), or the Multi-Angle Snowflake Camera (MASC Garrett et al., 2012). These systems are based on optical methods to capture particle size and terminal velocity. Unlike in the early studies, particle property relations (Barthazy & Schefold, 2006; Brandes et al., 2008; Garrett & Yuter, 2014; Zawadzki et al., 2010) are now based on a very large number of particles which are classified by automated algorithms rather than visual selection (Bernauer et al., 2016; von Lerber et al., 2017). All optical disdrometers have a smallest detectable size limit (e.g., 0.1–0.2 mm for 2DVD), which implies that measurements close to this limit should be interpreted with care. A general behavior, which is revealed by all instruments, is a “saturation” of aggregate terminal velocities (i.e., terminal velocities asymptotically approaching a limiting value) at approximately 1 m s⁻¹ for unrimed particles and sizes larger than a few millimeters (Figure 1a).

Most ice microphysics schemes use two categories for unrimed ice particles, which are commonly denoted as cloud ice and snow/aggregates. Relations between particle properties, such as size (e.g., the maximum dimension D_max), mass m, projected area A, or v_term, are defined for each category. Examples of the v_term dependence on size which are implemented in widely used two-moment schemes are shown in Figure 1b. When comparing these relations with observations (Figure 1a), we miss the saturation behavior of v_term for larger sizes in most relations. This discrepancy is expected as most schemes use power laws, which are unable to represent a saturation behavior. Alternative “Atlas-type” three-parameter fits have been suggested (Seifert et al., 2014), but so far, they have not been tested thoroughly. The recent Predicted Particle Properties (P3) scheme (Morrison & Milbrandt, 2015) uses only one ice category and a look-up table approach for v_term, which better matches the saturation at large sizes. At the smaller size range, the snow category is found for all schemes to fall significantly faster than the ice category with the same size. Considering that v_term depends strongly on m and A of the particle, it might sound plausible, that for example, an aggregate of a few plates should fall faster than a single plate of the same size. Unfortunately, most observations do not provide sufficiently detailed information about monomer number and type which would be needed to answer the question of whether there exists a “jump” in v_term for the number of monomers exceeding a certain threshold. Fairly direct observations of the particles' m and A are only available from manual, particle-based observations (e.g., Locatelli & Hobbs, 1974).

An interesting new tool to better understand the underlying principles of aggregation and its effects on particle properties are aggregation models (Hashino & Tripoli, 2011; Leinonen & Moisseev, 2015; Ori et al., 2014; Przybylo et al., 2019; Westbrook et al., 2004a). Those models use idealized monomer shapes (e.g., dendrites, needles, plates, and columns) with particle properties matched to in situ observations. Aggregates simulated with the model by Westbrook et al. (2004a) helped to better understand theoretical scaling relations associated to aggregation such as the increase of aggregate mass with size by a power of two (Westbrook et al., 2004b), which was known from several previous in situ observations. This model has been extended by Leinonen and Moisseev (2015) providing a large number of monomer shapes and also provides an option to rime the aggregate (Leinonen & Szyrmer, 2015). This allowed to better understand the evolution of size and mass of a large number of aggregates which were increasingly rimed (Seifert et al., 2019).

To infer v_term from modeled ice particles or aggregates, computational fluid dynamics is an accurate but also computational costly method. It has been recently applied to idealized ice particle shapes (Bürgesser et al., 2019; Hashino et al., 2016; Nettesheim & Wang, 2018), and more computations with more complex shapes can be expected shortly. Hydrodynamic theory is a computational cheaper alternative to calculate v_term based on a number of bulk particle characteristic, rather than the complex 3D shape (e.g., Böhm, 1992; Heymsfield & Westbrook, 2010; Khvorostyanov & Curry, 2005). The accuracy of hydrodynamic theories has recently been evaluated by ice particle analogs falling in an oil tank (Westbrook & Sephton, 2017). The experimental results show deviations smaller than 20% for the Heymsfield and Westbrook (2010) theory. A problematic aspect of these theories is still the formulation of the scaling toward higher Reynolds number (i.e., large particles) and the simulation of more complex particle shapes (Westbrook & Sephton, 2017).

Aggregation models in combination with hydrodynamic theory have recently been used to study v_term of aggregates (Hashino & Tripoli, 2011; Schmitt et al., 2019). Hashino and Tripoli (2011) identified a dependency of the aggregation rate and aggregate mass on the mean size and type of the monomers. Schmitt et al. (2019) analyzed v_term and its variability of simulated aggregates composed of hexagonal prisms taken from a monodisperse monomer size distribution. They found that the variability of v_term is caused by the variability of the number of monomers N_mono and the monomers' aspect ratio.

In this study, we aim to study the dependency of m, A, and v_term on size, monomer number, and type. For this, we create a large number of aggregates with various monomer types including also mixtures of different monomer types. The monomer size is sampled from a size distribution rather than a constant size to better represent real ensembles of aggregates. Central questions of this study are, how important is the monomer number and type information for parameterizing aggregate properties and how well can they be parameterized by different functional relations?

To answer these questions, we describe in section 2 the aggregation model and the created data set of unrimed aggregates as well as the hydrodynamic theory to calculate v_term based on m and A of these particles. The simulated particle properties are compared to in situ observations in section 3. Section 4 presents several parameterizations of the particle properties. Finally, in section 5, we use a 1D Lagrangian particle model to test the impact of including different complexity of particle properties for aggregation

Methods Aggregation Model

We use the aggregation model developed by Leinonen and Moisseev (2015) which includes a large number of realistic monomers (hexagonal plates, dendrites, columns, and needle). Originally, the aggregation model was designed to produce realistic snow particle structures which can then be used to calculate their scattering properties (Leinonen & Moisseev, 2015; Leinonen et al., 2018). The model has also been used to systematically investigate microphysical processes, such as riming (Seifert et al., 2019).

The shape characteristics (length, thickness, etc.) of the monomers are predefined by geometric relations based on in situ observations (Leinonen & Moisseev, 2015). The aggregation process starts with generating N_mono monomers with sizes following a predefined inverse exponential probability density function p_d(D_max), [Image Omitted. See PDF]where λ⁻¹ is the size parameter of the monomer distribution and D_max is the maximum size of the monomer. The higher λ⁻¹, the larger are the sizes of the monomers.

The monomers sizes are sampled from the monomer distribution and assembled until an aggregate consisting of N_mono monomers is build up. In each aggregation step, pairs of particles are selected according to a simplified gravitational collection kernel. The probability distribution of collision among each possible particle pair is calculated as being proportional to the particle geometric cross sections and differential fall speed (Westbrook et al., 2004a). The two colliding particles form an aggregate which then becomes one of the candidates for the next aggregation step. This process includes the collision between aggregates. The aggregation code is publicly available at https://github.com/jleinonen/aggregation, and more details on the implementation can be found in Leinonen and Moisseev (2015). During the aggregation process, the collecting particles are partially aligned with the principal axis in the x-y plane. Rotations around the principal axis are performed randomly with a standard deviation of 40°. The collected particles are randomly aligned, which mimics the complex flow in the vicinity of other particles (Leinonen & Moisseev, 2015).

The aggregation simulations performed in this study differ from previous studies in two main aspects. The first aspect is the resolution of the particle structure. The particle is internally represented by a three-dimensional lattice with a predefined distance of the volume elements of typically 40 μm. This distance was found to be sufficiently small for scattering computations, while being coarse enough in order to keep the numerical costs for the scattering computations in a reasonable range. However, we discovered that for small particle sizes, the theoretical relations for certain particle properties (see Fig. 1 in Leinonen and Moisseev (2015)) are not exactly matched by the discretized particle. This discrepancy can be easily explained when considering for example that plate monomers with D_max < 3.03 mm consist of only one layer of volume elements if the default resolution of 40 μm is used. This does not necessarily affect the aggregate properties of those monomers as shown in Leinonen and Moisseev (2015); however, in our study, the focus is to investigate the transition from small to larger sizes particles. Hence, we need to refine the resolution especially for small particles.

As a compromise between computational feasibility and having fine enough resolved particles, aggregates with N_mono ≤ 100 are simulated with a resolution of 5 μm, while aggregates with N_mono ≥ 100 are simulated with 10-μm resolution. With a resolution of 5 μm (10 μm), a plate monomer with D_max= 3 mm has a thickness of 4 (8) volume element layers. It should be noted that the sensitivity to resolution is smaller for monomer types with less extreme aspect ratios (e.g., columns).

The second major difference to previous aggregation studies using the model by Leinonen and Moisseev (2015) is that we extended the code in a way that we can also generate aggregates composed of monomers with different habits. The motivation for this new feature was based on observations that larger snowflakes often consist of a mixture of dendrites and needles (Lawson et al., 1998). The modified code extends Equation 1 to be the joint distribution of multiple mono-dispersed distributions. Each monomer distribution is defined by its own settings (e.g., monomer type, mean size, and truncation). The joint distribution is defined by the relative weights of each mono-dispersed distribution. These modifications have been merged to the main aggregation code and are also publicly available.

In order to account for a large variability of naturally observed particle shapes (Bailey & Hallett, 2009), we simulated a large suite of aggregates consisting of plates, columns, dendrites, needles, and mixtures of dendrites and columns. The m−D_max and A−D_max relations for the monomers are given in Table 1. Two sets of aggregates with mixed monomer types were created. For the first mixture, the selection of the monomer type is random with the same probability density function for both monomer types (“Mix1”). This would represent a scenario, where dendrites and needles coexist with similar PSD and likelihood of aggregation. For the second mixture, the monomers with D_max < 1 mm are columns, while dendrites are taken for larger monomers (“Mix2”). This choice is motivated by the fact that at temperatures below −20°C, the particle shape is less distinct but mostly described by polycrystals, while at temperatures between −20°C and −10°C, one finds more planar and dendritic crystals (Bailey & Hallett, 2009). Considering a thick cloud, we could assume that the small polycrystal or columnar crystals forming in the upper part of the cloud begin to form the first aggregate and then further grow by collection of larger dendrites at lower layers. Of course, both scenarios are quite ad hoc, and more detailed studies are needed to better understand the real properties of mixed-monomer aggregates. Our mixtures are thus rather intended to qualitatively analyze the differences of mixed monomer aggregates compared to single-monomer type aggregates (as done in another recent study by Dunnavan et al., 2019).

Table 1 Mass-Size ( $m(D_{\max },N_{\text{mono}}=1)=a_{m,1}{D}_{\max }^{b_{m,1}}$) and Projected Area-Size ( $A(D_{\max },N_{\text{mono}}=1)=a_{A,1}{D}_{\max }^{b_{A,1}}$) Relationships for Monomers (N_mono=1) used in the Aggregation Model

Monomer
type	am,1 (kg m ${}^{-b_m}$)	bm,1	aA,1 (m2 m ${}^{-b_A}$)	bA,1
Plate	0.788	2.48	0.631	1.99
Needle	0.005	1.89	0.002	1.42
Dendrite	0.074	2.33	0.142	1.94
Column	0.046	2.07	0.008	1.54

Note. All monomers have a grid resolution of 5 μm. The shapes are predefined in the aggregation model and mostly based on Pruppacher and Klett (1998) (see Fig. 1 in Leinonen and Moisseev, 2015).

The aggregation process strongly depends on the number concentration of particles and their relative terminal velocity differences. In conditions which are less favorable for aggregation (e.g., low number concentration), the particles can grow by depositional growth to relatively large sizes before aggregation becomes the dominant process. It is therefore possible that aggregation involves very different monomer sizes. In order to account for this variability, we vary λ⁻¹ in a large range from 50 μm to 10 mm with 500 different values of λ⁻¹, spaced evenly in the logarithmic space. The monomer distribution is limited to sizes of 100 μm up to 3 mm following Leinonen and Moisseev (2015) in order to be consistent with the typical size range of observed ice particles. Due to this truncation of the inverse exponential distribution, the mean monomer size differs from λ⁻¹ and ranges from 150 μm to 1.48 mm.

The spacing of the monomer number (Table 2) is finer at low N_mono and becomes more coarse at larger numbers. In this way, we can investigate the changes at small monomer numbers with greater detail. In fact, we expect the largest changes in snow properties at the transition from single monomers to aggregates composed of few pristine crystals as shown in earlier studies (Dunnavan et al., 2019; Schmitt & Heymsfield, 2010). The coarser spacing of N_mono also limits computational costs. With our settings, we obtain maximum aggregates sizes ranging from 3 to 5 cm which means that we include also the typically observed large snowflakes during intense snowfall on the ground (Lawson et al., 1998).

Table 2 Grid Resolution, Size Parameter λ⁻¹ of the Monomer Distribution, and Number of Monomers N_mono used to Create the Aggregate Data Set

Resolution	λ−1	Nmono	Dmax of the aggregate
5 μm	50 μm–10 mm	1, 2, 3, …, 10, 20, 30, …, 100	≈1–2 cm
10 μm	50 μm–10 mm	200, 300, …, 1,000	≈3–5 cm

Note. D_max denotes the maximum size range of the generated aggregates in the data set.

In Figure 2, several examples of similar sized aggregates simulated with different combinations of λ⁻¹, N_mono, and monomer types are shown. In total, 105,000 particles were simulated. Apart from the visual differences of shapes and structure, also the particle properties such as mass, area, or terminal velocity show a wide range of values although all aggregates have maximum sizes ranging between 3 and 5 mm.

Figure 2. Examples of simulated aggregates with various size parameters (λ−1), number of monomers Nmono, and monomer types. All aggregates have a comparable maximum size (in the range between 3 and 5 mm). The terminal velocity vterm is calculated using the hydrodynamic model by Böhm (1992; see section 2.2 ).

Hydrodynamic Models

Hydrodynamic models are needed in order to derive the terminal velocity v_term from the particle's mass m, projected area A, and maximum size D_max. The most commonly used hydrodynamic models are Böhm (1992, hereafter B92), Khvorostyanov and Curry (2005, hereafter KC05), and Heymsfield and Westbrook (2010, hereafter HW10). All models are based on particle boundary layer theory and rely on the Best number (X) approach (Abraham, 1970). v_term is calculated via, [Image Omitted. See PDF]where η is the dynamic viscosity, Re the Reynolds number (parameterized as a function of X), and ρ_a is the air density. X is defined as [Image Omitted. See PDF]where C_d is the drag coefficient. The proportionality of X to the particle properties is given by [Image Omitted. See PDF]for B92.

For this study, we decided to use B92 because it best represents the saturation of v_term for our simulated particles at larger aggregate sizes (Figure A2) in accordance with observations (Figure 1). B92 includes an empirical correction of X due to wake turbulence which increases the drag of large particles. X depends on the aspect ratio α, which is larger than one for prolate and smaller than one for oblate particles. For this study, we set α to 1.0, because aggregates with small values of N_mono are not easily classifiable as either prolate or oblate and show in general a large variability of α (Jiang et al., 2019).

To be able to interpret the dependency of v_term on N_mono in section 4.3, we sketch here how v_term scales with D_max in the simplified case of Re≪1 (Stokes drag) and Re≫1 (Newtonian drag). For Re≪1, C_D is approximately proportional to 1/Re. Inserting this approximation and Equations 3 and 4 into Equation 2 yields [Image Omitted. See PDF]

If we approximate m and A by the power laws $m=a_m{D}_{\max }^{b_m}$ and $A=a_A{D}_{\max }^{b_A}$, we can express v_term solely as a function of D_max: [Image Omitted. See PDF]

For Re≫1, C_D is approximately constant. In this case, Equation 3 gives us Re∼X^0.5, and by using again the Equations 2 and 4, we get [Image Omitted. See PDF]

In both extreme cases of Re, v_term increases the faster with size the higher b_m−0.25b_A is, and we expect this also to be in between these cases where Re transitions from Re∼X to Re∼X^0.5. This has certain implications for the dependency of v_term on N_mono (section 4.3).

The differences between the three hydrodynamic models as well as an analysis of the potential impact of changing to different hydrodynamic models is discussed in the Appendix A2.

Comparison of the Simulated Particle Properties to In Situ Observations Mass- and Area-Size Relations

Particle properties, such as m, A, and D_max, are used in hydrodynamic models to calculate v_term (section 2.2). We compare our relations of these particle properties and v_term with frequently used relations that are based on in situ measurements from Locatelli 306 and Hobbs (1974, LH74) and Mitchell (1996, M96). LH74 defined an equivalent diameter that is equal to “the diameter of the smallest circle into which the aggregate as photographed will fit without changing its density.” M96 collected observations as a function of D_max without specifying the exact definition. The definitions of particle size used in these studies are limited by the observation equipment used, and the conversion from one to the other is not trivial. In our simulation study, we can access the full 3D structure of the particles and use the true maximum size (i.e., the maximum distance between any two points of the particle) as size definition.

Except for the aggregates of dendrites, which have a considerably lower density than LH74 aggregates of dendrites, the absolute value of m of the simulated aggregates is similar to the observations, where the same monomer type is available (Figure 3). The slope of the m−D_max relation from this study is comparable to the slope from M96, while LH74 report lower slopes for the aggregates of dendrites. The m−D_max relation of the mixed aggregates (“aggregates of unrimed radiating assemblages of plates, side planes, bullets, and columns,” LH74 mix), however, has a similar slope to the simulated Mix2 aggregates. The mixS3 and sideplane aggregates from M96 are similar to many simulated aggregates (composed of different monomers).

Figure 3. Particle properties of simulated aggregates from this study (green and black), from previous studies (Kajikawa, 1972; Locatelli & Hobbs, 1974; Mitchell, 1996) (M96, LH74, and K72) and measurements of ice particle observed by PIP at the CARE site (brown, see text). (a) m versus Dmax; (b) A versus Dmax; (c) median (and 25th and 75th percentile for PIP CARE) of vterm versus side projected maximum dimension Dmax,side for data from this study and versus the size definition of the respective study (vterm is directly observed in K72 and LH74 and calculated with B92 from the m−Dmax and A−Dmax relations of M96) (d) same as (c) but for larger sizes. Note that K72 observations are for single monomers.

M96 derived A−D_max relations for “assemblages of planar polycrystals in cirrus clouds” (M96 polycrystal in Figure 3) based on observations in a relatively small size range and applied them to other aggregate types. This A−D_max relation is also used in several microphysics schemes (Brdar & Seifert, 2018; Morrison & Milbrandt, 2015). The absolute value of A given in M96 is slightly higher than A of the simulated particles from this study (except for the aggregates of plates). The slope of the A−D_max relations is slightly higher (b_A = 1.88) in M96 observations compared to the relations from this study (1.79 < b_A < 1.88). Observations of aggregates composed of the same monomer types than the one used in these studies are not available.

Terminal Velocity-Size Relations

Observations of v_term versus size have been reported using several different definitions of the diameter (Szyrmer & Zawadzki, 2010). To facilitate a consistent comparison between the observations from the PIP instrument (which are described in Section A0.1) and v_term of the simulated aggregates, we use the same bin sizes as the PIP instrument to derive the median v_term. Moreover, we derive the maximum dimension from a side projection of the modeled particle in the same way as in the observations from the PIP instrument described by (von Lerber et al., 2017) (D_max,side; Figures 3c and 3d). Displayed are the median and the 25th and 75th percentiles of v_term of the detected particles. Bins with fewer than 1,000 particles are excluded from the statistics. Although LH74, M96, and Kajikawa (1972, K72) did not use the same definition as the PIP-CARE data set, fits from this study are also shown in Figures 3c and 3d because they can ease the comparison with other studies.

At small sizes (D_max< 1 mm), v_term of the simulated aggregates of dendrites is close to v_term of the monomers from Kajikawa (1972, K72, Figure 3c). The plate monomers in K72 are reported with a similar v_term as the aggregates of plates, needles, and Mix1 (which all have similar values). Note that v_term of plates and dendrites from K72 and v_term of all aggregates simulated in this study (except for the aggregates of columns and “Mix2”) are considerably smaller than v_term of the aggregates from the PIP-CARE data set and LH74. The observations from LH74 are within the 25th and 75th percentile of the PIP-CARE data set. The median of v_term of the simulated aggregates of this study increases faster with size compared to the in situ observations at sizes of several millimeters (Figure 3d). Only v_term of the mixture of small columns and large dendrites (“Mix2”) has a comparably low slope. Potential reasons for this mismatch are limitations of the observations at these sizes (Brandes et al., 2008), turbulence affecting the observations (Garrett & Yuter, 2014), missing processes in the aggregation model (e.g., depositional growth on aggregates), imperfect parameterizations in the hydrodynamic model, or the dominance of monomer type mixtures in the aggregates.

Figures 3c and 3d also show v_term calculated with B92 and the m−D_max and A−D_max relations from M96 (which did not measure v_term directly). The simulated slope of v_term from M96 observed aggregates is similar to the one simulated in this study, while the absolute value is slightly higher.

At sizes larger than about 5 mm, the simulated and the observed v_term reach a saturation value close to 1 m s⁻¹. The median of v_term of most simulated aggregates lies within the 25th and 75th percentile in the subcentimeter range, except the aggregates with the most extreme density (aggregate of dendrites and aggregates of columns). Thus, based on this comparison, these aggregates can be considered most representative for many aggregates found in the atmosphere.

Parameterization of Particle Properties

The relationships between hydrometeor properties such as mass, size, projected area, and velocity are key components in any ice microphysics scheme, and they strongly influence various microphysical processes (e.g., sedimentation, depositional growth, aggregation, or riming). Different microphysics schemes require a more or less simplified parameterization of particle properties. To address these different needs, we derive in this section fits for m and A as a function of D_max and N_mono that can be used in microphysics schemes, which can predict m and N_mono given a certain D_max (section 4.2). Of course, most bulk schemes require less detailed fits, and hence, we also derive fits of m, A, and v_term as a function of D_max or the mass-equivalent diameter D_eq. This also allows us to assess the potential error of the less detailed fits (section 4.5), while their impact on modeled processes is studied later in section 5.

Fitting Approach for Monomer Number-Dependent Particle Properties

The particle properties of the monomers are defined a priori in the aggregation model and based on well-established observations. In contrast, the aggregate properties are determined by the aggregation process and change with increasing N_mono. As we are particularly interested in quantifying how key particle properties of aggregates differ from the properties of the same-sized monomers, we normalize the aggregate properties by the property of a monomer with the same D_max [Image Omitted. See PDF]

p represents the particle properties (mass or area), p(D_max,N_mono=1) is the property of single monomers (given in Table 1), and f_p is the normalizing function. A normalizing function which is larger (smaller) than 1 indicates that the aggregate properties are larger (smaller) than its composing monomer with the same size (Figure 4).

Figure 4. Schematic illustration of how compactness of aggregates can cause them to be heavier or lighter compared to a monomer of the same size. For simplicity a monodisperse monomer size distribution with monomer sizes of Dmax = 0.2 mm is used. The red line indicates the maximum theoretical compactness of mass of an ice sphere. The black lines shows the m−Dmax relation of the monomer (plate). The green line represents the m−Dmax relation of the least compact configuration of the plate monomers in an aggregate by aligning the plates along their maximum dimension. Particles have lower mass (fm [less than] 1) in the green shaded area and larger mass (fm > 1) in the red shaded region compared to an equal-size plate.

To fit f_p to various monomer types, we parameterize f_p by a power law and express the coefficients by rational functions to fit the dependency on N_mono similar to the approach presented in Frick et al. (2013). [Image Omitted. See PDF]

The coefficients of f_p for all monomer types can be found in Table 3. Note that we excluded the mixture of monomer types from the monomer-dependent analysis because our normalization approach cannot be applied to monomer mixtures.

Table 3 Coefficients in the Normalizing Functions f_m and f_A (Notation as in Equation 9) for Different Monomer Types

Monomer type	af,m	a′f,m	bf,m	b′f,m	af, A	a′f,A	bf, A	b′f,m
Plate	−0.673	0.364	−0.092	0.091	−0.473	0.322	−0.021	−0.166
Needle	0.162	−0.008	0.018	0.102	0.349	0.005	0.060	0.013
Dendrite	−0.288	0.215	−0.042	−0.056	−0.100	0.131	−0.019	−0.059
Column	0.079	−0.006	0.033	0.086	0.273	0.025	0.058	0.034

Dependence of Aggregate Mass and Area on Monomer Number

Motivated by the common classification of unrimed ice hydrometeors in cloud ice and snow in many bulk schemes, we will investigate in this section how mass and area change when building up an aggregate with an increasing number of monomers. In particular, we want to explore whether the properties change smoothly with monomer number or whether they show any sharp transition at certain monomer numbers.

When we compare the mass of an aggregate with the mass of its monomer of the same size, we find in some conditions the aggregate to be heavier or lighter than the monomer. The relevant mechanisms which explain this behavior are illustrated in Figure 4 for aggregates of plates. Note that we assume for simplicity a monodisperse monomer distribution in Figure 4. When we consider pure depositional growth, we obtain a specific m−D_max relation for each monomer type (Table 1; black line in Figure 4). One extreme aggregation scenario, which leads to the maximal size of an aggregate with a given number of monomers (which in this simplified case of a monodisperse distribution also determines its mass), would be if we assume that all monomers align along their maximum dimension. Clearly, the resulting aggregate would have a smaller m than a monomer of the same size. Of course, this maximal elongated assemblage of monomers is rather unlikely, and thus, the aggregate will have a more compact structure. If we imagine rearranging the monomers inside the aggregate in a progressively more packed configuration (indicated by the horizontal arrow in Figure 4), we might be able to reach the point where the size of the aggregate equals the one of the equal-mass monomer. At this point, it might be even possible to pack the monomers in a way that their size is smaller than an equal-mass monomer. A simple example of such an extreme packing would be to stack a number of plates on top of each other, that is, along their smallest axis. Whether an aggregate can be smaller than an equal-mass monomer is of course also dependent on how close the monomer m−D_max relation is to the theoretical maximum packing of an equal-mass sphere.

The dependency of A on N_mono can be understood analogously. Also, for A, the maximal elongated assemblage of the monomers leads to a lower A of the aggregate compared to the monomer of the same size, but in reality, the monomers will assemble in a more compact way. In addition, we have to consider that A is not simply additive as it is the case for m. Overlap (in the horizontally projected plane) and nonhorizontal alignment of the constituting monomers lead to a smaller A than the sum of A of the constituting monomers. Based on these simplified considerations, it becomes clear that the dependency of m and A on N_mono is determined by the exponent of the monomer power laws and the overall “compactness” of the aggregates.

When considering the monomer dependence of all simulated aggregates, we find the most different behavior for plate and needle aggregates. For plate aggregates, m and A steadily decrease for a given D_max with an increasing number of monomers (Figures 5b and 5d). From the principal considerations discussed in Figure 4, this behavior can be well understood. The plate monomers have the largest exponent (b_m,1 = 2.48) of all monomers (Table 1), while the monomers itself show relatively loose connections within the aggregate (Figures 2a–2c). Interestingly, the aggregate mass for very small N_mono can be slightly larger than the equal-size monomer, while A is immediately decreasing for N_mono > 1. This effect can be easily understood when considering, for example, two plates that connect in a 90° angle of their major axes.

Figure 5. (a and c) m and A of the simulated plate aggregates as a function of Dmax. (b and d) The normalizing functions fm and fA (defined in Equation 8) quantify the deviation of the aggregates' m or A from a monomer with same Dmax. The dots indicate the properties of individual particles with the color showing Nmono. Lines indicate m and A for constant Nmono as a result of the monomer number dependent fits and for all aggregates (Nmono > 1).

An opposite behavior is found for needle aggregates (Figures 6b and 6d). With increasing N_mono, both m and A of the aggregates become larger than the equal-size monomers. In contrast to plates, the needle monomers have the lowest exponents for the m and A power laws (Table 1). The aggregates of the more one-dimensional needles also show a more compact packing.

Figure 6. Same as Figure 5 but for aggregates of needles.

The deviation of the particle properties of the individual simulated particles from the fit is characterized by the mean absolute error (Table A2), which is smallest for plates (0.1190 for f_m and 0.0816 for f_A) and largest for needles (0.3737 for f_m and 0.3926 for f_A). The mean absolute error also shows that the monomer number-dependent fit is superior to the more simple power law fit (section 4.4) when there is a substantial dependence of the particle property on N_mono.

Dendrite and column aggregates have been analyzed similarly (according figures can be found in Supporting Information S1). The dendrites are similar to plates, while the columns are similar to needles. However, for all aggregate types, we find on average a relatively smooth transition of m and A when changing from single monomers to aggregates. For these two particle properties, we are unable to identify a “jump” due to the onset of aggregation. The next sections will show whether this behavior will change when deriving terminal velocity from m and A.

Dependence of Terminal Velocity on Monomer Number

The terminal velocity for all aggregates was calculated with the hydrodynamic model of B92 (section 2.2). In Figure 7a, v_term is shown as a function of D_max for plate aggregates. Note that the fits have been derived by applying B92 to the m−D_max and A−D_max fits (Table 3) rather than fitting them directly to the cloud of individual v_term. In this way, we are consistent with the way how v_term relations are usually connected to m−D_max in bulk schemes. The terminal velocity of plate aggregates steadily decreases with increasing N_mono. This dependency is much less pronounced at small D_max as compared to the largest sizes. However, it should be noted that the fits for very small monomer numbers are probably unrealistic for large D_max as we do not expect aggregates of centimeter sizes to be composed of only a few large plates. In fact, the here used geometrical relations for the plate monomers are only valid up to a size of 3 mm (Pruppacher & Klett, 1998).

Figure 7. vterm versus Dmax for the simulated aggregates of plates and needles. The dots indicate the properties of individual particles with the color showing Nmono. Lines indicate vterm for constant Nmono as a result of the monomer number dependent fits and for all aggregates (Nmono > 1). Note that the fits have been derived by applying B92 to the m−Dmax and A−Dmax (Table 3) fits rather than fitting them directly to the cloud of individual vterm.

We find a similar decreasing v_term with increasing N_mono for dendrites (see supporting information S1). As we might expect from the different change of m and A with N_mono seen in Figure 7a, also the behavior of v_term with increasing N_mono is different for needles (Figure 7). Needle aggregates seem to fall slightly faster when their monomer number increases. Interestingly, all aggregates reveal a very low dependence of v_term on monomer number at small sizes which is in contrast to assumptions in some microphysics schemes that distinguish between monomers and aggregates (e.g., Seifert & Beheng, 2006; Tsai & Chen, 2020). Besides, all aggregates reveal a saturation of v_term at large (centimeter) sizes which is in good agreement with observations (Figure 1). However, the absolute value of the saturation v_term ranges from 0.8 to 1.6 m s⁻¹ depending on the monomer type.

Because v_term of monomers and aggregates is converging toward the same value at small sizes (Figure 7), we can use the previously derived scaling relation (Equations 6 and 7) to relate the dependency of v_term on N_mono to the exponents b_m and b_A of the monomers (b_m,1 and b_A,1) and aggregates (b_m,agg and b_A,agg) in the m−D_max relation. Starting from a similar value of v_term at small sizes, v_term of an average aggregate increases slower than v_term of a monomer if s_monodep=b_m,agg−b_m,1−0.25(b_A,agg−b_A,1)<0 (cf. Equations 6 and 7). As a result, at larger sizes, v_term of the aggregate is lower than v_term of the monomer. In an analog way, v_term of an aggregate is larger than v_term of the monomer if s_monodep>0. As b_m,agg and b_A,agg are similar for all aggregates (Table 4), the sign of v_term with increasing N_mono depends mainly on b_m,1 and b_A,1. For plates and needles, s_monodep equals −0.21 and 0.12, respectively.

Table 4 Mass-Size ( $m(D_{\max })=a_{m,\text{agg}}{D}_{\max }^{b_{m,\text{agg}}}$) and Projected Area-Size ( $A(D_{\max })=a_{A,\text{agg}}{D}_{\max }^{b_{A,\text{agg}}}$) Relationships for Aggregates (N_mono > 1) in the Aggregate Model

Monomer
type	am,agg (kg m ${}^{-b_m}$)	bm,agg	aA,agg (m2 m ${}^{-b_A}$)	bA,agg
Plate	0.076	2.22	0.083	1.79
Needle	0.028	2.11	0.045	1.79
Dendrite	0.027	2.22	0.090	1.88
Column	0.074	2.15	0.060	1.79
Mix1	0.045	2.16	0.070	1.83
Mix2	0.017	1.94	0.066	1.79

How the particle properties change with increasing N_mono as well as the absolute values of calculated v_term depends on the choice of the hydrodynamic model. Finding the optimal formulation of hydrodynamic models for ice and snow particles is still an active field of research (Nettesheim & Wang, 2018; Westbrook & Sephton, 2017) and outside the scope of this study. In Appendix A2, we tested the sensitivity of the results to the choice of the hydrodynamic model for plate aggregates. HW10 seems to yield overall similar results to B92 except for the saturation at large diameters. Due to the absence of the turbulence correction in HW10, v_term increases also at large diameters. For KC05, the monomer dependence is much weaker. However, all hydrodynamic models show an overall small monomer dependence at small particle sizes.

It has also been observed (e.g., Garrett & Yuter, 2014) that tumbling of particles caused for example by turbulence might decrease the effective projected area and therefore increase v_term. We also tested the sensitivity of our results to different degrees of tumbling (section A0.2.2). As expected, the effect of tumbling is largest for single crystals (due to their more extreme aspect ratio) but strongly decreases for aggregates. Certainly, for aggregates, the choice of the hydrodynamic model has a larger effect of v_term than different assumptions on particle tumbling.

Mean Particle Properties of Monomers and Aggregates of Different Monomer Types

The relatively continuous change of particle properties with N_mono found in the last section justifies a simplified fit, which is also necessary for implementing the results into common bulk microphysics schemes. These schemes often only contain two classes for unrimed ice particles, usually denoted as cloud ice (monomers) and snow (aggregates).

Figures 8a and 8b show the derived m−D_max relations for single monomers (N_mono=1) and the derived v_term based on the m−D_max and A−D_max relations summarized in Table 1. Similar fits of m and v_term to aggregates of any monomer number larger than 1 are shown in Figures 8c and 8d; the fit coefficients can be found in Table 4.

Figure 8. Particle (a and c) m and (b and d) vterm as a function of Dmax calculated with B92 using the derived m/A−Dmax relations (Tables 1 and 4). Particles are separated into (a and b) single monomers and (c and d) aggregates composed of various monomer types (see legend).

The m−D_max relations for monomers show a larger spread especially for larger sizes as compared to the aggregates. This is expected considering that the exponents for monomers range between 1.89 and 2.48 (Table 1), while the exponents for aggregates are between 1.95 and 2.22 (Table 4). The values for aggregates agree well with theoretical aggregation studies (Westbrook et al., 2004b) as well as in situ observations (section 3.1). Despite the similar exponent, the effective density of the aggregates varies considerably (compare m at a given size in Figure 8c), which is in agreement with previous studies (Dunnavan et al., 2019; Hashino & Tripoli, 2011), even though their approaches to simulate aggregates are very different from the approach used in this study. Aggregates of columns exhibit the highest density, while aggregates of dendrites show the lowest density.

The differences in the m−D_max relation are linked to the resulting v_term−D_max relation (Figures 8c and 8d). At D_max = 5 mm, the v_term of different monomers spread nearly 1 m s⁻¹. The differences are in general smaller for aggregates. Interestingly, most aggregate types reveal very similar v_term. The main exceptions are dendrite aggregates with the slowest and column aggregates with the fastest v_term. v_term of the Mix2 aggregates increases slower with increasing D_max compared to the other aggregates.

Similar to the previous monomer number dependent fits, also the “two-category” fits show similar v_term at small sizes. The monomer type appears to have in general a much larger impact on v_term than the classification into certain N_mono regimes.

Power-Law and Atlas-Type Fits for Terminal Velocity

Power-law fits for m, A, and v_term are commonly used in bulk schemes. Especially for v_term, the power law introduces inconsistencies with observations because a saturation value for v_term as observed for raindrops or snowflakes cannot be represented. Instead of using standard power laws in the form, [Image Omitted. See PDF]and the two fit parameters $a_{v{D}_{\max }}$ and $b_{v_{D_{\max }}}$, Atlas et al. (1973) proposed a three-parameter ( ${\alpha }_{D_{\text{eq}}}$, ${\beta }_{D_{\text{eq}}}$, and ${\gamma }_{D_{\text{eq}}}$) formulation [Image Omitted. See PDF]

Formulating this “Atlas-type” fit with the mass equivalent diameter D_eq instead of D_max has been found to provide optimal fit quality for snow aggregates (Seifert et al., 2014). For small (large) values of D_eq, v_term approaches ${\alpha }_{D_{\text{eq}}}-{\beta }_{D_{\text{eq}}}$ ( ${\alpha }_{D_{\text{eq}}}$). With increasing values of γ, the transition from small to larger values of v_term is shifted toward larger values of D_eq. Approximations for bulk collision rates based on Atlas-type fits can be found in Seifert et al. (2014) which makes them usable in bulk microphysics schemes without the necessity of look-up tables.

Power-law and Atlas-type relations have been applied to the various aggregates and the fit coefficients are summarized in Table 5. For the fitting, we did not use v_term of the individual particles but directly applied to fit to v_term derived with B92 and the existing m−D_max and A−D_max relations.

Table 5 Derived Coefficients of the Power-Law and Atlas-Type Fits (Equations 10 and  11) for Monomers and Aggregates of Different Monomer Types

[Table omitted. See PDF]

In Figure 9, the different fits are compared for plate monomers and their aggregates. Note that the saturation region (D_max > 1 cm) has been excluded for the power-law fits. It can be seen in Figure 9b that the Atlas-type fit is very close to the theoretical line calculated with B92 and the m−D_max and A−D_max relations. The power-law fits (Figure 9a) provide only a close fit to the theoretical values at the smaller size range. Between 300 μm and 4 mm, they cause a slight underestimation, while at larger sizes, they increasingly overestimate v_term. Similar fits have been derived for all aggregate types (Table 5, figures for other monomer types similar to Figure 9 can be found in the supporting information S1).

Figure 9. vterm of individual plate aggregates (gray scale, a–c) and vterm derived with B92 and the m/A−Dmax of plate monomers (Table 1, solid blue line in a and b) and aggregates (Table 4, solid green line in a and b). Power-law (dashed-dotted, a) and Atlas-type fits (dashed-dotted, b) have been applied to the directly calculated vterm (solid lines) rather than the individual points. (c) vterm used in microphysics schemes (same as in Figure 1b).

When we compare the calculated v_term with some widely used microphysics schemes (Figure 9c), we find most schemes to overestimate v_term at small sizes (except of the cloud ice category in Morrison et al., 2005). The absolute values for v_term at small sizes are strongly dependent on monomer type, and hence, additional constraints should be provided by additional observations. However, the aggregation model shows independent on the monomer type that at submillimeter sizes, there should be no strong “jump” in v_term between ice particles and small aggregates. Also, in the centimeter-size range, models using a power-law formulation are strongly overestimating v_term for all aggregate types.

Application and Sensitivity Tests in the Lagrangian Particle Model McSnow

In this section, we will test the possible impact of implementing particle properties with different amount of complexity (monomer number dependence) or different fitting functions (power law vs. Atlas type) on the simulation of sedimentation, aggregation, and depositional growth. For this, we use a one-dimensional setup of the Lagrangian particle model McSnow (Brdar & Seifert, 2018), which provides the flexibility to implement the different particle property formulations.

For simplicity, only sedimentation, depositional growth, and aggregation are considered in our simulations. Aggregation is calculated with a Monte-Carlo algorithm following Shima et al. (2009), and the sticking efficiency of Connolly et al. (2012) is used. McSnow is based on the Lagrangian super particle approach (Shima et al., 2009), which allows deriving not only the particle mass and its multiplicity X_mult but also predicts the number of monomers the particle is composed of. This information is key to test the N_mono dependent particle relations. The thermodynamic profiles and the overall setup is similar to previous simulation studies with McSnow in Brdar and Seifert (2018) and Seifert et al. (2019). Particles are initialized at the upper boundary of the 5-km thick domain with a mass flux of F_m=2·10⁻⁵ kg s⁻¹ and a mean mass of the particle size distribution of m_mean= 2·10⁻¹⁰ kg. The initial ice particles follow a generalized gamma distribution of particle mass with a shape parameter of 0 and a dispersion parameter of 1/3 (following Eq. 9 in Khain et al., 2015). The temperature decreases linearly from 273.1 K at z = 0 km to 242.2 K at z = 5 km. The supersaturation over ice is held constant at 5% with respect to ice in the whole column and is not consumed by the growth of the particle. The simulations are performed with 250 vertical levels, which result in a vertical resolution of 20 m. The model time step is set to 5 s, and the initial multiplicity is chosen to be 1,000. The simulations are run for 10 hr, from which the last 5 hr are averaged in 10-min intervals to reduce noise in the analyzed profiles.

In the following, we will focus the comparison on particle number flux (F_N), mass flux (F_M), and mean mass m_mean (which is the ratio between the integrated mass density q_m and the integrated number density q_N).

In the first simulation experiment shown in Figure 10, we include particle properties for which the full N_mono dependence is taken into account (Table 6). In this setup, we call hereafter the control simulation (“CTRL”). Profiles are separated into single monomers (N_mono = 1) and aggregates (N_mono > 1) to better distinguish the effects on what we might define as “cloud ice“ and “snow” category in a bulk scheme. This separation might be important considering that there can be cases of weak aggregation. With weak aggregation, most of the particles will remain monomers, and thus, it is especially important to match profiles of these particles accurately.

Figure 10. Idealized McSnow simulation using the Nmono-dependent fit for plates (“monodep”; Table 3), the separation between Nmono=1 and Nmono > 1 (“binary”; Tables 1 and 4) and single relation (the one fitted to all aggregates) for all Nmono (“constant”; Table 4) for plates. For each individual super particle, B92 is used directly to calculate vterm. Shown are height profiles of (a and b) number flux FN, (c and d) mass flux Fm, and (e and f) mean mass mmean. The particles are categorized into Nmono = 1 (left) and Nmono > 1 (right).

Table 6 Settings of the McSnow Control (CTRL) and Sensitivity Runs

[Table omitted. See PDF]

Note. The second column specifies the monomer type from which the m−D_max and A−D_max (and v_term−D_max for the Atlas and power law run) fit is taken. The third column denotes how the N_mono dependency is represented. f_p(N_mono,D_max) is the normalizing function with full N_mono dependence (section 4.1), f_p(N_mono=1;N_mono>1,D_max) denotes only a binary seperation in N_mono=1 and N_mono>1, and f_p=f(N_mono>1,D_max) indicates that the fit for all aggregates N_mono>1 is taken for all particles (section 4.4). The fourth column indicates whether v_term is calculated using B92 or with a parameterized relation of v_term−D_max (section 4.5). The fifth column shows the precipitation rate (F_m(z = 0 m)) and in brackets its deviation from the “CTRL” run. The last column lists the mean mass m_mean at the surface and the ratio of m_mean between the simulation and its “CTRL” run (in brackets).

In general, aggregation decreases the number flux (F_N), while the increase in the mass flux (F_m) is due to depositional growth. The mass flux of aggregates increases also due to conversion from monomers to aggregates by aggregation. The combination of both processes is causing m_mean to continuously increase toward the surface. Aggregation rates in McSnow are proportional to the gravitational collection kernel (Eq. 21 in Brdar & Seifert, 2018). Thus, the probability of collision for two particles is high if they have strongly different v_term and if the sum of their cross-sectional areas is large. F_N of the monomers (N_mono = 1) decreases monotonically with decreasing height because the monomers are converted into aggregates (N_mono > 1) by the aggregation process and there is no source of monomers like nucleation considered (Figure 10a). This decrease of F_N (and increase of m_mean) is especially strong at heights between 2 and 3 km. This region of enhanced aggregation is found at heights where the temperature is close to −15°C where the sticking efficiency has a local maximum. As a result of the conversion of monomers to aggregates, F_N of the aggregates increases at heights higher than about 3.5 km (Figure 10b). At lower heights, the number of aggregate-aggregate collisions outweighs the number of monomer-monomer collisions, and thus, F_N of the aggregates decreases.

The signature of the conversion from monomers to aggregates is also seen in F_m of the monomers (Figure 10c). Especially in the region of enhanced aggregation, this leads to a strong decrease of F_m. In the heights above this region, depositional growth outweighs the loss of mass of the monomers to the aggregates, and thus, there is an increase of F_m with decreasing height. F_m of the aggregates increases monotonously due to both depositional growth of the aggregates and conversion from monomers to aggregates (Figure 10d). In this setup, the change of F_m with height is governed by v_term and q_N at a given height. For example, a combination of low v_term and high q_N at upper layers leads to a large increase in F_m. Simply speaking, a large number of slow falling ice crystals can grow efficiently by deposition which increases F_m.

Sensitivity to Representation of Monomer Number Dependency

The “CTRL” simulation is now compared to simulations with a binary separation into single-monomer particles and aggregates of any monomer number larger than 1 (binary). An additional simulation is performed with no monomer number dependence (constant). Here, the particle properties that were fitted to the mean of all aggregates are used for all particles. All simulations are done for plate and needle monomers and aggregates because we found the monomer dependence to be most pronounced for these monomer types. For the other monomer types, the effect of N_mono can be expected to be smaller.

The most apparent difference between the simulations with different representations of the N_mono dependencies for plate monomers and aggregates of plates is the faster decrease of F_N and F_m and faster increase of m_mean of the monomers (N_mono = 1) in the “constant” simulation (Figure 10). A slightly faster decrease of F_N (faster increase of m_mean) for aggregates (N_mono > 1) with decreasing height can be seen for both the “binary“ and the “CTRL” simulation. However, all of the simulations show very similar profiles.

Figure 11 shows the same experiment as Figure 10 but using the parameterizations for needles instead of plates. Also, for needles, the most remarkable difference between the simulations is the difference between the “constant“ and the “CTRL” run (Figures 11a and 11e). Also, aggregate-aggregate collections are less effective in the “CTRL” run (Figures 11b and 11f). Note that all monomers have been depleted by aggregation at a height of about 1,000 m, and thereby, m_mean is not defined below.

Figure 11. Same as Figure 10 but for needle monomers and aggregates.

Overall, the differences of m_mean at the ground of the total ice particle population are small (factor of 1.2 and 1.4 higher m_mean for the “binary“ and “constant” simulation for plates and factor of 0.8 lower m_mean for the “binary“ and “constant” simulation for needles, Table 6).

Also, the differences in the precipitation rates (F_m) are small (less than 5%; see Table 6). These small differences are due to the small difference of the absolute value of v_term at small sizes (Figure 7) and q_N at upper heights, which lead to a similar mass uptake (Figure 10). However, the precipitation rate between the “Plate CTRL” simulation and the “Needle CTRL” simulation is relatively large (Table 6), which might be due to the strongly different v_term of the monomers.

The N_mono dependency is even weaker for aggregates composed of other monomer types (sections 4.2 and 4.3). In summary, the simulation experiments with different monomer dependency indicate that a binary separation between single monomers and aggregates performs similarly well as relations which take into account a more detailed monomer dependency. Some but still small differences are found if no monomer dependency is taken into account, that is, a single ice class for all monomer numbers is assumed. In our simulation, particles with low N_mono are only prevalent at cold temperatures, where aggregation is less important due to the low sticking efficiency. Additional simulations (shown in the supporting information S1) with lower F_m and therefore weaker aggregation show that the “binary” simulations stay very close to the “CTRL” simulation, while the “monodep” simulations show considerably larger deviations. Hence, we find that the classical separation in cloud ice (monomers) and snow (aggregates) is sufficient for the aspects of monomer number-dependent particle properties.

Sensitivity to the Parameterization of Terminal Velocity

In this section, we test the sensitivity of the simulations to different implementations of the v_term−D_max formulation. In Figure 12, v_term of plate monomers and aggregates is parameterized either as power-law or Atlas-type fit.

Figure 12. Idealized McSnow simulation using m−Dmax and A−Dmax for plate monomers and aggregates of plates (see Tables 1 and 4) and power law (without (“powerlaw”) and with imposing an upper limit on vterm (“powerlawlimit”), which is consistent with the saturation value of the Atlas-type relation) and Atlas-type vterm−Dmax relations for plate monomers and aggregates of plates (see Table 5). Overlayed is the CTRL/monodep simulation in gray (see also Figure 10). Shown are height profiles of (a and b) number flux FN, (c and d) mass flux Fm, and (e and f) mean mass mmean. The particles are categorized into Nmono = 1 (left) and Nmono > 1 (right).

As we saw in Figure 9, the power-law and Atlas-type fits match very closely at small particle sizes. This explains the very close matching of the three simulations in the upper part of the simulated profiles (Figure 12) where the PSD is dominated by small particles. As soon as the aggregation becomes stronger (below ca. 3 km), F_N in the simulations using the power law (Figure 12b) is much lower than for Atlas-type. The decreasing number flux of aggregates with lower height (Figure 12b) also indicates that especially the self-collection of aggregates is stronger than for Atlas type. In the same height region, the mean mass of the aggregates (Figure 12f) is strongly increased for the power law (factor of 5). Instead of using an Atlas-type fit to consider the saturation of the terminal velocity, one can also think of imposing an upper limit on v_term in the power law relation. In the “Powerlawlimit” simulation, we chose the saturation value of the Atlas-type fit ( ${\alpha }_{D_{\text{eq}}}$) as an upper limit. This limit does not only affect the sedimentation but also all processes which depend on v_term (e.g., aggregation). In this way, the overestimation of m_mean, caused by an unlimited increase of v_term, can indeed be prevented, but the height profile of F_N and m_mean is not as well matched as with the Atlas-type approximation. As expected, the continuously increasing v_term in the unlimited power law leads to much stronger growth of aggregates as compared to relations which include the saturation velocity at large particle sizes. This is an interesting finding and could be one reason for the overestimation of radar reflectivities found at lower layers in ice clouds simulated with the Seifert-Beheng scheme (Heinze et al., 2017).

Although m_mean of the aggregates is much larger for the power law, the difference to the Atlas type in precipitation rates is very small (smaller than 5%; Figure 12d and Table 6). Note that in more realistic cases, as for example in presence of stronger sublimation layers, the difference in m_mean can induce larger differences in the precipitation rate because larger particles can fall through a thicker layer of subsaturated air before they sublimate completely.

Summary and Conclusions

In this study, we generated a large ensemble of ice aggregates (ca. 105,000 particles) using an aggregation model and hydrodynamic theory to study the change of particle properties such as mass m, projected area A, and terminal velocity v_term as a function of monomer number N_mono and size. The aggregates were composed of various monomers types (plates, dendrites, needles, and columns), monomer sizes, and monomer numbers. In order to test the impact of habit mixtures, we also included in our analysis two different mixtures of dendrites and columns. The choice of mixing specifically dendrites and columns was motivated by in situ observations of the composing monomers in large aggregates sampled on the ground (Lawson et al., 1998).

When comparing our aggregate properties with in situ observations, we find m and A to be very similar to the results presented in Mitchell (1996) but the slope of our m−D_max relations is larger than the slope given in Locatelli and Hobbs (1974). A better agreement with Locatelli and Hobbs (1974) and also with theoretical considerations in Westbrook et al. (2004b) is reached for mixtures of small columns and larger dendrites (Mix2). Interestingly, this monomer mixture also achieves the best agreement with observed v_term−D_max relations. Considering the large spread in the observations (Figure 3), we can overall conclude that our aggregate ensemble matches the observed range of variability and does not show any substantial bias.

Our synthetic aggregate ensemble allowed us to investigate the transition of particle properties from single crystals to aggregates with increasing number of monomers in a level of detail which is currently unavailable from in situ observations. For m and A as a function of size, we find the relations to change rather smoothly with increasing N_mono. The differences introduced by the choice of the monomer type are found to be overall larger than due to the number of monomers. We find the exponents in the A−D_max and m−D_max relations of the monomers to be closely connected to the resulting change with N_mono.

The derived A−D_max and m−D_max relations including the monomer type and number dependence were then used to calculate v_term−D_max relations. Again, we find a rather smooth transition from single crystals to aggregates rather than a “jump” as found in several microphysics schemes (Figure 1b). For small sizes below a few millimeters, our results suggests that the “ice” and “snow” category of microphysics schemes should have similar properties. At larger sizes, the aggregates v_term are found to deviate more from the monomers. Again, the monomer type is found to have a larger impact than the monomer number. Aggregates of plates tend to be faster, while aggregates of needles are slower than the equal-size monomer. In accordance to in situ observations, our simulations reveal for all aggregate types a saturation of v_term at centimeter sizes. However, the saturation value varies for the different aggregate types from 0.8 to 1.6 m s⁻¹.

In order to potentially implement our results in microphysics schemes, we derived two-parameter power-law fits and three-parameter Atlas-type fits for single monomers (N_mono = 1) and aggregates (N_mono > 1) representing the commonly used ice and snow classes in models. The new power-law fits match the small sizes well and avoid unrealistic “jumps” found in current schemes. However, the power laws are unable to represent the saturation of v_term at larger sizes. The Atlas-type fits are found to match the entire size range well and should thus be considered to be implemented in ice microphysics schemes as they do not substantially increase the computational costs while strongly improving the realism of the relations.

We finally tested the impact of implementing monomer dependence, habit type, and velocity fitting method on idealized aggregation simulations. For this, we used a new 1D Lagrangian Monte Carlo model which allowed us to implement the derived relations with different degree of complexity. The simulations experiments revealed that there is only a very small impact of using a relation of only two monomer categories (single particle and aggregate) as compared to a continuous monomer number dependence. A single category which does not take any monomer number into account shows slightly larger deviations, but the variability due to monomer type is in general larger than the impact of monomer number.

In a second simulation experiment, we investigated the impact of using a power law or an Atlas-type fit for v_term. The simulations show very small differences in the upper part of the cloud where the profiles are dominated by small particles which are fitted similarly well with the two relations. Once aggregation becomes more dominant and the spread of particles sizes shifts to larger sizes, the simulations using the power law lead to a much stronger aggregation and in particular stronger self-aggregation of particles as compared to the Atlas-type fit. The impact of the widely used power-law relations for v_term should thus be further studied for bulk schemes as it seems to be likely that they might cause an overestimation of aggregation and snow particle sizes.

We also shortly investigated the sensitivity of our derived relations to particle tumbling and the choice of the hydrodynamic theory. While tumbling can significantly affect the properties of single monomers, it has a surprisingly small effect on our results for the aggregates. The choice of the hydrodynamic theory is a larger source of uncertainty which should be further investigated in future studies. It seems to be important in the future to better constraint the composition of aggregates regarding the monomer type. This question could be approached by improved in situ techniques but also with detailed models that allow to predict the particle habit such as presented in (e.g., Jensen et al., 2017; Shima et al., 2019; Woods et al., 2007).

A Appendix A.1. Video Disdrometer Data Set

The terminal velocity v_term of the simulated aggregates from this study is compared to recent observations of falling ice particle properties and frequently used literature in section 3.2. These surface observations are from the Centre for Atmospheric Research Experiments (CARE), Canada. It is a research facility of the Air Quality Research Branch of the Meteorological Service of Canada, located about 80 km north of Toronto, Ontario (lat = 44 13′58^′′N, lon = 79 46′53^′′W). The instrumentation includes a video disdrometer, PIP, precipitation weighing gauge, and meteorological measurements of, for example, wind velocity.

More detail about PIP can be found in von Lerber et al. (2017) and references therein. The particle sizes are recorded in the range of 0.2–26 mm (disk equivalent diameter) with a resolution of 0.2 mm, which is converted to the side projected D_max. In practice, the minimum reliable size with measurement of v_term is approximately 0.5 mm. Observations of the side projected maximum dimension D_max,side can be conducted from the grayscale video images. The velocity v_term is obtained from the observations of the consecutive frames. The observed v_term utilized in the Figures 1a, and 3c, and 3d are separated from the whole data set by limiting the exponent of the “5-min m−D relation” between 1.7 and 2.2 to exclude rimed particles (von Lerber et al., 2017). To apply this m−D threshold, the mass of the single particle and D_max has to be retrieved. The mass estimate of a single particle is calculated from the observed v_term, corrected D_max, and area ratio using different parametrizations of the hydrodynamic theory (Bohm, 1989; Heymsfield & Westbrook, 2010; Mitchell & Heymsfield, 2005). For each snowfall event, each of these parameterizations is calculated and the one which minimizes the error in the estimate of the liquid water equivalent precipitation with respect to the precipitation gauge is selected for that event. This procedure and the related uncertainties are described more in detail in von Lerber et al. (2017). Additionally, observations during 5-min intervals, where the mean horizontal wind speed exceeds 4 m s⁻¹, are excluded to reduce turbulence effects (similar to Brandes et al., 2008).

After applying these filters, the data set, which covers the winters from 2014 to 2017 with 48 snowfall events, contains about 4.3 million ice particles. It should be noted that PIP is providing a measurement of the ensemble of particles and no particle by particle-based classification is performed. Hence, the measurement volume includes mixtures of different habits.

Table A1 Proportionality of the Best Number X on the Particle Properties (Mass m and Projected Area A), Scaling Relations of the v_term−D_max Relations and s_monodep in Different Hydrodynamic Models (Böhm, 1992, B92; Heymsfield & Westbrook, 2010, HW10; Khvorostyanov & Curry, 2005, KC05)

	B92	HW10	KC05
X∼	$m{D}_{\max }^{0.5}{A}^{-0.25}$	mDmaxA−0.5	$m{D}_{\max }^{2}A$
vterm,Re≪1∼	$D_{\max }^{b_m-0.25b_A-0.5}$	$D_{\max }^{b_m-0.5b_A}$	$D_{\max }^{b_m-b_A+1}$
vterm,Re≫1∼	${\left(D_{\max }^{b_m-0.25b_A-1.5}\right)}^{0.5}$	${\left(D_{\max }^{b_m-0.5b_A-1}\right)}^{0.5}$	${\left(D_{\max }^{b_m-b_A}\right)}^{0.5}$
smonodep=	bm,agg−bm,1	bm,agg−bm,1	bm,agg−bm,1
	−0.25(bA,agg−bm,1)	−0.5(bA,agg−bA,1)	−(bA,agg−bA,1)

Note. The derivation of the scaling relations is shown exemplary for B92 in section 2.2. s_monodep, which gives an estimate of the sign and strength of the dependency of v_term on N_mono is defined in section 4.3.

Figure A1. Same as Figure 7a (aggregates of plates) but using HW10 in (a) and KC05 in (b).

A.2. Sensitivity of the Terminal Velocity to the Hydrodynamic Model and Tumbling A.2.1. Hydrodynamic Models

As mentioned in section 2.2, the hydrodynamic models of B92, KC05, and HW10 differ in several aspects. The Re(X) relation requires assumptions about particle surface roughness, which are differently implemented in the models. Also, the definition of X is different (Table A1). While in B92 X is proportional to $m{D}_{\max }^{0.5}{A}^{-0.25}$, X is proportional to mD_maxA^−0.5 in HW10 and $m{D}_{\max }^2{A}^{-1}$ in KC05. As a result in B92 and HW10, v_term increases slower with decreasing area ratio (A_r=4Aπ⁻¹D⁻²) than in the formulation of KC05. The empirical correction of X due to wake turbulence is also applied in KC05 but not in HW10.

These differences affect the behavior of v_term at large sizes and the monomer number dependency (which we quantify by s_monodep). Without the empirical correction of X (which considers wake turbulence), v_term only saturates if v_term,Re≫1∼D⁰. For example with HW10, the saturation would be reached for b_m−0.5b_A−1=0 (Table A1). This is, for example, not the case for aggregates of plates simulated in this study, and therefore, HW10 does not predict a saturation of v_term at larger sizes (Figure A1a).

Also, the sign and the strength of the increase/decrease of v_term with increasing N_mono depends on the formulation of X. In section 4.3, we introduced s_monodep as a measure for this monomer number dependency. Applying this measure to the aggregates of plates yields s_monodep = −0.21 for HW10 and s_monodep = −0.06 for KC05. Both HW10 and KC05 show the decrease of v_term with increasing N_mono which we saw when using B92, but this decrease is very weak for KC05.

Figure A2. vterm based on m/A−D fits (Tables 1 and 4) and different hydrodynamic models. The particles are horizontally aligned (“no tumbling”) rotated by 20° or 40° around the principal axis to mimic different strength of tumbling. (a) plate monomers; (b) aggregates of plates.

A.2.2. Tumbling

To investigate the effect of the tumbling of the aggregates (as reported, e.g., by Garrett & Yuter, 2014) on the projected area A and v_term, the particles are tilted with a standard deviation of 0° , 20° , 40°, and 60° , around the principal axis (Figure A2). This is done only after the final aggregate is assembled and thereby does not influence the structure of the aggregates. This rotation reduces A, and in turn, v_term increases.

The monomers (top panel in Figure A2) are stronger effected by tumbling (especially at large D_max) due to their lower aspect ratio (not shown). The largest increase in v_term with increasing tumbling is found for KC05 due to the largest increase in the Best number with decreasing A (see section 2.2). B92 shows the least influence of tumbling, which increases v_term at maximum by about 0.1 m s⁻¹ and has a negligible effect on v_term for the aggregates.

Table A2 Mean Absolute Error of the Normalizing Function f_m and f_A for the Monomer Number Dependent “monodep” Fits (Section 4.2) and the Power Law Fits (Section 4.4) for the Particle Mass m and Projected Area A

Monomer type	fm monodep	fm power law	fA monodep	fA power law
Plate	0.1190	0.1261	0.0816	0.0950
Needle	0.3737	0.4539	0.3926	0.5916
Dendrite	0.1698	0.1721	0.1575	0.1566
Column	0.2436	0.2456	0.2728	0.3343

A.3. Mean Absolute Error of the Mass and Area-Size Relations

In sections  4.2 and  4.4, we provided fit relations for mass and area with and without taking into account the monomer number dependency of the simulated aggregates. The mean absolute error of the fits shown in Tables 3 and 4 (normalized by the properties of the monomers; e.g., shown for plates with the green dotted lines in Figure 5) is shown in Table A2.

Acknowledgments

Contributions by M.K., S.K. and D.O. were funded by the German Research Foundation (DFG) under grant KN 1112/2-1 as part of the Emmy-Noether Group Optimal combination of Polarimetric and Triple Frequency radar techniques for Improving Microphysical process understanding of cold clouds (OPTIMIce). A.v.L. gratefully acknowledges the funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project-ID 268020496 TRR 172, within the Transregional Collaborative Research Center “ArctiC Amplification: Climate Relevant Atmospheric and SurfaCe Processes, and Feedback Mechanisms (AC)³”. M.K. also acknowledges support by the Graduate School of Geosciences of the University of Cologne. The source code of McSnow used in this study is hosted at GitLab, and access can be granted by A. S. on request. We also thank Jussi Leinonen for fruitful discussions and for making his aggregation model freely available at GitHub (https://github.com/jleinonen/aggregation). We thank David Hudak and Peter Rodriguez from Environment Canada for providing the PIP measurement from the CARE site. The CARE site is part of the NASA GPM (Global Precipitation Measurement) Ground Validation (GV) program. Model output data of the aggregation model and McSnow are accessible at the ZENODO platform (https://doi.org/10.5281/zenodo.3606668). We also acknowledge the constructive comments of two anonymous reviewers, which helped to improve the manuscript.