Improving the representation of aggregation in a two-moment microphysical scheme with statistics of multi-frequency Doppler radar observations

Aggregation is a key microphysical process for the formation of precipitable ice particles. Its theoretical description involves many parameters and dependencies among different variables that are either insufficiently understood or difficult to accurately represent in bulk microphysics schemes. Previous studies have demonstrated the valuable information content of multi-frequency Doppler radar observations to characterize aggregation with respect to environmental parameters such as temperature. Comparisons with model simulations can reveal discrepancies, but the main challenge is to identify the most critical parameters in the aggregation parameterization, which can then be improved by using the observations as constraints. In this study, we systematically investigate the sensitivity of physical variables, such as number and mass density, as well as the forward-simulated multi-frequency and Doppler radar observables, to different parameters in a two-moment microphysics scheme. Our approach includes modifying key aggregation parameters such as the sticking efficiency or the shape of the size distribution. We also revise and test the impact of changing functional relationships (e.g., the terminal velocity–size relation) and underlying assumptions (e.g., the definition of the aggregation kernel). We test the sensitivity of the various components first in a single-column “snowshaft” model, which allows fast and efficient identification of the parameter combination optimally matching the observations. We find that particle properties, definition of the aggregation kernel, and size distribution width prove to be most important, while the sticking efficiency and the cloud ice habit have less influence. The setting which optimally matches the observations is then implemented in a 3D model using the identical scheme setup. Rerunning the 3D model with the new scheme setup for a multi-week period revealed that the large overestimation of aggregate size and terminal velocity in the model could be substantially reduced. The method presented is expected to be applicable to constrain other ice microphysical processes or to evaluate and improve other schemes.

1 Introduction

Ice growth processes which lead to precipitable particles are essential to understand because more than 60 % of the global precipitation reaching the surface is formed in the ice phase (Heymsfield et al., 2020). Besides depositional growth and riming, aggregation is one of the key growth mechanisms in ice clouds. Aggregation is found to be active in a large temperature range (Hobbs et al., 1974; Kajikawa and Heymsfield, 1989; Field, 2000). As revealed, for example, by radar observations (e.g., Barrett et al., 2019), aggregation can cause a rapid increase in the particle size in favorable conditions, such as the dendritic growth zone close to −15 ^∘C or close to the melting layer (Lamb and Verlinde, 2011). Unlike depositional growth, sublimation, or riming, aggregation does not directly modify the ice and snow water content. However, its strong influence on particle shape, particle size distribution, and terminal velocity v_t links aggregation to other processes, such as depositional growth, sublimation, and riming, that alter the mass flux considerably. Therefore, it is important to accurately represent aggregation in microphysical schemes.

A central component of the theoretical description of aggregation (see also Sect. 3.1) is the aggregation kernel. Therefore, many challenges in accurately simulating aggregation can be discussed by considering the various components of this kernel. The aggregation kernel is defined analogously to collision–coalescence of droplets in liquid clouds.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& K(D_i,D_j)={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}(D_i+D_j{)}^{\mathrm{2}}|v_{\text{t}}\left(D_i\right)-v_{\text{t}}\left(D_j\right)|\\ & E_{\text{coll}}(D_i,D_j)E_{\text{stick}}\left(T\right)E_{i,j}\left(D_i{D}_j\right)\end{array}\end{array}$$

The aggregation kernel is proportional to the probability K of two particles i and j colliding (Gillespie, 1975) and sticking together after collision. This probability increases with increasing size D and relative v_t of the particles, as well as the collision E_coll and sticking efficiency E_stick. Obviously, the size D is well-defined for spherical particles by their diameter, but this is already much more complex for ice and snow particles which have a nonspherical shape. How large v_t of ice and snow particles is also strongly depends on their size, shape, and orientation (Böhm, 1992; Mitchell and Heymsfield, 2005; Heymsfield and Westbrook, 2010). For smaller particles, v_t increases strongly, but the increase in v_t flattens with size and finally v_t approaches a constant value of 1 m s⁻¹ for centimeter-sized aggregates (Lohmann et al., 2016). The size ranges in which v_t increases most rapidly (i.e., has the largest slope) are highly shape-dependent (Barthazy and Schefold, 2006; Hashino and Tripoli, 2011; Karrer et al., 2020). Consequently, the slope of the v_t–size relation is uncertain but at the same time crucial for aggregation. Two remaining parameters, E_coll and E_stick, are also multiplicative in the kernel. E_coll describes the ratio between the actual collision cross section and the geometric cross section. E_coll is smaller than 1 for most particle pairs because typically the smaller and slower particle is deflected due to hydrodynamics as the larger particle approaches. For ice collisions, E_coll is generally poorly constrained (Wang, 2010). This can be easily understood given the enormous variety of particle shapes and orientation, leading to very complex flow fields. In the temperature region most relevant for aggregation ($T>-\mathrm{20}$ ^∘C), the number of activated ice-nucleating particles (INPs) and hence also the concentration of small ice particles rapidly decrease with increasing temperature (Kanji et al., 2017). Except for small ice particles generated by secondary ice production, the effect of small ice particles being deflected around larger particles might therefore be less important for aggregation. In fact, many bulk microphysical schemes (e.g., Seifert and Beheng, 2006) assume the bulk E_coll to be 1.

E_stick is the probability of two ice particles sticking after the collision. Although laboratory (Hosler and Hallgren, 1960; Connolly et al., 2012; Phillips et al., 2015) and in situ (Mitchell, 1988; Kajikawa and Heymsfield, 1989) estimates, as well as multi-frequency radar retrievals (Barrett et al., 2019) of E_stick exist, the exact value of E_stick and its dependency on parameters such as temperature or supersaturation are very uncertain. However, there is widespread agreement in the literature on two main temperature ranges in which E_stick is enhanced: around −15 ^∘C, the mechanical interlock of dendrites increases E_stick compared to the surrounding temperature regions (Pruppacher et al., 1998). In addition, sintering of ice particles due to an increasingly thick quasi-liquid layer (Slater and Michaelides, 2019) on the ice surface causes a general increase in E_stick when temperature rises up to 0 ^∘C. In addition to the aggregation kernel, the aggregation rate is also influenced by the particle size distribution. Simply put, the particles that have a high probability of aggregation, given by the aggregation kernel, must be present in the cloud to have efficient aggregation.

Bulk microphysics schemes cannot simulate aggregation on an individual particle level but require the calculation of bulk aggregation rates. Analytic solutions for the bulk aggregation rates are in principle possible (Verlinde et al., 1990). However, these solutions are computationally expensive and require the usage of power-law relationship for v_t and size, which cannot represent the asymptotic behavior known from observations for large sizes. Approximations of the bulk aggregation rates consider characteristic velocity differences (Wisner et al., 1972; Seifert and Beheng, 2006) and allow the use of more complex v_t–size relations, which consider the asymptotic behavior of v_t at large sizes and nonspherical particle shapes (Seifert et al., 2014).

In general, we need to distinguish between three different aspects of uncertainty in aggregation simulations: (1) a general lack of understanding or quantification of parameters, such as the absolute values of E_stick; (2) formulation of functional relationships, which cannot adequately represent the whole relevant range (e.g., v_t–size relationship); and (3) simplifications that must be made to keep the computational cost affordable, e.g., considering only bulk properties of the particle population. Because of these uncertainties, it is important to constrain the model by observations of aggregation in clouds.

In situ and remote sensing observations have provided valuable information on the general characteristics of aggregation and have allowed estimation of the relative importance of aggregation with respect to other processes. Decades ago, observations had already reported that the largest aggregates are found around −15 ^∘C, which is considered to be a consequence of mechanical interlocking of dendrites, and at temperatures a few degrees below 0 ^∘C, which is related to the quasi-liquid layer (Lamb and Verlinde, 2011).

Radar observations contain valuable information about the aggregation process, which also is the reason we rely on them in this study. The strong temperature dependence of aggregation observed in early studies could be confirmed by radar observations, especially in profiles of absolute and differential reflectivity (Kennedy and Rutledge, 2011; Andrić et al., 2013; Schrom and Kumjian, 2016; Moisseev et al., 2015). By additionally considering the mean Doppler velocity, the relative importance of aggregation and riming can be estimated (Mosimann, 1995; Mason et al., 2018; Kneifel et al., 2020). Furthermore, using radars of different frequencies allows for the estimation of mean particle size (Matrosov, 1998; Hogan et al., 2000; Liao et al., 2005; Szyrmer and Zawadzki, 2014; Kneifel et al., 2015) and therefore better characterization of temperatures at which aggregation is dominant.

Ori et al. (2020, O20) evaluated ice particle growth in simulations of the Icosahedral Nonhydrostatic Model (Zängl et al., 2015, ICON) using the Seifert–Beheng two-moment microphysics scheme (Seifert and Beheng, 2006, SB06) by comparing it with measurements in observational space. To this end, O20 used the multi-month cloud radar dataset from Dias Neto et al. (2019, D18). This quality-controlled dataset is particularly suitable because it contains multi-frequency and Doppler-measured data and thus fingerprints of aggregation and sedimentation. While model–observation comparisons based on a single or few cases can be difficult to interpret due to the specific conditions (specific water vapor field, synoptic situation, etc.) of the case, the statistical comparison of O20 could reveal model-inherent mean biases. The comparison of models and observations in the radar space using a radar forward operator simplifies the assessment of uncertainties because the deviations between models and observations can be directly compared to the variability of the observations. The alternative approach of applying a retrieval to the observations might seem more intuitive because microphysical variables, such as number density, can be compared directly. However, ensuring consistency between a model and retrieval as well as tracing the propagation of uncertainties, for example in the observables or the forward model, is often more complicated (e.g., Reitter et al., 2011).

O20 found an overall correct temperature dependency of aggregation but also revealed an overestimation of the snow size and v_t at temperatures above −7 ^∘C. O20 suggested that inaccurate E_stick and v_t–size parameterization might cause this overestimation. However, direct attribution of the observed biases (e.g., snow that is too large) to a specific component of the aggregation process (e.g., E_stick) requires simultaneous investigation of the influence of all parameters relevant for the aggregation process in a suitable modeling setup.

Microphysics schemes are usually tuned to improve the prediction of key variables, such as precipitation, the energy balance at the top of the atmosphere, or the near-surface temperature (Schmidt et al., 2017; Morrison et al., 2020). Only a small subset of variables (e.g., v_t of cloud ice) are varied during the tuning process, and tuning might be ad hoc rather than evidence-based. As the models simulate complex interacting processes, several parameter combinations can improve the predicting skill of modeled variables such as precipitation. Therefore, it is likely that tuning introduces compensating errors. For example, if two parameters are not accurately implemented, adjusting one of them might improve the model performance even when the adjustment leads away from the true value of the parameter. Detailed remote sensing observations can be used to adjust parameters and make improvements on the process level rather than improving the performance of the entire modeling system. However, because remote sensing observations are sensitive to a limited number of parameters and within a limited range of variability, there is a risk that model parameters may be adjusted to match observations well but still be inaccurate in regimes wherein these observations have low sensitivity. To reduce this risk, new methods for model improvement and development have been proposed whose parameter selection is still based on physical constraints, namely theory and laboratory measurements, but can be optimized by Bayesian inference of observations (Morrison et al., 2020). The advantage of this approach is that uncertainty of both laboratory measurements and remote sensing observations can be considered, and new knowledge about parameters can be continuously incorporated. The combination of several radar observables, such as multiple frequencies, Doppler spectra, and polarimetry, allows the observed signatures to be assigned to a specific microphysical process under some conditions (Kneifel et al., 2015; Kalesse et al., 2016; Pfitzenmaier et al., 2018; Barrett et al., 2019). For example, Barrett et al. (2019) focused their multi-frequency study on the dendritic growth zone, where aggregation is known to be particularly efficient. Hence, the rapidly changing size-dependent, multi-frequency variables could be clearly related to aggregation and a retrieval of E_stick could be obtained.

In addition, novel cloud radar techniques, e.g., multi-frequency Doppler observations, enable the identification of key growth mechanisms (Kneifel et al., 2015; Kalesse et al., 2016; Pfitzenmaier et al., 2018; Barrett et al., 2019). Barrett et al. (2019) identified a temperature range in which aggregation rapidly increases particle size and estimated E_stick from a retrieval using multi-frequency Doppler spectra. Identifying a dominant growth mechanism allows focusing on a single process, which simplifies the inverse problem by reducing the number of parameters and observables to be considered simultaneously.

In this study, we constrain the parameters that influence aggregation by confronting idealized and realistic simulations with the multi-frequency Doppler radar observations from D18. The methods used are described in Sect. 2. We revise all main parameters and functional relationships regarding the aggregation formulation in SB06 by incorporating recently published parameters and revising the bulk aggregation equations. We describe these parameters and formulations in Sect. 3.1 and compare them with the choices in the default SB06 scheme. In Sect. 3.1.5 the selection of the snow particle properties, which is a critical component of both aggregation and radar simulations, is described. The sensitivity of the aggregation and associated radar variables to individual parameters of the revised formulation is evaluated with an ensemble of 1D model simulations (Sect. 3.2). The optimal combination of these simulations is chosen and tested in sensitivity studies in ICON-LEM simulations (Sect. 3.3). Finally, we perform ICON-LEM simulations of several weeks, which we evaluate against the default simulations from O20 and the observations from D18 (Sect. 3.4). This approach allows testing many different parameters against observed statistics of several weeks in a numerically efficient way. Section 4 summarizes the approach and draws conclusions regarding the following questions: how can we investigate the sensitivity of aggregation to the components of its parameterization? How can we improve the representation of aggregation in a two-moment microphysical scheme? Which microphysical parameters influence the simulation of aggregation the most?

2 Methods

The Icosahedral Nonhydrostatic Model (ICON; Zängl et al., 2015) has numerous applications due to its different configurations. ICON-NWP (ICON numerical weather prediction) is used by the Deutscher Wetterdienst (DWD) for operational weather forecast in a global and recently also in all regional setups. ICON's large-eddy mode is called ICON-LEM (Dipankar et al., 2015; Heinze et al., 2017). We use the SB06 two-moment microphysics scheme instead of the single-moment scheme currently used in operational weather forecasting, as do most studies that perform simulations with ICON-LEM. Since simulations with ICON-LEM are relatively computationally expensive, we also use a simple 1D model to efficiently test different parameterizations and their influence on the simulation.

Since we want to further investigate the causes and reduce the discrepancies between modeled and simulated observables, we use the same simulation setup of ICON-LEM as in O20. We only briefly describe the setup here, since an extensive description can be found in O20. The modifications we make to the SB06 microphysics scheme are described in detail in Sect. 3.1.

2.1 “Snowshaft” model

Simple 1D models have been used to assess the influence of several parameters or processes on microphysical or observed quantities (e.g., precipitation rates, polarimetric variables) and to test new parameterizations (Seifert, 2008; Kumjian and Ryzhkov, 2010; Milbrandt and Morrison, 2016; Paukert et al., 2019). These models are much simpler than full 3D models (like ICON-LEM) and are therefore also referred to as rain-shaft models. Because we apply such a simple model to ice microphysics we use the term “snowshaft” model. In these simple models, the atmospheric variables (e.g., temperature gradient, relative humidity) are predefined and feedback mechanisms from microphysics to thermodynamic and thus dynamic variables are neglected. These simplifications allow the analysis of selected processes and their sensitivity to a range of parameters without having to consider the full range of complexity. Another advantage of the snowshaft model is the low computational effort, which allows testing a large number of parameter combinations and process formulations.

The snowshaft model has 250 layers and the temperature spans the range from 0 to −40 ^∘C, which covers the most relevant range for precipitating ice clouds. The temperature profile is linear with a gradient of 0.0062 K m⁻¹. Consequently, the top of the model is at 6450 m. The relative humidity with reference to ice (RH_i) is constant for h>3000 m and increases linearly until it reaches RH = 100 % (RH is the relative humidity with reference to water; Fig. B5). The thermodynamic variables are constant in time and there is no air motion. These simplifications can be justified by the nearly stationary nature of many clouds and the low vertical velocity seen in the dataset of D18.

At the top of the model, a gamma distribution (following the size distribution parameter as described in Table 3) is initialized for cloud ice and snow. Together with the size distribution parameter, the mass density Q and the number density N completely define the size distribution at the model top. Below the model top, the size distribution evolves through the following microphysical processes: sedimentation, depositional growth, and aggregation. These processes are considered dominant below the cloud top (where nucleation is especially important) and above temperatures near the melt layer, where riming rates increase sharply (Kneifel and Moisseev, 2020). The values of RH_i, Q, and N are chosen in Sect. 3.2.1 to match profiles of observables with substantial precipitation.

2.2 ICON-LEM setup

In our simulations, we use a small domain setup of ICON-LEM. This setup has been shown to be both computationally efficient and to represent clouds well in various conditions (Marke et al., 2018; Schemann and Ebell, 2020; Schemann et al., 2020). The domain is circular with a radius of 110 km, and the observational site Jülich Observatory for Cloud Evolution Core Facility (JOYCE-CF; Löhnert et al., 2015) is in the center. At JOYCE the TRIple-frequency and Polarimetric radar Experiment for improving process observation of winter precipitation (Tripex, D18) took place, which we use in the model–observation comparison. The horizontal grid spacing of the simulations is ca. 400 m, and the vertical grid spacing ranges from 20 m at the surface to 380 m at the model top. With a total of 150 vertical layers, the atmosphere is simulated up to a height of 21 km. Initial and lateral boundary conditions are taken from the ECMWF Integrated Forecasting System (IFS). Initialization is carried out each day at 00:00 UTC. IFS data are incorporated as forcing on the lateral boundary every hour.

2.3 SB06 scheme

The SB06 scheme is used in the snowshaft simulations (Sect. 3.2) and as the microphysics scheme in the ICON-LEM simulations (Sect. 3.3 and 3.4). The SB06 scheme is a two-moment scheme that simulates the evolution of the number (N=M⁽⁰⁾) and mass density (L=M⁽¹⁾) from which the mixing ratio ($Q=L\cdot {\mathit{\rho }}_{\text{air}}^{-\mathrm{1}}$) can be easily derived. ρ_air is the air density and Mⁿ (Eq. 2) represents the moments of the mass distribution (Eq. 5).

$$\begin{array}{}\text{(2)}& M^{\left(n\right)}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{m}^{n}f\left(m\right)dm\end{array}$$

The scheme simulates six different hydrometeor classes (cloud water, cloud ice, rain, snow, graupel, and hail). The conversion from one class to another is in general associated with a specific microphysical process. For example, if cloud ice forms aggregates, Q and N of cloud ice are converted to snow (Sect. 3.1). Therefore, it is consistent to assume properties of monomers for cloud ice and properties of aggregates for snow. The predefined particle properties of the default setting of the scheme are listed in Table 2 for each hydrometeor, along with the properties of the cloud ice and snow class alternatives proposed in Sect. 3.1.

In the SB06 scheme, aggregation rates are the product of collision rates and E_stick because E_coll is assumed to be 1. In the scheme, the variance approximation (SB06), based on the work of Wisner et al. (1972), provides a computationally feasible analytical solution of bulk collision rates. The variance approximation of Seifert and Beheng (2006) avoids the usage of pre-calculated lookup tables (Seifert et al., 2014) and, unlike Wisner et al. (1972), is able to estimate collision rates of self-collection, i.e., aggregation within a particle class. The latter is made possible by considering the square root of the second moment of the velocity differences, which also has the advantage over the approximation by Wisner et al. (1972) that the collision rates between different particles are nonzero even if their bulk velocities are equal. The default SB06 scheme assumes power-law relations for the v_t–size relation in the calculation of the collision rates. The extension of the variance approximation of Seifert et al. (2014), which allows using Atlas-type v_t–size relations (Sect. 3.1.3), is applied in the SB06 scheme for the first time in this study.

Details of the components of the aggregation process considered in the SB06 scheme can be found in Sect. 3.1 and Appendix A.

2.4 Passive and Active Microwave radiative TRAnsfer tool (PAMTRA)

The Passive and Active Microwave radiative TRAnsfer tool (PAMTRA; Mech et al., 2020) is used to simulate synthetic radar observations. Microphysical properties are represented consistently in the SB06 scheme and PAMTRA (Table 2).

Throughout the study, we adopt the same scattering assumptions for each of the hydrometeor classes in the SB06 default scheme (“SB cloud ice”, “SB snow”, “cloud droplet”, “rain”, “graupel”, and “hail” in Table 2). As in O20, we apply the self-similar Rayleigh–Gans approximation (SSRGA; Hogan and Westbrook, 2014; Hogan et al., 2017) and coefficients derived from 3D models of aggregates of plates for cloud ice and aggregates of needles for snow. In O20, the coefficients used for the snow class were slightly adjusted to closely match the observed triple-frequency signature. The SSRGA parameters of aggregates of plates are also used for the new cloud ice categories (“column” and “needle” in Table 2). For Mix2, SSRGA parameters derived from the same 3D models used for the determination of particle properties (Karrer et al., 2020, K20) are available (Ori et al., 2021). Since we find little influence of SSRGA parameters in Sect. 3.1.5, we use the adjusted SSRGA properties of the aggregates of needles from O20 for the Mix2 aggregates throughout the study to be consistent with 020, although using the SSRGA parameters derived from the same 3D aggregate models would be most physically consistent.

2.5 Multi-frequency radar approach

Like O20, we use multi-frequency observations to derive information about the aggregation process. Multi-frequency observations are useful to distinguish the size of particles, since the ratio of wavelength and particle size along with the particle density are the main factors that determine their scattering properties. The scattering of particles much smaller than the wavelength can be approximated well by the Rayleigh approximation. For larger particles, however, the interference of waves scattered from different parts of the particles must be considered (Kneifel et al., 2020), which leads to differential scattering among the various frequencies.

The ratio between the reflectivities of two radars with operating wavelengths λ₁ and λ₂ (λ₁<λ₂),

$$\begin{array}{}\text{(3)}& {\text{DWR}}_{{\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}}}={\displaystyle \frac{\text{Ze}\left({\mathit{\lambda }}_{\mathrm{1}}\right)}{\text{Ze}\left({\mathit{\lambda }}_{\mathrm{2}}\right)}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{4}}}{{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}}}{\displaystyle \frac{\int {\mathit{\sigma }}_b(m,{\mathit{\lambda }}_{\mathrm{1}})f\left(m\right)dm}{\int {\mathit{\sigma }}_b(m,{\mathit{\lambda }}_{\mathrm{2}})f\left(m\right)dm}},\end{array}$$

quantifies the amount of differential scattering. DWR is called the dual-wavelength ratio, Ze is the reflectivity, and σ_b is the backscattering cross section (all variables in linear units). Although differential attenuation can also contribute to DWR (Battaglia et al., 2020), we did not include this effect in Eq. (3) because the processing of D18 already corrects for the impact of differential attenuation on DWR. D18 evaluated the absolute calibration of the observed Ze values from the Ka-band radar using disdrometer measurements during rainfall. After correcting differential attenuation due to gases at all three frequencies, the Ka-band radar was then used as a reference for estimating calibration biases and differential attenuation effects due to hydrometeors by comparing the three Ze values at cloud top. The DWRs caused by differential scattering are usually close to 0 dB for small ice particles present at the cloud top. Calibration biases can be identified as DWR biases which are relatively constant over time; differential attenuation effects due to supercooled liquid water, rainfall, or the melting layer vary more strongly on shorter timescales (minutes to hours). The path-integrated differential attenuation estimated at cloud top was then used to correct the DWRs in the entire profile. A more in-depth discussion of various correction methods for multi-frequency radar observations is provided in Tridon et al. (2020). If differential scattering effects are the only contributor to DWR, it correlates well with the mean mass of the distribution f(m) (Sect. 3.1.1), as can be seen from Eq. (3). For small particles, the Rayleigh approximation is valid for all frequencies and σ_b scales with the mass squared. However, for larger particles and shorter wavelengths, σ_b is smaller than predicted by the Rayleigh approximation and σ_b(m,λ₂) is smaller than σ_b(m,λ₁). As a result, particle populations that contain larger particles, e.g., due to their large mean mass, have larger DWRs than particle populations with smaller mean masses. Mason et al. (2019) and others have shown that not only the mean mass, but also the shape of the distribution, the particle density, and the internal structure of the particles (through σ_b) can substantially affect the DWRs. Given the radars available in D18, we investigate the sensitivity of aggregation by analyzing DWR_X,Ka and DWR_Ka,W. The subscripts W, Ka, and X denote the radar bands and, more specifically, the wavelengths of 3.3, 8.6, and 31 mm. Each combination of wavelengths is sensitive to a different range of particle sizes. For example, DWR_Ka,W is most sensitive to mean particle sizes of unrimed cloud ice and snow between 0.5 and 3 mm, and DWR_X,Ka is sensitive between 1.5 and 15 mm (O20). Outside this sensitivity range, DWRs are zero (small mean size) or asymptotically approach (saturate) a DWR value (large mean sizes) that depends on the scattering properties of the particles present. More detailed information on the approach and its sensitivities can be found in O20.

Moreover, D18 reported that strong riming is rare in their dataset, so aggregation is the main contributor to particle growth and thus the increasing DWRs from cloud top to cloud bottom.

3 Results and discussion

3.1 Ice microphysical parameters influencing aggregation

To interpret the following sensitivity experiments, we describe which parameters need to be considered in the simulation of aggregation in a bulk scheme, which parameters and process formulations are currently used in the SB06 scheme, and how the assumptions could be updated with recently published parameterizations.

The stochastic collection equation (SCE) describes how the particle distribution (PSD_m) changes with time under the action of aggregation (Khain et al., 2015).

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\displaystyle \frac{\partial f\left(m_i\right)}{\partial t}}=\underset{\mathrm{0}}{\overset{m_i/\mathrm{2}}{\int }}f\left(m_j\right)f(m_i-m_j)K(m_i-m_j,m_j)d{m}_j\\ & -\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}f\left(m_i\right)f\left(m_j\right)K(m_i,m_j)d{m}_j\end{array}\end{array}$$

Here, f(m) is the particle distribution as a function of mass and K is the aggregation kernel described in Sect. 3.1.2. The first term of Eq. (4) describes the gain of particles of mass m_i by aggregation of particles with masses m_j and m_i−m_j. The second term considers the loss of particles of mass m_i by aggregation with particles of mass m_j (illustrated in Fig. 1a and b). In general, PSDs cannot be perfectly described by simple functional relationships (e.g., gamma distribution) but can have complex shapes (Fig. 1a). Thus, explicit prediction of the evolution of PSDs must take into account the full SCE.

Bulk schemes, however, can only account for the evolution of the PSD in a simplified form. The tendencies of the moments in the SB06 scheme (mass density: $\partial Q/\partial t$, number density: $\partial N/\partial t$) can be calculated by considering only the loss term. The reason for this can be further explained by Fig. 1c–h, where the collision events are separated among the ice (monomers) and snow (aggregates) classes. In fact, because of the mass conservation, the total mass of particles gained (integral of the first term) has to match the total mass of particles lost (integral of the second term). Since it is assumed that within one time step a particle can participate only in one collision event, only one snow particle results from the collision of two ice particles (number of arrows in Fig. 1c and d). The same applies for the ice–snow and snow–snow collisions, but here there is no conversion of N from one category to another but only a loss of N_i or N_s. Thus, it is sufficient to calculate only one collision rate for each of the three considered collision scenarios (ice–ice, ice–snow, snow–snow) and moments (N and Q).

Figure 1 Illustration of the SCE (Eq. 4) for an explicitly resolved PSD_m (a, b) and when applied to the cloud ice and snow classes of the SB06 scheme (c, h). The left column depicts the loss term (second term in Eq. 4) and the middle column the gain term (first term in Eq. 4). The right column shows the sign and connection of the tendency of the bulk moments. Arrows indicate whether the number density is rising or falling at the specified mass. Red lines indicate the ice distribution and blue the snow distribution. The arrows are red if ice particles are collected and blue if snow particles are collected or are created as a result of the collision.

3.1.1 Size distribution

In most bulk schemes, the PSD is described by the generalized gamma distribution or simplifications thereof. With the mass m as a primary variable, the generalized gamma distribution can be written as

$$\begin{array}{}\text{(5)}& f_m\left(m\right)=N_{\mathrm{0},m}{m}^{{\mathit{\nu }}_m}\exp (-{\mathit{\lambda }}_m{m}^{{\mathit{\mu }}_m}).\end{array}$$

For some applications, using the mass-equivalent diameter

$$\begin{array}{}\text{(6)}& D_{\text{eq}}={\left({\displaystyle \frac{\mathrm{6}m}{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}}\right)}^{\mathrm{1}/\mathrm{3}}\end{array}$$

as a primary variable and the ordinary gamma distribution is more convenient:

$$\begin{array}{}\text{(7)}& f\left(D_{\text{eq}}\right)=N_{\mathrm{0},\mathrm{eq}}{D}_{\text{eq}}^{{\mathit{\mu }}_{\text{eq}}}\exp (-{\mathit{\lambda }}_{\text{eq}}{D}_{\text{eq}}),\end{array}$$

where D_eq is the mass-equivalent diameter. One such application is the use of the Atlas-type v_t–size relationship (Eq. 11) in the calculation of collision rates in Appendix A. Size distributions derived from in situ observations are usually presented as a function of the maximum dimension D_max, which is often derived by circumscribing a sphere or spheroid to the projected particle image.

$$\begin{array}{}\text{(8)}& f\left(D_{\text{max}}\right)=N_{\mathrm{0},max}{D}_{\text{max}}^{{\mathit{\mu }}_{\text{max}}}\exp (-{\mathit{\lambda }}_{\text{max}}{D}_{\text{max}})\end{array}$$

In general, a distribution described by Eq. (5) cannot be expressed by Eq. (7) or (8). Only when ${\mathit{\mu }}_m=\mathrm{1}/\mathrm{3}$ can Eq. (5) be expressed by Eq. (7). To allow conversion of Eq. (5) to (8), μ_m must be set to $b_m^{-\mathrm{1}}$ (exponent in the m–D_max relation; Eq. 12). As we calculate the collision rates of particles following an Atlas-type v_t–size relation (Appendix A), we need to set ${\mathit{\mu }}_m=\mathrm{1}/\mathrm{3}$. Since b_m≠3 for cloud ice and snow, μ_max in Eq. (8) can only be approximated.

The PSD shape can vary strongly, e.g., for nonstationary events (Seifert, 2008). Furthermore, ν_m, or equivalent parameters in distributions that use a different primary variable, is often described as a function of other parameters (e.g., the mean size; Heymsfield, 2003). Nevertheless, in the current version of the SB06 scheme, we must choose a single value of ν_m in each simulation. Therefore, we test two different values of ν_m in the simulations and later (Sect. 3.2) select the one with which the simulations can reproduce the observations the best. The SB06 standard configuration (ν_m=0) corresponds to μ_eq=2.0 and ${\mathit{\mu }}_{\text{max}}=-\mathrm{0.11}$. If we use ν_m=2 instead of ν_m=0, we obtain a narrower distribution with much fewer particles at small masses and a peak near the mean mass. Considering Heymsfield (2003), ν_m=0 is representative of a mean mass diameter D_mean of about 0.5 mm, and ν_m=2 is representative for D_mean of about 0.2 mm. Many studies have shown that the size distribution parameters are correlated (e.g., Field et al., 2005; Mcfarquhar et al., 2015), further complicating the selection of ν_m. Moreover, PSDs can exhibit bimodalities, e.g., due to secondary ice generation (Korolev and Leisner, 2020), which can be accounted for by the two classes of cloud ice and snow in the SB06 scheme.

Table 1 Size distribution parameter for ${\mathit{\mu }}_m=\mathrm{1}/\mathrm{3}$, the mass–size relationship of the Mix2 particles (Table 2), $q=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ kg m⁻³, and N=10⁴ m⁴ (same as Fig. 2). Ze_Ka, MDV_Ka, and DWR_X,Ka are calculated using the self-similar Rayleigh–Gans approximation (SSRGA) and the SSRGA parameters of Mix2 as provided by Ori et al. (2021). μ_max is estimated by the zeroth, third, and sixth moment of the distribution.

Figure 2 Particle distribution as a function of mass (PSD_m) with a mass density of $q=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ kg m⁻³ and a number density of qn=10⁴ m⁻³ illustrating PSDs with a different shape parameter μ_m.

The PSD width affects the aggregation rates and the radar variables. The narrower the distribution is, the lower the aggregation rates are. This is obvious from the bulk collision rates in Appendix A and can be explained by the small v_t difference of similarly sized particles (Sect. 3.1.2). The PSD width also affects the radar observables. The reflectivity in the Rayleigh regime is proportional to the second moment of the PSD. A narrower distribution reduces the number of large particles (above 10⁻⁷ kg in Fig. 2). Therefore, the reflectivity (Ze) and mean Doppler velocity (MDV) are slightly lower for a narrower distribution compared to a broader distribution with the same Q and N. This effect is even stronger for DWRs, as the large particles contribute the most to the differential scattering signal (Table 1).

3.1.2 Collision kernel

The D-kernel (Eq. 1), defined analogously to the collisional coalescence of droplets in liquid clouds, is often used not only for particles that can be approximated well by spheres (e.g., cloud droplets, hail), but for all particles. However, the collision cross section of nonspherical particles is smaller than the one of spheres with the same D_max because of the presence of voids in their circumscribing sphere. This deviation was previously considered, e.g., a part of E_coll (Keith and Saunders, 1989; Pruppacher and Klett, 2010) by using the equivalent circular radii $r_i=(A_i/\mathit{\pi }{)}^{\mathrm{0.5}}$ as a characteristic length. Using the D-kernel with a constant E_coll that does not depend on particle size (as done, e.g., in SB06), the D-kernel approximation cannot account for the decrease in A_r with increasing size (Fig. 3d). Therefore, we test whether an alternative formulation of the collision kernel that takes the projected areas into account (A-kernel; Connolly et al., 2012) provides a better approximation.

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& K(D_i,D_j)={\left(A_i(D_i{)}^{\mathrm{0.5}}+A_j(D_j{)}^{\mathrm{0.5}}\right)}^{\mathrm{2}}\\ & \left|v_i\right(D_i)-v_j(D_j\left)\right|E_{\text{stick}}\left(T\right)E_{\text{coll}}(D_i,D_j)\end{array}\end{array}$$

The A-kernel approximation has been used previously in the same or similar formulation (Kienast-Sjögren et al., 2013; Morrison and Milbrandt, 2015; Dunnavan, 2021). In these studies, the aggregation rates are calculated numerically and, in the case of the scheme proposed by Morrison and Milbrandt (2015), stored in lookup tables that are used at the model run time. Lookup tables can accurately store precomputed process rates and might be numerically more efficient than analytical solutions, depending on the computer architecture, size of the lookup table, and complexity of the analytical solution. However, Seifert et al. (2014) argue that the use of lookup tables also has disadvantages, like increasing complexity during preprocessing, additional memory access, and difficult reproducibility for follow-up studies. To avoid these disadvantages, the SB06 scheme uses analytical solutions of the variance approximation introduced by Seifert and Beheng (2006). To use the A-kernel we have to generalize the collision rates. For brevity, we moved the lengthy derivations to Appendix A. To our knowledge, this is the first application of an A-kernel in a bulk microphysics scheme that uses an analytical formulation of the aggregation rates. How large the difference is between the D- and the A-kernel depends on the particle properties (e.g., area–size and v_t–size relation).

3.1.3 Particle properties

Particle properties influence aggregation because they are an essential part of the aggregation kernel. According to Eqs. (1) and (9) collection is enhanced if the product of the difference in v_t and the joint cross section is large. Thus, a particle population will aggregate rapidly if the mean mass is relatively large and particles with largely different v_t are present. The coefficients of area–size and v_t–size relations of the SB06 default scheme and the particle from K20 are included in Table 2.

While the particle properties of the SB06 default scheme particle classes are taken from in situ observations, K20 used an aggregation model and hydrodynamic theory to simulate the particle properties. The advantages of this approach are that particle properties can be studied over a large size range, are physically consistent, and can be studied in great detail. Particle property relations from in situ observations have a comparably small sample size. Thus, extrapolation to small and large sizes is unavoidable because microphysics schemes need information about particle properties in a large size range. This extrapolation might lead to inaccuracies, such as the overestimation of v_t at large sizes (K20). Since we take all snow particle properties (m–size, A–size, v_t–size; Table 2) from the same aggregate type within the dataset, all properties are physically consistent. By comparing with in situ observations, K20 found that their mixed aggregates consisting of small columns and large dendrites (Mix2) can approximate mean aggregate properties well. Besides aggregates (including aggregates of columns and aggregates of dendrites; Sect. 3.1.5), K20 also summarized different monomer particle properties, e.g., the columns and needles shown in Fig. 3.

v_t of the default cloud ice and snow class increases continuously with increasing size (Fig. 3) due to the power-law relation used.

$$\begin{array}{}\text{(10)}& v_{\text{t}}=a_{\text{vel}}{m}^{b_{\text{vel}}}\end{array}$$

Due to this continuous increase, the self-collection rates of these hydrometeor classes stay relatively large at large sizes (Figs. B3 and B4). In contrast, the asymptotic approach to a limit of v_t in the new relations leads to rapidly decreasing collision rates at large sizes. The asymptotic approach is evident from in situ observations and can be accounted for by using an Atlas-type v_t–size relation.

$$\begin{array}{}\text{(11)}& v_{\text{t}}={\mathit{\alpha }}_v-{\mathit{\beta }}_v\exp (-{\mathit{\gamma }}_v{D}_{\text{eq}})\end{array}$$

The relative v_t of cloud ice and snow particles also plays a role in ice–snow collection rates. In the SB06 default scheme, v_t of cloud ice and v_t of snow differ greatly. However, K20 showed that v_t of cloud ice and snow should have similar values. The difference between cloud ice and snow v_t determines the location and magnitude of the minimum of the collection rates.

The projected area A is derived differently in the D- and the A-kernel. In the D-kernel, the m–D_max relation,

$$\begin{array}{}\text{(12)}& m=a_m{D}_{\text{max}}^{b_m},\end{array}$$

determines the relation between A and size. Since m is the primary variable in the SB06 scheme, it is most useful to consider the differences between the kernels and the particle classes as a function of D_eq (which is directly related to the mass).

$$\begin{array}{}\text{(13)}& A_{\text{sphere}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{D}_{\text{max}}^{\mathrm{2}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{\left({\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{w}}{D}_{\text{eq}}^{\mathrm{3}}}{\mathrm{6}{a}_m}}\right)}^{\frac{\mathrm{2}}{b_m}}\end{array}$$

Thus, the particles which have the lowest effective density,

$$\begin{array}{}\text{(14)}& {\mathit{\rho }}_{\text{eff}}={\displaystyle \frac{\mathrm{6}m}{\mathit{\pi }{\mathit{\rho }}_{\text{ice}}{D}_{\text{max}}^{\mathrm{3}}}},\end{array}$$

have the largest A for a given D_eq (e.g., needles of K20 in Fig. 3b). The other particles have similar A. In the A-kernel, the actual projected area A_act derived from the particle shape is relevant.

$$\begin{array}{}\text{(15)}& A_{\text{act}}={\mathit{\gamma }}_A{D}_{\text{eq}}^{{\mathit{\sigma }}_A}\end{array}$$

The particle shapes and thus A_act are not defined for the SB06 default classes because this property is not required. The area ratio A_r is commonly defined as the ratio of A_act to the area of a sphere with diameter D_max.

$$\begin{array}{}\text{(16)}& A_{\text{r}}={\displaystyle \frac{\mathrm{4}{\mathit{\gamma }}_A{D}_{\text{eq}}^{{\mathit{\sigma }}_A}}{\mathit{\pi }{D}_{\text{max}}^{\mathrm{2}}}}\end{array}$$

At small sizes, A_r is close to 1, indicating compact particles and small differences between the D- and the A-kernel (Fig. 3d). With increasing size, A_r decreases down to 0.2 at D_eq=5 mm for the Mix2 class and lower for the cloud ice classes needle and column. A_act is similar to observations of Mitchell (1996) as shown in K20.

However, the low values of A_r of the cloud ice classes are less important because such large sizes of cloud ice are rarely reached in the model. The decrease in A_r leads to a decrease in collision rates, especially at large sizes, similar to the Atlas-type v_t–size relations. Thus, combining the new v_t–size relations with the A-kernel substantially decreases collision rates at large sizes.

While the properties of snow can be validated well against mean observed quantities (as done in K20 and in Sect. 3.1.5 of this study), selecting a single habit for cloud ice is a strong simplification that is necessary for a simplified microphysics scheme.

Figure 3 Particle properties from the default (SB06 default cloud ice, SB06 default snow) and modified version (column, needle, Mix2) of the scheme. (a) Terminal velocity v_t, (b) projected area A of a circumscribing sphere (as assumed in D-kernel), (c) “real” projected area A considering the voids in the circumscribing sphere (as assumed in A-kernel), and (d) area ratio (Eq. 16). The default scheme does not assume an A–D relation explicitly, and therefore the real projected area and the area ratio are not given.

Table 2 Parameterizations used in ICON-LEM, the snowshaft model, and radar forward simulations of hydrometeor properties in PAMTRA. D represents the particle maximum dimension and $D_{\text{eq}}={\left(\frac{\mathrm{6}\text{m}}{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}\right)}^{\mathrm{1}/\mathrm{3}}$ the mass-equivalent diameter; m is the particle mass and ρ_w the density of water. The mass–size (m–D), terminal velocity v_t–size, and projected area–size (A–D) relations are reported in their full mathematical form. For the SSRGA scattering model, the four parameters (κ, β, γ, ζ₀) are given in parentheses. SB indicates that the properties are exclusively used in the default setup. Cloud droplets, rain, graupel, and hail (which are only relevant for the 3D simulations) follow the same properties in all simulations. The aspect ratio is 1.0 for all classes except for the snow classes (SB snow, Mix2, and Mix2; O20 scat), for which an aspect ratio of 0.6 is assumed. All variables are in SI units.

3.1.4 Sticking efficiency

The parameters discussed so far determine how often collisions occur. The percentage of the colliding particles that stick together after a collision is defined by the sticking efficiency E_stick.

Figure 4 The sticking efficiency (E_stick) in the SB06 scheme for collisions among ice particles (ice self-collection) follows L83; for other collisions (ice–snow collection, snow self-collection) it applies the C86 parameterization. Our new relation (red) combines the relations from M88 and C12 with a characteristic maximum around −15^∘ C and values quickly approaching unity for temperature larger than −5^∘ C.

E_stick is mostly only described as a function of the temperature (Mitchell, 1988; Connolly et al., 2012, M88, C12). To stick to each other, ice particles must form ice bonds (Lamb and Verlinde, 2011), which is highly unlikely for colliding solid-ice particles when the temperature is well below the melting temperature and the particles only touch for a short time. There are two main mechanisms that increase the likelihood of adhesion after a collision and explain the temperature dependence. The first mechanism is explained by the quasi-liquid layer (QLL) on the ice particle surface. The phenomenon of QLL has been studied since the mid-19th century (Slater and Michaelides, 2019). QLL thickens with increasing temperature and consists of weakly bound molecules on the particle surface (Slater and Michaelides, 2019). When two particles touch, the molecules form a solid bond at the point of contact. The second mechanism is the mechanical interlocking of relatively large particles with dendritic features (Pruppacher et al., 1998). These dendritic features occur at temperatures between −17 and −12 ^∘C.

The SB06 default scheme uses the E_stick parameterization of Cotton et al. (1982) for ice–ice collisions and Lin et al. (1983) for ice–snow and snow–snow collisions (Fig. 4). The exponential shape of both parameterizations can be justified by the approximately exponentially increasing QLL thickness. These relations, however, miss the maximum of E_stick suggested by studies (M88, C12) that consider the mechanical interlocking mechanism.

We combine M88 and C12 to propose a new parametrization. For $T<-\mathrm{20}$ ^∘C we follow C12, then linearly approach the plateau proposed by M88 with E_stick=1 between −17 and −12 ^∘C. As discussed in the Introduction, there is ample evidence from both in situ and remote sensing observations that E_stick is likely to be present at temperatures near −15 ^∘C (at which particles with dendritic features are present) and near the melt boundary. At −10 ^∘C the new parameterization again follows C12 but increases towards 1 at higher temperatures, at which C12 does not provide an estimate of E_stick. One might prefer to follow C12 rather than M88, since C12 derived E_stick directly from laboratory measurements and M88 provided only an ad hoc parameterization. However, C12 analyzed only the initial stage of aggregation, during which few monomers compose the aggregates. The interlocking mechanism might be more efficient for more complex aggregates compared to early aggregates as discussed in C12. Even considering only the initial stage of aggregation, the confidence interval of E_stick at −15 ^∘C ranges from 0.35 to 0.85 (C12).

3.1.5 Selecting a particle type representative for a large aggregate ensemble

Figure 5 Comparison of the modeled and observed relationship between MDV and DWR: (a) DWR_Ka,W, (b) DWR_X,Ka. The histogram shows the observations from the Tripex campaign (D18). The lines show the theoretical MDV at a given DWR for the v_t size relations of snow particles as assumed in the SB default (black) and as modeled in K20. For the dashed lines, the SSRGA parameters have been directly derived from the corresponding aggregate ensemble properties (as found in Ori et al., 2021). The solid lines use SSRGA parameters as used in O20 in order to illustrate the uncertainty due to the scattering parameters. The lines are calculated using PAMTRA and the properties of the US standard winter atmosphere at 700 hPa.

After discussing the various components of the aggregation process formulation, we need to decide which aggregate type to use to best represent the physical particle properties (e.g., v_t) and scattering properties. In O20, the particle properties were defined by the assumptions in the standard SB06 scheme. The best-fitting aggregate model and associated SSRGA parameters were selected based on the best fit in the triple-frequency DWR space. In this section, we ask whether there is an aggregate type in the database of K20 and Ori et al. (2021) that reproduces both the physical and scattering properties well compared to the observations.

O20 already noted that the representation of MDV as a function of DWR resembles to some extent the underlying v_t–size relation. In contrast to the triple-frequency DWR-DWR, the MDV-DWR space is rather insensitive to the PSD width. Different aggregate types composed of different monomer types generated and studied in K20 are used to simulate their corresponding MDV-DWR signatures (Fig. 5). The underlying distribution shows the observed values from D18, which naturally contain larger scatter and even negative DWR values, mainly due to imperfect radar volume matching (for a detailed discussion, see D18). Fortunately, as shown in D18, the dataset contains only very short and weak riming events. This scarcity of substantial riming is important because the increased MDV due to riming would bias our comparison. Moderately or strongly rimed particles would exceed 1.5 m s⁻¹ upon reaching a size that results in a nonzero DWR_Ka,W (Mason et al., 2018). The MDV-DWR space is also well-suited to evaluate our aggregate choice, as it combines the two radar variables that showed the largest discrepancies with the model simulations in O20.

O20 already recognized the overestimation of v_t at large sizes, which is also evident in Fig. 5. For example, at DWR_X,Ka=5 dB the observed MDV scatters around 1 m s⁻¹, while the snow falls at 1.7 m s⁻¹ in the SB06 default scheme. From the aggregate dataset of K20 the aggregates of dendrites fall the slowest and the aggregates of columns fall the fastest. A mixture composed of small columns and large dendrites (Mix2), which fit in situ observations (K20) best, also matches the observations in the MDV-DWR space well. Therefore, we utilize the Mix2 aggregate properties as an improved description for the snow class in the following.

Interestingly, the use of the SSRGA coefficients of the aggregate type O20 does not lead to a strong change in the curves in the MDV-DWR space. Although it would be most consistent to use the SSRGA coefficients of Mix2 directly, we will use the scattering properties of O20 in the following analysis to allow a fair comparison of our new results with the discrepancies found in O20.

3.2 Exploring sensitivity to microphysical parameters in the snowshaft model

The snowshaft model (Sect. 2.1) allows us to test the influence of the particle properties, the formulation of the collision kernel, E_stick, and the size distribution on the aggregation rates with low computational effort and with reduced complexity. In Sect. 3.1 we showed how these parameters affect aggregation. We not only examine the influence of the parameters on the predicted model variables but also on the radar observables. After carefully setting up the model, the comparison in radar space enables us to directly contrast the statistics of the simulation and the observations, as given in O20 and D18. Since we compare the statistics of the model and observations over a relatively long time range this analysis already attempts to select a combination of parameters that can reproduce the observational statistics well. The optimal parameter combinations found in the snowshaft simulations will then be applied in the 3D model to simulate a case study (Sect. 3.3) before we use it to rerun simulations for the whole time period of the Tripex campaign (Sect. 3.4).

This comparison between the model and observation benefits from the simultaneous consideration of multiple model parameters and multiple observables. When looking at a single observable only, one might reduce a bias by an adjustment of a single process or parameter, even though this might just compensate for an inaccurate choice in another parameter, introducing compensating errors. As the number of independent observables increases, this problem is reduced as the inaccurate choice of a parameter might be detectable in one of the remaining observables. In other words, the larger degree of freedom in the observations helps to better constrain the parameters by comparison with the model when several observables are considered. We focus our comparison on the DWRs (as a measure of particle size) and the MDV (as a measure of v_t). These two quantities constrain the strength of aggregation and the assumed v_t–size relationship, and the statistical comparison in O20 also revealed the largest differences between observations and the model in these variables.

3.2.1 Optimizing the snowshaft model and selecting microphysical parameters for new setup

O20 pointed out that the inconsistencies between observed and synthetic MDV and DWRs are especially evident for raining periods. As we attempt to remove these inconsistencies, the atmospheric variables and the hydrometeor contents at the top of the simulation are chosen so that the hydrometeor profiles in the snowshaft simulation roughly follow the profiles of the ICON-LEM simulations from O20 wherein RR is larger than 1 mm h⁻¹; compare “default” and the histogram in Fig. 6. To match the profiles the RH_i has to be set to 1 % above about −18 ^∘C with increasing values up to about 6 % at about −7 ^∘C. These values of RH_i, which are relatively high compared to those from the ICON-LEM simulations (Fig. B5), might be necessary because of the absence of nucleation and advection in the snowshaft simulations. Also, the values of Q_i, N_i, Q_s, and N_s at the model top are chosen so that the hydrometeor profiles of the CTRL simulation (performed with the SB06 default setup) match those of the profiles of the ICON-LEM simulations of O20 with RR > 1 mm (Fig. 6). After this optimization of the snowshaft model, the simulated profiles from ICON-LEM (O20 and Figs. 11 and 12) and the snowshaft model (Figs. 6) reveal that a simple initialization (nucleation) of the profiles at cloud top is sufficient at least for testing the sensitivities of aggregation to our set of parameters and various formulations.

After iterating over many parameter combinations, we found one particular setup (which we refer to as colMix2_Akernel or simply as NEW) to match the observed profiles particularly well. In these iterations, we varied mostly the less-known components, e.g., the size distribution width, while parameters that we were already better able to constrain (Sect. 3.1.5), e.g., the v_t size relation, were not varied. Our approach can hence be seen as a combination of a purely physically based approach, incorporating current knowledge of parameters obtained, e.g., through laboratory studies, and an empirical correction based on observations.

3.2.2 Sensitivity of aggregation to individual ice microphysical parameters in the snowshaft model

The hydrometeor profiles (Fig. 6) and radar observables (Fig. 7) of the NEW setup exhibit many interesting differences from the profiles of the CTRL run. In the following, we discuss where the differences originate from by looking at the different sensitivity runs. In each sensitivity run only one set of parameters is different from the NEW run (Table 3).

The cloud ice mixing ratio Q_i and the cloud ice number density N_i are lower in the NEW run than in the CTRL run for $T<-\mathrm{10}$ ^∘C (Fig. 6). At the same time, the snow mixing ratio Q_s and number density N_s are slightly larger in the NEW run at temperatures below −17 ^∘C. These differences can be explained by the higher E_stick at lower temperatures in the NEW setup (Fig. 4), which leads to more collisions among cloud ice particles, and therefore more particles are converted from the cloud ice to the snow category. When using the E_stick parameterization of Cotton et al. (1982) and Lin et al. (1983) (colMix2_Akernel_LinCot;), Q_i and N_i are larger at lower temperatures (and Q_s and N_s are smaller).

The smaller values of E_stick in colMix2_Akernel_LinCot compared to NEW at lower temperatures (compare L83 and C86 with “new” in Fig. 4) lead to further differences; colMix2_Akernel_LinCot has a smaller mean mass $\stackrel{\mathrm{‾}}{x}$, which is the mean mass of the sum of the cloud ice and snow class ($(Q_i+Q_{\text{s}})(N_i+N_{\text{s}}{)}^{-\mathrm{1}}$), and correspondingly lower DWRs for $T<-\mathrm{7}$ ^∘C (Fig. 7c and d). The smaller mean size also leads to slower-falling particles (visible in MDV; Fig. 7b). For $T>-\mathrm{7}$ ^∘C the strong increase in E_stick in colMix2_Akernel_LinCot triggers a strong increase in $\stackrel{\mathrm{‾}}{x}$ and DWR_X,Ka. A similar increase in the mean and median of the investigated statistics of DWR_X,Ka was already discussed in O20. As in O20 the strong increase is not visible in DWR_Ka,W, since this observable already reaches saturation for mass median diameters of about 3 mm (Sect. 2.5). The local maximum of the new E_stick parameterization at temperatures from −17 to −12.5 ^∘C leads in the NEW run to a rapid increase in the $\stackrel{\mathrm{‾}}{x}$, DWRs, and MDVs in the same temperature range and therefore matches the observed profile of DWR_Ka,W better than the CTRL run.

O20 speculated that the overestimation of the particle sizes at high temperatures and the mismatch in the profiles of the DWRs might be mainly due to the E_stick parameterization and the v_t–size relation. However, Figs. 6 and 7 as well as the aggregation rates (Appendix A) reveal that the v_t relation at smaller sizes and the aggregation kernel formulation also strongly affect the aggregation rates. Both $\stackrel{\mathrm{‾}}{x}$ (Fig. 6i) and DWR_X,Ka are lower in colMix2_Akernel_LinCot than in the CTRL run. If E_stick were the dominating driver, these two simulations would be very similar. The differences in the $\stackrel{\mathrm{‾}}{x}$ profiles of these two simulations can only be explained by relevant influences of other parameters on the aggregation rates. The simulations with the D-kernel (colMix2_Dkernel) exhibit a strong influence on aggregation. This is evident in the rapid decrease in Q_i and N_i and a rapid increase in ${\stackrel{\mathrm{‾}}{x}}_i$, ${\stackrel{\mathrm{‾}}{x}}_{\text{s}}$, and $\stackrel{\mathrm{‾}}{x}$ caused by high aggregation rates (supported by Appendix A). From this simulation, it is evident that the use of the new particle properties (including the Atlas-type v_t–size relation) together with the D-kernel results in even larger particles than in the default run, and thus DWRs are strongly overestimated (Fig. 7d). This overestimation can only be reduced by using the A-kernel.

The vertical gradients of Q result from mass uptake by depositional growth and divergence of v_t (Fig. 6h). First, Q increases from the cloud top to the cloud bottom due to depositional growth. Second, deposition growth and aggregation increase particle size and thus v_t increases. If there were no mass uptake (no deposit growth) but only aggregation, Q could only decrease because the product of v_t and Q would be conserved. The v_t–size relation plays an important role in these processes: on the one hand, smaller v_t for a given particle size, e.g., as in NEW vs. CTRL, means more time for mass uptake, leading to a faster increase in Q per height. On the other hand, smaller v_t could also lead to less ventilation and thus less mass uptake due to depositional growth. The v_t–size relationship, which defines the slope of v_t with increasing size, influences the divergence of v_t with height and the aggregation rates (Sect. 3.1.3). These multiple effects also interact, which further complicates the interpretation of the profiles of Q. Nevertheless, we attempt to interpret the most obvious features of the profiles of Q.

At about −17 ^∘C, MDV increases sharply in the NEW run (Fig. 7b), causing a decrease in Q at these temperatures (Fig. 6f), while Q increases continuously in the CTRL run. The differences in the profiles of Q between the sensitivity runs are relatively large. These large differences are likely due to the different conversion rates of cloud ice to snow and differently strong increasing $\stackrel{\mathrm{‾}}{x}$ near the model top. For example, in colMix2_Dkernel the cloud ice converts rapidly to larger snow particles. As a result, particles near the model top fall faster and therefore have less time to grow by depositional growth (the increase in Q is weaker compared to the NEW run); colMix2_Akernel_LinCot shows a weaker increase in $\stackrel{\mathrm{‾}}{x}$ for $T>-$15 ^∘C compared to the NEW run (Fig. 6i). This weaker increase in $\stackrel{\mathrm{‾}}{x}$ leads to a weaker increase in MDV (Fig. 7b) and thus to a stronger increase in Q (Fig. 6h). The reflectivity Ze_Ka is closely related to Q so that colMix2_Dkernel (colMix2_Akernel_LinCot) has the lowest (highest) reflectivity. However, the CTRL run has the highest Ze_Ka, although Q is lower than in some sensitivity runs. The large Ze_Ka here could be caused by the relatively dense snow particles assumed in CTRL (Fig. 3). Overall, Q and Ze_Ka show relatively large sensitivity to the varied parameters in these snowshaft simulations. However, this observation must be interpreted with caution. The simulations assume a relatively large humidity in order to match the hydrometeor profiles and compensate for processes not considered (Sect. 3.2.1). This high humidity could lead to an overestimation of mass uptake due to depositional growth. Additionally, considering that supersaturation is not consumed by depositional growth but is held constant in our snowshaft simulations, one could hypothesize that Q and Ze might be more similar among the sensitivity runs in the ICON-LEM simulations.

The new particle properties reduce the bias of the scheme regarding MDV to a large extent (Fig. 7b). While all simulations with the new particle properties are within the deciles of the observations, the standard run is already outside the deciles at −35 ^∘C and is more than 0.5 m s⁻¹ larger than the median at some temperatures (e.g., at T=5 ^∘C). The other parameters change the profile of the MDV to a much lesser extent. At temperatures from −18 to −12 ^∘C, all simulations show an increase in MDV, while all quantiles of the observed MDV decrease. This discrepancy could be due to the lack of habit prediction, underestimated or missing upwinds, or the lack of collisional fragmentation (Korolev and Leisner, 2020) in the model. At these temperatures dendritic growth occurs, which could lead to decreasing particle density and thus decreasing v_t and/or updrafts as a result of strong latent heat release. Collisional fragmentation could furthermore lead to the formation of new small particles with low v_t, which also reduces the MDV.

In addition to the particle properties, the width of the size distribution changes the MDV the most. The simulation with the wider size distribution (colMix2broad_Akernel) has a larger MDV (Fig. 7b) than the NEW run, which is due to the increasing number of large particles at the larger end of the distribution (Sect. 3.1.1). These large particles contribute more to the MDV than the smaller particles; to calculate MDV, each particle must be weighted by reflectivity, which for Rayleigh scatterers scales approximately with mass to the power of 2. The higher weight of the large participants also explains why the DWRs in colMix2broad_Akernel are significantly higher compared to the NEW run, even though the mean size of the hydrometeors is relatively similar. This sensitivity illustrates that the DWRs can only to some extent be used to infer $\stackrel{\mathrm{‾}}{x}$ and the size distribution width has to be additionally considered.

Despite the various simplifications in the snowshaft model (no nucleation, no advection, constant humidity) the mean profile of the radar profiles from the ICON-LEM simulations of O20 could be well-matched. This allowed us to investigate the sensitivity of aggregation to the individual model components and to select a model setup that best matches the observed radar profiles. The particle properties of the snow, the aggregation kernel formulation, and E_stick have a strong influence on the hydrometeor contents and the simulated radar observables. Interestingly, the choice of particle size distribution has little effect on the hydrometeor profiles but a large effect on the DWR values. The choice of cloud ice properties (needle or column) is less important than the choice of the other parameters in this cloud regime. However, the choice of cloud ice properties might be more important for clouds with smaller aggregation rates, e.g., cirrus. If we combine the A-kernel, the particle properties of Mix2 from K20, the newly proposed E_stick parameterization, and a relatively narrow size distribution the observed profiles of MDV and DWRs could be better matched. To test whether these sensitivities and improvements in NEW are also found persistently in more realistic simulations, in the next section we test whether these observations occur similarly in the ICON-LEM simulations.

Table 3 Overview of parameters and settings varied in the microphysical sensitivity experiments. The sensitivity runs have the same settings as colMix2_Akernel unless otherwise noted. K is the collision kernel, D the maximum dimension, and A the particle's projected area; μ and ν are parameters in the generalized gamma function describing the mass distribution in the microphysics scheme (Eq. 5).

.

Figure 6 Profiles of model variables in the snowshaft simulations. Number density N (a, d, g), mass mixing ratio Q (b, e, h), and mean mass $\stackrel{\mathrm{‾}}{x}$ (c, f, i) of the cloud ice (a, b, c), snow (d, e, f), and the sum of cloud ice and snow (g, h, i). Lines: simulations using different model settings as described in Table 3. Greyscale: histogram of the hydrometeor contents vs. temperature from the ICON-LEM simulations of the Tripex campaign (O20) filtered to include only profiles for which the precipitation rate exceeds 1 mm h⁻¹. The simulations in O20 used the default model settings.

Figure 7 (a) Reflectivity Ze_Ka, (b) mean Doppler velocity MDV_Ka, (c) DWR_Ka,W, and (d) DWR_X,Ka. Lines: simulated profiles based on snowshaft simulations (Fig. 6) as well as the median and quartiles of the observations. Greyscale: histogram of observations from the Tripex campaign (O20).

3.3 ICON-LEM case study simulation using the new parameterizations

In the snowshaft simulations (Sect. 3.2) we had to use several idealized assumptions. ICON-LEM (Sect. 2.2) contains additional processes (e.g., advection, nucleation, varying humidity field) and therefore simulates a more realistic representation of the atmosphere. In this section, we investigate the impact of the various parameters studied in the sensitivity analysis in a more complex case study with an ICON-LEM simulation. Furthermore, the ICON-LEM simulations provide an opportunity to extend the analysis to various conditions (e.g., nonstationary regime during the frontal passage, sublimation layers).

The case study of interest was 3 January 2016, when a low-pressure area over the British Isles and an accompanying frontal system over western and central Europe determined the synoptic situation over the modeled domain. Shallow mixed-phase clouds are present in the morning and dissipate around noon (Fig. 8a). The passage of a warm front manifests itself at 10:00 UTC, first in high clouds and then in sinking cloud bases. These frontal clouds start to precipitate at 18:00 UTC. The selected case is particularly interesting because it contains clouds in different regimes and precipitation of weak to moderate intensity.

The observed and simulated Ze_Ka fields match relatively well for all simulations in terms of cloud structure and precipitation (Fig. 8). Both the shallow mixed-phase clouds and the frontal cloud are very well-captured in terms of temporal and spatial structure.

Ze_Ka exhibits strong differences between the observations and the simulation only in the rain and ice slightly above the melting temperature in the period from 19:00 to 23:00 UTC. The sharp decrease in the observed Ze_Ka indicates strong sublimation. The presence of sublimation is also revealed by the model showing subsaturated air in this time range (Fig. B6). There are three main reasons that explain why the model does not accurately represent the sharp decrease in Ze_Ka in this sublimation scenario. First, the humidity could be overestimated in the model, e.g., due to inaccurate forcing data. Second, particle sizes could be overestimated due to processes in microphysics that weaken the effect of sublimation. We cannot completely rule out the humidity mismatch, but we found good agreement between the model and radiosonde data when available. Unfortunately, there was no radiosonde launched on the considered day. Thus, we are confident in the general ability of the model to accurately simulate the humidity field, but we cannot rule out the possibility that inaccuracies in the simulated humidity field contribute to the bias in Ze_Ka. Lastly, the parameterization of sublimation could also be an error source. For example, the evolution of the PSD during sublimation is challenging to represent in a two-moment scheme (Seifert, 2008). Since all of these reasons might be able to explain the mismatch in Ze_Ka, we should be cautious in assessing the validity of the assumptions of the individual model settings based on this sublimation feature. However, regardless of the accuracy of the model in predicting the humidity or simulating sublimation, the following differences in Ze_Ka of the model simulations underscore the importance of accurate prediction.

Figure 8 Time–height profile of Ze_Ka from 3 January 2016 as observed (a) and simulated (b–g) with various model settings (Table 3). Selected temperature isolines from CloudNet (Illingworth et al., 2007) for the observations (a) and the corresponding ICON-LEM output (b–g) are also shown.

While the NEW (Fig. 8c) and most sensitivity runs show a slight decrease in Ze_Ka due to sublimation in the time period when the air is subsaturated, the sublimation is barely seen in Ze_Ka of some other simulations (e.g., CTRL, colMix2_Dkernel; Fig. 8b and g). The differences between the simulations are caused by the differences in the particle size indicated by DWR_X,Ka (Fig. 9). Similar to the snowshaft simulations, DWR_X,Ka is strongly overestimated in colMix2_Dkernel and the CTRL run. In contrast, DWR_X,Ka is well-matched closely above the melting temperature in the NEW simulation. The hydrometeor populations with realistic particle sizes are more strongly affected by the subsaturated air and sublimate quickly, whereas particles that are too large sublimate less and therefore retain more mass. Thus, the overestimated particle size leads to overestimated precipitation. Between 18:00 and 24:00 UTC, 1.40 mm of accumulated rain was observed, 8.91 mm simulated by the default simulation and 2.29 mm by colMix2_Akernel. While this represents an overestimates of 536 % by the CTRL run during this time period, we emphasize the overall good agreement between modeled and observed precipitation reported by O20 for the entire campaign. While E_stick appeared to be important for the simulated DWR_Ka,W in the snowshaft simulations (Fig. 7), the differences between the simulation with the old (colMix2_Akernel_LinCot; Fig. 9d) and the new E_stick parameterization (colMix2_Akernel; Fig. 9c) are relatively small. In the ICON-LEM simulation, the weaker growth of the particles in colMix2_Akernel_LinCot at lower temperatures might be partly compensated for by advection or nucleation.

Figure 9 Same as Fig. 8 but displaying DWR_X,Ka, which is sensitive to mean mass diameters of 1 to 20 mm.

Besides the DWRs, MDV provides valuable information about the microphysical properties. As also reported by O20, MDV is overestimated in the SB06 default simulation, especially in regions where the particle sizes are overestimated (Fig. 10). MDV is often used to distinguish rimed from unrimed particles (e.g., Mosimann, 1995). Using this method, we detect some smaller episodes in which rimed particles dominate at about 04:00, 18:00, and 22:00 UTC. At other times, the observations indicate unrimed or only slightly rimed particles. In the SB06 default simulation, high MDVs are obtained in the whole time range after 18:00. Since the profiles of the hydrometeors show only very little mass of rimed particles during this period, the larger predicted MDV can be attributed to the overestimation of the unrimed snow particle v_t.

The new simulations, all using the new particle properties, have significantly lower values of MDV at all temperatures. This reduction of MDV compared to the SB06 default setup constitutes a significant reduction of the bias in MDV at temperatures below −10 ^∘C. For $T>-\mathrm{10}$ ^∘C, MDV is even slightly underestimated. Considering that Fig. 5 shows good agreement of MDV between the observations and the v_t–size relation of Mix2, we assume that the underestimation of MDV is not caused by the underestimation of the v_t–size relation of aggregates. Since DWR_X,Ka also matches well at these temperatures, processes other than aggregation and sedimentation of unrimed aggregates most probably cause this underestimation of MDV. One could speculate that riming rates are underestimated or that the vertical air motion is not well-simulated.

Figure 10 Same as Fig. 8 but displaying MDV_Ka, which is strongly linked to v_t.

Most of the findings from the snowshaft simulations (e.g., the strong reduction of MDV and DWR at temperatures close to the melting temperature) are confirmed by the ICON-LEM simulation of this case study. However, the ICON-LEM simulations reveal that the influence of E_stick seems to be overestimated in the snowshaft simulations. Moreover, accurate modeling of particle sizes and v_t in the presence of a sublimating layer is critical. The simulations with the new particle properties showed a slight underestimation of the MDV. This underestimation most likely does not arise from an inaccurate representation of the particle properties or the aggregation rates but is caused by another process (e.g., riming, vertical air motion). In previous analyses of the SB06 default setup, this underestimation could not have been detected because it was masked by the overestimation of the aggregate's v_t. Because errors can be specific to the chosen day, such as a particular mismatch of the relative humidity, relying on only one case to detect a discrepancy in the microphysical properties is prone to error. Therefore, we analyze the statistics of a multi-month simulation in the next section.

3.4 Statistical comparison

After evaluating the choices of the new scheme in the snowshaft model and in a case study with ICON-LEM, we perform ICON-LEM simulations for the entire Tripex time period. By comparing observed and modeled histograms of DWR and MDV as a function of temperature, we can evaluate the new scheme. Since we additionally contrast the histograms of the NEW and CTRL simulations, we can test whether the reduction in the bias of DWRs and MDV found in Sect. 3.3 is specific to the selected case or rather a consistent feature of the model changes. As DWRs are related to the mean particle size, we can assess whether the chosen parameter combination can accurately simulate aggregation in various weather situations present in the simulated days. The same applies to MDV profiles, which are especially valuable in evaluating the suitability of the assumed v_t–size relationship.

The observed and synthetic radar profiles are filtered in the same way for comparability. For example, the first 6 h of simulation and observation are not considered because the model output could contain artifacts during this spin-up time. Moreover, we include only profiles in which the rain rate RR exceeds 1 mm h⁻¹. The latter filter enables us to focus on the most relevant cases for precipitation. Interestingly, O20 found the discrepancy between the model and observations to be especially obvious for these profiles. For a detailed description of the processing, we refer to O20.

To quantify the agreement between the histograms of the simulations and the observation, the Hellinger distance H is used. H can be defined for two distributions P =(p₁, …, p_k) and Q =(q₁, …, q_k) as

$$\begin{array}{}\text{(17)}& H(P,Q)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\sqrt{\sum _{i=\mathrm{1}}^k{\left(\sqrt{p_i}-\sqrt{q_i}\right)}^{\mathrm{2}}}.\end{array}$$

H is zero for two identical distributions and 1 if the distributions do not overlap at all.

The medians and larger quantiles of the observed distributions of DWRs indicate a strong increase in particle size around −15 ^∘C (most evident in DWR_Ka,W; Fig. 11a) and just above the melting temperature (most evident in DWR_X,Ka; Fig. 11e). Both of these characteristic increases in the particle sizes are found to some extent in CTRL (panels b and f in Fig. 11) and NEW (panels c and g in Fig. 11). The increase in particle sizes between −15 and −10 ^∘C happens in the new simulations at slightly lower temperatures, and the different profiles reveal a greater variability (visible, e.g., in the difference of DWR_Ka,W between the lower and upper decile). H indicates a slightly better match by CTRL in this temperature range. For $T>-\mathrm{10}$ ^∘C the mean and higher quantiles of DWR_X,Ka increase very rapidly in CTRL and more slowly in NEW and the observation. The increase in particle sizes as simulated by NEW is in much better agreement with the observed profiles. The upper quartile of DWR_X,Ka only slightly exceeds 5 dB in the observations and NEW but is higher than 10 dB in CTRL for $T>-\mathrm{1}$ ^∘C. This better match is also indicated by H (Fig. 11h), which is about 5 times larger for CTRL compared to NEW.

Figure 11 Contour frequency by temperature diagrams (CFTDs) for all profiles with RR > 1 mm h⁻¹ of the dual-wavelength ratios between the X and Ka band (DWR_X,Ka, top) and between the Ka and W band (DWR_X,Ka, bottom) from the default simulation (a, e), the new simulation (colMix2_Akernel; b and f), and measured (c, g). The black lines represent the statistical measures (median, mean, quartiles, and deciles) at different temperatures. Panels (d) and (h) show the Hellinger distance between the simulated and observed distributions for all temperatures.

Besides the overestimation of DWR_X,Ka closely above the melting temperature, O20 also highlighted the overestimation of MDV by CTRL. This overestimation is present at all temperatures (compare panels a and b in Fig. 12) and can be attributed to the overestimation of the v_t–size relationship of the snow class as reported in Karrer et al. (2020) and the overestimated particle sizes for the higher temperatures. The overestimation of MDV by CTRL is most pronounced for $T>-\mathrm{15}{}^{\circ }C$. In this temperature range, CTRL cannot reproduce the asymptotic approach because of the power-law v_t–size relationship (Sect. 3.1.3). For example, the median of MDV_Ka at −5 ^∘C is 1 m s⁻¹ in the observations and 1.3 m s⁻¹ in CTRL. In contrast, the new simulations agree better with the observations and H is about half as large as for CTRL. The new scheme setup is more accurate in this temperature range because the Atlas-type v_t–size relationship of the Mix2 particles (Fig. 3) correctly considers the asymptotic approach to 1 m s⁻¹ at large sizes. However, MDV is slightly underestimated by NEW for $T>-\mathrm{10}$ ^∘C. Values substantially above 1 m s⁻¹ occur in the observations and the new simulations only closely above the melting temperature, where rain is present. At temperatures below −15 ^∘C, both simulations perform similarly, with H ranging from 0.2 to 0.5. While CTRL exhibits a continuous overestimation of MDV, the new simulations lack the observed increase in MDV for $T<-\mathrm{20}$ ^∘C. At these temperatures, the selected PSD width (Sect.3.1.1) and cloud ice particle properties (Sect. 3.1.3) may not be ideal.

The statistical comparison shows that the changes we made to the model could eliminate the most striking biases, namely the overestimation of DWR_X,Ka and MDV closely above the melting temperatures. The match of these quantities is important for accurate simulation of precipitation, as exemplified in the case study in Sect. 3.3. Some discrepancies remain, namely the overly strong increase in the DWRs at temperatures between −15 and −10 ^∘C and the overestimation (underestimation) of MDV temperatures below −25 ^∘C (above −10 ^∘C). These discrepancies can be caused by several model errors (inaccurate simulation of, e.g., PSD shape, E_stick, degree of riming, variability in cloud ice properties) that cannot be fully deciphered by this observational setup and could benefit from advances in laboratory measurements, observational setup, and representation of cloud ice habits and riming degree in the model.

Figure 12 Contour frequency by temperature diagrams (CFTDs) of the mean Doppler velocity of the Ka band (MDV_Ka) from the default simulation (a, e), the new simulation (colMix2_Akernel, b and f), and measured (c, g). The black lines represent the statistical measures (median, mean, quartiles, and deciles) at different temperatures. The histograms on top are calculated including all data and on the bottom only data from profiles for which the precipitation rate RR exceeds 1 mm h⁻¹. A vertical line at 1 m s⁻¹ eases the comparison of the different distributions. Panels (d) and (h) show the Hellinger distance between the simulated and observed distributions for all temperatures.

4 Conclusions

Aggregation is a key ice microphysical growth process for the formation of precipitating ice particles, which are the precursor of raindrops in cold rain formation. Recent studies using statistics from multi-frequency Doppler radar observations provided observational constraints on how critical radar quantities, such as DWRs or MDV, change with temperature. In this study, we aimed at a deeper analysis of the underlying causes for the observed discrepancies between radar statistics and a state-of-the-art two-moment microphysical scheme, and we improved its simulation of aggregation.

To this end, as a first step, we revisited all relevant components of aggregation as considered in the two-moment scheme to see how well they represent current knowledge of physics. These components are the size distribution width, the temperature dependence of E_stick, the particle properties (with a focus on the v_t–size relation of aggregates), and the representation of nonspherical particles in the aggregation kernel formulation.

To systematically test the sensitivities of various parameter combinations, we performed 1D simulations with the snowshaft model, which uses simple profiles of thermodynamic variables and a simple initialization of particles at the model top. Moreover, the model only accounts for a subset of all the microphysical processes that occur in real clouds. Nevertheless, by adjusting the model setup we could match the average profiles of radar observables obtained by the 3D simulations of O20, which used the SB06 default scheme setup.

The snowshaft simulations revealed high sensitivity of aggregation to particle properties, the aggregation kernel formulation, and E_stick. Surprisingly, the size distribution width had a relatively small effect on the modeled mean mass but a considerable influence on the simulated DWRs. The influence of the cloud ice properties was small in both the model and radar variables.

By comparing the profiles from the snowshaft simulation with the average observed profiles, we were able to select a set of parameters that provided the best agreement with the observations. In this selection process, we mainly varied the less-known components, e.g., the size distribution width, and held other parameters constant that we could better constrain, e.g., the v_t–size relationship. The size distribution width proved to be a critical component in linking modeled $\stackrel{\mathrm{‾}}{x}$ to observed DWRs and at the same time is difficult to constrain with the given observational setup. Therefore, using a microphysical scheme that explicitly simulates the width of the size distribution (e.g., a three-moment scheme) would provide a more consistent link between the model and observation. However, additional observational constraints from radar (e.g., Doppler spectrum width) and in situ observations should be considered in this case. In particular, we find that the v_t–size relationship, which accounts for the asymptotic behavior of v_t at large sizes, leads to better agreement with the observations. Moreover, the A-kernel appears to be a better approximation of the aggregation kernel when combined with a constant E_coll.

We implemented this improved scheme setup in the ICON-LEM and also tested the individual model modifications in a case study. These more realistic ICON-LEM simulations allowed us to derive potential differences in the analysis of sensitivities compared to the snowshaft simulations, possibly caused by effects such as dynamics and advection. Overall, the ICON-LEM simulations yielded similar sensitivities as the snowshaft simulations, but slight differences were apparent with respect to sensitivity to E_stick. The differences between simulations with different E_stick parameterization were less pronounced in the ICON-LEM simulations. This discrepancy between the simulation frameworks could result from accounting for feedback from microphysics to model humidity in the ICON-LEM simulation.

On the day considered in the case study, relatively dry low-level air resulted in strong sublimation of ice particles. This sublimation feature demonstrated the relevancy of accurately simulating $\stackrel{\mathrm{‾}}{x}$. The SB06 default scheme with its largely overestimated aggregate sizes strongly overestimated the rainfall rate on the ground because the large snowflakes could not sublimate fast enough. In contrast, the more realistic aggregate sizes obtained with the new scheme were able to fit the observations much better.

Finally, the entire period of the campaign dataset (46 d) was simulated again with ICON-LEM using the best-matching parameter combination from the previous tests. This allowed us to directly compare the new statistics with the previous analysis of the default scheme provided in O20. The new aggregation formulation is clearly able to reduce the observed overestimation of MDV. This improvement can be attributed to the new Atlas-type v_t–size relationship. The overestimation of the mean particle size at high temperatures revealed in the DWRs was also substantially reduced by the new aggregation parameterization.

Remaining discrepancies are found for DWRs at temperatures of about −12 ^∘C and for MDV at low and high temperatures. The overestimated DWRs by the new simulations could result in an overestimated $\stackrel{\mathrm{‾}}{x}$ or overly broad size distribution in the model. Inclusion of a higher-frequency radar (Battaglia et al., 2014) may help to infer the particle growth above −12 ^∘C. The analysis of Doppler spectra (e.g., similar to Barrett et al., 2019) or observational techniques, e.g., in situ probing of the particle size distribution, would provide additional constraints on the size distribution and ease the interpretation of the MDV and DWR. The mismatch of the MDV at lower temperatures could be caused by an inaccurate size distribution width, as well as E_stick or cloud ice properties. Future studies could focus on this temperature region, which is highly relevant for cloud radiative effects. The slight underestimation of MDV at high temperatures could be due to underestimated riming rates, the representation of partially rimed particles, or other effects such as vertical air motion. Further insight could be gained, e.g., from the analysis of the Doppler spectra or comparison with other microphysical schemes with a different representation of the riming process (Morrison and Milbrandt, 2015; Tsai and Chen, 2020).

In addition to the results obtained in this study for aggregation in the SB06 scheme, we think that our approach for how to utilize state-of-the-art radar datasets to improve parameterizations may also serve as a blueprint for future studies focusing on other processes or microphysical schemes. Therefore, we shortly summarize the approach in general terms with the following points.

1. Revisit components of the physical parameterization.

2. Set up single-column simulations which match the average profiles of simulated observables obtained from long-term 3D simulations with the default scheme setup.

3. Systematically test the sensitivities of various parameter combinations in 1D simulations.

4. Select the model configuration that best matches the observations.

5. Implement model modifications in the 3D model and infer possible differences in sensitivities between a 3D simulation and 1D simulations in a case study.

6. Rerun the long-term 3D simulation using the best-matching parameter combination and investigate the improvements by comparing observations with simulations using the default and the new scheme setup.

Appendix A: Bulk aggregation rates

We summarize the bulk aggregation formulas for all aggregation processes: ice–snow collection, ice self-collection, and snow self-collection. While the formulations using the D-kernel were already given by Seifert et al. (2014), the formulas using the A-kernel were newly derived in this study.

Combining the definition of the moments,

$$\begin{array}{}\text{(A1)}& M_n=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{m}^{n}f\left(m\right)dm,\end{array}$$

the SCE (Eq. 4) and its simplifications in the SB06 scheme (Sect. 3.1), an equation can be derived that allows for the calculation of all relevant aggregation rates between particles of the classes i and j:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& {\displaystyle \frac{\partial M_{i,n}}{\partial t}}{|}_{\text{coll,ij}}=\mathrm{\Phi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{f}_i\\ & \left(D_i\right)f_j\left(D_j\right)K_{i,j}(D_i,D_j)m_i^{n}d{D}_{j}d{D}_i,\end{array}\end{array}$$

where M_j,n is the nth moment of the hydrometeor class j, f is the particle size distribution for a selected size variable (D_max, D_eq, or m), K is the aggregation kernel, and m is the particle mass.

Seifert et al. (2014) use the variance approach proposed in Seifert and Beheng (2006), which parameterizes the bulk velocity difference by the square root of the second moment of the velocity differences. In this way, the integral is separated into a term containing the geometrical properties (𝒞_n,ij) and a part which contains the velocity difference (${\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ij}$) to enable the analytical integration.

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\partial M_{i,n}}{\partial t}}{|}_{\text{coll},ij}={\stackrel{\mathrm{‾}}{E}}_{i,j}{\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ij}{\mathcal{C}}_{n,ij}\end{array}$$

The expressions of 𝒞_n,ij and ${\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ij}$ depend on the expression of the PSD_m (Sect. 3.1.1), the formulation of the aggregation kernel (Sect. 3.1.2), and the particle properties (Sect. 3.1.3). The SB06 scheme assumes a modified gamma distribution as a function of mass (Eq. 5), which can be easily converted to a gamma distribution as a function of D_eq if ${\mathit{\mu }}_m=\mathrm{1}/\mathrm{3}$ (Eq. 7). The particle properties are characterized by power-law relations of m (Eq. 12) and A_act (Eq. 15) vs. D_max and D_eq. In the new scheme, v_t of cloud ice and snow is parameterized by an Atlas-type relation as a function of D_eq (Eq. 11). Coefficients of the relations can be found in Table 2.

A1 D-kernel

Inserting the D-kernel (Eq. 1) into Eq. (A2), the 𝒞_n,ij and ${\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ij}$ can be written as

$$\begin{array}{}\text{(A4)}& {\displaystyle }\begin{array}{rl}& {\mathcal{C}}_{n,ij}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(D_{max,i}+D_{max,j}\right)}^{\mathrm{2}}{f}_i\left(m_i\right)\\ & f_j\left(m_j\right)m_j^{n}d{m}_{i}d{m}_j,\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ij}=\mathit{\{}{\displaystyle \frac{\mathrm{1}}{{\mathcal{N}}_{n,ij}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left[v_i\left(D_{\text{eq},i}\right)-v_j\left(D_{\text{eq},j}\right)\right]}^{\mathrm{2}}\\ & \times D_{\text{eq},i}^{\mathrm{2}}{D}_{\text{eq},j}^{\mathrm{2}}\\ & f_{\text{eq},i}\left(D_{\text{eq},i}\right)f_{\text{eq},j}\left(D_{\text{eq},j}\right)m_i^{n}d{D}_{\text{eq},i}d{D}_{\text{eq},j}{\mathit{\}}}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

where 𝒩 is the normalization factor given by

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& {\mathcal{N}}_{n,ij}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{D}_{\text{eq},i}^{\mathrm{2}}{D}_{\text{eq},j}^{\mathrm{2}}{f}_{\text{eq},i}\left(D_{\text{eq},i}\right)f_{\text{eq},j}\left(D_{\text{eq},j}\right)\\ & m_i^{n}d{D}_{i}d{D}_j.\end{array}\end{array}$$

Inserting the D_max−m relation,

$$\begin{array}{}\text{(A7)}& D_{max,i}=a_i{m}_i^{b_i}={\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{w}}{a}_i}{\mathrm{6}}}{D}_{\text{eq},i}^{\mathrm{3}{b}_i},\end{array}$$

and the PSD_m (Eq. 5) into 𝒞_n,ij (Eq. A4) and solving the integral, we obtain

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}& {\mathcal{C}}_{n,ij}={\left({\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}{\mathrm{6}}}\right)}^n{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{N}_i{N}_j\\ & [{\mathit{\delta }}_{D,i}^{\mathrm{0}}{\stackrel{\mathrm{‾}}{D}}_i^{\mathrm{2}}+{\mathit{\delta }}_{D,ij}^n{\stackrel{\mathrm{‾}}{D}}_i{\stackrel{\mathrm{‾}}{D}}_j+{\mathit{\delta }}_j^n{\stackrel{\mathrm{‾}}{D}}_j^{\mathrm{2}}],\end{array}\end{array}$$

where ${\mathit{\delta }}_i^n$ and ${\mathit{\delta }}_j^n$ are equal to ${\mathit{\delta }}_p^{\mathrm{0}}$ of Eq. (90) of SB2006 and ${\mathit{\delta }}_{ij}^n$ is equal to ${\mathit{\delta }}_g^{\mathrm{0}}$ of Eq. (91) of SB2006.

$$\begin{array}{}\text{(A9)}& {\displaystyle }{\mathit{\delta }}_{D,i}^n={\displaystyle \frac{\mathrm{\Gamma }\left(\right(\mathrm{2}{b}_i+{\mathit{\nu }}_{m,i}+\mathrm{1}+n)/{\mathit{\mu }}_{m,i})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i})}}{\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{2})/{\mathit{\mu }}_{m,i})}}\right]}^{\mathrm{2}{b}_i+n}\text{(A10)}& {\displaystyle }\begin{array}{rl}& {\mathit{\delta }}_{D,ij}^n=\mathrm{2}{\displaystyle \frac{\mathrm{\Gamma }\left(\right(b_i+{\mathit{\nu }}_{m,i}+\mathrm{1}+n)/{\mathit{\mu }}_{m,i})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i})}}{\displaystyle \frac{\mathrm{\Gamma }\left(\right(b_j+{\mathit{\nu }}_{m,j}+\mathrm{1})/{\mathit{\mu }}_{m,j})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,j}+\mathrm{1})/{\mathit{\mu }}_{m,j})}}\\ & \times {\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{2})/{\mathit{\mu }}_{m,i})}}\right]}^{b_i+n}{\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,j}+\mathrm{1})/{\mathit{\mu }}_{m,j})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,j}+\mathrm{2})/{\mathit{\mu }}_{m,j})}}\right]}^{b_j}\end{array}\end{array}$$

Inserting the velocity relation (Eq. 11) and the size distribution using D_eq (Eq. 7) into the velocity variance (Eq. A4) and solving the integral, we obtain

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ij}=\left[\right({\mathit{\alpha }}_{v,j}-{\mathit{\alpha }}_{v,i}{)}^{\mathrm{2}}-\mathrm{2}{\mathit{\beta }}_{v,j}({\mathit{\alpha }}_{v,j}-{\mathit{\alpha }}_{v,i}){\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,j}}{{\mathit{\lambda }}_{\text{eq},j}}}\right)}^{-{\mathit{\xi }}_{D,i}^n}\\ & -\mathrm{2}{\mathit{\beta }}_{v,i}({\mathit{\alpha }}_{v,i}-{\mathit{\alpha }}_{v,j}){\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-{\mathit{\xi }}_{D,i}^n}+{\mathit{\beta }}_{v,j}^{\mathrm{2}}{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_{v,j}}{{\mathit{\lambda }}_{\text{eq},j}}}\right)}^{-{\mathit{\xi }}_{D,j}}\\ & +{\mathit{\beta }}_{v,i}^{\mathrm{2}}{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-{\mathit{\xi }}_{D,i}^n}-\mathrm{2}{\mathit{\beta }}_{v,j}{\mathit{\beta }}_{v,i}{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,j}}{{\mathit{\lambda }}_{\text{eq},j}}}\right)}^{-{\mathit{\xi }}_{D,j}}\\ & \times {\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-{\mathit{\xi }}_{D,i}^n}{]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

with

$$\begin{array}{}\text{(A12)}& {\mathit{\xi }}_{D,i}^n={\mathit{\mu }}_{\text{eq},i}+\mathrm{3}+\mathrm{3}n,{\mathit{\xi }}_{D,j}={\mathit{\mu }}_{\text{eq},j}+\mathrm{3}.\end{array}$$

A2 A-kernel

Inserting the A-kernel (Eq. 9) into Eq. (A2), the velocity variance and the geometric part of the bulk collision rates can be written as

$$\begin{array}{}\text{(A13)}& {\displaystyle }{\mathcal{C}}_{n,ij}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left(A_i^{\mathrm{0.5}}+A_j^{\mathrm{0.5}}\right)}^{\mathrm{2}}{f}_i\left(D_i\right)f_j\left(D_j\right)m_j^{n}d{D}_{i}d{D}_j\text{(A14)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ij}=\mathit{\{}{\displaystyle \frac{\mathrm{1}}{{\mathcal{N}}_{n,ij}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\left[v_i\left(D_{\text{eq},i}\right)-v_j\left(D_{\text{eq},j}\right)\right]}^{\mathrm{2}}\\ & \times D_{\text{eq},i}^{{\mathit{\sigma }}_{A,i}}{D}_{\text{eq},j}^{{\mathit{\sigma }}_{A,j}},\\ & f_{\text{eq},i}\left(D_{\text{eq},i}\right)f_{\text{eq},j}\left(D_{\text{eq},j}\right)m_i^{n}d{D}_{\text{eq},i}d{D}_{\text{eq},j}{\mathit{\}}}^{\frac{\mathrm{1}}{\mathrm{2}}}\end{array}\text{(A15)}& {\displaystyle }\begin{array}{rl}& {\mathcal{N}}_{n,ij}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{D}_{\text{eq},i}^{{\mathit{\sigma }}_{A,i}}{D}_{\text{eq},j}^{{\mathit{\sigma }}_{A,j}}{f}_{\text{eq},i}\left(D_{\text{eq},i}\right)f_{\text{eq},j}\left(D_{\text{eq},j}\right)\\ & m_i^{n}d{D}_{\text{eq},i}d{D}_{\text{eq},j}.\end{array}\end{array}$$

Inserting the A–D_eq relation (Eq. 15) and the size distribution as a function of D_eq (Eq. 7) into the geometric part (Eq. A13) and solving the integral leads to

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}& {\mathcal{C}}_{n,ij}={\left({\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}{\mathrm{6}}}\right)}^n{N}_i{N}_j\\ & [{\mathit{\delta }}_{A,i}^n{\stackrel{\mathrm{‾}}{D}}_{max,i}^{{\mathit{\sigma }}_{A,i}^{*}}+{\mathit{\delta }}_{A,ij}^n{\stackrel{\mathrm{‾}}{D}}_{max,i}^{{\mathit{\sigma }}_{A,i}^{*}/\mathrm{2}}{\stackrel{\mathrm{‾}}{D}}_{max,j}^{{\mathit{\sigma }}_{A,j}^{*}/\mathrm{2}}+{\mathit{\delta }}_{A,j}^n{\stackrel{\mathrm{‾}}{D}}_{max,j}^{{\mathit{\sigma }}_{A,j}^{*}}],\end{array}\end{array}$$

with the following.

$$\begin{array}{}\text{(A17)}& {\displaystyle }{\mathit{\delta }}_{A,i}^n={\mathit{\gamma }}_{A,i}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+{\mathit{\sigma }}_{A,i}+\mathrm{1}+\mathrm{3}n)}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1})}}{c}_{\mathit{\lambda },i}^{{\mathit{\sigma }}_{A,i}+\mathrm{3}n}\text{(A18)}& {\displaystyle }\begin{array}{rl}& {\mathit{\delta }}_{A,ij}^n=\mathrm{2}({\mathit{\gamma }}_{A,i}{\mathit{\gamma }}_{A,j}{)}^{\mathrm{0.5}}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+{\mathit{\sigma }}_{A,i}/\mathrm{2}+\mathrm{1}+\mathrm{3}n)}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1})}}{c}_{\mathit{\lambda },i}^{{\mathit{\sigma }}_{A,i}/\mathrm{2}+\mathrm{3}n}\\ & \times {\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},j}+{\mathit{\sigma }}_{A,j}/\mathrm{2}+\mathrm{1})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},j}+\mathrm{1})}}{c}_{\mathit{\lambda },j}^{{\mathit{\sigma }}_{A,j}/\mathrm{2}}\end{array}\text{(A19)}& {\displaystyle }\begin{array}{rl}& {\mathit{\delta }}_{A,j}^n={\mathit{\gamma }}_{A,j}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},j}+{\mathit{\sigma }}_{A,j}+\mathrm{1})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},j}+\mathrm{1})}}{c}_{\mathit{\lambda },j}^{{\mathit{\sigma }}_{A,j}}\\ & \times {\left[{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{4})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1})}}\right]}^n{c}_{\mathit{\lambda },i}^{\mathrm{3}n}\end{array}\text{(A20)}& {\displaystyle }{\mathit{\sigma }}_{A,i}^{*}={\displaystyle \frac{b_{m,i}{\mathit{\sigma }}_{A,i}}{\mathrm{3}}}\text{(A21)}& {\displaystyle }{c}_{\mathit{\lambda },i}={\left[{\displaystyle \frac{\mathrm{6}{a}_{m,i}}{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{4})}}\right]}^{\mathrm{1}/\mathrm{3}}\end{array}$$

Inserting the velocity relation (Eq. 11) and the size distribution as a function of D_eq (Eq. 7) into the velocity variance (Eq. A14) and solving the integral, we obtain

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ij}=\left[\right({\mathit{\alpha }}_{v,j}-{\mathit{\alpha }}_{v,i}{)}^{\mathrm{2}}-\mathrm{2}{\mathit{\beta }}_{v,j}({\mathit{\alpha }}_{v,j}-{\mathit{\alpha }}_{v,i}){\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,j}}{{\mathit{\lambda }}_{\text{eq},j}}}\right)}^{-{\mathit{\xi }}_{A,j}}\\ & -\mathrm{2}{\mathit{\beta }}_{v,i}({\mathit{\alpha }}_{v,i}-{\mathit{\alpha }}_{v,j}){\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-{\mathit{\xi }}_{A,i}^n}+{\mathit{\beta }}_{v,j}^{\mathrm{2}}{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_{v,j}}{{\mathit{\lambda }}_{\text{eq},j}}}\right)}^{-{\mathit{\xi }}_{A,j}}\\ & +{\mathit{\beta }}_{v,i}^{\mathrm{2}}{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-{\mathit{\xi }}_{A,i}^n}-\mathrm{2}{\mathit{\beta }}_{v,j}{\mathit{\beta }}_{v,i}{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,j}}{{\mathit{\lambda }}_{\text{eq},j}}}\right)}^{-{\mathit{\xi }}_{A,j}}\\ & \times {\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-{\mathit{\xi }}_{A,i}^n}{]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

with

$$\begin{array}{}\text{(A23)}& {\displaystyle }{\mathit{\xi }}_{i,A}^n={\mathit{\mu }}_{\text{eq},i}+{\mathit{\sigma }}_{A,i}+\mathrm{1}+\mathrm{3}n,\text{(A24)}& {\displaystyle }{\mathit{\xi }}_{j,A}={\mathit{\mu }}_{\text{eq},j}+{\mathit{\sigma }}_{A,j}+\mathrm{1}.\end{array}$$

A3 Ice self-collection

A3.1 D-kernel

For ice self-collection the geometry part (Eq. A13) simplifies to

$$\begin{array}{}\text{(A25)}& {\mathcal{C}}_{n,ii}={\left({\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}{\mathrm{6}}}\right)}^n{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{N}_i^{\mathrm{2}}[\mathrm{2}{\mathit{\delta }}_{D,i}^{\mathrm{0}}+{\mathit{\delta }}_{D,ii}^n]{\stackrel{\mathrm{‾}}{D}}_i^{\mathrm{2}},\end{array}$$

where ${\mathit{\delta }}_i^n$ is equal to ${\mathit{\delta }}_p^{\mathrm{0}}$ of Eq. (90) of SB2006 and ${\mathit{\delta }}_{ii}^n$ is equal to ${\mathit{\delta }}_g^{\mathrm{0}}$ of Eq. (91) of SB2006.

$$\begin{array}{}\text{(A26)}& {\displaystyle }\begin{array}{rl}& {\mathit{\delta }}_{D,i}^n={\displaystyle \frac{\mathrm{\Gamma }\left(\right(\mathrm{2}{b}_i+{\mathit{\nu }}_{m,i}+\mathrm{1}+n)/{\mathit{\mu }}_{m,i})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i})}}\\ & {\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{2})/{\mathit{\mu }}_{m,i})}}\right]}^{\mathrm{2}{b}_i+n}\end{array}\text{(A27)}& {\displaystyle }\begin{array}{rl}& {\mathit{\delta }}_{D,ii}^n=\mathrm{2}{\displaystyle \frac{\mathrm{\Gamma }\left(\right(b_i+{\mathit{\nu }}_{m,i}+\mathrm{1}+n)/{\mathit{\mu }}_{m,i})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i}{)}^{\mathrm{2}}}}\\ & \mathrm{\Gamma }\left(\right(b_i+{\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i})\\ & \times {\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{1})/{\mathit{\mu }}_{m,i})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,i}+\mathrm{2})/{\mathit{\mu }}_{m,i})}}\right]}^{\mathrm{2}{b}_i+n}\end{array}\end{array}$$

The velocity variance simplifies to

$$\begin{array}{}\text{(A28)}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ii}={\mathit{\beta }}_{v,i}\sqrt{\mathrm{2}}[{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-{\mathit{\xi }}_{D,i}^n}-{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-\mathrm{2}{\mathit{\xi }}_{D,i}^n}{]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

with

$$\begin{array}{}\text{(A29)}& {\mathit{\xi }}_{D,i}^n={\mathit{\mu }}_{\text{eq},i}+\mathrm{3}+\mathrm{3}n.\end{array}$$

A3.2 A-kernel

For ice self-collection 𝒞 (Eq. A13) simplifies to

$$\begin{array}{}\text{(A30)}& \begin{array}{rl}& {\mathcal{C}}_{n,ii}={\left({\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}{\mathrm{6}}}\right)}^n{N}_i^{\mathrm{2}}[{\mathit{\delta }}_{A,i}^n+{\mathit{\delta }}_{A,ii}^n+{\mathit{\delta }}_{A,i\mathrm{2}}^n]{\left({\displaystyle \frac{\mathrm{6}{a}_{m,i}}{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}}\right)}^{\frac{{\mathit{\sigma }}_{A,i}}{\mathrm{3}}}\\ & min\left({\mathit{\gamma }}_{A,i}{\left({\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}{\mathrm{6}{a}_{m,i}}}\right)}^{\frac{{\mathit{\sigma }}_{A,i}}{\mathrm{3}}}{\stackrel{\mathrm{‾}}{D}}_{max,i}^{{\mathit{\sigma }}_{A,i}^{*}},{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{\stackrel{\mathrm{‾}}{D}}_{max,i}^{\mathrm{2}}\right),\end{array}\end{array}$$

with the following.

$$\begin{array}{}\text{(A31)}& {\displaystyle }{\mathit{\delta }}_{A,i}^n={\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+{\mathit{\sigma }}_{A,i}+\mathrm{1}+\mathrm{3}n)}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1})}}{c}_{\mathit{\lambda },i}^{{\mathit{\sigma }}_{A,i}+\mathrm{3}n}\text{(A32)}& {\displaystyle }\begin{array}{rl}& {\mathit{\delta }}_{A,ii}^n=\mathrm{2}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+{\mathit{\sigma }}_{A,i}/\mathrm{2}+\mathrm{1}+\mathrm{3}n)}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1}{)}^{\mathrm{2}}}}\mathrm{\Gamma }\\ & ({\mathit{\mu }}_{\text{eq},i}+{\mathit{\sigma }}_{A,i}/\mathrm{2}+\mathrm{1})c_{\mathit{\lambda },i}^{{\mathit{\sigma }}_{A,i}+\mathrm{3}n}\end{array}\text{(A33)}& {\displaystyle }{\mathit{\delta }}_{A,i\mathrm{2}}^n={\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+{\mathit{\sigma }}_{A,i}+\mathrm{1})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1})}}{c}_{\mathit{\lambda },i}^{{\mathit{\sigma }}_{A,i}+\mathrm{3}n}{\left[{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{4})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1})}}\right]}^n\text{(A34)}& {\displaystyle }{\mathit{\sigma }}_{A,i}^{*}={\displaystyle \frac{b_{m,i}{\mathit{\sigma }}_{A,i}}{\mathrm{3}}}\text{(A35)}& {\displaystyle }{c}_{\mathit{\lambda },i}={\left[{\displaystyle \frac{\mathrm{6}{a}_{m,i}}{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{1})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq},i}+\mathrm{4})}}\right]}^{\mathrm{1}/\mathrm{3}}\end{array}$$

For small sizes, the parametrization of A_act yields values of A_r larger than 1 (e.g., columns smaller than $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{5}}$; Fig. 3d). For small mean sizes, these particles with unphysical A_r can substantially contribute to 𝒞_n,ii. Therefore, we limit A_act to A_sphere in Eq. (A30). The effect of this limiter can be seen in the kink of the bulk collision rates (Fig. B3c and d).

Inserting the velocity relation (Eq. 11) and the size distribution using D_eq (Eq. 7) into the velocity variance (Eq. A14) and solving the integral, we find

$$\begin{array}{}\text{(A36)}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{n,ii}=\sqrt{\mathrm{2}}{\mathit{\beta }}_{v,i}[{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-{\mathit{\xi }}_{A,i}^n}-{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,i}}{{\mathit{\lambda }}_{\text{eq},i}}}\right)}^{-\mathrm{2}{\mathit{\xi }}_{A,i}^n}{]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

with

$$\begin{array}{}\text{(A37)}& {\mathit{\xi }}_{i,A}^n={\mathit{\mu }}_{\text{eq},i}+{\mathit{\sigma }}_{A,i}+\mathrm{1}+\mathrm{3}n.\end{array}$$

A4 Snow self-collection

A4.1 D-kernel

For snow self-collection only the first moment is relevant and 𝒞 (Eq. A13) simplifies to

$$\begin{array}{}\text{(A38)}& {\mathcal{C}}_{\mathrm{0},\mathrm{ss}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{N}_{\text{s}}^{\mathrm{2}}[\mathrm{2}{\mathit{\delta }}_{D,s}^{\mathrm{0}}+{\mathit{\delta }}_{D,ss}^{\mathrm{0}}]{\stackrel{\mathrm{‾}}{D}}_{\text{s}}^{\mathrm{2}},\end{array}$$

where ${\mathit{\delta }}_{\text{s}}^n$ is equal to ${\mathit{\delta }}_p^{\mathrm{0}}$ of Eq. (90) of SB2006 and ${\mathit{\delta }}_{\text{ss}}^n$ is equal to ${\mathit{\delta }}_g^{\mathrm{0}}$ of Eq. (91) of SB2006.

$$\begin{array}{}\text{(A39)}& {\displaystyle }{\mathit{\delta }}_{D,s}^{\mathrm{0}}={\displaystyle \frac{\mathrm{\Gamma }\left(\right(\mathrm{2}{b}_{\text{s}}+{\mathit{\nu }}_{m,\text{s}}+\mathrm{1})/{\mathit{\mu }}_{m,\text{s}})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,\text{s}}+\mathrm{1})/{\mathit{\mu }}_{m,\text{s}})}}{\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,\text{s}}+\mathrm{1})/{\mathit{\mu }}_{m,\text{s}})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,\text{s}}+\mathrm{2})/{\mathit{\mu }}_{m,\text{s}})}}\right]}^{\mathrm{2}{b}_{\text{s}}}\text{(A40)}& {\displaystyle }{\mathit{\delta }}_{D,ss}^{\mathrm{0}}=\mathrm{2}{\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\right(b_{\text{s}}+{\mathit{\nu }}_{m,\text{s}}+\mathrm{1})/{\mathit{\mu }}_{m,\text{s}})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,\text{s}}+\mathrm{1})/{\mathit{\mu }}_{m,\text{s}})}}\right]}^{\mathrm{2}}{\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,\text{s}}+\mathrm{1})/{\mathit{\mu }}_{m,\text{s}})}{\mathrm{\Gamma }\left(\right({\mathit{\nu }}_{m,\text{s}}+\mathrm{2})/{\mathit{\mu }}_{m,\text{s}})}}\right]}^{\mathrm{2}{b}_{\text{s}}}\end{array}$$

The velocity variance simplifies to

$$\begin{array}{}\text{(A41)}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{\mathrm{0},\mathrm{ss}}=\sqrt{\mathrm{2}}{\mathit{\beta }}_{v,\text{s}}[{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_{v,\text{s}}}{{\mathit{\lambda }}_{\text{eq, s}}}}\right)}^{-{\mathit{\xi }}_{D,s}}-{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,\text{s}}}{{\mathit{\lambda }}_{\text{eq, s}}}}\right)}^{-\mathrm{2}{\mathit{\xi }}_{D,s}}{]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

with

$$\begin{array}{}\text{(A42)}& {\mathit{\xi }}_{D,s}={\mathit{\mu }}_{\text{eq, s}}+\mathrm{3}.\end{array}$$

A4.2 A-kernel

𝒞 of the A-kernel for snow self-collection simplifies to

$$\begin{array}{}\text{(A43)}& \begin{array}{rl}& {\mathcal{C}}_{\mathrm{0},\mathrm{ss}}=N_{\text{s}}^{\mathrm{2}}[\mathrm{2}{\mathit{\delta }}_{A,\text{s}}^{\mathrm{0}}+{\mathit{\delta }}_{\text{A,ss}}^{\mathrm{0}}]{\left({\displaystyle \frac{\mathrm{6}{a}_{m,\text{s}}}{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}}\right)}^{\frac{{\mathit{\sigma }}_{A,\text{s}}}{\mathrm{3}}}\\ & min\left({\mathit{\gamma }}_{A,\text{s}}{\left({\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}{\mathrm{6}{a}_{m,\text{s}}}}\right)}^{\frac{{\mathit{\sigma }}_{A,\text{s}}}{\mathrm{3}}}{\stackrel{\mathrm{‾}}{D}}_{\text{max,s}}^{{\mathit{\sigma }}_{A,\text{s}}^{*}},{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{\stackrel{\mathrm{‾}}{D}}_{\text{max,s}}^{\mathrm{2}}\right),\end{array}\end{array}$$

with the following.

$$\begin{array}{}\text{(A44)}& {\displaystyle }{\mathit{\delta }}_{A,\text{s}}^{\mathrm{0}}={\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq, s}}+{\mathit{\sigma }}_{A,\text{s}}+\mathrm{1})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq, s}}+\mathrm{1})}}{c}_{\mathit{\lambda },\text{s}}^{{\mathit{\sigma }}_{A,\text{s}}}\text{(A45)}& {\displaystyle }{\mathit{\delta }}_{\text{A,ss}}^{\mathrm{0}}=\mathrm{2}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\mathrm{eq},\mathrm{s}}+{\mathit{\sigma }}_{A,\text{s}}/\mathrm{2}+\mathrm{1}{)}^{\mathrm{2}}}{\mathrm{\Gamma }({\mathit{\mu }}_{\mathrm{eq},\mathrm{s}}+\mathrm{1}{)}^{\mathrm{2}}}}{c}_{\mathit{\lambda },\mathrm{s}}^{{\mathit{\sigma }}_{A,\text{s}}}\text{(A46)}& {\displaystyle }{\mathit{\sigma }}_{A,\text{s}}^{*}={\displaystyle \frac{b_{m,\text{s}}{\mathit{\sigma }}_{A,\text{s}}}{\mathrm{3}}}\text{(A47)}& {\displaystyle }{c}_{\mathit{\lambda },\text{s}}={\left[{\displaystyle \frac{\mathrm{6}{a}_{m,\text{s}}}{\mathit{\pi }{\mathit{\rho }}_{\text{w}}}}{\displaystyle \frac{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq, s}}+\mathrm{1})}{\mathrm{\Gamma }({\mathit{\mu }}_{\text{eq, s}}+\mathrm{4})}}\right]}^{\mathrm{1}/\mathrm{3}}\end{array}$$

The area ratios are limited in the same way as for ice self-collection.

The velocity variance simplifies to

$$\begin{array}{}\text{(A48)}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }v}}_{\mathrm{0},\mathrm{s}}=\sqrt{\mathrm{2}}{\mathit{\beta }}_{v,\text{s}}{\left[{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{\mathit{\gamma }}_{v,\text{s}}}{{\mathit{\lambda }}_{\text{eq, s}}}}\right)}^{-{\mathit{\xi }}_{A,\text{s}}}-{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\gamma }}_{v,\text{s}}}{{\mathit{\lambda }}_{\text{eq, s}}}}\right)}^{-\mathrm{2}{\mathit{\xi }}_{A,\text{s}}}\right]}^{\mathrm{0.5}},\end{array}$$

with

$$\begin{array}{}\text{(A49)}& {\mathit{\xi }}_{A,\text{s}}={\mathit{\mu }}_{\text{eq, s}}+{\mathit{\sigma }}_{A,\text{s}}+\mathrm{1}.\end{array}$$

Table A1 Prefactor Φ of the aggregation rates (Eq. A2) for different aggregation processes and the predicted moments of the cloud ice and snow distribution.

Appendix B: Atmospheric setup for 1D simulation and atmospheric fields of the case study predicted by ICON-LEM

Figure B5 shows the atmospheric variables from O20 simulations and the setup for the snowshaft simulations. Figure B6 shows the atmospheric variables of the case study.

Figure B1 Numeric and analytic solution of the bulk collision rates for ice–snow (column and Mix2, respectively) collisions for Atlas-type and power-law velocity size relations. The shape parameter is μ_eq=2 (Eq. 7), which is equal to μ_m=0 (Eq. 5) for cloud ice and snow. Left: number density, right: mass density; top: D-kernel; bottom: A-kernel. (a, c) normalized number collision rate, (b, d) normalized mass collision rate, (a, b) D-kernel, (c, d) A-kernel.

Figure B2 Same as Fig. B1 but with μ_eq=8 (Eq. 7), which is equal to μ_m=2 (Eq. 5) for snow.

Figure B3 Numeric and analytic solution of the bulk collision rates for ice–ice (both column) collisions: (a, c) normalized number collision rate, (b, d) normalized mass collision rate, (a, b) D-kernel, (c, d) A-kernel.

Figure B4 Numeric and analytic solution of the bulk collision rates for snow–snow (both Mix2) collisions: (a) D-kernel; (b) A-kernel.

Figure B5 Setup of atmospheric variables in the 1D simulations (Sect. 3.2) (black line), which was chosen based on the histograms from the ICON-LEM simulation (the histogram is shown in the background, O20). The histogram is filtered to include only profiles for which the rain rate exceeds 1 mm h⁻¹. (a) Temperature, (b) relative humidity with respect to water, and (c) relative humidity with respect to ice. The height of the melting temperature 0^∘C is set to 0 m, and other heights are calculated assuming a temperature gradient of 0.0062 K m⁻¹. Counts in the ICON-LEM simulations from O20 are color-coded.

Figure B6 Temperature (a), vertical velocity (b), and relative humidity with respect to water (c) and ice (d) over Jülich in the SB06 default simulation on 3 January 2016. Temperature isolines are shown in each plot.