Investigating ice formation pathways using a novel two-moment multi-class cloud microphysics scheme

For pure ice clouds in the cold-temperature regime (T<235 K), two major formation pathways are possible. Liquid origin ice clouds stem from droplets that freeze close to water saturation. In-situ-formed ice clouds form directly from the vapor phase below water saturation. For a better investigation of these pathways, we developed a novel microphysics scheme. The new two-moment scheme distinguishes between five ice classes (“ice modes”) each with their own unique formation mechanism: homogeneous freezing of solution droplets, deposition nucleation, homogeneous freezing of cloud droplets and raindrops, immersion freezing, and secondary ice from rime splintering. The ice modes interact with each other, e.g., in competition for growth by deposition of water vapor and aggregation, but also with the other cloud particle classes, i.e., cloud droplets, rain, snow, graupel, and hail.

This scheme was employed to investigate the liquid origin vs. in situ formation in the fully glaciated parts of an idealized convective cloud. The majority of the cloud ice in the deep convection cloud consisted of frozen droplets (liquid origin). This was caused by the high number concentration of cloud droplets available for freezing. In-situ-formed ice was only relevant for the overshoot where ice from both formation pathways mixed.

The new scheme is also useful for investigation of the ice formation in the mixed-phase parts of the convective cloud. We find a vertical layering of ice modes in the cloud. The lowermost layer consists of secondary ice from rime splintering and occurs near the updraft core at temperatures around the Hallett–Mossop zone. At altitudes between 6 and 9 km, ice mostly stems from immersion freezing. We find a correlation between the abundance of ice from immersion freezing and snow. The majority of ice crystals above 9 km stems from homogeneously frozen cloud droplets since ice-nucleating particles (INPs) required for immersion freezing were quickly depleted.

1 Introduction

Clouds cover a large part of the planet and constitute an important component of the Earth–atmosphere system. We can discriminate between different thermodynamic regimes for clouds. For temperatures above melting point (T>273 K), clouds consist of liquid water droplets. In the $\mathrm{235}\mathrm{K}\le T\le \mathrm{273}\mathrm{K}$ temperature regime, supercooled liquid droplets and ice particles can exist; however, clouds containing both kinds of hydrometeors are called mixed-phase clouds. At lower temperatures (T<235 K), only the solid phase exists; i.e., we find pure ice clouds.

Clouds influence the hydrological cycle, e.g., by forming precipitation. Moreover, clouds influence general circulation by diabatic processes induced by phase changes (Reed et al., 1992; Wernli and Gray, 2024). Clouds also affect the system's energy budget by interaction with solar and thermal radiation. Incoming solar radiation is partly scattered and reflected back to space (albedo effect), whereas the Earth’s outgoing infrared radiation is partly absorbed and re-emitted at a different temperature (greenhouse effect). The IPCC report of 2021 states that the overall cloud feedback on climate is positive, hence enhancing global warming. However, there are still many cloud regimes with large uncertainties regarding their impact on the global energy budget (Forster et al., 2021). One of these cloud regimes is high-altitude (i.e., cold-temperature regime) pure ice clouds, referred to as cirrus. Since the albedo and greenhouse effects of cirrus are of the same magnitude but different sign, the combined effect can lead to a net cooling or net warming depending on microphysical properties (see, e.g., Zhang et al., 1999; Fusina et al., 2007; Joos et al., 2014). Overall, our understanding of clouds and the relevant and dominant processes is still quite limited, especially for pure ice clouds (Krämer et al., 2020).

For pure ice clouds at low temperatures, the pathway of formation remains unclear a priori. The different formation mechanisms have recently been summarized and classified in the following ways (see, e.g., Krämer et al., 2016; Wernli et al., 2016):

• Formation pathways related to freezing of pre-existing cloud droplets (i.e., in the mixed-phase temperature regime and close to water saturation) are termed liquid origin formation.

• Formation pathways related to direct formation of ice crystals from vapor at (liquid or solid) aerosols at low temperatures and below water saturation are termed in situ formation.

It is quite difficult to determine the formation pathway of ice crystals from observations. For ice particles measured with imaging techniques, there are some indications that complex shapes and large particle sizes are most probable for liquid origin ice crystals, whereas in-situ-formed ice crystals remain small, and their shapes remain simple, e.g., as quasi-spherical shapes (Wolf et al., 2018). This can be explained by the amount of available water vapor, which can be used for diffusional growth. In the mixed-phase regime, more water vapor is available, leading to larger and more complex shapes, whereas at low temperatures the amount of water vapor is very limited, leading to small-sized and simply shaped ice crystals. It also seems to be probable that liquid origin ice clouds consist of much more ice particles than in-situ-formed ice clouds (Krämer et al., 2016, 2020). This can be explained by the formation process. Liquid origin ice particles stem from pre-existing cloud droplets; the number concentration of water droplets in liquid clouds is usually of the order of some hundreds of particles per cubic centimeter and thus much larger than the amount of available ice nuclei at low temperatures.

There are some indications that we might be able to use the particles' properties themselves for a classification of the formation pathways (see, e.g., Luebke et al., 2016; Wolf et al., 2018). However, such evaluations are just based on single cases of measurements together with trajectory calculations and have an inherent uncertainty, which cannot be quantified. Another approach to identify the ice formation pathway is to employ retrievals from remote sensing to calculate cloud properties and then to search for prescribed characteristics (see Urbanek et al., 2017; Huo et al., 2020).

The formation pathway can also be estimated by using model simulations. Gasparini and Lohmann (2016) utilized a temperature-based criterion to classify cirrus cloud types. Wernli et al. (2016) utilized Lagrangian trajectories to identify the origin of the cirrus cloud. If cirrus originated from the liquid phase (via a mixed-phase cloud) or from clear sky, the resulting ice cloud was classified as liquid origin or in situ, respectively. However, both approaches do not assess the formation pathway directly and do not consider transport relative to the air parcel, i.e., neither sedimentation nor turbulent mixing. In addition, the evaluations were carried out on coarse grid resolutions (climate model or reanalysis data).

A complementary viewpoint from the (even small-scale) modeling perspective is necessary in order to give more rigorous insights into the different formation mechanisms in cloud systems, especially for vertically extending systems such as convective clouds or warm conveyor belts that show both the mixed-phase and the pure ice thermodynamic regimes.

A particle‐by‐particle direct numerical simulation of each cloud particle is not feasible for realistic clouds (Morrison et al., 2020). Instead, bulk schemes are a commonly used model approach to parameterize cloud physics in large-scale atmospheric models with spatial resolutions of several kilometers. They have also been utilized in high-resolution large-eddy simulations and, on the other end of the spatial and temporal scale, in climate models. Instead of describing the evolution of individual cloud particles, bulk schemes describe the evolution of mean quantities of an ensemble of particles. The cloud particles are sorted into preset classes and assigned a (semi-empirical) type of size distribution. Integrating the evolution equation over the size distribution, we obtain the temporal evolution of general moments (Hulburt and Katz, 1964). Instead of size distributions, the particle mass is often used as a variable of distributions, thus leading to general moments of mass; in this work, we use this approach, which is more physical since mass is the conserved quantity (see, e.g., Seifert and Beheng, 2006). Bulk schemes were first introduced by Kessler (1969) but only considered cloud droplets and raindrops. Since the work of Lin et al. (1983), ice particles have also been commonly included in bulk schemes. For practical reasons, we use a finite number of moments in order to describe the system. The order of the bulk scheme is set by the number of predicted moments. One-moment bulk schemes typically predict the mass content (first moment) of a cloud particle class. Two-moment schemes predict the number concentration (zeroth moment) in addition to the first moment. For a meaningful closure of the system, we usually fix the type of the underlying mass distribution.

Some common choices for bulk schemes in numerical weather prediction and atmospheric research models are Seifert and Beheng (2006) for the Icosahedral Nonhydrostatic (ICON) model (Zängl et al., 2015; Heinze et al., 2017), Thompson et al. (2008) and Morrison et al. (2009) for the WRF model (Skamarock et al., 2019), Field et al. (2023) for the Unified Model (Walters et al., 2019), and the bulk scheme employed in the IFS model (ECMWF, 2023). The predicted particle properties (P3) scheme (Morrison and Milbrandt, 2015) is also commonly employed in cloud research. Instead of multiple distinct ice particle classes (cloud ice, snow, graupel, hail), they employ free ice-phase classes that can evolve into any ice particle type (or mixture thereof).

In this paper we introduce a novel bulk scheme, the ice mode scheme. This new scheme is based on the standard two-moment scheme in ICON (Seifert and Beheng, 2006) (SB hereafter). However, instead of a single ice crystal (or “cloud ice”) class, it features several classes, each equipped with their own unique formation process. Apart from the unique formation processes, these ice modes use the same parameterizations for shape, terminal velocity, and all other microphysical processes. They are subject to advection, turbulent mixing, and sedimentation. Ice modes (like all microphysics classes) can also coexist; i.e., like in the original SB scheme, multiple ice formation processes can still be active at the same time if the relevant thermodynamic requirements are met.

The purpose of the new ice mode scheme is not to perform better than the base scheme of SB concerning quantifiable metrics like precipitation rates and patterns or even scores of operational numerical weather prediction models. In that regard the aim is that the ice mode scheme performs similarly to and ideally the same as the original SB scheme. The benefit of employing multiple cloud ice classes with unique formation processes is that the ice mode scheme retains the information of how cloud ice originally formed. This allows for a novel way of addressing research questions for glaciated clouds. Liquid origin vs. in situ formation can be quantified by comparing the amount of cloud ice stemming from frozen droplets and directly from the vapor-phase-formed ice. The relevance of secondary ice processes can be assessed as the mass content and number concentrations of primary and secondary ice are separate tracers. Also note that the separation into these ice modes itself is not artificial. Conceptually, we could attribute a distinct formation process to each ice crystal measured in real clouds or laboratory experiments as it is a unique property assigned to the particle at its formation.

Hence, the new ice mode scheme is a useful tool for the investigation of cloud evolution. Additionally, it can evaluate the sensitivity of different ice formation pathways to environmental conditions, e.g., thermodynamic variables and domain properties like spatial and temporal resolution. It is also useful for testing new ice nucleation parameterizations as the individual impact of each routine on the model cloud can be directly assessed.

There are two other approaches to investigating ice formation pathways in models. The first is to track the (usually temporally or spatially accumulated) process rates of ice nucleation as additional model output variables. This approach has been used in studies that evaluate the impact of individual cloud processes including nucleation (see, e.g., Barrett and Hoose, 2023; Han et al., 2023; Oertel et al., 2023; Han et al., 2024; Schwenk and Miltenberger, 2024). However, this approach only provides the amount of ice that formed from each individual formation process. It does not track what happens to this ice after formation as the signal gets blurred by advection; turbulent mixing; sedimentation; and removal by collisions, melting, and evaporation. Hence, this approach provides information on how much ice was formed by each process, and the ice mode scheme provides information on how much ice from each formation process is left. These are not derivative but complementary viewpoints.

One other approach that can (in theory) track all microphysical processes, including formation, and retain that information is (Lagrangian) particle models. These models simulate the properties of “super particles” that represent an ensemble of usually 10¹ to 10⁴ members. These particles move freely through the model domain in a Lagrangian sense and in general do not belong to a certain preset class but instead are able to evolve dynamically depending on environmental conditions and their interaction with other particles. There are only a few particle models that include the ice phase: Sölch and Kärcher (2010), Unterstrasser and Sölch (2010), Brdar and Seifert (2018), and Shima et al. (2020). However, there is no particle model yet that (directly) models all ice formation processes relevant to both the liquid origin and the in situ formation pathways. Another issue is that particle models are computationally expensive (Morrison et al., 2020). Hence, they are only employed in parcel models or high-resolution (cloud-resolving) models to investigate atmospheric systems on the micro- to mesoscale. The ice mode bulk scheme introduced in this work can be used to investigate synoptic-scale systems like warm conveyor belts and can even be used in global simulations.

In a first application of the new model, we investigate the formation of ice clouds for an idealized test case of convection and determine the formation pathways of (anvil) ice clouds. Although it has always been assumed that the anvil of a convective thunderstorm should mainly be formed by frozen cloud droplets transported upwards, it has not been shown before (Gasparini et al., 2018). Using our newly developed ice mode scheme, we can investigate this statement quantitatively.

In Sect. 2 we describe the model equations and relevant microphysical parameterizations affecting ice in the model and then present an idealized simulation of a convective cloud in Sect. 4. We compare our results with a reference simulation that uses the standard SB scheme and also demonstrate the capabilities of the new scheme to characterize the origin of cloud ice and the impact of different nucleation parameterizations.

2 Model description

The basis for our model development is the two-moment scheme of SB. The original SB scheme distinguished between six classes of cloud particles: cloud droplets, rain, (cloud) ice, snow, graupel, and hail. While the formulation for the liquid droplet's treatment in the model remains unchanged, we reformulate the ice microphysics and also the interaction of ice and water particles in terms of collision processes. In the following the new scheme is described.

2.1 General settings

The model formulation relies on the usual approach of a bulk model with two (general) moments. Instead of computing the temporal and spatial evolution of a (maybe multivariate) mass distribution with a Boltzmann-type evolution equation, we use the integrated version, i.e., the evolution equations for moments of the underlying mass distribution of the respective class of hydrometeors. The moments are defined in the usual way; i.e.,

$$\begin{array}{}\text{(1)}& M_i^k=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{x}^k{f}_i\left(x\right)\mathrm{d}x\end{array}$$

with the mass of particles x in kilograms and the mass distribution f_i in inverse kilograms for the respective class i of hydrometeors. The mass distribution describes the number and is normalized by the total number concentration of particles n_i. For a meaningful closure of the systems of equations, we have to choose a suitable type of mass distribution. We generally assume that the particle mass distributions can be represented by generalized gamma distributions (see also Seifert and Beheng, 2006) of the form

$$\begin{array}{}\text{(2)}& f_i\left(x\right)=A_i{x}^{{\mathit{\nu }}_i}\exp \left(-{\mathit{\lambda }}_i{x}^{{\mathit{\mu }}_i}\right),\end{array}$$

where the shape parameters ν_i and μ_i are prescribed, and λ_i and A_i are linked to the zeroth and first distribution moments, the number concentration n_i, and the mass content q_i:

$$\begin{array}{}\text{(3)}& {\displaystyle }{n}_i=M_i^{\mathrm{0}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}f\left(x\right)\mathrm{d}x={\displaystyle \frac{A}{\mathit{\mu }{\mathit{\lambda }}^{\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}}}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}}\right),\text{(4)}& {\displaystyle }{q}_i=M_i^{\mathrm{1}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}f\left(x\right)x\mathrm{d}x={\displaystyle \frac{A}{\mathit{\mu }{\mathit{\lambda }}^{\frac{\mathit{\nu }+\mathrm{2}}{\mathit{\mu }}}}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathit{\nu }+\mathrm{2}}{\mathit{\mu }}}\right).\end{array}$$

See Appendix B for details on the properties of the gamma distribution and the analytical solution of the integrals.

For the formulation of single particle processes (e.g., growth or sedimentation), we have to introduce size–mass relations and velocity–mass relations of the form

$$\begin{array}{}\text{(5)}& D_i\left(x\right)=a_i{x}^{b_i},v_i\left(x\right)={\mathit{\alpha }}_i{x}^{{\mathit{\beta }}_i}\end{array}$$

with the (generalized) size D_i and the terminal velocity v_i for a particle of mass x within a class i. These relations represent the different shapes of particles from different classes. The remaining classes of hydrometeors are labeled by indices, i.e., c for cloud droplets, r for raindrops, s for snow, g for graupel, and h for hail. The old class cloud ice is now split into five classes. For details of the formulation and the determination of the shapes for the hydrometeor classes used, we refer to Seifert and Beheng (2006).

For our newly introduced ice modes we use the same parameters of the distribution (Eq. 2) and the relations (Eq. 5) for all classes. While the formation pathway of an ice crystal might have an impact on its shape, the morphology of the particle is mainly determined by the environmental conditions encountered during its growth, e.g., temperature and humidity (see, e.g., Magono and Lee, 1966; Kobayashi, 1967; Libbrecht, 2005; Pruppacher et al., 2010), which are not accounted for in the SB scheme. Using the same coefficients and parameterizations for each ice mode, apart from its source, also has the advantage that it allows for a more concise interpretation of the impact of the individual ice formation pathways on the cloud evolution.

The time evolution of the kth-moment $M_i^k$ of an ice mode i (or another class of hydrometeors) is governed by

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\partial M_i^k}{\partial t}}& +\mathrm{\nabla }\cdot \left[v{M}_i^k\right]+{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{v}}_{i,k}{M}_i^k}{\partial z}}-\mathrm{\nabla }\cdot \left[K_{\mathrm{h}}\mathrm{\nabla }{M}_i^k\right]\\ & =\text{source/sink}.\end{array}\end{array}$$

The terms on the left-hand side describe the effects of advection with the mean wind velocity v, sedimentation with the weighted mean fall speed ${\stackrel{\mathrm{‾}}{v}}_{i,k}$, and turbulent mixing with the mean turbulent diffusivity of heat K_h. On the right-hand side are the source and sink terms for the particle formation, which is unique for each ice mode; deposition of water vapor and evaporation; and a number of collision processes, most importantly aggregation and riming. These equations must be coupled with a model for atmospheric flows, such as, e.g., suitable approximations of the Navier–Stokes equations within the ICON model.

In the following we focus on the description of ice-related physics and refer to Seifert and Beheng (2006) for an in-depth description of warm- and mixed-phase microphysics.

2.2 Ice formation pathways

For the treatment of ice particles, we introduce new classes of ice particles, discriminated by their formation mechanisms. The ice mode scheme features five independent ice classes instead of a single one (former class cloud ice), called ice modes in this work, each with their own unique formation pathway. For each new ice mode, we introduce number concentrations and mass content. The other classes (graupel, hail, and snow) remain unchanged; however, we have to reformulate the processes of interactions between the former class cloud ice and other variables (e.g., see below for collision processes). The new classes are as follows:

• homogeneous freezing of cloud droplets n_frz, q_frz

• immersion freezing of cloud and rain droplets n_imm, q_imm

• freezing of solution droplets (homogeneous nucleation) n_hom, q_hom

• deposition nucleation n_dep, q_dep

• secondary ice from rime splintering n_sec, q_sec.

The sum of all ice modes represents all ice crystals present (n_tot, q_tot).

2.2.1 Homogeneous freezing of cloud droplets

Cloud droplets are considered to be pure water droplets, thus freezing homogeneously at temperatures below the triple point T_m but at water saturation. Homogeneous freezing of cloud droplets is the source of the FRZ ice mode with the bulk quantities q_frz and n_frz. The stochastic process of homogeneous freezing is described by a nucleation rate J_hom in ${\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$, which depends on temperature only (since existence of cloud droplets requires water saturation); thus the change in the number and mass concentrations of the ice mode can be described by

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial n_{\mathrm{frz}}}{\partial t}}=J_{\mathrm{hom}}{q}_{\mathrm{c}},{\displaystyle \frac{\partial q_{\mathrm{frz}}}{\partial t}}=J_{\mathrm{hom}}{q}_{\mathrm{c}}{x}_{\mathrm{c}}{\displaystyle \frac{{\mathit{\nu }}_{\mathrm{c}}+\mathrm{2}}{{\mathit{\nu }}_{\mathrm{c}}+\mathrm{1}}}\end{array}$$

using the mass content of cloud droplets q_c, the mean mass of cloud droplets x_c, and the distribution shape parameter ν_c of the cloud mass distribution (generalized gamma distribution). For the formulation of the homogeneous freezing coefficient of cloud droplets in $\mathrm{\mu }{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$, we use the fit of Cotton and Field (2002) in the formulation of Jeffery and Austin (1997):

$$\begin{array}{}\text{(8)}& \log \left({\displaystyle \frac{J_{\mathrm{hom}}}{{\mathit{\rho }}_{\mathrm{w}}}}\right)=\left\{\begin{array}{l}-\mathrm{243.4}-\mathrm{14.75}T-\mathrm{0.307}{T}^{\mathrm{2}}\\ -\mathrm{0.00287}{T}^{\mathrm{3}}-\mathrm{1.02}\times {\mathrm{10}}^{-\mathrm{5}}{T}^{\mathrm{4}}& \mathrm{for}T\le \mathrm{243}\mathrm{K}\\ -\mathrm{7.63}-\mathrm{2.996}(T+\mathrm{30})& \mathrm{for}T>\mathrm{243}\mathrm{K},\end{array}\right.\end{array}$$

where ρ_w denotes the liquid water mass density. The parameterization for the nucleation rate is only valid at water saturation. However, due to saturation adjustment (see SB), the system stays at water saturation permanently as long as cloud droplets are present.

Homogeneous freezing of raindrops in the SB scheme is omitted, since raindrops freeze rapidly by heterogeneous freezing before reaching temperature levels close to the homogeneous freezing temperature.

Figure 1 Homogeneous nucleation following Kärcher et al. (2006) for no pre-existing ice. (a) Critical supersaturation S_cr serves as a nucleation threshold. (b) The number of ice crystals after a homogeneous nucleation event at S_cr.

2.2.2 Homogeneous freezing of solution droplets

Aqueous solution droplets, i.e., liquid aerosol particles, can be supercooled to lower temperatures than pure liquid droplets. The solute obviously impedes the establishment of a critical cluster. Koop et al. (2000) showed that the effect of solutes on the freezing temperature is driven by their thermodynamic (equilibrium) quantities, which can be expressed in terms of water activity, defined as the water saturation vapor pressure ratio between the solution and pure water $a_{\mathrm{w}}=e_{\mathrm{sol}}/e_{\mathrm{liq}}$. When we assume that the solution droplet is in a thermodynamic equilibrium with its environment, then the freezing temperature is independent of the choice of chemical droplet composition, at least for inorganic compounds. This translates into a nucleation rate for solution droplets, which solely depends on the environmental conditions, i.e., temperature and supersaturation with respect to ice, but not on the substance in the water drops. The nucleation rate can be formulated using a threshold of critical supersaturation, as could be shown by Spichtinger et al. (2023), e.g., using the formulation by Ren and MacKenzie (2005):

$$\begin{array}{}\text{(9)}& S_{\text{cr}}=\mathrm{2.349}-T/\mathrm{259}.\end{array}$$

Panel (a) in Fig. 1 shows that quite high supersaturations (S_i>1.4) are necessary for homogeneous nucleation events to occur, especially in low-temperature regimes (T<230 K) where freezing of pure water droplets is no longer possible and homogeneous nucleation is commonly observed. This ice formation process is called homogeneous nucleation and is the source term for the HOM ice mode with the bulk properties n_hom and q_hom. Generally, we would obtain similar equations for the change in n_hom,q_hom to the freezing of cloud droplets (see Sect. 2.2.1). However, since the system is in a non-equilibrium state and there is no boundary condition (like the assumption of water saturation as above), we would have to represent the evolution of the saturation ratio, also changed by diffusional growth; this would require a small time step for the numerical scheme and is thus not feasible.

Therefore, we use the parameterization by Kärcher and Lohmann (2002), which describes the homogeneous nucleation event in an ascending air parcel. When an air parcel ascends adiabatically, supersaturation is generated by adiabatic cooling. When the critical supersaturation S_cr is reached, homogeneous nucleation is triggered. The newly nucleated ice crystals deplete supersaturation by depositional growth. The competition between generating supersaturation by adiabatic cooling and depleting supersaturation by depositional growth in an adiabatically ascending air parcel, driven by a constant wind velocity w, can be described as

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=a_{\mathrm{1}}Sw-\left(a_{\mathrm{2}}+a_{\mathrm{3}}S\right)R_i\left(t\right)\end{array}$$

with the parameters

$$\begin{array}{}\text{(11)}& a_{\mathrm{1}}={\displaystyle \frac{L_{\mathrm{s}}{M}_{\mathrm{w}}g}{c_{\mathrm{p}}R{T}^{\mathrm{2}}}}-{\displaystyle \frac{M_{a}g}{RT}},a_{\mathrm{2}}={\displaystyle \frac{T{k}_{\mathrm{b}}}{e_{\mathrm{sat},i}}},a_{\mathrm{3}}={\displaystyle \frac{L_{\mathrm{s}}^{\mathrm{2}}{M}_{\mathrm{w}}{m}_{\mathrm{w}}}{c_{\mathrm{p}}pT{M}_{\mathrm{a}}}}.\end{array}$$

Here, a₁ describes the effect of adiabatic cooling, where a₂, a₃, and R_i describe the depletion of supersaturation by depositional growth as a function of time. The parameters and process of depositional growth are further explained in Sect. 2.3.

If the air parcels continue to rise, aqueous solution droplets continue to freeze. When the newly nucleated ice crystals deplete more supersaturation than is being generated by adiabatic cooling, the maximum value of the supersaturation S_* is reached. The supersaturation is depleted, and once it reaches S_cr, the freezing of solution droplets stops. The integrated number of solution droplets frozen within the freezing time interval t_* describes the number of ice crystals from the homogeneous nucleation event. Numerical simulations suggest that the freezing time interval t_* is short, and the maximum value of supersaturation S_* is approximately equal to S_cr (Kärcher et al., 2006). Thus the number concentration of ice particles in a homogeneous nucleation event can then be estimated as

$$\begin{array}{}\text{(12)}& n_{\mathrm{hom}}=w{\displaystyle \frac{a_{\mathrm{1}}{S}_{\text{cr}}}{a_{\mathrm{2}}+a_{\mathrm{3}}{S}_{\text{cr}}}}{\displaystyle \frac{\mathrm{1}}{R_{im}\left(r_{\mathrm{0}}\right)}},\end{array}$$

where R_im is the analytical approximation of the integral that describes the growth of the solution droplet within the freezing time interval where we refer to Kärcher et al. (2006) for a detailed derivation. Here, r₀=27.2 nm is the radius of the monodisperse aqueous solution droplets. We obtain the bulk mass of the newly nucleated homogeneous ice mode as

$$\begin{array}{}\text{(13)}& q_{\mathrm{hom}}={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\pi }{\mathit{\rho }}_i{\left({\displaystyle \frac{r_{\mathrm{0}}\left(\mathrm{1}+b_{\mathrm{KH}}{r}_{\mathrm{0}}\right)-\mathrm{1}}{b_{\mathrm{KH}}}}\right)}^{\mathrm{3}}{n}_{\mathrm{hom}}\end{array}$$

with the parameter $b_{\mathrm{KH}}=\frac{{\mathit{\alpha }}_{\mathrm{d}}{v}_{\mathrm{th}}}{\mathrm{4}{D}_{\mathrm{v}}}$. Here, α_d is the deposition coefficient, v_th is the mean molecular velocity of water vapor, and D_v denotes the diffusivity.

Panel (b) in Fig. 1 shows the number concentration of ice n_i after a homogeneous nucleation event at critical supersaturation S_cr for vertical velocities ranging from synoptic velocities (up to 10 cm s⁻¹) to gravity waves (up to 100 cm s⁻¹) and convection (up to 1000 cm s⁻¹). Homogeneous nucleation is strongly sensitive to the model-level vertical velocity, and large nucleation events can only be represented if vertical velocity is resolved and no other nucleation mechanisms disturb the effect of homogeneous nucleation. Nucleation events at synoptic vertical velocity will be more common in the model and nucleate ice crystals of up to 10 cm⁻³. Note that the parameterization of Kärcher et al. (2006) does not explicitly scale with the microphysics time step as it describes an entire nucleation event and not a nucleation rate. Shorter time steps might lead to higher number concentrations because newly nucleated ice crystals have less time to deplete the supersaturation until a new event is allowed to be triggered. However, the homogeneous ice mass content is still constrained by the availability of supersaturation, i.e., by thermodynamics.

The ice mode scheme uses the extension of Kärcher et al. (2006), where the effect of pre-existing ice depleting supersaturation is taken into account as equivalent to a fictitious downdraft velocity w_pre:

$$\begin{array}{}\text{(14)}& w^{\prime }=w-w_{\mathrm{pre}}=w-{\displaystyle \frac{a_{\mathrm{2}}+a_{\mathrm{3}}{S}_{\mathrm{i}}}{a_{\mathrm{1}}{S}_{\mathrm{i}}}}{R}_{i,\mathrm{pre}},\end{array}$$

where w is the model-level vertical velocity, and R_i,pre is the mean radius of the pre-existing ice. Pre-existing ice usually originates from a heterogeneous ice mode or a prior homogeneous nucleation event, as shown in the simulation results.

2.2.3 Heterogeneous nucleation

In the new scheme we use the two nucleation pathways, i.e., immersion freezing and deposition nucleation. Both ice formation pathways depend on temperature, supersaturation, and properties of the INPs, like effectiveness and the number of sites causing nucleation per unit surface area (Vali et al., 2015). Immersion freezing of both raindrops and cloud droplets is considered in this model.

We obtain the number concentration for the immersion freezing ice mode as the sum of freezing cloud and rain droplets as

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial n_{\mathrm{imm}}}{\partial t}}=C_{\mathrm{IMM}}{|}_{\mathrm{c}}-n_{\mathrm{inact}}+{\displaystyle \frac{\partial n_{\mathrm{imm}}}{\partial t}}{|}_{\mathrm{r}}\end{array}$$

and the mass content as

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial q_{\mathrm{imm}}}{\partial t}}={\displaystyle \frac{\partial n_{\mathrm{imm}}}{\partial t}}{|}_{\mathrm{c}}{x}_{\mathrm{c}}+{\displaystyle \frac{\partial q_{\mathrm{imm}}}{\partial t}}{|}_{\mathrm{r}},\end{array}$$

with the mean cloud droplet mass x_c. Similarly, we obtain the number concentration of the deposition nucleation ice mode as

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\partial n_{\mathrm{dep}}}{\partial t}}=C_{\mathrm{DEP}}-n_{\mathrm{inact}}\end{array}$$

and the mass content

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\partial q_{\mathrm{dep}}}{\partial t}}={\displaystyle \frac{\partial n_{\mathrm{dep}}}{\partial t}}{x}_{i,\min }\end{array}$$

with the minimum mass of ice crystals in the model x_i,min. C_IMM and C_DEP describe the number concentration in m⁻³ of activated INPs for immersion freezing (of cloud droplets) and deposition nucleation, respectively. The model does not feature an explicit model for INPs or aerosols in general. Instead, the number of activated INPs n_inact in m⁻³ is tracked as an additional tracer and advected along with the other ice bulk properties. In further nucleation events, n_inact is subtracted from the INP number. n_inact relaxes back to zero in an ice-free environment:

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial n_{\mathrm{inact}}}{\partial t}}=\left\{\begin{array}{ll}-\frac{n_{\mathrm{inact}}}{{\mathit{\tau }}_{\mathrm{inact}}}& \text{if}{q}_{\mathrm{tot}}=\frac{\partial n_{\mathrm{imm}}}{\partial t}=\frac{\partial n_{\mathrm{dep}}}{\partial t}=\mathrm{0}\\ \frac{\partial n_{\mathrm{imm}}}{\partial t}+\frac{\partial n_{\mathrm{dep}}}{\partial t}& \text{else},\end{array}\right.\end{array}$$

with the relaxation timescale τ_inact=600 s and the total ice content (as the sum of all ice modes) q_tot.

First, we describe immersion freezing for rain droplets, i.e., for large water droplets. Bigg (1953) performed laboratory studies investigating the freezing of supercooled purified water in the presence of INPs and estimated the probability of a raindrop freezing depending on temperature and droplet volume. Using these results, the freezing rate of rain can be expressed as

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\displaystyle \frac{\partial n_{\mathrm{imm}}}{\partial t}}{|}_{\mathrm{r}}& =-J_{\mathrm{bigg}}{q}_r\\ & =-A_{\mathrm{imm}}{\mathit{\rho }}_w^{-\mathrm{1}}\exp \left(B_{\mathrm{imm}}\left(T_{\mathrm{m}}-T\right)-\mathrm{1}\right)q_{\mathrm{r}},\end{array}\end{array}$$

where we use the coefficients $A_{\mathrm{imm}}=\mathrm{200}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ and $B_{\mathrm{imm}}=\mathrm{0.65}{\mathrm{K}}^{-\mathrm{1}}$ for rainwater following Barklie and Gokhale (1959). We obtain the mass of frozen raindrops as

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial q_{\mathrm{imm}}}{\partial t}}{|}_r=-A_{\mathrm{imm}}{\mathit{\rho }}_{\mathrm{w}}^{-\mathrm{1}}\exp \left(B_{\mathrm{imm}}\left(T_{\mathrm{m}}-T\right)-\mathrm{1}\right)q_{\mathrm{r}}{\stackrel{\mathrm{‾}}{x}}_r{Z}_{\mathrm{r}},\end{array}$$

with the second moment of the rain mass distribution $Z_{\mathrm{r}}=M_{\mathrm{r}}^{\mathrm{2}}$ and the melting temperature T_m. The mass and number concentration of frozen raindrops is partitioned into the immersion freezing ice mode (IMM), graupel, and hail depending on its diameter. Raindrops smaller than 0.5 mm freeze into ice, raindrops with sizes between 0.5 and 1.25 mm are shifted into the graupel class, and raindrops larger than 1.25 mm are considered hail.

Second, we present the nucleation schemes for immersion freezing of cloud droplets and deposition nucleation, resulting in changes in the ice modes n_imm,q_imm (immersion freezing) and n_dep,q_dep (deposition nucleation). The ice mode scheme offers a choice between the three heterogeneous nucleation schemes, i.e., Hande et al. (2015), Ullrich et al. (2017), and Phillips et al. (2008). In the following, we introduce all three options. The choice of a heterogeneous nucleation scheme has a large impact on ice formation. We compare the impact of the heterogeneous nucleation choice on the ice modes for an idealized convective case in Sect. 4.3.

Hande et al. (2015) (HA15) regarded dust as the main source of INPs over Europe. They used the COnsortium for Small-scale MOdelling (COSMO) meteorological model coupled to the MUlti-Scale Chemistry Aerosol Transport (MUSCAT) model to simulate Sahara dust outbreaks for the year 2008. From the statistics of simulated dust concentrations, HA15 calculated atmospheric profiles of potential INPs. Since Sahara dust outbreaks show a strong seasonal variability, HA15 provided mean profiles for each season. The number of active INPs for immersion freezing was then parameterized using the laboratory results of Niemand et al. (2012) for dust particles,

$$\begin{array}{}\text{(22)}& C_{\mathrm{IMM}}\left(T\right)=A_{\mathrm{H}}\cdot \exp \left[-B_{\mathrm{H}}{\left(T-T_{\mathrm{H},\min }\right)}^{C_{\mathrm{H}}}\right],\end{array}$$

valid for the temperature range from $T_{\mathrm{H},\min }=\mathrm{237.15}\mathrm{K}$ to 261.15 K with the set of coefficients $A_{\mathrm{H}},B_{\mathrm{H}},C_{\mathrm{H}}$, and T_H,min being chosen depending on the season (see Table 1 in Hande et al., 2015). The number concentration of active INPs for deposition nucleation was estimated as follows using the parameterization of Steinke et al. (2015):

$$\begin{array}{}\text{(23)}& C_{\mathrm{DEP}}\left(T,S_{\mathrm{i}}\right)=C_{\mathrm{IMM}}\left(T\right)\cdot \left(a_{\mathrm{H}}\arctan \left(b_{\mathrm{H}}\left(S_{\mathrm{i}}-\mathrm{1}\right)+c_{\mathrm{H}}\right)+d_{\mathrm{H}}\right)\end{array}$$

for temperatures between $T_{\mathrm{H},\min }=\mathrm{220}\mathrm{K}$ and 253 K. Again, the coefficients depend on the season (see Table 1 in Hande et al., 2015).

Ullrich et al. (2017) (UL17) used 11 years of data from ice nucleation experiments in the Aerosol Interaction and Dynamics in the Atmosphere (AIDA) cloud chamber to develop an empirical parameterization for both immersion freezing and deposition nucleation. Dust and soot samples of different types, collected in various locations of the world, were analyzed in AIDA. The dataset includes the results published in Niemand et al. (2012), which were also used in the HA15 parameterization. However, there are some differences in their approach to calculating the density of active surface sites (see UL17), resulting in a slight shift to higher densities. In this study we only consider the dust mode. We obtain the number of active INPs for immersion freezing and deposition nucleation for dust by

$$\begin{array}{}\text{(24)}& C_{\mathrm{IMM}/\mathrm{DEP}}\left(T,S_{\mathrm{i}}\right)=n_{\mathrm{a}}\left(\mathrm{1}-\exp \left[-n_{\mathrm{s},\mathrm{IMM}/\mathrm{DEP}}\left(T,S_{\mathrm{i}}\right){\mathrm{SA}}_{\mathrm{a}}\right]\right),\end{array}$$

with the aerosol number concentration n_a in m⁻³ and the aerosol surface area concentration SA_a and the ice-nucleating active surface site density $n_{\mathrm{s},\mathrm{IMM}/\mathrm{DEP}}$ for immersion freezing and deposition nucleation, respectively, in m⁻². The ice-nucleating active surface site density for immersion freezing follows a simple exponential temperature profile,

$$\begin{array}{}\text{(25)}& n_{\mathrm{s},\mathrm{IMM}}=\exp [\mathrm{150.577}-\mathrm{0.517}T],\end{array}$$

and for deposition nucleation,

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{n}_{\mathrm{s},\mathrm{DEP}}=& \exp [\mathit{\alpha }{\left(S_{\mathrm{i}}-\mathrm{1}\right)}^{\mathrm{1}/\mathrm{4}}\cos {\left(\mathit{\beta }\left(T-\mathit{\gamma }\right)\right)}^{\mathrm{2}}\\ & \mathrm{arccot}\left(\mathit{\kappa }(T-\mathit{\lambda })\right)],\end{array}\end{array}$$

where the fixed coefficients α, β, γ, κ, and λ can be found in UL17. For dust size distribution we use a sum of three lognormal distribution modes. The distribution parameters were chosen such that the immersion freezing temperature profile for UL17 is similar to that of HA15.

Phillips et al. (2008) (PH08) developed an empirically derived framework for heterogeneous nucleation of multiple aerosol species.

INPs are grouped into three basic types: dust and/or metallic aerosols, inorganic black carbon, and insoluble biological aerosols like bacteria and pollen. The basic assumption is that the ice-nucleating activity of insoluble aerosol depends on its number of sites causing nucleation (active sites) and is hence proportional to the total aerosol surface area. The number concentration of active INPs for an aerosol group X was parameterized as

$$\begin{array}{}\text{(27)}& C_{\mathrm{INP},X}=\underset{D_{X,\min }}{\overset{\mathrm{\infty }}{\int }}\left(\mathrm{1}-\exp \left[-{\mathit{\mu }}_X\left(D_X,T,S_{\mathrm{i}}\right)\right]\right)f_X\left(D_X\right)\mathrm{d}{D}_X,\end{array}$$

where μ_X is the Poisson-distributed average of the number of activated INPs for each aerosol particle of size D_X. Here, μ_X was empirically determined by using observational and laboratory data (see PH08). The ice mode scheme uses a simple lookup table and 2D interpolation to determine the fraction of activated INPs as a function of temperature and supersaturation for each aerosol type. Lognormal size distributions are assumed for all aerosol types. The number of INPs is the sum over all three aerosol types with the associated number concentration n_X and active fraction C_INP,X:

$$\begin{array}{}\text{(28)}& C_{\mathrm{INP}}=\sum _X{C}_{\mathrm{INP},X}(T,S_{\mathrm{i}})n_X,\end{array}$$

where we obtain the number of INPs for immersion freezing at and for deposition nucleation below saturation with respect to water. The initial number concentration of soot particles and biological aerosol is fixed, and the number concentration of dust varies with altitude z following a prescribed background profile:

$$\begin{array}{}\text{(29)}& n_{\text{dust}}\left(z\right)=n_{\text{dust},\mathrm{0}}\exp \left[{\displaystyle \frac{z_{\text{dust},\mathrm{0}}-z}{z_{\text{dust},e}}}\right].\end{array}$$

Unless otherwise noted, we only consider the dust mode for PH08 in this work. We chose z_dust,0 such that the maximum number of INPs activated for immersion freezing is the same as for HA15.

Panel (a) in Fig. 2 shows the temperature profiles for immersion freezing parameterizations (IMM) and contour plots of deposition nucleation (c–d). PH08 uses geometric altitude as a coordinate for dust concentration. We transform the altitude into a temperature coordinate using the ICAO Standard Atmosphere. The profile labeled PH08+ also considers soot and biological aerosols in addition to dust.

Figure 2 Number concentration of activated INPs for (a) immersion freezing with profiles for the HA15, UL17, PH08, and PH08+ schemes and (b–d) deposition nucleation as a function of T and S_i for (b) HA15, (c) UL17, and (d) PH08.

The exponential profiles of immersion freezing for HA15 and UL17 are very similar, both starting close to 260 K and capped at the homogeneous freezing temperature threshold. It is not surprising that we only observe minor differences since active INPs for immersion freezing in both parameterizations are based on mostly the same laboratory datasets. On the other hand, immersion freezing for PH08 starts at higher temperatures but increases less steeply past 260 K. For convectively driven clouds in particular, the earlier onset of immersion freezing potentially has a strong impact on the evolution of the ice phase. PH08+ shows activation of INPs at even higher temperatures caused by soot and biological aerosol. However, since we did not tune this version of the PH08 scheme to HA15, the maximum number of active INPs is lower.

2.2.4 Secondary ice production

The only secondary ice process considered in this study is rime splintering (RS). Ice crystals stemming from secondary ice processes are referred to as secondary ice particles (SIPs). Supercooled droplets colliding with ice particles (riming), especially graupel and hail, can throw off small ice splinters, which can grow into SIPs. This phenomenon was first investigated by Macklin (1960). Various physical mechanisms have since been suggested: Macklin (1960) reported that riming causes the growth of fine ice structures that would break to create SIPs, while Choularton et al. (1978) suggested that shattering of freezing droplets produces SIPs. Hallett and Mossop (1974) counted ice crystals in a light beam beneath a metal rod, which moved through the cloud chamber and swept up water droplets, causing riming. By comparing the enhanced number of ice crystals to the background number of the cloud chamber, they derived a profile depending on cloud temperature. The standard parameterization for RS is based on this study and is often referred to as the Hallett–Mossop process.

Emersic and Connolly (2017) investigated riming events using high-speed cameras and concluded that even small ice crystals and not only large rimers like graupel could potentially produce a significant amount of SIPs. Seidel et al. (2024) investigated rime splintering with high-speed video microscopy, IR thermography, and a custom-built ice counter. They could not reproduce the results of Hallett and Mossop (1974). In general, they found only insignificant amounts of SIP production during riming, which cannot explain the amount of SIPs expected in convective and frontal clouds.

Overall, Korolev and Leisner (2020) found no consistency in measured and estimated rime splintering rates between various groups and attributed the discrepancy to different laboratory setups and techniques. Despite RS being the most commonly included secondary ice mechanism in numerical cloud models, the physical understanding is severely lacking, and thus the development of a physically based parameterization seems unfeasible at the moment (Korolev and Leisner, 2020). Still, we use the parameterization of RS based on the dataset of Hallett–Mossop in this study as it is part of the standard SB scheme.

In Hallett and Mossop (1974), rime splintering occurred within a narrow temperature range of $T_{\mathrm{rs},\min }=\mathrm{265}$ and $T_{\mathrm{rs},\max }=\mathrm{270}\mathrm{K}$. A fit onto their dataset is used as a function of temperature and rimed mass q_rime in kilograms as it is implemented in SB:

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\displaystyle \frac{\partial n_{\sec }}{\partial t}}& =q_{\mathrm{rime}}{C}_{\mathrm{RS}}{\displaystyle \frac{T-T_{\mathrm{rs},\min }}{T_{\mathrm{rs},\mathrm{opt}}-T_{\mathrm{rs},\min }}}{\displaystyle \frac{T-T_{\mathrm{rs},\max }}{T_{\mathrm{rs},\mathrm{opt}}-T_{\mathrm{rs},\max }}}\\ {\displaystyle \frac{\partial q_{\sec }}{\partial t}}& ={\displaystyle \frac{\partial n_{\sec }}{\partial t}}{x}_{\mathrm{i},\min },\end{array}\end{array}$$

with the multiplication factor $C_{\mathrm{RS}}=\mathrm{3.5}\times {\mathrm{10}}^{\mathrm{8}}{\mathrm{kg}}^{-\mathrm{1}}$ and the optimal temperature $T_{\mathrm{rs},\mathrm{opt}}=\mathrm{268}\mathrm{K}$, minimum temperature $T_{\mathrm{rs},\min }=\mathrm{265}\mathrm{K}$, and maximum temperature $T_{\mathrm{rs},\max }=\mathrm{270}\mathrm{K}$ for rime splintering. The rimed mass q_rime represents the mass content of cloud droplets q_c converted into graupel q_g and hail q_h. For riming, the model considers the collision of supercooled droplets with ice crystals, snow, graupel, and hail. The triangular profile is centered at T_rs,opt and ejects up to 350 splinters per milligram rimed ice mass. The fragments are initialized as SIPs of the minimum ice mass $x_{i,\min }={\mathrm{10}}^{-\mathrm{12}}\mathrm{kg}$. Additional sources for SIPs from raindrop freezing and shattering as well as collisional breakup will be presented in another study.

2.3 Depositional growth of ice particles

Deposition of water vapor onto a cloud particle in a supersaturated environment is the most fundamental growth process for all cloud particles. The SB scheme uses saturation adjustment for the water phase; hence whenever the environment is supersaturated with respect to water (S_w>1) at the beginning of a microphysical time step, S_w is relaxed to water saturation (i.e., S_w≡1), and the change in water vapor q_v is converted to cloud water q_c, considering latent heat release (using a Newton–Raphson scheme). However, supersaturation with respect to ice (i.e., S_i>1) is explicitly resolved. The mass growth of an ice particle by deposition of water vapor can be described by formulating the flux of vapor and heat between the particle and the environment considering mass and energy conservation and then integrating these fluxes over the particle surface and evaluating the effect of latent heat of sublimation on the particle surface temperature. Assuming that the temperature difference between the particle surface and the environment is small, we can derive the general growth equation of a single ice particle with mass x using the Clausius–Clapeyron relation (Pruppacher and Klett, 1998):

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\partial x}{\partial t}}{|}_{\mathrm{dep}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }C\left(x\right)(S_{\mathrm{i}}-\mathrm{1})F_{\mathrm{v}}\left(x\right)}{\frac{T{R}_{\mathrm{v}}}{D_{\mathrm{v}}{e}_{\mathrm{si}}}+\frac{L_{\mathrm{s}}\left(\frac{L_{\mathrm{s}}}{R_{\mathrm{v}}T}-\mathrm{1}\right)}{K_{\mathrm{a}}T}}},\end{array}$$

with the saturation ratio with respect to ice $S_{\mathrm{i}}=\frac{e}{e_{\mathrm{si}}\left(T\right)}$ and the ventilation coefficient F_v and capacitance C accounting for the enhancement of the depositional growth by the flow field and the non-spherical shape, respectively. The left term in the denominator represents the mass flux relation with the diffusion coefficient D_v and water vapor saturation pressure with respect to ice e_si; the right term describes the heat flux with the thermal conductivity of air K_a and the latent heat of sublimation L_s. The capacitance can be related to the maximum diameter $C\left(x\right)=D\left(x\right)/c$, with c_k depending on the class of the particle k.

Integrating (Eq. 31) over the entire particle distribution for an ice mode k, we obtain

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{g}_{\mathrm{dep},k}& ={\displaystyle \frac{\partial q_k}{\partial t}}{|}_{\mathrm{dep}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\partial x_k}{\partial t}}{|}_{\mathrm{dep}}\mathrm{d}{x}_k\\ & ={\displaystyle \frac{\mathrm{4}\mathit{\pi }\left(S_{\mathrm{i}}-\mathrm{1}\right)c_k^{-\mathrm{1}}}{\frac{T{R}_{\mathrm{v}}}{D_{\mathrm{v}}{e}_{\mathrm{si}}}+\frac{L_{\mathrm{s}}\left(\frac{L_{\mathrm{s}}}{R_{\mathrm{v}}T}-\mathrm{1}\right)}{K_{\mathrm{a}}T}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}D\left(x_k\right)F_{\mathrm{v}}\left(x_k\right)f\left(x_k\right)\mathrm{d}{x}_k\\ & ={\displaystyle \frac{\mathrm{4}\mathit{\pi }\left(S_{\mathrm{i}}-\mathrm{1}\right)c_k^{-\mathrm{1}}}{\frac{T{R}_{\mathrm{v}}}{D_{\mathrm{v}}{e}_{\mathrm{si}}}+\frac{L_{\mathrm{s}}\left(\frac{L_{\mathrm{s}}}{R_{\mathrm{v}}T}-\mathrm{1}\right)}{K_{\mathrm{a}}T}}}D\left(\stackrel{\mathrm{‾}}{x_k}\right)\stackrel{\mathrm{‾}}{F_{\mathrm{v}}},\end{array}\end{array}$$

with the mean particle mass $\stackrel{\mathrm{‾}}{x_k}$ and averaged ventilation coefficient $\stackrel{\mathrm{‾}}{F_{\mathrm{v}}}$; for details, see calculations in the appendix of SB.

However, we must consider that all ice particle classes compete for water vapor. Thus we follow the semi-analytic approach of Morrison et al. (2005) to estimate the depositional growth rate with an exponential relaxation towards equilibrium. The change in the mass content of ice particle class k within a physical time step is

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{\mathrm{\Delta }{q}_k}{\mathrm{\Delta }t}}{|}_{\mathrm{dep}}={\displaystyle \frac{{\mathit{\delta }}_{i}X}{{\mathit{\tau }}_k}}\left(\mathrm{1}-\exp \left(-{\displaystyle \frac{\mathrm{\Delta }t}{X}}\right)\right),\end{array}$$

with ${\mathit{\delta }}_i=q_{\mathrm{v}}-q_{\mathrm{si}}$ expressing the supersaturation with respect to ice and the microphysical time step Δt. τ_k is defined as the depositional timescale

$$\begin{array}{}\text{(34)}& {\mathit{\tau }}_k^{-\mathrm{1}}={\displaystyle \frac{g_{\mathrm{dep},k}}{{\mathit{\delta }}_i}}.\end{array}$$

X is the conjoined relaxation time needed to describe the competition between the ice particle classes as a sum of the individual relaxation timescales:

$$\begin{array}{}\text{(35)}& X={\left[\sum _j^N{\mathit{\tau }}_j^{-\mathrm{1}}\right]}^{-\mathrm{1}},\end{array}$$

where we consider graupel, hail, and snow, as well as all five ice modes. For a derivation and application of this method to resolve supersaturation, see, e.g., Khvorostyanov (1995), Morrison et al. (2005), and Köhler and Seifert (2015).

2.4 Sedimentation of cloud particles

For the sedimentation of ice particles, we use the mass–velocity relation (Eq. 5) for ice crystals; these relations were used for all new ice modes. For use in the evolution equation of moments (Eq. 6), an averaged sedimentation velocity must be calculated via a weighted integration; i.e.,

$$\begin{array}{}\text{(36)}& {\stackrel{\mathrm{‾}}{v}}_{i,k}={\displaystyle \frac{\mathrm{1}}{M_i^k}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{x}^k{f}_i\left(x\right)v_i\left(x\right)\mathrm{d}x.\end{array}$$

Using the analytical results for generalized gamma distributions, this results in the following expression:

$$\begin{array}{}\text{(37)}& {\stackrel{\mathrm{‾}}{v}}_{i,k}={\mathit{\alpha }}_i{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{k+{\mathit{\nu }}_i+{\mathit{\beta }}_i+\mathrm{1}}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{k+{\mathit{\nu }}_i+\mathrm{1}}{{\mathit{\mu }}_i}\right)}}{\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{1}}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{2}}{{\mathit{\mu }}_i}\right)}}\right]}^{{\mathit{\beta }}_i}{\stackrel{\mathrm{‾}}{x}}_i^{{\mathit{\beta }}_i}.\end{array}$$

See also Sect. 3.7 in SB.

2.5 Collision processes of cloud particles

Collision is an important process for producing large cloud particles. For the formation of large raindrops, collision is essential. However, considering the formation of rain via the ice phase, collisions between water and ice particles again play an important role. In the presented scheme, the rates for collisions between liquid particles (cloud droplets and raindrops) remain unchanged. However, the collision rates for ice particles must be reconsidered.

As formulated in the SB scheme, we have the following conceptual treatment for sorting new particles stemming from the collision of ice particles with others:

• Collision between a cloud ice particle and a cloud droplet leads to a larger ice particle or to a graupel particle, depending on the size of the ice particle (riming).

• Collision between two cloud ice particles leads to a snow particle (self-aggregation).

• Collision between a cloud ice particle and a snow particle leads to a snow particle (aggregation).

This concept must be adapted for ice particles from the five different ice modes. The collision of ice particles from different ice modes in particular must be taken into account as a new process, extending the existing formulation of self-aggregation of ice particles. In the following the changes are documented.

2.5.1 Collision of liquid and solid particles

We adapt the formulation by SB of riming (ice particle collides with cloud droplet or raindrop) for all new ice modes. Analogous to the implementation for a single ice class in SB, all ice modes can be collected by graupel, hail, or snow and contribute to riming.

2.5.2 Aggregation of ice particles

There is no separate treatment for each ice mode regarding aggregation; all collisions between ice modes and self-collection within a single ice mode contribute to the same snow class. The collision processes between two ice classes i and j leads to the formation of snow “s” and can be described as

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{\partial f_i\left(x\right)}{\partial t}}{|}_{\mathrm{coll},ij}=-\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{f}_i\left(x\right)f_j\left(x^{\prime }\right)K_{ij}(x,x^{\prime })\mathrm{d}{x}^{\prime },\text{(39)}& {\displaystyle \frac{\partial f_j\left(x\right)}{\partial t}}{|}_{\mathrm{coll},ij}=-\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{f}_i\left(x^{\prime }\right)f_j\left(x\right)K_{ij}(x,x^{\prime })\mathrm{d}{x}^{\prime },\text{(40)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial f_{\mathrm{s}}\left(x\right)}{\partial t}}{|}_{\mathrm{coll},ij}=& \underset{\mathrm{0}}{\overset{x}{\int }}{f}_i\left(x^{\prime }\right)f_j(x-x^{\prime })\\ & K_{ij}(x^{\prime },x-x^{\prime })(x-x^{\prime })\mathrm{d}{x}^{\prime }.\end{array}\end{array}$$

The collection kernel is defined as

$$\begin{array}{}\text{(41)}& K_{ij}\left(x_i,x_j\right)=A_{ij}{E}_{\mathrm{coll}}\left|v_i\right(x_i)-v_j(x_j\left)\right|,\end{array}$$

where A_ij is the cross section of the sweep-out volume, E_coll is the mean sticking efficiency, and v_i and v_j are the terminal velocities. Let D_i and D_j denote the maximum diameters of ice crystals, then we obtain the following for A_ij:

$$\begin{array}{}\text{(42)}& A_{ij}\left(x_i,x_j\right)={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\left(D_i(x_i\right)+D_j\left(x_j\right){)}^{\mathrm{2}},\end{array}$$

where we simplified the expression by assuming spherical particles. The sticking efficiency is parameterized following Cotton et al. (1986):

$$\begin{array}{}\text{(43)}& E_{\mathrm{coll}}=\min ({\mathrm{10}}^{\mathrm{0.035}\cdot (T-T_m)-\mathrm{0.7}},\mathrm{0.2}),\end{array}$$

with a maximum efficiency of 0.2 for T=273 K. The mass relation provided in Eq. (5) is not sufficient to characterize the terminal velocity of ice particles since complex shapes and atmospheric turbulence affect the flow field around the particle (Seifert, 2002). The velocity distribution function describes the probability of an ice particle of mass x to have the terminal velocity v^′:

$$\begin{array}{}\text{(44)}& P\left(v^{\prime }\right|x)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{vel}}\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left[{\displaystyle \frac{v^{\prime }-v\left(x\right)}{{\mathit{\sigma }}_{\mathrm{vel}}}}\right]}^{\mathrm{2}}\right],\end{array}$$

with the variance of terminal velocity σ_vel.

The absolute velocity difference in the collision kernel would split the integral and impose incomplete gamma functions, thus complicating the result. We use the Wisner approximation, which assumes characteristic mean values for the terminal velocities that are constant and can hence be separated from the integrand. As proposed in Seifert (2002), the integral cross section is used as a weight:

$$\begin{array}{}\text{(45)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{v_{ij,k}}}^{\mathrm{2}}=& {\displaystyle \frac{\mathrm{1}}{N_k}}\underset{x_i=\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{x_j=\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{D}_i^{\mathrm{2}}\left(x_i\right)D_j^{\mathrm{2}}\left(x_j\right)f_i\left(x_i\right)f_j\left(x_j\right)x_i^k\\ & \underset{v_j=-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{v_i=-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\left(v_i^{\prime }\left(x_i\right)-v_j^{\prime }\left(x_j\right)\right)}^{\mathrm{2}}\\ & P\left(v_i^{\prime }|x_i\right)P\left(v_j^{\prime }|x_j\right)\mathrm{d}{x}_i\mathrm{d}{x}_j\mathrm{d}{v}_i^{\prime }\mathrm{d}{v}_j^{\prime }\\ =& \stackrel{\mathrm{‾}}{v_i^{\mathrm{2}}}{\mathit{\vartheta }}_i^k+\stackrel{\mathrm{‾}}{v_i}\stackrel{\mathrm{‾}}{v_j}{\mathit{\vartheta }}_{ij}^k+{\stackrel{\mathrm{‾}}{v_j}}^{\mathrm{2}}{\mathit{\vartheta }}_j^{\mathrm{0}}+\mathrm{2}{\mathit{\sigma }}_{\mathrm{vel}}^{\mathrm{2}},\end{array}\end{array}$$

where we used the properties of the gamma function and the scaling factor,

$$\begin{array}{}\text{(46)}& N_k=\underset{x_i=\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{x_j=\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{D}_i^{\mathrm{2}}\left(x_i\right)D_j^{\mathrm{2}}\left(x_j\right)f_i\left(x_i\right)f_j\left(x_j\right)x_i^k\mathrm{d}{x}_i\mathrm{d}{x}_j.\end{array}$$

Integrating Eq. (38) into Eq. (40) and again using the properties of the gamma function, we finally obtain an analytical solution for the collision rate:

$$\begin{array}{}\text{(47)}& \begin{array}{rl}{\displaystyle \frac{\partial n_{\mathrm{s}}}{\partial t}}{|}_{\mathrm{coll},ij}& =-{\displaystyle \frac{\partial n_i}{\partial t}}{|}_{\mathrm{coll},ij}=-{\displaystyle \frac{\partial n_j}{\partial t}}{|}_{\mathrm{coll},ij}\\ & ={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{E}_{\mathrm{coll}}{n}_i{n}_j[{\stackrel{\mathrm{‾}}{D_i}}^{\mathrm{2}}{\mathit{\delta }}_i^{\mathrm{0}}+\stackrel{\mathrm{‾}}{D_i}\stackrel{\mathrm{‾}}{D_j}{\mathit{\delta }}_{ji}^{\mathrm{0}}\\ & +{\stackrel{\mathrm{‾}}{D_j}}^{\mathrm{2}}{\mathit{\delta }}_j^{\mathrm{0}}]/\stackrel{\mathrm{‾}}{v_{ij,\mathrm{0}}},\end{array}\end{array}$$

and the change in mass contents due to aggregation,

$$\begin{array}{}\text{(48)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial q_i}{\partial t}}{|}_{\mathrm{coll},ij}& =-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{E}_{\mathrm{coll}}{q}_i{n}_j[{\stackrel{\mathrm{‾}}{D_i}}^{\mathrm{2}}{\mathit{\delta }}_i^{\mathrm{1}}+\stackrel{\mathrm{‾}}{D_i}\stackrel{\mathrm{‾}}{D_j}{\mathit{\delta }}_{ji}^{\mathrm{1}}\\ & +{\stackrel{\mathrm{‾}}{D_j}}^{\mathrm{2}}{\mathit{\delta }}_j^{\mathrm{0}}]/\stackrel{\mathrm{‾}}{v_{ij,\mathrm{1}}},\end{array}\text{(49)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial q_j}{\partial t}}{|}_{\mathrm{coll},ij}& =-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{E}_{\mathrm{coll}}{q}_j{n}_i[{\stackrel{\mathrm{‾}}{D_i}}^{\mathrm{2}}{\mathit{\delta }}_i^{\mathrm{0}}+\stackrel{\mathrm{‾}}{D_i}\stackrel{\mathrm{‾}}{D_j}{\mathit{\delta }}_{ij}^{\mathrm{1}}\\ & +{\stackrel{\mathrm{‾}}{D_j}}^{\mathrm{2}}{\mathit{\delta }}_j^{\mathrm{1}}]/\stackrel{\mathrm{‾}}{v_{ji,\mathrm{1}}},\end{array}\text{(50)}& {\displaystyle \frac{\partial q_{\mathrm{s}}}{\partial t}}{|}_{\mathrm{coll},ij}=-{\displaystyle \frac{\partial q_i}{\partial t}}{|}_{\mathrm{coll},ij}-{\displaystyle \frac{\partial q_j}{\partial t}}{|}_{\mathrm{coll},ij}.\end{array}$$

See Appendix C for the notation of the coefficients. Note that with the inclusion of four additional cloud ice classes, the number of collision processes contributing to snow has increased. There are 4 additional calls to ice self collection and 10 new calls to collisions between ice modes. The equations of the collision rates introduced in this section are in general non-linear. Hence, the increased number of collision processes might affect the amount of snow in the model. We analyze this bias (compared to the standard SB scheme) in Sect. 4.2.

3 Methods

Idealized simulations of a deep convective cloud were conducted with the Icosahedral Nonhydrostatic (ICON) weather and climate model version 2.6.

3.1 Model setup

The setup of the ICON model largely follows the setup described in Heinze et al. (2017). ICON is run in limited-area mode with constant boundary conditions derived from the initial conditions. The initial conditions are provided by the Weisman–Klemp setup described in Sect. 3.2. The domain extends over 2° in longitude and latitude without topography.

The model grid is an unstructured triangular grid, R02B11, with 2 initial divisions and 11 bisecting iterations (Wan et al., 2013). The effective horizontal resolution of the grid is ∼1.23 km. The vertical grid contains 128 levels. The physical time step is 2 s. The simulation time for each experiment is 240 min. Output is written every 10 min and interpolated on a regular longitude–latitude grid with a spacing of 0.005 °× 0.005 °.

The model uses parameterizations for subgrid-scale orographic drag (Lott and Miller, 1997) and non-orographic gravity wave drag (Orr et al., 2010), and the turbulence and vertical diffusion scheme is from Raschendorfer (2001). The setup does not use a surface scheme. Radiation transfer physics are disabled. Convection is treated explicitly; i.e., all convection parameterizations (shallow and deep) are switched off.

3.2 Weisman–Klemp setup

The Weisman–Klemp setup (Weisman and Klemp, 1982, WK82 hereafter) is a suitable test environment for evaluating the ice mode scheme microphysics, as all microphysics classes and associated model routines are active, and both mixed-phase and pure ice clouds – the latter mostly located in the anvil – are observed. The WK82 setup is frequently used to achieve confidence in newly developed models (Zängl et al., 2015). It has also been utilized in many studies to investigate convection or cloud microphysics (see, e.g., Bluestein and Weisman, 2000; Takemi, 2006; Miglietta and Rotunno, 2014; Rousseau-Rizzi et al., 2017; Miltenberger et al., 2020; Barrett and Hoose, 2023). The WK82 setup is available in ICON version 2.6 as one of the standard test cases.

The vertical profiles of temperature and relative humidity are chosen following WK82 to represent a typical sounding of a deep convection event with a maximum water vapor mixing ratio near the surface of 0.014 kg kg⁻¹ and a surface temperature of 300 K. The horizontal wind u_z varies with altitude from 0 to $U_{\mathrm{s}}=\mathrm{5}{\mathrm{ms}}^{-\mathrm{1}}$. Convection was initiated with a symmetric thermal perturbation of 1400 m vertical extent and 10 000 m diameter with a temperature amplitude of 2 K at the center, which decays to 0 K at the edge of the bubble. Temperature, wind, and moisture profiles of the setup can be found in WK82.

Table 1 List of sensitivity experiments for idealized simulations.

3.3 Experiments

Experiments with the ice mode scheme in several configurations were performed to evaluate the impact of different heterogeneous nucleation schemes on the distribution of the ice modes and the liquid origin vs. in situ formation pathway. Table 1 lists all experiments. HA, UL, and PH label experiments using the heterogeneous ice nucleation schemes of Hande et al. (2015), Ullrich et al. (2017), and Phillips et al. (2008), respectively. The choice of parameters for the dust aerosol profiles is explained in Sect. 2.2.3. PH3 refers to a simulation with the Phillips et al. (2008) parameterization including the soot and biological aerosol modes. The HA experiment denotes the default setup for the ice mode scheme, as the Hande et al. (2015) scheme is the most widely used heterogeneous ice nucleation scheme for ICON simulations with the SB microphysics scheme.

The reference simulation (REF) uses the standard SB scheme with only a single class for cloud ice. However, the same microphysical parameterizations and assumptions as the ice mode microphysics scheme are employed. In particular, REF uses the same setup for the heterogeneous ice nucleation as the HA experiment. Thus we can assess the ramifications of using multiple ice classes instead of a single one for the overall cloud evolution.

3.4 Ice cloud origin

We investigate the origin of cloud ice by introducing the ice mode mass fractions f_x for ice mode x:

$$\begin{array}{}\text{(51)}& f_x={\displaystyle \frac{q_x}{q_{\mathrm{tot}}}}.\end{array}$$

Thus, f_x describes the relative contribution of ice mode x to the total cloud ice at each grid point. Insignificant trace amounts of cloud ice can distort the statistics of ice mode mass fractions. Hence, we only calculate mass fractions for grid points where the (total) cloud ice mass content exceeds 0.1 mg kg⁻¹. This threshold value was used by previous studies that investigated ice cloud origins to identify regions where ice clouds of significant thickness are present (Wernli et al., 2016).

f_liq is the liquid origin mass fraction describing the ratio of cloud ice stemming from liquid origin formation processes:

$$\begin{array}{}\text{(52)}& f_{\mathrm{liq}}={\displaystyle \frac{q_{\mathrm{frz}}+q_{\mathrm{imm}}+q_{\sec }}{q_{\mathrm{frz}}+q_{\mathrm{imm}}+q_{\mathrm{hom}}+q_{\mathrm{dep}}+q_{\sec }}},\end{array}$$

where HOM and DEP contribute to in-situ-formed cirrus and FRZ and IMM to liquid origin cirrus.

4 Results

In this section we present the results of the idealized simulations of deep convection. In Sect. 4.1 we discuss the evolution and spatial distribution of hydrometeor classes and the ice modes, in particular for the default ice mode experiment (see Sect. 3.3). In Sect. 4.2 we compare the HA and REF simulations to validate the new ice mode scheme against the established SB scheme. In Sect. 4.3 we investigate the impact of different heterogeneous ice nucleation schemes on the ice modes by comparing simulations HA, UL, PH, and PH3. In Sect. 4.4 we address the research question of liquid vs. in situ ice formation in the convective cloud.

Figure 3 Temporal evolution of the horizontally averaged ice water path for each hydrometeor class for (a) ice modes (HA) and (b) the reference SB simulation (REF). Temporal evolution of average integrated number concentration for (c) HA and (d) REF. “Ice” is the sum of all ice modes for HA.

Figure 4 Temporal evolution of the horizontally averaged (ice) water path for each ice mode for the heterogeneous ice nucleation schemes (a) HA, (b) PH, (c) PH3, and (d) UL. Average integrated number concentration for each ice mode for (e) HA, (f) PH, (g) PH3, and (h) UL.

4.1 Ice mode simulation

Figure 3 shows the temporal evolution of the different hydrometeor classes as the domain-averaged (ice) water path, (I)WP (a, b), and vertically integrated number concentrations, (I)WPN (c, d), beneath the cloud system. Ice is displayed as the sum of all ice modes (TOT). Note that the output interval of the model is 10 min, and we only display domain-averaged water path values above ${\mathrm{10}}^{-\mathrm{6}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$. Hence, the plots only show the time where hydrometeor classes appear in a “significant” number.

First, we focus on the experiment with the default setup for the ice mode scheme (HA) in Fig. 3a and c. Panels (b) and (d) are discussed in Sect. 4.2. In panel (a) we note that cloud droplets condense within the first 10 min of the simulation, followed by the initiation of the warm-rain process. Ice and graupel appear at the same time around 18 min and quickly become, within 5 min, the two major contributors to the water path. Snow becomes relevant with a short delay of up to 5 min after the occurrence of graupel. Strong riming is the last microphysical process to become active, with hail first occurring at 30 min. At a 60 min simulation time the cell matures. Hence, all hydrometeor classes are present, and their avg water path increases only slightly after that. Note that the cell shows no signs of dissipation within the simulation time. The ratio of the avg water path between the classes also remains constant. Graupel, followed by cloud ice, is the dominate class with regard to IWP, while hail, rain, and cloud droplets show mostly the same avg IWP values. Snow is distinctly the weakest hydrometeor class with regard to avg IWP. The mass content of cloud particles in general is primarily governed by the availability of water vapor and thus thermodynamics. The number concentrations of the classes, on the other hand, are more varied since they directly depend on the formation pathways. Panel (c) in Fig. 3 shows the avg integrated number concentrations (IWPNs) of all hydrometeor classes. We observe a large gap in avg IWPNs between cloud droplets and ice compared to the other classes, with an increase of 3.2 and 4 orders of magnitude for simulation times after 130 min, respectively. Cloud droplet number concentrations are tied to the cloud condensation nuclei (CCN) activation scheme of Hande et al. (2016), which depends mainly on upward vertical velocity and thus produces a large droplet number concentration in a strong convective case like the WK82 setup. However, the scheme is still limited by the number of available CCN. We discuss the connection between high cloud ice and cloud droplet number concentrations in the following paragraphs.

Figure 4 provides insight into the development of each ice mode as domain-averaged IWP and IWPN. We focus on the results of the default setup of the ice mode schemes (HA) first (panels a and e). At the beginning, ice consists of frozen cloud droplets with both immersion and homogeneous freezing, represented by the IMM and FRZ modes with avg IWP values of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{1}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$, respectively. FRZ is the most dominant ice mode by several orders of magnitude, especially with regard to avg IWPN. Secondary (SEC) ice from rime splintering occurs with a 20 min offset to the first occurrence of graupel and becomes almost as important as immersion freezing in regard to avg IWP with values up to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ at the end of the simulation. This is consistent with the findings of Miltenberger et al. (2020), who also performed simulations with the ice mode schemes for the same test case but utilized Lagrangian trajectory analysis to investigate rime splintering and spread of secondary ice through the cloud system. The first occurrence of secondary ice of significant number concentrations ($n_{\sec }>\mathrm{0.1}{\mathrm{l}}^{-\mathrm{1}}$) was observed after 30 min, even for simulations of higher wind shear. They found riming with graupel that sedimented from higher levels into the Hallett–Mossop temperature zone ($\mathrm{265}<T<\mathrm{270}\mathrm{K}$) to be the dominant source of secondary ice.

Ice from homogeneous (HOM) and deposition (DEP) nucleation occurs first from nucleation events in the overshoot starting at 30 and 40 min for DEP and HOM, respectively. Both types of nucleation events occur multiple times during the simulation, and the avg IWP of both modes increases over time. Initially, DEP is overall stronger in terms of avg IWP, although with significant IWP ($\ge {\mathrm{10}}^{-\mathrm{6}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$) occurring at 95 min. However, at the end of the simulation, HOM shows almost the same avg IWP and IWPN values as DEP, with $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ and $\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{2}}$, respectively.

In summary, FRZ dominates with regard to avg IWP with up to 2.5 orders of magnitude and with regards to avg IWPN with up to 4 orders of magnitude. These results are strongly tied to the CCN activation scheme. Droplets freeze preferably heterogeneously (here by immersion freezing) and then homogeneously once the active INPs for immersion freezing are depleted. Thus, if more CCN than INPs are available, the FRZ mode becomes more dominant than the IMM mode as more (activated) droplets will freeze homogeneously once the INPs are depleted. The integrated total ice number concentration (as the sum of all ice modes) is overall dictated by FRZ in this case and thus tied to the high number concentrations of droplets that the CCN scheme activates.

Figure 5 Vertical slice through domain center at 200 min simulation time. Mass content of (a) total ice, (b) snow, (c) graupel, (d) cloud droplets, (e) raindrops, and (f) hail. The black line indicates the melting temperature level. The red and blue contours show the critical LWC and IWC value of 0.1 mg m⁻³, respectively.

Figure 6 Vertical slice through domain center at 200 min simulation time. Mass content of (a) sum of all ice modes, (b) secondary ice, (c) heterogeneously frozen droplets, (d) homogeneously frozen droplets, (e) homogeneous nucleation, and (f) deposition nucleation. The black line indicates the melting temperature level. The red and blue contours show the critical LWC and IWC value of 0.1 mg m⁻³, respectively.

We now investigate the vertical distribution of cloud particles. Figure 5 shows a vertical slice of the cloud through the domain center after 200 min simulation time and the mass content of all hydrometeor classes, where ice is represented as the sum of all ice modes. The red and blue contours show the critical LWC and IWC value of 0.1 mg m⁻³, respectively. The matured convective cell has a vertical extent of 15 km and developed an anvil between 9 and 12 km. The cell produces strong precipitation beneath the cell, as can be seen from panel (e) for rain. Most of the cloud is glaciated, with rimed particles, mainly graupel, being the most prevalent in terms of mass, with mass contents of up to ${\mathrm{10}}^{-\mathrm{3}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. Hail is mostly found in the narrow updraft core at 130 km, where riming rates and resident times are the highest. Ice is present in most parts of the cloud above the melting temperature and forms an anvil of (pure) ice with an overshooting top above 12 km. Snow is a product of aggregation and thus depends on the mass content and number concentrations of ice. Cloud regions with high cloud ice number concentrations and mass content have high aggregation rates. However, snow has a lower mean size than cloud ice and thus higher sedimentation rates. Thus, snow is mostly located below the maximum cloud ice mass content.

Figure 6 shows the same vertical slice as in Fig. 5 but for the mass content of each ice mode after 200 min simulation time. The lower part of the cloud from 3 to 9 km around the updraft core is dominated by secondary ice from rime splintering (SEC) with mass contents of up to ${\mathrm{10}}^{-\mathrm{5}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. This coincides with the region where riming is active, which is evident from the high mass content of graupel and hail in this region, as observed in Fig. 5. The prevalence of secondary ice from rime splintering in this region was also observed by Miltenberger et al. (2020).

Ice between 6 and 10 km partially stems from heterogeneously frozen droplets (IMM), while in the upper part of the cloud, homogeneously frozen droplets (FRZ) are prevalent. The latter shows the highest mass contents for cloud ice with ${\mathrm{10}}^{-\mathrm{4}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ within a wide region extending from altitudes of 8 to 12 km and the length of the convective cloud (150 km). When cloud droplets are transported upwards into colder regions, they freeze first through immersion freezing due to the lower temperature threshold. Once all the INPs are activated, there is still a significant number of cloud droplets left or provided by the constant supply of humidity (and then saturation adjustment) from convection. At this point homogeneous freezing begins to dominate over immersion freezing. Moreover, the ice mode scheme only shows the ice formation pathway of the ice observed at this time; it does not track sources and sinks individually. This means that ice from IMM could have been converted to snow by aggregation or lost in collisions with graupel and is for this reason not as prevalent as FRZ.

The part of the cloud consisting of ice from FRZ stays structurally mostly the same during the simulation run time, only widening as the anvil expands. However, the cloud parts where IMM and SEC are present fluctuate with the evolution of the cell. SEC follows the liquid core where riming occurs. IMM is at the beginning located above the liquid core between the sections where SEC and FRZ are prevalent and is later mostly present in the flanks.

Ice from homogeneous (HOM) nucleation is limited to the overshoot, and deposition (DEP) nucleation also occurs in small amounts outside the liquid core in the glaciated parts of the cloud between 6 and 9 km.

Overall, we notice a mostly clear spatial separation of the ice modes with some overlapping areas where ice of a different origin is found. Ice in the upper part of the cloud above 9 km is mostly of liquid origin, hence consisting of frozen droplets and stemming from rime splintering below them. This general spatial distribution of the ice modes is also true for the other sensitivity experiments.

Figure 7 Vertical slice through domain center at 200 min simulation time. Mass content of snow for the (a) ice mode scheme (HA) and (b) standard SB (REF) as well as the sum of all cloud ice and snow for (c) HA and (d) REF. The black line indicates the melting temperature level. The red and blue contours show the critical LWC and IWC value of 0.1 mg m⁻³, respectively.

4.2 Comparison – reference simulation

In this section we compare the reference simulation (REF) performed with the (unmodified) SB scheme to the ice mode simulation (HA). As elaborated on in Sect. 3.3, both experiments use the same setup for heterogeneous nucleation (HA17). Thus the only difference between them is the split of the cloud ice class into the five ice modes.

The overall cloud evolution of the reference (REF) simulation is very close to the ice mode scheme (HA). As can be seen in Fig. 3b, the evolution of microphysics classes directly related to the liquid phase, which is cloud droplets, rain, graupel, and hail, shows the same temporal evolution. Even the onset times for cloud ice, where the governing parameterizations of both schemes differ, are mostly the same. We observe, however, a large difference in the average water path of snow, which shows an increase of an order of magnitude for the simulation with the ice mode scheme (HA), although the difference decreases with time to a factor of ∼2.5.

This is also evident in the vertical slice of the convective cloud. Figure 7 shows the vertical slice of the convective cloud at 200 min through the domain center for snow with (a) the ice mode simulation and (b) the reference simulation. There is significantly more snow present in the ice mode simulation than in the reference simulation. While the maximum snow mass content is still at 10⁴ kg m⁻³, snow is distributed over a wider area of the cloud, especially towards the flanks and higher altitudes (up to 11 km). However, comparing the sum of ice and snow for both simulations in panels (c) and (d) shows that there is in general no increase in combined ice particle mass content. Rather, we notice there is a tendency for the ice mode scheme (represented by the HA experiment) to shift mass from the ice to the snow class(es); hence the aggregation rates have to be higher. This is linked to the conceptional differences between both schemes: the SB standard scheme only produces snow by self-collection of a single cloud ice class, where in the ice mode schemes, there are five independent cloud ice classes. They aggregate snow not only by self-collection, but also by collisions with each other. While the physics of self-collection and ice mode–ice mode collision are the same, the increased number of collision processes leads to higher aggregation rates (see also Sect. 2.5.2).

While the shift of ice to snow should affect overall sedimentation rates of ice particles, there seems to be little impact on the cloud evolution and dynamics of the convective cell since the IWC and LWC outlines as well as the general evolution of the cloud observed in Fig. 3 remain the same. It is also important to note that pure ice and especially cirrus clouds will be mostly unaffected by this effect since collision efficiencies are small at low temperatures (see Sect. 2.5.2). Generally, ice number concentrations are also insufficient to produce significant amounts of snow for pure ice clouds.

Figure 8 Vertical slice through domain center at 200 min simulation time comparing the mass content of IMM (a–d), FRZ (f–h), and snow (i–l) for the experiments HA (a, e, i), PH (b, f, j), PH3 (c, g, k), and UL (d, h, l). The black line indicates the melting temperature level. The red and blue contours show the critical LWC and IWC value of 0.1 mg m⁻³, respectively.

4.3 Comparison – heterogeneous ice nucleation

A change in assumptions regarding the distribution of INPs or a switch to a different INP activation scheme entirely can have a large impact on the evolution of the ice phase. We compare the experiments HA, PH, PH3, and UL here. For the general evolution of an average cloud water path and ice water path, we do not observe major differences (not shown here). Hence, we discuss the general evolution of the convective cloud only briefly here. Since the same CCN activation scheme is used for all experiments, we do not observe significant change in the evolution of the liquid phase. The rainwater path is slightly lower for PH after the initial warm-rain formation but catches up later. Hail formation starts sooner for PH but also reaches the same average value as for the other experiments. Further, we observe a shift towards more snow and slightly less (total) cloud ice for the PH simulation. Thus, the overall evolution of the convective clouds is similar for all heterogeneous ice nucleation schemes tested here.

We now investigate the spatial and temporal distribution of the ice modes more closely for HA, PH, PH3, and UL. When studying the evolution of the ice mode average IWP in Fig. 4a to d, we observe that FRZ and SEC stay mostly the same for all runs at an avg IWP of 10⁻³ and ${\mathrm{10}}^{-\mathrm{1}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$, respectively. SEC stems from rime splintering and thus depends on the evolution of the graupel and hail class, which again is mostly sensitive to CCN activation and the resulting number of supercooled droplets. We observe in Sect. 4.1 that SEC occupies a distinct part of the lower cloud around the updraft core. A change in the heterogeneous nucleation scheme seems not to significantly affect the dominance of ice from rime splintering in this region or change the avg IWP from secondary ice.

FRZ is also strongly sensitive to the CCN activation scheme, which remains unchanged between the simulations. A weaker IMM mode should affect FRZ since droplets that do not freeze heterogeneously freeze homogeneously instead, unless they are removed by evaporation or collision. Hence, both modes are in direct competition for unfrozen cloud droplets. However, FRZ is overall so dominant in this case study that this effect is not noticeable with regard to avg mass content and number concentrations.

Compared to HA, IMM is weaker in the simulation with the PH scheme for both avg IWP and IWPN (Fig. 4b, f), with values of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ and 10⁶ m⁻², respectively. However, for UL (panel d, f), the avg IWP and IWPN of IMM stays the same as for HA. The latter is not surprising since the assumptions about immersion freezing and the number concentrations of INPs are almost the same. In Sect. 2.2.3 we choose the dust concentration parameters for the PH parameterization such that the maximum number of INPs activated for immersion freezing was the same as for the UL and HA schemes. The important difference between the PH scheme and the HA and UL schemes is that immersion freezing occurs at higher temperatures (close to 263 instead of 258 K), thus, in lower and warmer parts of the cloud, where hydrometeor classes other than cloud ice are abundant, e.g., graupel. Collisions with those classes are an efficient sink for the IMM mode. Hence, the average IWP and IWPN of IMM are smaller for PH compared to HA and UL.

For PH3 (Fig. 4c, h) we observe a higher avg IWP of IMM than for all other simulations, with values of up to $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{2}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$. Even though we stated in Sect. 2.2.3 that the maximum number of active INPs is lower with this scheme than for the others, the inclusion of soot and biological aerosols triggers immersion freezing at higher temperatures even close to the melting temperature level T_m=273 K. While the shift to higher activation temperatures lowered the IMM content for the PH scheme with the dust-only mode, for PH3, with the inclusion of the soot and biological aerosol mode, IMM content increased. This even “earlier” activation of INPs likely changes the cloud evolution such that the IMM mode becomes more important. Indeed, we note that the cloud shape is different for PH3 when we later present a vertical slice of the cloud.

DEP is in general weak in this convective case since deposition nucleation events usually only trigger outside mixed-phase clouds (see Sect. 2.2.3). In Fig. 4 we observe different triggering times and strengths of deposition nucleation events. However, these differences are mostly caused by non-linear realizations of the model dynamics and are not sensitive to the choice of nucleation scheme.

Figure 9 Average liquid origin and ice mode fractions for (a) HA, (b) PH, (c) PH3, and (d) UL.

While DEP is not continuously persistent in the overshoot, there are smaller nucleation events being triggered in the anvil and the flanks of the convective cloud. Also note that output is written only every 10 min; hence not all nucleation events are sampled in the output data, and ice from DEP and HOM could be removed by aggregation or evaporation before it is sampled. When using the PH or PH3 scheme DEP is almost not present at all. This is consistent with our description of the scheme in Sect. 4.1, where even weak nucleation events are only triggered at high S_i compared to UL and HA. Consequently, the HOM mode is strengthened for PH since homogeneous nucleation events are not suppressed by pre-existing ice from DEP.

Figure 8 shows vertical slices of the cloud at 200 min through the domain center. Plotted are the mass contents of ice from immersion freezing (IMM) and homogeneously frozen droplets (FRZ) as well as snow for HA, PH, PH3, and UL (first, second, third, and fourth columns, respectively). As already discussed, FRZ dominates the upper part of the convective cloud with values of up to ${\mathrm{10}}^{-\mathrm{4}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ for all heterogeneous nucleation schemes. The structure of vertically layered ice modes described in Sect. 4.1 holds true for all simulations. As immersion freezing is mostly the same for HA and UL, the distribution of ice modes and snow shows only marginal differences between both experiments (first and fourth columns).

For PH, however, we found that IMM is much weaker in terms of mass content (and number densities). This also limits its horizontal distribution, which is confined more towards the core of the cloud compared to HA and UL. This directly affects the snow class, which shows lower number densities and less spreading throughout the cloud. This underlines again that IMM is important for aggregation because it produces ice crystals large enough for efficient collection kernels (see Sect. 2.5.2).

Using PH3 affects not only the IMM mode but also the overall structure of the cloud. As discussed above, PH3 favors immersion freezing, and its biological and soot INP modes show an onset of higher temperatures compared to other heterogeneous ice schemes. We observe the result of INP activation at higher temperatures with a shift of ice mass towards the lower model levels. This can also be observed in the shape of the convective cloud. For example, the PH3 cloud shows a less developed left flank. Additionally, there is a pocket of IMM ice above the melting temperature line which is not present for HA and UL, for which immersion freezing is not possible in this temperature region. While the general structure and layering of ice modes are the same, the shape and location of liquid water zones change (as can be seen from the red LWC outlines), which has a direct effect on the SEC mode that in general follows the liquid core (where riming occurs) closely. The PH3 simulation also shows an increase in snow, underlining the importance of aggregation as a sink for the IMM mode (stronger in PH3).

4.4 Liquid origin vs. in situ formation

Figure 9 shows the average mass fraction for all ice modes at 200 min simulation time. The mass fraction describes the ratio of ice from a particular mode x to the total cloud ice q_tot (see Sect. 3.4).

The majority of the fully glaciated parts of the cloud above the secondary ice region between 6 and 12 km is dominated by IMM (for HA, PH, and UL) or FRZ (for PH3) classifying the parts of (pure) ice cloud as liquid origin (hence, LIQ is close to 1). This includes the anvil of the deep convective cloud. The overshoot is located above 12 km, and the liquid origin fraction is determined by the amount of FRZ ice being mixed into this region and the strength of DEP nucleation events. PH and PH3 show a weak DEP mode with incursions from HOM, resulting in liquid origin fractions between 0.3 and 0.7. For HA and UL the liquid origin fraction is below 0.4. This makes the overshoot a region where ice stemming from different formation pathways mixes. The dynamic forcing transporting FRZ and IMM ice into the overshoot, the occurrence and strength of HOM and DEP nucleation events, and the sedimentation of newly formed ice from the overshoot into the cloud below all affect and change the liquid origin fraction. That the majority of the fully glaciated parts of the deep convection cloud is of liquid origin is in accordance with previous studies (Gasparini et al., 2018).

5 Discussion and summary

We presented a new microphysics bulk model adapted from the Seifert and Beheng (2006) two-moment scheme by introducing multiple ice classes each with their own unique particle formation mechanism: homogeneous freezing of solution and pure water droplets, immersion freezing, deposition nucleation, and secondary ice from rime splintering. The microphysical processes governing these ice modes have been described with an emphasis on the particle formation mechanisms.

Idealized simulations of a convective cloud, using the Weisman–Klemp setup, were performed to validate the model with a comparison to a reference simulation using the standard SB scheme. The Weisman–Klemp test is a highly idealized description of the typical life cycle of a deep convective scenario. However, it contains the major features of the development of a convective storm, such as, e.g., the development of liquid clouds at lower levels, which are transformed into mixed-phase clouds, and also the typical development of the anvil cloud. Although the scenario cannot be compared on a one-on-one basis to measurements, it contains the main features of the convective life cycle. Thus, we can study the formation and evolution of ice clouds exemplarily and can derive insights into the general distribution of ice clouds and their formation pathways. Therefore, the derived statistics can serve as a first, but still idealized, evaluation of ice clouds in convective storms, but it nevertheless gives the first quantitative results about the distribution of formation pathways.

It was shown that the ice mode scheme reproduces the same cloud evolution for the dynamics and all particle classes except for snow, where we observed a shift of mass content from ice crystals to snow. This is linked to an increase in collision processes due to the introduction of multiple ice classes.

We found that the scheme showed a reasonable distribution of the ice modes with liquid origin ice, formed from homogeneously and immersion frozen droplets, constituting the majority of ice present in the matured cell. The occurrence of in situ ice, formed from homogeneous and deposition nucleation origins, was low compared to other modes, as expected for a convective cloud. Most in situ ice modes were only present in the overshoot. But even there they mixed with liquid origin ice modes that were transported into the overshoot.

Simulations with four different heterogeneous nucleation parameterizations showed that the ice mode scheme provides the basis for a nuanced analysis to evaluate the impact of parameterization choice not only on the total ice content and number concentrations but also on the competition between ice particle formation pathways.

We investigated rime splintering as a secondary ice mechanism. In general we observed that it did not significantly enhance ice content in parts of the convective clouds where primary ice formation was active. Rather, secondary ice was relevant in thermodynamic regions where primary ice formation was insufficient. Thus, it helped to expand cloud ice, especially towards higher temperature levels. However, we note that the importance and high number densities of secondary ice from RS in the lower parts of the cloud (below 7 km) are concerning since the underlying Hallett–Mossop parameterization could not be confirmed by recent laboratory studies (Seidel et al., 2024). In future work we will implement additional secondary ice processes, droplet shattering and collisional breakup, and compare them to rime splintering.

In the future we will consider combining the ice mode schemes with tracking (accumulated) microphysical process rates. As discussed in the Introduction, both viewpoints are complementary as they allow us to combine information on regions where ice nucleation is active with the distribution of ice modes later. The ice mode scheme might “miss” cloud ice that formed, e.g., by immersion freezing but is quickly (within the output writing steps) removed due to collision or evaporation. Another future model extension is implementing the ice mode schemes into an ICON version that includes the Aerosol and Reactive Trace gases (ART) model component, which allows for explicit modeling of aerosol serving as ice-nucleating particles.

In the next study we are going to investigate the origin of cloud ice in the cloud band and outflow region of warm conveyor belts, especially addressing the open research question of in situ vs. liquid origin cirrus. We will also compare the classification of the ice mode scheme with other cirrus classification methods, e.g., Wernli et al. (2016).

Appendix A: List of abbreviations

Hande et al. (2015)Phillips et al. (2008)Ullrich et al. (2017)Seifert and Beheng (2006)

Table A1 List of abbreviations.

Appendix B: Generalized gamma distribution

The generalized gamma distribution is used to describe the mass or size distribution of a cloud particle class:

$$\begin{array}{}\text{(B1)}& f\left(x\right)=A{x}^{\mathit{\nu }}\exp (-\mathit{\lambda }{x}^{\mathit{\mu }}),\end{array}$$

where we use the particle mass x as the variable, and the shape and scale parameters ν and μ are prescribed. λ and A are linked to the prognostic distribution moments, the zeroth and first moments, corresponding to the number concentration M⁰=n and mass mixing ratio M¹=q, respectively. To obtain the nth mass-weighed moment of the generalized gamma distribution Mⁿ, we multiply Eq. (B1) by xⁿ and integrate over the entire domain, which yields

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}{M}^n& =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}f\left(x\right)x^n\mathrm{d}x=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}A{x}^{\mathit{\nu }+n}\exp (-\mathit{\lambda }{x}^{\mathit{\mu }})\mathrm{d}x\\ & ={\displaystyle \frac{A}{\mathit{\mu }}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{y}^{\frac{\mathit{\nu }+n+\mathrm{1}}{\mathit{\mu }}-\mathrm{1}}\exp (-\mathit{\lambda }y)\mathrm{d}y\\ & ={\displaystyle \frac{A}{\mathit{\mu }{\mathit{\lambda }}^{\frac{\mathit{\nu }+n+\mathrm{1}}{\mathit{\mu }}}}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathit{\nu }+n+\mathrm{1}}{\mathit{\mu }}}\right),\end{array}\end{array}$$

where we used the substitution y=x^μ and the relation of the gamma function ${\int }_{\mathrm{0}}^{\mathrm{\infty }}{y}^{\mathit{\zeta }}\exp (-\mathit{\eta }y)\mathrm{d}y=\mathrm{\Gamma }(\mathit{\zeta }+\mathrm{1}){\mathit{\eta }}^{-(\mathit{\zeta }+\mathrm{1})}$ for $\mathit{\zeta }>-\mathrm{1}$, which can be reduced to Euler's definition of the gamma function with the use of a suitable substitution.

Prescribing constant values for the shape parameters ν and μ, we can now solve Eq. (B2) for A and λ:

$$\begin{array}{}\text{(B3)}& {\displaystyle }\mathit{\lambda }={\left({\displaystyle \frac{n}{q}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{2}}{\mathit{\mu }}\right)}{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}\right)}}\right)}^{\mathit{\mu }},\text{(B4)}& {\displaystyle }A={\displaystyle \frac{\mathit{\mu }n}{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}\right)}}{\mathit{\lambda }}^{\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}}={\displaystyle \frac{\mathit{\mu }n}{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}\right)}}{\left({\displaystyle \frac{n}{q}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{2}}{\mathit{\mu }}\right)}{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}\right)}}\right)}^{\mathit{\nu }+\mathrm{1}},\end{array}$$

and so we can finally describe the distribution function as a function of particle mass and the two prognostic moments:

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}f(x,q,n)=& {\left[{\displaystyle \frac{x}{\stackrel{\mathrm{‾}}{x}}}\right]}^{\mathit{\nu }}{\displaystyle \frac{n}{\mathit{\mu }\stackrel{\mathrm{‾}}{x}\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}\right)}}{\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{2}}{\mathit{\mu }}\right)}{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}\right)}}\right]}^{\mathit{\nu }+\mathrm{1}}\\ & \exp \left(-{\left[{\displaystyle \frac{x}{\stackrel{\mathrm{‾}}{x}}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{2}}{\mathit{\mu }}\right)}{\mathrm{\Gamma }\left(\frac{\mathit{\nu }+\mathrm{1}}{\mathit{\mu }}\right)}}\right]}^{\mathit{\mu }}\right)\end{array}\end{array}$$

with the mean particle mass $\stackrel{\mathrm{‾}}{x}=\frac{q}{n}$.

Appendix C: Notation collision integrals

These shorthand notations are used for the coefficients in the analytical results of the collision integrals in Sect. 2.5.2 and can be calculated beforehand since they only contain fixed parameters:

$$\begin{array}{}\text{(C1)}& {\displaystyle }{\mathit{\delta }}_i^k={\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{2}{b}_i+\mathrm{1}+{\mathit{\nu }}_i+k}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{1}}{{\mathit{\mu }}_i}\right)}}{\left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{1}}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{2}}{{\mathit{\mu }}_i}\right)}}\right)}^{\mathrm{2}{b}_i+k},\text{(C2)}& {\displaystyle }\begin{array}{rl}{\mathit{\delta }}_{ji}^k=& {\displaystyle \frac{\mathrm{\Gamma }\left(\frac{b_i+\mathrm{1}+{\mathit{\nu }}_i+k}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+k}{{\mathit{\mu }}_i}\right)}}{\left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{1}}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{2}}{{\mathit{\mu }}_i}\right)}}\right)}^{b_i+\mathrm{1}}\\ & {\displaystyle \frac{\mathrm{\Gamma }\left(\frac{b_j+\mathrm{1}+{\mathit{\nu }}_j}{{\mathit{\mu }}_j}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_j+\mathrm{1}}{{\mathit{\mu }}_j}\right)}}{\left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_j+\mathrm{1}}{{\mathit{\mu }}_j}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_j+\mathrm{2}}{{\mathit{\mu }}_j}\right)}}\right)}^{b_j},\end{array}\text{(C3)}& {\displaystyle }{\mathit{\vartheta }}_i^k={\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{2}{b}_i+\mathrm{1}+\mathrm{2}{\mathit{\beta }}_i+k}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{\mathrm{2}{b}_{\mathrm{1}}+\mathrm{1}+{\mathit{\nu }}_i+k}{{\mathit{\mu }}_i}\right)}}{\left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{1}}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{2}}{{\mathit{\mu }}_i}\right)}}\right)}^{\mathrm{2}{\mathit{\beta }}_i},\text{(C4)}& {\displaystyle }\begin{array}{rl}{\mathit{\vartheta }}_{ji}^k=& {\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_j+\mathrm{2}{b}_j+{\mathit{\beta }}_j+\mathrm{1}}{{\mathit{\mu }}_j}\right)}{\mathrm{\Gamma }\left(\frac{\mathrm{2}{b}_j+\mathrm{1}+{\mathit{\nu }}_j}{{\mathit{\mu }}_j}\right)}}{\left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_j+\mathrm{1}}{{\mathit{\mu }}_j}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_j+\mathrm{2}}{{\mathit{\mu }}_j}\right)}}\right)}^{{\mathit{\beta }}_j},\\ & {\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{2}{b}_i+\mathrm{1}+{\mathit{\beta }}_i+k}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{\mathrm{2}{b}_i+\mathrm{1}+{\mathit{\nu }}_i+k}{{\mathit{\mu }}_i}\right)}}{\left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{1}}{{\mathit{\mu }}_i}\right)}{\mathrm{\Gamma }\left(\frac{{\mathit{\nu }}_i+\mathrm{2}}{{\mathit{\mu }}_i}\right)}}\right)}^{{\mathit{\beta }}_i}.\end{array}\end{array}$$