the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Predicting the morphology of ice particles in deep convection using the super-droplet method: development and evaluation of SCALE-SDM 0.2.5-2.2.0, -2.2.1, and -2.2.2

### Shin-ichiro Shima

### Yousuke Sato

### Akihiro Hashimoto

### Ryohei Misumi

The super-droplet method (SDM) is a particle-based numerical scheme that enables accurate cloud microphysics simulation with lower computational demand than multi-dimensional bin schemes. Using SDM, a detailed numerical model of mixed-phase clouds is developed in which ice morphologies are explicitly predicted without assuming ice categories or mass–dimension relationships. Ice particles are approximated using porous spheroids. The elementary cloud microphysics processes considered are advection and sedimentation; immersion/condensation and homogeneous freezing; melting; condensation and evaporation including cloud condensation nuclei activation and deactivation; deposition and sublimation; and coalescence, riming, and aggregation. To evaluate the model's performance, a 2-D large-eddy simulation of a cumulonimbus was conducted, and the life cycle of a cumulonimbus typically observed in nature was successfully reproduced. The mass–dimension and velocity–dimension relationships the model predicted show a reasonable agreement with existing formulas. Numerical convergence is achieved at a super-particle number concentration as low as 128 per cell, which consumes 30 times more computational time than a two-moment bulk model. Although the model still has room for improvement, these results strongly support the efficacy of the particle-based modeling methodology to simulate mixed-phase clouds.

Mixed-phase clouds, which are clouds comprising droplets and ice particles, appear under multiple atmospheric conditions, from the tropics to the poles, and throughout the year (Shupe et al., 2008). Accurately simulating the evolution of droplets and ice particles in mixed-phase clouds is crucial to understanding cloud dynamics, precipitation formation, water transport, radiative properties, aerosol–cloud interaction, cloud electrification, and lightning. These features are all crucial to many environmental and societal issues, such as climate change and variability, numerical weather prediction, weather modification, and icing on infrastructure (e.g., wind turbines and power lines) and aircraft (e.g., Korolev et al., 2017).

Through their 70-year history, numerical models of cloud microphysics have become increasingly sophisticated (e.g., Khain et al., 2015; Khain and Pinsky, 2018; Grabowski et al., 2019; Morrison et al., 2020). However, recent model intercomparison studies revealed that the models do not show any sign of converging toward the truth. Even the most sophisticated models do not correspond well, and the divergence in model results is as large in sophisticated models as it is in simple models (VanZanten et al., 2011; Xue et al., 2017). Mixed-phase cloud microphysics modeling is particularly challenging because we still lack a sufficient scientific understanding of mixed-phase cloud microphysics, and an algorithm appropriate for mixed-phase cloud microphysics does not exist. This study aims to address the second problem.

Every numerical model is an approximation of a phenomenon's mathematical model, which is a theoretical description that should express the system's behavior accurately. We apply a numerical scheme to construct a numerical model, which we use to produce an approximate solution of the phenomenon's underlying mathematical model for given spatiotemporal boundary conditions. This general philosophy of simulation is well documented, e.g., in Stevens and Lenschow (2001).

There are several types of cloud microphysics numerical models that are based on different levels of theoretical descriptions.

The first of these is the bulk model, which is the most widely used cloud microphysics model type (see, e.g., Khain et al., 2015; Morrison and Milbrandt, 2015; Khain and Pinsky, 2018; Grabowski et al., 2019; Morrison et al., 2020, for a review). Bulk models consider only the particle population's statistical features and are thus based on macroscopic descriptions of cloud microphysics. They solve a mathematical model that is closed in the lower moments of the distribution function of cloud droplets, rain droplets, and ice particle categories (e.g., mass and number mixing ratios). The basic premise of bulk models is that the distribution function can be determined by the lower moments, but such a universal relationship is unknown. In other words, in bulk models, to predict the time evolution of a chosen set of moments, their time derivatives are approximated by some functions of the moments being predicted, but this is not generally possible (see, e.g., Beheng, 2010). It would be also informative to note the analogy and difference between the Navier–Stokes equation and bulk models (Morrison et al., 2020), which highlights the difficulty in deriving bulk models. Therefore, for cloud microphysics, a more bottom-up approach to construct more accurate and reliable numerical models would be desired.

Kinetic description provides a more detailed microscopic mathematical model of cloud microphysics, with the evolution and motion of individual aerosol, cloud, and precipitation particles being explicitly considered. Assuming that particles are locally well mixed, particle collisions are regarded as a stochastic process. Each particle is characterized by its position and internal state, the latter of which is specified by variables known as attributes, such as size, mass, ratio of the ice crystal's minor axis to the major axis (hereafter called “aspect ratio”), velocity, and chemical composition.

Mixed-phase cloud microphysics are far more complicated than those of liquid-phase clouds, with various ice crystal formation mechanisms, diffusional growth by deposition/sublimation, diverse ice particle morphologies, ice melting and shedding, riming and wet growth, aggregation, spontaneous/collisional breakup of ice particles, and rime splintering at play (e.g., Pruppacher and Klett, 1997; Hashino and Tripoli, 2007, 2008, 2011a, b; Khvorostyanov and Curry, 2014; Khain and Pinsky, 2018). Although our scientific understanding is not yet sufficient, it is plausible that mixed-phase cloud microphysics could be accurately described under a kinetic description framework. Indeed, direct comparison with laboratory data suggests that a kinetic description could express ice particle morphology evolution accurately (Jensen and Harrington, 2015). This is crucial because ice particle morphology significantly influences the fall speed, growth by diffusion and collision, and radiative properties of ice particles. Because of their direct correspondence to elementary processes, it should also be easier to refine kinetic descriptions using laboratory measurements.

Two numerical scheme types exist for kinetic descriptions, namely bin schemes and particle-based schemes.

The development of bin schemes started independently of bulk models in the 1950s (e.g., Mason and Ramanadham, 1954; Hardy, 1963; Srivastava, 1967). For a review, see, e.g., Khain et al. (2015), Khain and Pinsky (2018), Grabowski et al. (2019), and Morrison et al. (2020).

Particle-based cloud microphysics modeling is a new approach that has emerged since the mid-2000s (e.g., Paoli et al., 2004; Jensen and Pfister, 2004; Shirgaonkar and Lele, 2006; Andrejczuk et al., 2008, 2010; Shima et al., 2009; Sölch and Kärcher, 2010; Riechelmann et al., 2012; Brdar and Seifert, 2018; Seifert et al., 2019; Jaruga and Pawlowska, 2018; Grabowski and Abade, 2017; Abade et al., 2018; Grabowski et al., 2018; Hoffmann et al., 2019). During particle-based modeling's early development, calculating the coalescence process was a numerical challenge. Shima et al. (2009), Andrejczuk et al. (2010), Sölch and Kärcher (2010), and Riechelmann et al. (2012) proposed different algorithms, and among those four schemes, the super-droplet method (SDM) developed by Shima et al. (2009) provides a computationally efficient Monte Carlo algorithm (Unterstrasser et al., 2017; Dziekan and Pawlowska, 2017). Several other coalescence algorithms were proposed in different research areas such as the weighted flow algorithm for aerosol dynamics (DeVille et al., 2011); O'Rourke's method (1981), and the no-time counter method (Schmidt and Rutland, 2000) for spray combustion; and Ormel and Spaans's method (2008) and Johansen et al.'s method (2012) for astrophysics. Li et al. (2017) confirmed that the performance of SDM is better than Johansen et al.'s method (2012), but direct comparison with other algorithms remains to be assessed.

The essential difference between bin schemes and particle-based schemes lies in the representation of particles. Bin schemes adopt an Eulerian approach and the particle distribution function is approximated using a finite number of control volumes (histogram). The time evolution is solved using a finite volume method or a finite difference method. In contrast, particle-based schemes rely on a Lagrangian approach and the population of real particles is approximated by using a population of weighted samples, sometimes referred to as super-droplets or super-particles. As discussed in Grabowski et al. (2019), bin schemes face problems that are challenging to overcome such as numerical diffusion, computational cost, and the breakdown of the Smoluchowski equation (Smoluchowski, 1916; Alfonso and Raga, 2017; Dziekan and Pawlowska, 2017). However, SDM could resolve, or at least mitigate, those problems.

Therefore, SDM and similar particle-based schemes should be more suitable for mixed-phase cloud microphysics simulations than bin schemes. Mainly because of computational costs, it is practically impossible to apply bin schemes to the most comprehensive form of kinetic description, which inevitably involves multiple attributes to express each particle's internal state. Instead, many existing bin models solve a simplified kinetic description that uses particle distribution functions with a one-dimensional attribute space approximation. For example, most rely on artificially separated categories of ice particles, with predefined mass–dimension and area–dimension relationships in each category. Another approach is adopted in the SHIPS model developed by Hashino and Tripoli (2007, 2008, 2011a, b), which is a bin model that solves sophisticated and comprehensive kinetic descriptions and does not use ice categories or mass–dimension relationships. However, to justify using the one-dimensional particle distribution function, they rely on the “implicit mass sorting assumption”, stating that different solid hydrometeor species do not belong to the same bin because they are naturally sorted by mass. Such simplifications could be a significant source of errors. SDM and similar particle-based schemes could directly simulate comprehensive kinetic descriptions with lower computational demand.

This study's primary objective is to assess particle-based modeling methodology's capability to simulate mixed-phase clouds. Therefore, we develop and evaluate the performance of a detailed numerical mixed-phase cloud model using SDM, wherein ice particle morphologies are explicitly predicted.

We first construct a mixed-phase cloud microphysics mathematical model, which is based on kinetic description. The fluid dynamics of moist air is described by the compressible Navier–Stokes equation, and aerosol, cloud, and precipitation particles are represented by point particles. Following Chen and Lamb (1994a, b) and Misumi et al. (2010), ice particles are approximated using porous spheroids. The elementary cloud microphysics processes considered in the model are advection and sedimentation; immersion/condensation and homogeneous freezing; melting; condensation and evaporation including the cloud condensation nuclei (CCN) activation and deactivation; deposition and sublimation; and coalescence, riming, and aggregation. We base the mathematical models used for those elementary processes on revised versions of existing formulas. Additionally, our model does not rely on ice categories or predefined mass–dimension relationships. For simplicity, and due to the lack of appropriate algorithms, we do not consider spontaneous/collisional breakup or rime splintering. We then develop a numerical model called SCALE-SDM to solve the mathematical model. Mixed-phase cloud microphysics is solved using the SDM. The fluid dynamics of moist air is solved by adopting a forward temporal integration scheme to both horizontal and vertical directions using a finite volume method with an Arakawa-C staggered grid. To evaluate our model's performance, we conduct a two-dimensional (2-D) simulation of an isolated cumulonimbus, and find that our model well reproduces the life cycle of a cumulonimbus typically observed in nature. The mass–dimension and velocity–dimension relationships our model predicts show a reasonable agreement with existing formulas based on laboratory measurements and field observations. We also investigate the simulation's numerical convergence and confirm that our model can produce an accurate approximate solution with lower computational demand than multi-dimensional bin schemes. We then explore the possibility of further refining and sophisticating the model; however, advancing our understanding of mixed-phase cloud microphysics is beyond the scope of this study.

Several previous works are closely relevant to this
study. Chen and Lamb (1994a, b) developed a detailed
multi-dimensional bin model, which Misumi et al. (2010) extended
and added ice volume as a new particle attribute. We follow that
strategy and approximate ice particles as porous spheroids; however,
their kinetic description is more detailed than ours because they also considered
spontaneous/collisional breakup, shedding, rime splintering, and
surface chemical reactions. They solved the
model using a multi-dimensional bin scheme; hence, their numerical model
carries a high computational cost. Hashino and Tripoli (2007, 2008, 2011a, b) further extended
Chen and Lamb (1994a, b)'s kinetic description to account for
polycrystals that can form below
−20 ^{∘}C. They solve the mathematical model using a one-dimensional bin
scheme; however, careful validation is needed to justify their implicit mass sorting assumption.
Paoli et al. (2004), Jensen and Pfister (2004), and Shirgaonkar and Lele (2006) separately developed
a particle-based model for ice-phase clouds, but neither
the evolution of ice particle morphologies nor the aggregation of ice particles were considered in their models.
Sölch and Kärcher (2010) also developed a
particle-based model for ice-phase clouds, but that model
relies on ice categories and mass–dimension relationships. Brdar and Seifert (2018)
developed McSnow, the first particle-based model for mixed-phase
clouds. McSnow is a multi-dimensional expansion of
the P3 bulk model (Morrison and Milbrandt, 2015; Milbrandt and Morrison, 2016)
and thus free from ice
categories; however, it still relies on mass–dimension relationships. Further,
a kinetic approach is applied to ice particles but not to droplets
or aerosol particles.

In this study, we demonstrate that a large-eddy simulation of a cumulonimbus that predicts ice particle morphologies without assuming ice categories or mass–dimension relationships is possible if we use SDM.

The organization of the remainder of this paper is as follows. In Sects. 2–4, our mixed-phase cloud mathematical model is described in detail. The cloud microphysics model is based on kinetic description and is coupled with moist air fluid dynamics. Note that this model is an expansion of Shima et al. (2009)'s warm cloud model. In Sect. 5, we develop a numerical model called SCALE-SDM by applying SDM. To evaluate SCALE-SDM's performance, we conduct a 2-D simulation of an isolated cumulonimbus. Section 6 presents the design of the numerical experiments, and in Sect. 7, the overall properties of the simulated cumulonimbus and ice particle morphologies are analyzed. The numerical convergence characteristics of the model are investigated in Sect. 8. In Sect. 9, possible improvements of the model are discussed, and a summary and conclusions are presented in Sect. 10. Lastly, lists of symbols and abbreviations are provided in Appendices A and B, respectively. Note that a comprehensive table of contents is provided as PDF bookmarks.

## 2.1 Notion of a particle

Let us represent aerosol, cloud, and precipitation particles as point
particles. The particle state is then characterized by two
types of variables: position ** x** and attributes

**. Attributes consist of several variables representing the particle's internal state, and the attributes considered in this study are $\mathit{a}=\mathit{\{}r,\mathit{\{}{m}_{\mathit{\alpha}}^{\mathrm{sol}}\mathit{\}},\mathit{\{}{m}_{\mathit{\beta}}^{\mathrm{insol}}\mathit{\}},{T}^{\mathrm{fz}},a,c,{\mathit{\rho}}^{\mathrm{i}},{m}^{\mathrm{rime}},{n}^{\mathrm{mono}},\mathit{v}\mathit{\}}$, i.e., liquid water amount, masses of soluble substances, masses of insoluble substances, freezing temperature, equatorial radius, polar radius, apparent density, rime mass, number of monomers, and velocity.**

*a*In this study, for simplicity, partially frozen/melted particles are not considered.
We assume that each particle completely freezes or melts
instantaneously (see Sects. 4.1.4 and 4.1.5).
Therefore, either the equivalent droplet radius *r* or
ice particle attributes $\mathit{\{}a,c,{\mathit{\rho}}^{\mathrm{i}}\mathit{\}}$ are always zero in our model.
Furthermore, we assume that all particles contain
soluble substances and are always deliquescent even when the
humidity is low (see Sect. 4.1.6).
Further, as a crude representation of “pre-activation”,
we do not allow the complete sublimation of an ice particle (see Sect. 4.1.7).
Therefore, *r* and $\mathit{\{}a,c,{\mathit{\rho}}^{\mathrm{i}}\mathit{\}}$ cannot be simultaneously zero.

In the remainder of this section, we provide a detailed explanation of each attribute.

## 2.2 Liquid water amount

The amount of liquid water contained in a particle is expressed by the
volume-equivalent sphere's radius *r*. That is, the volume of water in
a particle is (4∕3)*π**r*^{3}.

## 2.3 Masses of soluble and insoluble substances

Let ${m}_{\mathit{\alpha}}^{\mathrm{sol}}$, $\mathit{\alpha}=\mathrm{1},\mathrm{2},\mathrm{\dots},{N}^{\mathrm{sol}}$ be the masses of soluble substances contained in the particle, and let ${m}_{\mathit{\beta}}^{\mathrm{insol}}$, $\mathit{\beta}=\mathrm{1},\mathrm{2},\mathrm{\dots},{N}^{\mathrm{insol}}$ be the masses of insoluble substances.

## 2.4 Freezing temperature and ice nucleation active surface site

We only consider homogeneous freezing and condensation/immersion freezing in this study because these are dominant in mixed-phase clouds (e.g., Cui et al., 2006; De Boer et al., 2011; Murray et al., 2012).

Based on the “singular hypothesis” (Levine, 1950), we consider that each insoluble particle
has its own freezing temperature *T*^{fz}, and that a supercooled droplet freezes as soon as
the ambient temperature *T* decreases below *T*^{fz}.
The freezing process is described in detail in Sect. 4.1.4.

Each particle's *T*^{fz} is directly connected to the ice
nucleation active surface site (INAS) density concept
(e.g., Fletcher, 1969; Connolly et al., 2009; Niemand et al., 2012; Hoose and Möhler, 2012).

An INAS is a localized structure, such as lattice mismatches, cracks, and hydrophilic sites, on an insoluble substance's surface
that catalyzes ice formation at temperatures lower than a specific temperature. INAS
density *n*_{S}(*T*) gives the accumulated number of INAS per
unit surface area of the insoluble substance. Therefore,
*n*_{S}(*T*) is a function that increases as *T* decreases. The freezing temperature *T*^{fz} corresponds to the
highest temperature at which the first INAS appears on the insoluble
substance's surface. Let *A*^{insol} be the insoluble substance's surface area.
Then, the probability that *T*^{fz} is larger than *T* can be
calculated as $P({T}^{\mathrm{fz}}>T)=\mathrm{1}-\mathrm{exp}[-{A}^{\mathrm{insol}}{n}_{\mathrm{S}}(T\left)\right].$
The probability density function of *T*^{fz} then becomes

We can determine *T*^{fz} by selecting a random number that follows this probability distribution.

For mineral dust, biogenic substances, and soot, we can use the INAS density formulas of Niemand et al. (2012), Wex et al. (2015), and Ullrich et al. (2017), respectively.
If a particle consists of multiple insoluble substances,
we assume that *T*^{fz} is the highest of all.

It is possible that a single INAS does not appear until
−38 ^{∘}C, meaning that the particle is ice nucleation (IN) inactive and
will not freeze by immersion/condensation freezing but only by homogeneous
freezing. To account for this, we set
${T}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C. If a particle contains
only soluble substances, we also set
${T}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C.

There are various ice nucleation pathways (e.g., Kanji et al., 2017); however, in this study, we do not consider other ice nucleation pathways, such as deposition nucleation, deliquescent freezing, pore freezing, and contact freezing. The possibility of extending our model to incorporate these mechanisms is discussed in Sect. 9.3.1.

## 2.5 Porous spheroid approximation of ice particles

Ice particles have diverse morphologies such as columns, hexagonal
plates, dendrites, rimed crystals, graupel, hailstones, and
aggregates (e.g., Magono and Lee, 1966; Kikuchi et al., 2013).
Following the strategies of Chen and Lamb (1994a, b),
Misumi et al. (2010), and Jensen and Harrington (2015), let us approximate each
ice particle as a porous spheroid, which is characterized
by three variables, namely equatorial radius *a*, polar radius *c*, and
apparent density *ρ*^{i}. That is, the ice particle's apparent volume is $V=(\mathrm{4}\mathit{\pi}/\mathrm{3}){a}^{\mathrm{2}}c$, and its mass can be evaluated as
*m*=*ρ*^{i}*V*. The two radii *a* and *c* represent the ice particle's spatial extent
and *ρ*^{i} represents its
internal structure. Let us define the aspect ratio as $\mathit{\varphi}:=c/a$.
A spheroid is considered a prolate spheroid if *ϕ*>1, and
columns could be approximated by prolate
spheroids. In contrast, plates and
dendrites are approximated by oblate spheroids, i.e., *ϕ*<1. If
an ice particle is hollowed out or intricately branched, *ρ*^{i} becomes
smaller than the ice crystal's true density
${\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}\approx \mathrm{916.8}$ kg m^{−3}.

## 2.6 Rime mass and number of monomers

Following Brdar and Seifert (2018) we introduce two additional
ice particle attributes, namely rime mass *m*^{rime} and number
of monomers *n*^{mono}. Rime mass *m*^{rime}
records the mass of ice a particle has obtained through the
riming process. The number of monomers *n*^{mono} is an integer
representing the number of primary ice crystals in the particle. In
this study, *m*^{rime} and *n*^{mono} are
used only for analyzing the simulation results. Unlike the McSnow model of
Brdar and Seifert (2018), this study's time evolution equations
do not depend on *m*^{rime} or *n*^{mono}, as
will be detailed in Sect. 4.1.

## 2.7 Velocity

We approximate that each particle is always moving at its terminal velocity.
Therefore, a particle's velocity ** v** is a diagnostic attribute.

## 2.8 Effective number of attributes

In summary, particle attributes consist of
$\mathit{a}=\mathit{\{}r,\mathit{\{}{m}_{\mathit{\alpha}}^{\mathrm{sol}}\mathit{\}},\mathit{\{}{m}_{\mathit{\beta}}^{\mathrm{insol}}\mathit{\}},{T}^{\mathrm{fz}},a,c,{\mathit{\rho}}^{\mathrm{i}},{m}^{\mathrm{rime}},{n}^{\mathrm{mono}},\mathit{v}\mathit{\}}$. We
need the mass of insoluble substances
$\mathit{\{}{m}_{\mathit{\beta}}^{\mathrm{insol}},\mathit{\beta}=\mathrm{1},\mathrm{2},\mathrm{\dots},{N}^{\mathrm{insol}}\mathit{\}}$
(and corresponding INAS densities) to specify freezing temperature
*T*^{fz}. However, as described in Sect. 4.1,
time evolution equations do not depend on
$\mathit{\left\{}{m}_{\mathit{\beta}}^{\mathrm{insol}}\mathit{\right\}}$. Rime mass *m*^{rime} and
the number of monomers *n*^{mono} do not affect time
evolution either. Particle velocity ** v** is a diagnostic
attribute. Therefore, the attributes directly relevant to time
evolution are reduced to
$\mathit{\{}r,\mathit{\{}{m}_{\mathit{\alpha}}^{\mathrm{sol}}\mathit{\}},{T}^{\mathrm{fz}},a,c,{\mathit{\rho}}^{\mathrm{i}}\mathit{\}}$.
Compared to the warm cloud SDM model of Shima et al. (2009),
we have introduced four new attributes.

We only consider dry air and water vapor for the gas phase and ignore other trace gases.
In this section, we introduce several variables that describe
the state of moist air:
wind velocity $\mathit{U}=(U,V,W)$,
density of dry air *ρ*_{d},
density of water vapor *ρ*_{v},
density of moist air $\mathit{\rho}:={\mathit{\rho}}_{\mathrm{d}}+{\mathit{\rho}}_{\mathrm{v}}$,
specific humidity ${q}_{\mathrm{v}}:={\mathit{\rho}}_{\mathrm{v}}/\mathit{\rho}$,
mass of dry air per unit mass of moist air ${q}_{\mathrm{d}}:={\mathit{\rho}}_{\mathrm{d}}/\mathit{\rho}$,
temperature *T*,
pressure *P*,
and
potential temperature of moist air $\mathit{\theta}:=T/\mathrm{\Pi}:=T/(P/{P}_{\mathrm{0}}{)}^{R/{c}_{\mathrm{p}}}$.
Here, *P*_{0}=1000 hPa is a reference pressure;
*R*_{d}, *R*_{v}, and $R:={q}_{\mathrm{d}}{R}_{\mathrm{d}}+{q}_{\mathrm{v}}{R}_{\mathrm{v}}$
are the gas constants of dry air, water vapor, and moist air, respectively; and
*c*_{pd}, *c*_{pv}, and ${c}_{\mathrm{p}}:={q}_{\mathrm{d}}{c}_{\mathrm{pd}}+{q}_{\mathrm{v}}{c}_{\mathrm{pv}}$
are the isobaric specific heats of dry air, water vapor, and moist air, respectively.
To simplify notation, we introduce a variable representing the state of moist air:
$\mathit{G}:=\mathit{\{}\mathit{U},\mathit{\rho},{q}_{\mathrm{v}},\mathit{\theta},P,T\mathit{\}}$.

In this section, we describe our model's time evolution equations, first from cloud microphysics and then moist air fluid dynamics. Our model is detailed; however, it still falls short in completely describing mixed-phase cloud microphysics. To keep the model description concise, discussions on the shortcomings and how to overcome them are left for Sect. 9.

## 4.1 Cloud microphysics

Let us assign a unique index *i* to each particle.
This section explains the time evolution equations of particles
$\mathit{\left\{}\mathit{\right\{}{\mathit{x}}_{i}\left(t\right),{\mathit{a}}_{i}\left(t\right)\mathit{\}},i=\mathrm{1},\mathrm{2},\mathrm{\dots},{N}_{\mathrm{r}}^{\mathrm{wp}}\mathit{\}}$.
Here, ${N}_{\mathrm{r}}^{\mathrm{wp}}$ represents the total number of particles accumulated over the whole period.
However, because of coalescence, precipitation, and other processes, some particles might not exist all the time; thus, we
let *I*_{r}(*t*) be the set of particle indices existing in the domain at time *t*.

### 4.1.1 Advection and sedimentation

Particle *i*'s motion equation is

where *m*_{i} is the particle's mass,
${\mathit{F}}_{i}^{\mathrm{drg}}$ is the force of drag from moist air,
*g* is Earth's gravity, and
$\widehat{\mathit{z}}$ is the unit vector in the *z*-axis direction.
Note that $-{\mathit{F}}_{i}^{\mathrm{drg}}$ gives the reaction force acting on moist air.
The momentum of moist air changes as described in Eqs. (73) and (81).

If terminal velocity is reached, the motion equation becomes

where *U*_{i}:=** U**(

*x*_{i}) is the

*i*th particle's ambient wind velocity, and ${v}_{i}^{\mathrm{\infty}}$ is the terminal velocity, which is a function of attributes

*a*_{i}and the state of the ambient air

*G*_{i}.

In this study, we assume that terminal velocity is always achieved instantaneously;
however, this is a simplification. For example, the relaxation time of large droplets is a few seconds
(Fig. 3 of Wang and Pruppacher, 1977) though that of micrometer-sized droplets is approximately 10^{−5} s
(see, e.g., Eq. 1 of Chen et al., 2018, and the discussion that follows).
The acceleration of particles can be considered by explicitly solving the motion equation
(see, e.g., Naumann and Seifert, 2015), but extremely small time steps would be required for small particles.

The next two subsections explain the formulas used to calculate droplet and ice particle terminal velocities.

### 4.1.2 Droplet terminal velocity

To calculate droplet terminal velocity, we use the formula of
Beard (1976):
${v}_{i}^{\mathrm{\infty}}={v}_{\mathrm{Beard}}^{\mathrm{\infty}}(min({r}_{i},\mathrm{3.5}\phantom{\rule{0.125em}{0ex}}\mathrm{mm});{\mathit{\rho}}_{i},{P}_{i},{T}_{i})$,
where *ρ*_{i}:=*ρ*(*x*_{i}) and *P*_{i}:=*P*(*x*_{i}) are the density and pressure of ambient moist air,
respectively. This formula
applies to droplets with radii smaller than 3.5 mm.
If we use the formula for droplets larger than this,
the fall speed becomes unrealistically fast.
Therefore, we use the fall speed of
a droplet with a 3.5 mm radius for droplets larger than the size limit.

### 4.1.3 Ice particle terminal velocity

For ice particle terminal velocity, we use the formula of Böhm (1989, 1992c, 1999):
${v}_{i}^{\mathrm{\infty}}={v}_{\text{B\xf6hm}}^{\mathrm{\infty}}({m}_{i},{\mathit{\varphi}}_{i},{d}_{i},{q}_{i};{\mathit{\rho}}_{i},{T}_{i})$,
where *d*_{i} is the characteristic length, and *q*_{i} is the area ratio.

In Böhm's theory, *d*_{i} is defined by 2*a*_{i},
and *q*_{i} is defined by the area ratio regarding
circumscribed ellipse ${q}_{i}^{\mathrm{ce}}:={A}_{i}/{A}_{i}^{\mathrm{ce}}$,
where
*A*_{i} is the projected area perpendicular to the flow direction,
and ${A}_{i}^{\mathrm{ce}}$ is the area of the circumscribed ellipse of *A*_{i},
i.e., the area of the smallest ellipse that completely contains *A*_{i}.

However, in this study, we start from a slightly different definition of *d*_{i} and *q*_{i},
which we adopted mistakenly:

where
*D*_{i} is the maximum dimension,
${q}_{i}^{\mathrm{cc}}$ is the area ratio regarding circumcircle,
and ${A}_{i}^{\mathrm{cc}}$ is the area of the circumcircle of *A*_{i},
i.e., the area of the smallest circle that completely contains *A*_{i}.

Consequently, Eq. (4) underestimates the fall speeds of columnar ice particles. Nevertheless, based on the assessment detailed in Sect. 9.2, we will confirm that this difference does not change the results of our simulation significantly, and hence we conclude that this flaw causes only a minor impact on this study. We also note that in Sect. 9.2 we will develop and release a fixed version of the model, SCALE-SDM 0.2.5-2.2.2.

In our model,
we assume that ice particles are falling with their maximum dimension perpendicular to the flow direction. Therefore,
the circumcircle area becomes ${A}_{i}^{\mathrm{cc}}=\mathit{\pi}max({a}_{i},{c}_{i}{)}^{\mathrm{2}}$.
The projected area *A*_{i} can be roughly evaluated by the area of the circumscribed ellipse ${A}_{i}^{\mathrm{ce}}=\mathit{\pi}{a}_{i}max({a}_{i},{c}_{i})$; however,
we must subtract pores and indentations at boundaries from ${A}_{i}^{\mathrm{ce}}$.
We assume that the ratio ${A}_{i}/{A}_{i}^{\mathrm{ce}}$ is a power of the volume fraction ${\mathit{\rho}}_{i}^{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$, and that the exponent *κ* is a function of the aspect ratio *ϕ*_{i}:

Based on the following arguments, we propose a value *κ* of the form

Following Jensen and Harrington (2015),
we assume *κ*→1 as *ϕ*_{i}→0, and
*κ*→0 as *ϕ*_{i}→∞. *ϕ*_{i}≪1 means that
the ice particle is thin and extends horizontally.
Therefore, we can expect that the structure is
uniform along the vertical axis and that the ratio ${A}_{i}/{A}_{i}^{\mathrm{ce}}$ is equal to
the volume fraction ${\mathit{\rho}}_{i}^{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$.
Thus, $\mathit{\kappa}({\mathit{\varphi}}_{i}=\mathrm{0})=\mathrm{1}$. At the other extreme, *ϕ*_{i}≫1 indicates
that the ice particle is columnar.
Such ice crystals typically hollow inward along
their basal face; therefore, the volume fraction ${\mathit{\rho}}_{i}^{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$ will not affect the ratio ${A}_{i}/{A}_{i}^{\mathrm{ce}}$.
Thus, $\mathit{\kappa}({\mathit{\varphi}}_{i}\to \mathrm{\infty})=\mathrm{0}$.

For *ϕ*_{i}≈1, Jensen and Harrington (2015) argued that $({\mathit{\rho}}_{i}^{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}{)}^{\mathit{\kappa}}=\mathrm{1}$, i.e., *κ*=0.
However, this cannot be justified for aggregates with low apparent densities.
Thus, we estimate *κ* through a dimensional analysis.
We assume that the power laws ${m}_{i}\propto {D}_{i}^{\mathit{\beta}}$ and ${A}_{i}\propto {D}_{i}^{\mathit{\beta}/s}$ hold.
Thus, by the definition of apparent density, ${\mathit{\rho}}_{i}^{\mathrm{i}}={m}_{i}/\left(\right(\mathrm{4}/\mathrm{3}\left)\mathit{\pi}{a}_{i}^{\mathrm{2}}{c}_{i}\right)\propto {D}_{i}^{\mathit{\beta}-\mathrm{3}}$.
From Eq. (5), ${D}_{i}^{\mathit{\beta}/s}={D}_{i}^{\mathrm{2}}{D}_{i}^{(\mathit{\beta}-\mathrm{3})\mathit{\kappa}}$.
Hence, $\mathit{\kappa}=(\mathrm{2}s-\mathit{\beta})/\mathit{\left\{}s\right(\mathrm{3}-\mathit{\beta}\left)\mathit{\right\}}$ holds. Schmitt and Heymsfield (2010)
estimated that $(\mathit{\beta},s)=(\mathrm{2.22},\mathrm{1.30})$ for aggregates observed during
the Cirrus Regional Study of
Tropical Anvils and Cirrus Layers – Florida Area Cirrus
Experiment (CRYSTAL-FACE) field project. Therefore, *κ*=0.375 for CRYSTAL-FACE aggregates.
They also estimated that $(\mathit{\beta},s)=(\mathrm{2.20},\mathrm{1.25})$ for aggregates observed during
an Atmospheric Radiation Measurement (ARM) field project, which results in *κ*=0.300.

The *κ* given by Eq. (6) yields
*κ*(0)=1, *κ*(1)=0.368, and *κ*(∞)=0,
which agree with the aforementioned estimation.

### 4.1.4 Immersion/condensation and homogeneous freezing

As explained in Sect. 2.4, a supercooled droplet freezes when the ambient temperature drops below its freezing temperature. This section provides a more precise description of when and how freezing occurs in our model.

We consider that the *i*th particle freezes immediately when the
following three conditions are all satisfied: (1) the particle is a
droplet, i.e., *r*_{i}>0; (2) the ambient water vapor is supersaturated over liquid
water, i.e., ${e}_{i}>{e}_{\mathrm{s}}^{\mathrm{w}}\left({T}_{i}\right)$; and (3) the ambient
temperature is lower than the particle's freezing temperature, i.e.,
${T}_{i}<{T}_{i}^{\mathrm{fz}}$. Here, *e*_{i}:=*e*(*x*_{i}) and
*T*_{i}:=*T*(*x*_{i}) are the ambient vapor pressure and temperature of
the *i*th particle, respectively, and ${e}_{\mathrm{s}}^{\mathrm{w}}\left(T\right)$ is
the saturation vapor pressure over a planar liquid water surface at
temperature *T*.

We assume that the resulting ice crystal is spherical, with the true ice crystal density ${\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$. Therefore,
attributes are initiated as follows: ${r}_{i}^{\prime}=\mathrm{0}$,
${a}_{i}^{\prime}={c}_{i}^{\prime}={r}_{i}({\mathit{\rho}}^{\mathrm{w}}/{\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}{)}^{\mathrm{1}/\mathrm{3}}$,
${\mathit{\rho}}_{i}^{\mathrm{i}\prime}={\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$,
${n}_{i}^{\mathrm{mono}\prime}=\mathrm{1}$, and ${m}_{i}^{\mathrm{rime}\prime}=\mathrm{0}$. The
primed variables here denote values after the update, and *ρ*^{w}
is the density of liquid water. $\mathit{\left\{}{m}_{\mathit{\alpha}i}^{\mathrm{sol}}\mathit{\right\}}$,
$\mathit{\left\{}{m}_{\mathit{\beta}i}^{\mathrm{insol}}\mathit{\right\}}$, and ${T}_{i}^{\mathrm{fz}}$ remain unchanged.

When freezing occurs, each particle releases latent heat of fusion to the moist air, as described in Eqs. (74), (79), and (80).

### 4.1.5 Melting

When ambient temperature rises above
0 ^{∘}C, we consider that melting occurs
immediately. Thus, the attributes are updated as follows:
${r}_{i}^{\prime}=({a}_{i}^{\mathrm{2}}{c}_{i}{\mathit{\rho}}_{i}^{\mathrm{i}}/{\mathit{\rho}}^{\mathrm{w}}{)}^{\mathrm{1}/\mathrm{3}}$ and
${a}_{i}^{\prime}={c}_{i}^{\prime}={\mathit{\rho}}_{i}^{\mathrm{i}\prime}={n}_{i}^{\mathrm{mono}\prime}={m}_{i}^{\mathrm{rime}\prime}=\mathrm{0}$.
$\mathit{\left\{}{m}_{\mathit{\alpha}i}^{\mathrm{sol}}\mathit{\right\}}$, $\mathit{\left\{}{m}_{\mathit{\beta}i}^{\mathrm{insol}}\mathit{\right\}}$,
and ${T}_{i}^{\mathrm{fz}}$ remain unchanged. When melting occurs, each
particle absorbs latent heat of fusion from the moist air,
as indicated in Eqs. (74), (79), and (80).

### 4.1.6 Condensation and evaporation

Following, e.g., Rogers and Yau (1989), the time evolution equation describing droplet growth by condensation/evaporation can be derived as follows.

The growth rate is identical to vapor flux at the droplet surface. If the diffusion of vapor around the droplet is in a quasi-steady state, we obtain

Here, *D*_{v} is water vapor's diffusivity in air,
*ρ*_{vi}:=*ρ*_{v}(*x*_{i}) is
the ambient moist air's water vapor density, and ${\mathit{\rho}}_{\mathrm{v}i}^{\mathrm{sfc}}$ is
water vapor density at the surface of the droplet.

If we further assume that thermal diffusion is also in a quasi-steady state,
and that surface temperature ${T}_{i}^{\mathrm{sfc}}$ and ambient temperature *T*_{i} are close to each other,
i.e., $({T}_{i}^{\mathrm{sfc}}-{T}_{i})/{T}_{i}\ll \mathrm{1}$,
Eq. (7) can be reduced to

where ${S}_{i}^{\mathrm{w}}:={e}_{i}/{e}_{\mathrm{s}}^{\mathrm{w}}\left({T}_{i}\right)$ is the ambient saturation ratio over liquid water, and

where *L*_{v} is the latent heat of vaporization, *k* is the thermal conductivity of moist air,
and ${e}_{\mathrm{s}i}^{\mathrm{w},\mathrm{eff}}$ is the effective saturation vapor pressure regarding the *i*th droplet's surface.
Following Köhler's theory (Köhler, 1936),
an approximate formula of ${e}_{\mathrm{s}i}^{\mathrm{w},\mathrm{eff}}$ can be derived as

where $a\approx \mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{5}}\phantom{\rule{0.125em}{0ex}}\mathrm{cm}\phantom{\rule{0.125em}{0ex}}\mathrm{K}/{T}_{i}$,
$b\approx \mathrm{4.3}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{\mathrm{3}}{\sum}_{\mathit{\alpha}}{I}_{\mathit{\alpha}}{m}_{\mathit{\alpha}i}^{\mathrm{sol}}/{M}_{\mathit{\alpha}}^{\mathrm{sol}}$,
*I*_{α} is the van 't Hoff factor, which represents the degree of ionic dissociation, and
${M}_{\mathit{\alpha}}^{\mathrm{sol}}$ is the molecular weight of the solute *α*.
The second and third terms of Eq. (10) account for curvature and solute effects, respectively.

The growth of a droplet by condensation/evaporation is governed by Eqs. (8)–(10) in our model. When a droplet or an ice particle falls through the air, the flow around it enhances the diffusional growth, a phenomenon known as the ventilation effect. It does not essentially affect the growth of droplets smaller than 50 µm in radius (see Sect. 13.2.3 of Pruppacher and Klett, 1997). Therefore, for simplicity, we do not consider the ventilation effect on droplets in this study. Notably, Eqs. (8)–(10) also describe the respective activation and deactivation of cloud droplets from and to aerosol particles (see, e.g., Arabas and Shima, 2017; Hoffmann, 2017; Abade et al., 2018).

Vapor and latent heat couplings to moist air through condensation and evaporation are calculated by Eqs. (71), (72), (74), (76), (77), and (79).

### 4.1.7 Deposition and sublimation

The shapes of ice crystals formed by depositional growth exhibit strong dependencies on temperature and, to a lesser extent, supersaturation (e.g., Nakaya, 1954; Hallett and Mason, 1958; Kobayashi, 1961). The former is known as the primary growth habit and the latter as the secondary growth habit. The primary growth habit determines the preferred growth direction, i.e., columnar or planar, and the secondary growth habit determines the mode of growth, i.e., whether the columnar crystal becomes solid or hollow, and whether the planar crystal becomes plate-like, sectored, or dendritic. In this study, we use the model of Chen and Lamb (1994a) with various modifications.

The mass growth rate can be derived similarly to Eqs. (7) and (8):

where ${S}_{i}^{\mathrm{i}}:={e}_{i}/{e}_{\mathrm{s}}^{\mathrm{i}}\left({T}_{i}\right)$ is the ambient saturation ratio over ice, and
${e}_{\mathrm{s}}^{\mathrm{i}}\left(T\right)$ is the saturation vapor pressure over ice at temperature *T*,

where *L*_{s} is the latent heat of sublimation,
$C=C({a}_{i},{c}_{i})$ is the electric capacitance of the spheroid,
and $\stackrel{\mathrm{\u203e}}{{f}_{\mathrm{vnt}}}$ is the particle-averaged ventilation coefficient.

The exact form of capacitance *C*(*a*_{i},*c*_{i}) is given by Chen and Lamb (1994a).
$C\approx (\mathrm{2}{a}_{i}+{c}_{i})/\mathrm{3}$ gives a good approximation for *ϕ*_{i}≈1.

The coefficient $\stackrel{\mathrm{\u203e}}{{f}_{\mathrm{vnt}}}$ accounts for the ventilation effect, i.e., the enhancement of diffusional growth by air flow. Hall and Pruppacher (1976) suggested that $\stackrel{\mathrm{\u203e}}{{f}_{\mathrm{vnt}}}$ could be described by

where
$({b}_{\mathrm{1}},{b}_{\mathrm{2}},\mathit{\gamma})=(\mathrm{1.0},\mathrm{0.14},\mathrm{2})$ for *X*≤1,
$({b}_{\mathrm{1}},{b}_{\mathrm{2}},\mathit{\gamma})=(\mathrm{0.86},\mathrm{0.28},\mathrm{1})$ for *X*>1,
$X={N}_{\mathrm{Sc}}^{\mathrm{1}/\mathrm{3}}({N}_{\mathrm{Re}i}^{\mathrm{i}}{)}^{\mathrm{1}/\mathrm{2}}$,
${N}_{\mathrm{Sc}}=\mathit{\mu}/\left(\mathit{\rho}{D}_{\mathrm{v}}\right)$ is the Schmidt number,
${N}_{\mathrm{Re}i}^{\mathrm{i}}=\mathit{\rho}{v}_{i}^{\mathrm{\infty}}{D}_{i}/\mathit{\mu}$ is the Reynolds number of ice particle *i*,
and *μ* is the dynamic viscosity of moist air.

Note that *m*_{i} in Eq. (11) can become zero through sublimation over a finite time.
However, in this study, we prohibit complete sublimation, and instead, we impose a limiter to d*m*_{i} as follows:

where ${m}_{\mathrm{min}}^{\mathrm{i}}$ is an arbitrary small mass taken from the mass of a spherical ice particle with a radius of 1 nm
and the true ice density ${\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$.
This is a crude representation of pre-activation (see, e.g., Marcolli, 2017, for a review).
Each particle keeps the memory of ice activation until the ambient temperature rises above 0 ^{∘}C.
A particle with ${m}_{\mathrm{min}}^{\mathrm{i}}$ ice grows immediately after the ambient air is supersaturated over ice,
irrespective of its freezing temperature ${T}_{i}^{\mathrm{fz}}$.

In Chen and Lamb's (1994a) model, the primary growth habit is expressed by
an empirical function known as the inherent growth ratio Γ(*T*),
which modulates the *c*-axis to *a*-axis growth rate ratio:

where *f*_{vnt} is the primary growth habit's ventilation coefficient,
and Γ^{*} is the effective inherent growth ratio, including the ventilation effect.

For purely diffusional growth, $\text{d}{c}_{i}/\text{d}{a}_{i}={c}_{i}/{a}_{i}$ holds;
therefore, the aspect ratio does not change, i.e., *d**ϕ*_{i}=0.
Γ(*T*) represents the lateral
redistribution of vapor on the ice crystal surface through kinetic processes.
We use the Γ(*T*) proposed by Chen and Lamb (1994a)
but set Γ(*T*)=1 for *D*<10 µm,
as observations suggest that ice crystals are quasi-spherical if *D*<60 µm
(Baran, 2012; Korolev and Isaac, 2003; Lawson et al., 2008).
Additionally, the Γ(*T*) provided in Chen and Lamb (1994a) is for temperatures between −30
and 0 ^{∘}C.
For lower temperatures, we simply assume

The ventilation coefficient *f*_{vnt} represents the preferential enhancement of vapor flux
toward the ice crystal's major axis because of the air flow around it.
Chen and Lamb (1994a) derived a *f*_{vnt} of the form

The secondary growth habit is expressed by deposition density *ρ*_{dep},
which represents the apparent density of the ice fraction newly created by deposition.
Then, the change in ice particle volume d*V*_{i} is given by

Deposition density *ρ*_{dep} can be expressed as

Here, following Jensen and Harrington (2015),
we assume that planar crystal branching does not occur if the equatorial radius *a*_{i} is smaller than 100 µm.
${\mathit{\rho}}_{\mathrm{dep}}^{\mathrm{CL}\mathrm{94}}$ is an empirical formula of deposition density proposed by Chen and Lamb (1994a),

where $\mathrm{\Delta}{\mathit{\rho}}_{i}:={\mathit{\rho}}_{\mathrm{v}i}-{\mathit{\rho}}_{\mathrm{v}i}^{\mathrm{sfc}}$.
From Eq. (11), Δ*ρ*_{i} becomes

Here, following Miller and Young (1979), we limit *ρ*_{vi} by water saturation
and replace the Δ*ρ*_{i} in Eq. (20) with

For sublimation, the particle volume change d*V*_{i} is given by

where sublimation density *ρ*_{sbl} represents the apparent
density of the ice fraction removed by sublimation. For simplicity,
we assume that the ice particle's apparent density will not be
changed through sublimation, i.e.,

We can now calculate the attributes at time *t*+d*t*.
The apparent density becomes

where d*m*_{i} is given in Eqs. (11) and (14),
and d*V*_{i} is given in Eqs. (18) and (23).

From Eq. (15) and the definition of volume ${V}_{i}=(\mathrm{4}\mathit{\pi}/\mathrm{3}){a}_{i}^{\mathrm{2}}{c}_{i}$,
after d*t*, the two radii become

Applying those equations to a small ice particle's sublimation creates
an extremely small planar or columnar ice particle.
However,
observations suggest that ice crystals are quasi-spherical if *D*<60 µm
(Baran, 2012; Korolev and Isaac, 2003; Lawson et al., 2008).
Therefore, we regard the ice particle as spherical with the true ice density
if the minor axis predicted by Eqs. (26) and (27) is smaller than 1 µm.
That is, if $min\mathit{\left\{}{a}_{i}\right(t+\text{d}t),{c}_{i}(t+\text{d}t\left)\mathit{\right\}}<\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m}$,

where primed variables indicate values after correction.

For simplicity, we assume that the rime mass fraction ${m}_{i}^{\mathrm{rime}}/{m}_{i}$ does not change through sublimation, following Brdar and Seifert (2018):

Vapor and latent heat couplings to moist air through deposition and sublimation are calculated by Eqs. (71), (72), (74), (76), (78), and (79).

In this section, we detailed the deposition and sublimation model used in
SCALE-SDM; however, there is significant room for improvement.
For example, as we will discuss in Sect. 9.1.4,
using Γ(*T*) for sublimation is questionable.
Instead, we propose using Γ(*T*)=1 for sublimation (Eq. 110),
and validate this correction in Sect. 9.1.5.
Furthermore, in Sect. 9.2, to prohibit the creation of unnaturally slender ice particles,
we will propose to impose a limiter to the
effective inherent growth ratio Γ^{*} (Eq. 123).
Several other issues of our deposition/sublimation model,
such as the representation of polycrystals, will be
discussed in Sect. 9.3.5.

### 4.1.8 Stochastic description of coalescence, riming, and aggregation

Particle coalescence, riming, and aggregation can be considered a stochastic process. Following Gillespie (1972),
consider a region with volume Δ*V*. If Δ*V* is
sufficiently small, we can consider that particles within this
region are well mixed, e.g., by atmospheric turbulence
(see, e.g., Shima et al., 2009; Dziekan and Pawlowska, 2017).
Then, all particle pairs in the volume can collide and
coalesce/rime/aggregate during an infinitesimal time interval d*t*.
The probability that a particle pair *j* and *k* inside Δ*V* will collide and coalesce/rime/aggregate within an infinitesimal time interval
$(t,t+\text{d}t)$ is given by

where the function $K({\mathit{a}}_{j},{\mathit{a}}_{k};\mathit{G})$ is called the
collision–coalescence/–riming/–aggregation kernel, and
** G** denotes the state of the moist air in Δ

*V*.

In this study, we consider coalescence, riming, and aggregation induced by differential gravitational settling of particles because this mechanism is dominant in mixed-phase clouds.

### 4.1.9 Coalescence between two droplets

First, we consider droplet coalescence, which accounts for the formation of rain droplets from cloud droplets (autoconversion), the collection of cloud droplets by rain droplets (accretion), and the coalescence of two rain droplets (self-collection).

The collision–coalescence kernel is given by

where *E*_{coal}(*r*_{j},*r*_{k}) is the collection efficiency of collision–coalescence,
which can be decomposed into ${E}_{\mathrm{coal}}={E}_{\mathrm{coal}}^{\mathrm{collis}}{E}_{\mathrm{coal}}^{\mathrm{coal}}$.
Here, collision efficiency ${E}_{\mathrm{coal}}^{\mathrm{collis}}$
considers the effect that a smaller droplet is swept aside by
the flow around a larger droplet, or a droplet being caught in the wake of a similarly sized droplet collides on the downstream side.
We adopt the collision efficiency used in Seeßelberg et al. (1996) and Bott (1998).
Here, Davis (1972) and Jonas (1972) are used for small droplets,
and Hall (1980) for larger droplets, with modifications to the collector droplet radius range 70–300 µm
to incorporate the wake effect suggested by Lin and Lee (1975).
Not all the collisions end up with coalescence.
Rebound or breakup (fragmentation) could also occur.
Coalescence efficiency ${E}_{\mathrm{coal}}^{\mathrm{coal}}$
represents the fraction of collisions that result in permanent coalescence.
In this study, we assume ${E}_{\mathrm{coal}}^{\mathrm{coal}}=\mathrm{1}$ for simplicity.

If coalescence takes place, droplets *j* and *k* then merge into a single droplet. Thus, we keep *j* and remove *k* from the system.
The attributes of the new droplet *j* can be calculated as follows:

where primed values indicate the resultant droplet. Here, we assumed that the resultant particle's ${T}_{j}^{\mathrm{fz}\prime}$ is given by $max({T}_{j}^{\mathrm{fz}},{T}_{k}^{\mathrm{fz}})$, i.e., the higher freezing temperature of the two constituent particles. We also assume that the same applies to riming and aggregation.

Let us emphasize that the stochastic model introduced in this section describes the underlying mathematical model of the coalescence process, not the Monte Carlo algorithm of SDM that solves the stochastic process numerically. In the preceding paragraph, droplet *k* was removed from the system because both *j* and *k* are real particles.
On the contrary, in the SDM, the number of super-particles is (almost always) conserved through coalescence (Shima et al., 2009).

### 4.1.10 Riming between an ice particle and a droplet

Riming usually refers to the collection of small supercooled droplets by a larger ice particle, but we also include the collection of small ice particles by a larger droplet. The latter case could be regarded as a type of contact freezing. However, ice particles grow preferentially when ice particles and supercooled droplets coexist (Wegener–Bergeron–Findeisen mechanism). Therefore, we can expect that the latter case happens less frequently in mixed-phase clouds.

Hereafter we assume, without loss of generality, that particle *j* is an ice particle and particle *k* is a droplet. The collision–riming kernel is expressed as

where *E*_{rime} is the collision–riming collection efficiency
and *A*_{g} is the geometric cross-sectional area of *j* and *k*.

Figure 1 of Wang and Ji (2000) defines *A*_{g} for riming, but
calculating it rigorously for porous spheroid models is impossible.
Thus, we approximate *A*_{g} by

i.e., the indentation of the ice particle $({A}_{j}^{\mathrm{ce}}-{A}_{j})$ is subtracted from
the area of an ellipse with semi-axes (*a*_{j}+*r*_{k}) and $\mathit{\{}max({a}_{j},{c}_{j})+{r}_{k}\mathit{\}}$. Therefore, if ${r}_{k}\ll {a}_{j},{c}_{j}$, then *A*_{g}≈*A*_{j}. At the other extreme, if ${r}_{k}\gg {a}_{j},{c}_{j}$, then ${A}_{\mathrm{g}}\approx \mathit{\pi}({a}_{j}+{r}_{k})\mathit{\{}max({a}_{j},{c}_{j})+{r}_{k}\mathit{\}}$.

To evaluate collision–riming collection efficiency *E*_{rime}, we combine formulas proposed by Beard and Grover (1974) and Erfani and Mitchell (2017).

If ${v}_{j}^{\mathrm{\infty}}<{v}_{k}^{\mathrm{\infty}}$, we consider droplet *k* as the collector and adopt the formula of Beard and Grover (1974):

where
${p}^{\mathrm{i}/\mathrm{w}}:={r}_{j}^{\mathrm{i}}/{r}_{k}$, ${r}_{j}^{\mathrm{i}}:=({a}_{j}^{\mathrm{2}}{c}_{j}{)}^{\mathrm{1}/\mathrm{3}}$,
${N}_{\mathrm{Re}k}^{\mathrm{w}}=\mathit{\rho}{v}_{k}^{\mathrm{\infty}}\mathrm{2}{r}_{k}/\mathit{\mu}$ is the Reynolds number of droplet *k*, ${N}_{\mathrm{St}}^{\mathrm{i}/\mathrm{w}}=({p}^{\mathrm{i}/\mathrm{w}}{)}^{\mathrm{2}}{\mathit{\rho}}_{j}^{\mathrm{i}}{N}_{\mathrm{Re}k}^{\mathrm{w}}{C}_{\mathrm{SC}}/(\mathrm{9}\mathit{\rho})$ is the Stokes impaction parameter
when droplet *k* is collecting an ice particle, and
*C*_{SC} is the Cunningham slip correction factor.

If ${v}_{j}^{\mathrm{\infty}}\ge {v}_{k}^{\mathrm{\infty}}$, we consider ice particle *j* to be the collector. For spherical ice particle *ϕ*_{j}≈1, we again use the formula of Beard and Grover (1974) but replace the Stokes impaction parameter ${N}_{\mathrm{St}}^{\mathrm{i}/\mathrm{w}}$ with the mixed Froude number *N*_{mFr}
following Hall (1980), Rasmussen and Heymsfield (1985), and Heymsfield and Pflaum (1985).
For columnar and planar ice particles, we use formulas
${E}_{\mathrm{EM}\mathrm{17}}^{\mathrm{clm}}$ and ${E}_{\mathrm{EM}\mathrm{17}}^{\mathrm{pln}}$ from Erfani and Mitchell (2017),
which were obtained by fitting the numerical results of Wang and Ji (2000).
For the intermediate case, we calculate an average weighted by the aspect ratio *ϕ*_{j}.
For *ϕ*_{j}≤1 (planar),

For *ϕ*_{j}>1 (columnar),

Here, ${p}^{\mathrm{w}/\mathrm{i}}:=\mathrm{1}/{p}^{\mathrm{i}/\mathrm{w}}={r}_{k}/{r}_{j}^{\mathrm{i}}$,
${N}_{\mathrm{mFr}}=({v}_{j}^{\mathrm{\infty}}-{v}_{k}^{\mathrm{\infty}}){v}_{k}^{\mathrm{\infty}}/(g{D}_{j}/\mathrm{2})$,
and ${N}_{\mathrm{Re}j}^{\mathrm{clm}}=\mathit{\rho}{v}_{j}^{\mathrm{\infty}}\mathrm{2}{a}_{j}/\mathit{\mu}$ is the
Reynolds number based on the width of column 2*a*_{j}. Note that
there is a typo in Eq. (19) of Erfani and Mitchell (2017); i.e.,
the two case conditions are opposite.

If riming takes place, the ice particle *j* and droplet *k*
merge and instantaneously freeze into a single ice particle.
Thus, we keep *j* and remove *k* from the system.

If $max({a}_{j},{c}_{j})<{r}_{k}$, we assume that the resultant ice particle is spherical with the true ice density:

where primed values indicate the resultant ice particle.

If $max({a}_{j},{c}_{j})\ge {r}_{k}$,
we preserve the ice particle's maximum dimension, i.e., ${D}_{j}^{\prime}={D}_{j}$, until the ice particle becomes quasi-spherical.
This accounts for the gradual growth of an unrimed ice crystal
to a graupel particle with a quasi-spherical shape.
This filling-in simplification was introduced by Heymsfield (1982), and is used in various models
(e.g., Chen and Lamb, 1994b; Morrison and Grabowski, 2008, 2010; Jensen and Harrington, 2015; Morrison and Milbrandt, 2015). As graupels have an aspect
ratio of approximately 0.8 (Heymsfield, 1978),
we preserve the minor dimension if $\mathrm{0.8}<{\mathit{\varphi}}_{j}\le \mathrm{1}/\mathrm{0.8}=\mathrm{1.25}$, which mimics graupel's tumbling.
When an accreted droplet freezes, the air will be trapped inside.
Let rime density *ρ*_{rime} be the frozen droplet's apparent density.
Then, for *ϕ*_{j}≤0.8 (planar) and $\mathrm{1.0}<{\mathit{\varphi}}_{j}\le \mathrm{1.25}$ (columnar but quasi-spherical),

Other attributes are updated using Eqs. (44)–(48).
For *ϕ*_{j}>1.25 (columnar) and $\mathrm{0.8}<{\mathit{\varphi}}_{j}\le \mathrm{1.0}$ (planar but quasi-spherical),

Other attributes are updated using Eqs. (44)–(48).
Following Chen and Lamb (1994b), we use the formula of Heymsfield and Pflaum (1985)
to calculate rime density
*ρ*_{rime}:

where $Y:=(-{r}_{k}{v}_{\mathrm{imp}}/{T}_{j}^{\mathrm{sfc}})/\left(\mathrm{\mu}\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}/{}^{\circ}\mathrm{C}\right)$,
*v*_{imp} is impact velocity,
and ${T}_{j}^{\mathrm{sfc}}$ is the surface temperature of ice particle *j*.

where $A=\mathrm{0.30}\phantom{\rule{0.125em}{0ex}}\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$, *B*_{1}=0.44, ${B}_{\mathrm{2}}=-\mathrm{0.03115}$, ${B}_{\mathrm{3}}=-\mathrm{1.7030}$, *B*_{4}=0.9116, and ${B}_{\mathrm{5}}=-\mathrm{0.1224}$.

Impact velocity can be calculated using the formula of Rasmussen and Heymsfield (1985):
${v}_{\mathrm{imp}}=|{v}_{j}^{\mathrm{\infty}}-{v}_{k}^{\mathrm{\infty}}|max\mathit{\left\{}{f}_{\mathrm{RH}\mathrm{85}}\right({N}_{\mathrm{Re}j}^{\mathrm{i}},$ ${N}_{\mathrm{St}}^{\mathrm{w}/\mathrm{i}}),\mathrm{0}\mathit{\}}$,
where ${N}_{\mathrm{St}}^{\mathrm{w}/\mathrm{i}}=({p}^{\mathrm{w}/\mathrm{i}}{)}^{\mathrm{2}}{\mathit{\rho}}^{\mathrm{w}}{N}_{\mathrm{Re}j}^{\mathrm{i}}/(\mathrm{9}\mathit{\rho})$ is the Stokes impaction parameter
when an ice particle collects a droplet.
Because the *f*_{RH85} given in Rasmussen and Heymsfield (1985)
becomes slightly negative around $\mathrm{0.1}<{N}_{\mathrm{St}}^{\mathrm{w}/\mathrm{i}}<\mathrm{1.0}$, we impose a limiter to ensure it is positive.
Surface temperature ${T}_{j}^{\mathrm{sfc}}$ can be evaluated as

where Δ*ρ*_{j} is given in Eq. (21).
This equation is derived under an assumption of quasi-steady vapor and thermal diffusion.

When riming occurs, the frozen droplet releases the latent heat of fusion to the moist air as described in Eqs. (74), (79) and (80).

As we will discuss in Sect. 9.1.1, the rime density formula of Heymsfield and Pflaum (1985) must be revised slightly. We propose to replace the *Y* in Eq. (56) (not in Eq. 57)
with ${Y}^{\downarrow}=min(Y,\mathrm{3.5})$ (Eq. 107), because the rime density derived from Eq. (56) becomes too small for larger values of *Y*, which affects the shape of hailstones near the freezing level.

Another issue discussed in Sect. 9.1.2 is related to the filling-in model. Assuming that the diameter of the frozen droplet is preserved, if the diameter is larger than the ice particle's maximum dimension, we propose replacing Eq. (50) by Eq. (108) and Eq. (54) by Eq. (109).

We validate these two corrections in Sect. 9.1.5. More discussions to refine our riming model will be presented in Sect. 9.3.7.

### 4.1.11 Aggregation between two ice particles

Finally, we consider the aggregation of ice particles. Following Connolly et al. (2012), we use the projected area of particles to evaluate the geometric cross-sectional area. The collision–aggregation kernel is then given by

where *E*_{agg} is the collision–aggregation collection efficiency.
Following Morrison and Grabowski (2010),
we assume that the efficiency is given by a constant, *E*_{agg}=0.1, in this study. Field et al. (2006) confirmed that *E*_{agg}=0.09 produces a good agreement with aircraft observations.

If aggregation takes place, ice particles *j* and *k* merge into a single ice particle. Thus, we keep *j* and remove *k* from the system.
However, no reliable model exists for calculating the next porous spheroid.
Chen and Lamb (1994b) proposed a model, but it tends to create
snow aggregates with impossibly low apparent densities (lighter than vapor).
In this study, we propose another intuitive model by incorporating the
compaction of fluffy snowflakes to cope with the problem.

Snow aggregates have complicated fractal structures. However, if we circumscribe them using a spheroid, the growth by aggregation is in three dimensions, rather than one (columnar) or two (planar). Therefore, as in the case of riming, we assume that only the minor dimension grows by aggregation.

If the volume-weighted average density $\stackrel{\mathrm{\u203e}}{{\mathit{\rho}}_{jk}^{\mathrm{i}}}=({m}_{j}+{m}_{k})/({V}_{j}+{V}_{k})$ is closer to the true density of ice ${\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$, the two particles aggregate without changing their shapes. Hence, when we approximate the resultant aggregate with a spheroid, there are more empty spaces inside, thus reducing the apparent density. Let us denote the minimum possible apparent density as ${\mathit{\rho}}_{jk}^{\mathrm{i},\mathrm{min}}$, which can be evaluated using Eq. (61), which we will derive shortly.

In contrast, if $\stackrel{\mathrm{\u203e}}{{\mathit{\rho}}_{jk}^{\mathrm{i}}}$ is small, compaction of the fluffy snowflakes occurs, and the empty space of the larger ice particle could be filled with the smaller ice particle or the particles might deform because of the collision–aggregation impact. Because of this compaction mechanism, we assume there is a limiting value of the apparent density, and let it be ${\mathit{\rho}}_{\mathrm{crt}}^{\mathrm{i}}=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}$. This choice of value is roughly consistent with observations by Magono and Nakamura (1965). If $\stackrel{\mathrm{\u203e}}{{\mathit{\rho}}_{jk}^{\mathrm{i}}}$ is closer to ${\mathit{\rho}}_{\mathrm{crt}}^{\mathrm{i}}$, we consider that the apparent density of the resultant aggregate is closer to the maximum possible density ${\mathit{\rho}}_{jk}^{\mathrm{i},\mathrm{max}}$. Let us assume ${\mathit{\rho}}_{jk}^{\mathrm{i},\mathrm{max}}=\stackrel{\mathrm{\u203e}}{{\mathit{\rho}}_{jk}^{\mathrm{i}}}$.

In the following, we derive equations describing how to update the attributes.

Without loss of generality, assume that *D*_{j}≥*D*_{k}.
For *ϕ*_{j}≤1 (planar),

because we assumed the maximum dimension is preserved. The longest possible minor axis length is ${c}_{j}+min({a}_{k},{c}_{k})$; hence, the largest possible volume becomes ${V}_{\mathrm{max}}=(\mathrm{4}\mathit{\pi}/\mathrm{3}){a}_{j}^{\mathrm{2}}\mathit{\{}{c}_{j}+min({a}_{k},{c}_{k}\left)\mathit{\right\}}$. The minimum possible apparent density ${\mathit{\rho}}_{jk}^{\mathrm{i},\mathrm{min}}$ then becomes

The resultant particle's apparent density is given by a weighted average of ${\mathit{\rho}}_{jk}^{\mathrm{i},\mathrm{max}}=\stackrel{\mathrm{\u203e}}{{\mathit{\rho}}_{jk}^{\mathrm{i}}}$ and ${\mathit{\rho}}_{jk}^{\mathrm{i},\mathrm{min}}$:

where primed values indicate the resultant ice particle. All other attributes are updated as follows:

For *ϕ*_{j}>1 (columnar), the polar axis length is preserved

If approximating the largest possible particle using an ellipsoid, the largest possible volume becomes ${V}_{\mathrm{max}}=(\mathrm{4}\mathit{\pi}/\mathrm{3}){c}_{j}\mathit{\{}{a}_{j}+min({a}_{k},{c}_{k}\left)\mathit{\right\}}max({a}_{j},{a}_{k},{c}_{k})$. Then, the resultant ice particle's apparent density ${\mathit{\rho}}_{j}^{\mathrm{i}\prime}$ can be calculated using Eqs. (61) and (62). Then, the minor axis is updated by

and other attributes are updated by Eqs. (64)–(68).

Note that our aggregation outcome model does not produce particles lighter than ${\mathit{\rho}}_{\mathrm{crt}}^{\mathrm{i}}=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}$.

### 4.1.12 Limitations of our cloud microphysics model

Equations (2)–(70) provide time evolution equations for mixed-phase cloud microphysics. Our model is based on a detailed kinetic description, and all aerosol, cloud, and precipitation particles in the system are followed. The respective activation and deactivation of cloud droplets from and to CCN, and their growth by diffusion and collision are also explicitly predicted. Additionally, the formation of ice particles by condensation/immersion and homogeneous freezing, and gradual morphology changes in ice particles during their growth by diffusion and collision are also predicted explicitly without relying on artificial ice categories or predefined mass–dimension relationships. However, because our basic understanding of mixed-phase cloud microphysics is still insufficient, the introduced models have room for improvement. Further, several processes critical for mixed-phase clouds are ignored for simplicity. For example, collisional breakup of ice particles and rime-splintering are not considered, although they are thought to be responsible for secondary ice production (e.g., Field et al., 2017). In Sect. 9.3, we will discuss more on the limitations and possible future refinements of our model.

## 4.2 Fluid dynamics of moist air

Moist air fluid dynamics can be described by the compressible Navier–Stokes equation for moist air:

where the four terms with the form $\partial \cdot /\partial t{|}_{\mathrm{cm}}$ represent cloud microphysics coupling terms. $\partial \mathit{\rho}/\partial t{|}_{\mathrm{cm}}=\partial \mathit{\rho}{q}_{\mathrm{v}}/\partial t{|}_{\mathrm{cm}}$ is the source of vapor:

Here,
*s*_{v} and *s*_{s} are sources of vapor through condensation/evaporation and deposition/sublimation, respectively:

where *δ*^{3}(** x**) is the three-dimensional Dirac delta
function, and the time derivatives for condensation/evaporation and deposition/sublimation are given by Eqs. (7) and (11), respectively.

$\partial \mathit{\rho}\mathit{\theta}/\partial t{|}_{\mathrm{cm}}$ represents heating due to the phase transition of water:

where *L*_{f} is the latent heat of fusion,
and *s*_{f} is the production rate of liquid water through freezing, melting, or riming.
Let ${t}_{n}^{\mathrm{fz}}$ be the time of the *n*th freezing event and ${i}_{n}^{\mathrm{fz}}$ be the index of the frozen droplet.
Similarly, let ${t}_{n}^{\mathrm{mlt}}$ and ${i}_{n}^{\mathrm{mlt}}$ be the time and melted ice particle of the *n*th melting event, respectively.
Let ${t}_{n}^{\mathrm{rime}}$ and ${i}_{n}^{\mathrm{rime}}$ be the time and rimed droplet of the *n*th riming event, respectively.
Then, *s*_{f} is given by

$\partial \mathit{\rho}\mathit{U}/\partial t{|}_{\mathrm{cm}}$ is the drag force from the particles.
From Eq. (2), we can derive ${\mathit{F}}_{i}^{\mathrm{drg}}={m}_{i}g\widehat{\mathit{z}}+\text{d}\left({m}_{i}{\mathit{v}}_{i}\right)/\text{d}t$.
The terminal velocity assumption does not mean that the second term vanishes
because *m*_{i} and *v*_{i} are still time dependent.
However, even if a droplet accelerated from 0 to 10 m s^{−1} in 100 s
through rapid precipitation development,
the contribution of the second term is much smaller than that of the first term: $\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.25em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}/\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\ll g$.
Thus, we finally obtain

## 4.3 Summary of the section

Now, we have the complete set of the system's time evolution equations: Eqs. (2)–(70) for cloud microphysics (i.e., aerosol, cloud, and precipitation particles) and Eqs. (71)–(81) for cloud dynamics (i.e., moist air). With suitable initial and boundary conditions, our mathematical model can predict mixed-phase cloud behavior. In the next section, we explain how SCALE-SDM solves those time evolution equations numerically.

We develop a numerical model known as SCALE-SDM to solve the mathematical model of mixed-phase clouds presented in the preceding sections.

SCALE is a library of weather and climate models of the Earth and other planets (Nishizawa et al., 2015; Sato et al., 2015, https://scale.riken.jp/, last access: 26 August 2020). We implemented SDM into SCALE version 0.2.5, thus constructing a mixed-phase cloud model called SCALE-SDM 0.2.5-2.2.0.

In our model, we use SDM to solve cloud microphysics as defined by Eqs. (2)–(70). SDM is a particle-based scheme using an efficient Monte Carlo algorithm for coalescence, riming, and aggregation, which enables the accurate simulation of aerosol, cloud, and precipitation particles with lower computational demand (Shima et al., 2009).

Moist air fluid dynamics are solved using SCALE's dynamical core. We solve the compressible Navier–Stokes equation for moist air (Eqs. 71–81) using a forward temporal integration scheme using a finite volume method with an Arakawa-C staggered grid. In this study, we resolve only large eddies and do not use a subgrid-scale (SGS) turbulence model. To stabilize the calculation, we add an artificial fourth-order hyper-diffusion term. Numerical schemes and implementation are described in further detail.

## 5.1 Spatial discretization of moist air

We consider
the density of moist air *ρ*,
density of water vapor *ρ**q*_{v},
momentum of moist air *ρ*** U**,
and mass-weighted potential temperature

*ρ*

*θ*as prognostic variables for moist air. We employ the Arakawa-C staggered grid for discretization:

*ρ*,

*ρ*

*q*

_{v}, and

*ρ*

*θ*are defined at the center of each grid cell, and the three components of

*ρ*

**are defined on the faces of each grid cell. To simplify the notation, we use**

*U*

*G*_{lmn}to denote the status of moist air at each point on the center grid and the face grid. Let Δ

*x*, Δ

*y*, and Δ

*z*represent the grid sizes.

## 5.2 Super-particles and real particles

There are many particles in the atmosphere; thus, it is
practically impossible to follow all of them in a numerical model.
However, it is reasonable to assume that only the collective properties of the particle population
are relevant to predict the behavior of clouds,
because clouds are insensitive to each individual particle.
Therefore, let us approximate
the population of real particles $\mathit{\left\{}\mathit{\right\{}{\mathit{x}}_{i}\left(t\right),{\mathit{a}}_{i}\left(t\right)\mathit{\}},i=\mathrm{1},\mathrm{2},\mathrm{\dots},{N}_{\mathrm{r}}^{\mathrm{wp}}\mathit{\}}$
by a population of super-particles:
$\mathit{\left\{}\mathit{\right\{}{\mathit{\xi}}_{i}\left(t\right),{\mathit{x}}_{i}\left(t\right),{\mathit{a}}_{i}\left(t\right)\mathit{\}},i=\mathrm{1},\mathrm{2},\mathrm{\dots},{N}_{\mathrm{s}}^{\mathrm{wp}}\mathit{\}}$ (see, e.g., Fig. 4 of Grabowski et al., 2019).
A super-particle is characterized by multiplicity *ξ*_{i}, position *x*_{i}, and attributes *a*_{i}.
We consider that the *i*th super-particle
represents *ξ*_{i} real particles *{**x*_{i},*a*_{i}*}*.
Note that multiplicity *ξ*_{i} is an integer and is time dependent.
${N}_{\mathrm{s}}^{\mathrm{wp}}$ is the total number of super-particles accumulated over the whole period.

The relationship between super-particles and real particles can be expressed more precisely as follows.
Let $n(\mathit{a},\mathit{x},t)$ be the particle distribution function, i.e., the mean number density of particles with attributes ** a** at position

**and time**

*x**t*. The following relation then holds:

where 〈⋯〉 denotes the mean, and *δ*^{d}(** a**) is the

*d*-dimensional Dirac delta function. Super-particles reproduce the behavior of particles in expectation:

where $p(\mathit{\xi},\mathit{a},\mathit{x},t)$ is the probability density that a super-particle has
multiplicity *ξ*, attributes ** a**, and position

**at time**

*x**t*;

*I*

_{s}(

*t*) is the set of super-particle indices existing in the domain at time

*t*; and

*N*

_{s}(

*t*):=

*#*

*I*

_{s}(

*t*) is the number of super-particles existing at time

*t*.

## 5.3 Initialization of super-particles

There is an arbitrariness in how to initialize super-particles. In this study, we use the uniform sampling method.

Any probability density function $p(\mathit{\xi},\mathit{a},\mathit{x},t=\mathrm{0})$ that satisfies Eq. (83) can be used to initialize super-particles; however, Unterstrasser et al. (2017) showed that SDM's performance is sensitive to the choice of the probability density function.

Let us consider a specific type of procedure wherein we assign
** a** and

**based on the probability density function**

*x**p*(

**,**

*a***), and determine the super-particle's multiplicity**

*x**ξ*by using a deterministic function of

**and**

*a***, i.e., $\mathit{\xi}=\mathit{\xi}(\mathit{a},\mathit{x})$. Then, Eq. (83) at**

*x**t*=0 reduces to

If we set *ξ*(** a**,

**) as a constant, the probability density function must be proportional to the initial distribution function of real particles: $p(\mathit{a},\mathit{x})\propto n(\mathit{a},\mathit{x},\mathrm{0})$. This so-called constant multiplicity method was adopted in Shima et al. (2009). However, Unterstrasser et al. (2017) found that the numerical convergence of this method regarding the super-particle number is slow. Note that constant multiplicity method is referred to as**

*x**ν*

_{const}-init in Unterstrasser et al. (2017).

Instead, we can set *p*(** a**,

**) as a constant (i.e., uniform sampling). Multiplicity then becomes proportional to the initial distribution function of real particles:**

*x*Using the uniform sampling method, we can more frequently sample rare but important particles in the tail of the distribution, thus improving the numerical convergence. This uniform sampling method was used in various studies (e.g., Arabas and Shima, 2013; Shima et al., 2014; Sato et al., 2017, 2018).

Unterstrasser et al. (2017) proposed several other procedures using a grid, known as SingleSIP-init, multiSIP-init, and
*ν*_{random}-init to more uniformly distribute super-particles along the particle size axis. They confirmed that their methods had much better
performance than the constant
multiplicity method but did not try the uniform sampling
method. Dziekan and Pawlowska (2017) also proposed a similar
procedure. However, both works focused on
coalescence and their initialization procedures are tested only in a zero-dimensional
simulation (box model) with one particle attribute (size). It is questionable
whether their procedures would work efficiently for three-dimensional (3-D) simulations with
several particle attributes. The “discrepancy” of axis-aligned grid
decreases slowly in higher dimensions (e.g., Niederreiter, 1978).
Therefore, an axis-aligned grid is generally unsuitable for sampling high-dimensional spaces.
A uniform sampling method should be more efficient for such a purpose and using
quasi-random numbers would further improve performance.
Meanwhile, as indicated in Grabowski et al. (2018), we should also note that the unbalanced
mass of super-particles could cause larger statistical fluctuations
when super-particles are advected from one grid cell of moist air to another.

Overall, further investigation is required to determine an optimal method for initializing super-particles. In this study, we use the uniform sampling method given by Eq. (85). More details of our procedure will be specified in Sect. 6.1.7. As shown in Fig. 9, our model's numerical convergence regarding super-particle numbers is good for at least the 2-D cumulonimbus simulation that we will conduct to evaluate our model.

## 5.4 Operator splitting of the time integration

We separately evaluate each process using the first-order operator splitting scheme. Let Δ*t* be the common time step. Here, we explain
how $\mathit{\left\{}\mathit{\right\{}{\mathit{\xi}}_{i},{\mathit{x}}_{i},{\mathit{a}}_{i}\mathit{\left\}}\mathit{\right\}}$ and *G*_{lmn} are updated from time *t* to *t*+Δ*t*.

Let Δ*t*_{adv}, Δ*t*_{fz∕mlt}, Δ*t*_{cnd∕evp}, Δ*t*_{dep∕sbl},
and Δ*t*_{collis} be the time steps for the advection and sedimentation of particles, freezing and melting,
condensation and evaporation, deposition and sublimation, and collision–coalescence, –riming, and –aggregation, respectively.

Let Δ*t*_{dyn} be the time step for
moist air fluid dynamics.

These process time steps are all divisors of the common time step Δ*t*.

We first calculate fluid dynamics without the coupling terms from particles to moist air (Eqs. 76–81),
and update moist air from *G*_{lmn}(*t*) to ${\mathit{G}}_{lmn}^{\prime}\left(t\right)$.
Then, we update super-particles $\mathit{\left\{}\mathit{\right\{}{\mathit{\xi}}_{i},{\mathit{x}}_{i},{\mathit{a}}_{i}\mathit{\left\}}\mathit{\right\}}$ from *t* to *t*+Δ*t*.
We select one elementary cloud microphysics process, integrate it forward by one time step, and then move on to the next process.
Here, processes lagging in time are calculated preferentially.
Simultaneously, we evaluate feedback from the particles to moist air through the coupling terms (Eqs. 76–81),
and update the moist air from ${\mathit{G}}_{lmn}^{\prime}\left(t\right)$ to *G*_{lmn}(*t*+Δ*t*).
Table 1 shows an example of the calculation order.

## 5.5 Time integration of cloud microphysics

We use SDM to solve cloud microphysics.
We provide details of the numerical schemes used to calculate cloud microphysics in this section.
The state of ambient air *G*_{i}:=** G**(

*x*_{i}) around a super-particle

*i*is often needed. For scalar variables, we use the value at the center point of the grid cell in which the super-particle is located, whereas we interpolate wind velocities from face grids, as detailed in the next section.

### 5.5.1 Advection and sedimentation

For each super-particle, the motion equation (Eq. 3) is solved using a time step Δ*t*_{adv}.
We normally select a short enough Δ*t*_{adv} to satisfy the Courant–Friedrichs–Lewy (CFL) condition for wind velocity.
So that we can predict the particle number concentration accurately,
we use the predictor-corrector scheme with the “simple linear interpolation” of wind velocities from the face grid following Grabowski et al. (2018).
The momentum *ρ*** U** is defined on the face grid and
density

*ρ*is defined on the center grid. Therefore, we average the

*ρ*

_{lmn}on both sides of the face grid to calculate wind velocity

*U*_{lmn}on the face grid. We then interpolate

*U*_{lmn}to the super-particle position using the simple linear scheme of Grabowski et al. (2018), which ensures that the wind velocity divergence over any subgrid volume becomes exactly the same as that over the grid cell volume.

The reaction force acting on moist air is calculated using Eq. (81).
Feedback from each super-particle is imposed only on the (*ρ**W*)_{lmn} nearest to the super-particle.

### 5.5.2 Freezing and melting

Every Δ*t*_{fz∕mlt} interval,
for each super-particle,
freezing and melting are examined following the model detailed in Sects. 4.1.4 and 4.1.5.
The exchange of latent heat of fusion is calculated using Eqs. (74), (79), and (80).
Feedback from each super-particle is imposed only on the grid cell where the super-particle is located.

### 5.5.3 Condensation and evaporation

For each super-droplet, we solve the condensation and evaporation equation (Eq. 8)
with a time step of Δ*t*_{cnd∕evp}.
The activation/deactivation timescale is much shorter than that of other processes.
To eliminate stiffness, we convert the equation to the time evolution equation of *r*^{2} following Shima et al. (2009)
and adopt the backward Euler scheme.

The exchange of vapor and latent heat with moist air is calculated using Eqs. (71), (72), (74), (76), (77), and (79). Feedback from each super-droplet is imposed only on the grid cell where the super-droplet is located.

The growth of droplets is calculated implicitly;
however, the evolution of supersaturation through feedback is calculated explicitly.
Therefore, the length of Δ*t*_{cnd∕evp} is restricted mostly by supersaturation's phase relaxation time regarding condensation and evaporation,
which is the timescale on which a supersaturation fluctuation decays through condensation or evaporation.

### 5.5.4 Deposition and sublimation

For each ice super-particle, we solve the deposition and sublimation time evolution equations
detailed in Sect. 4.1.7 using the time step Δ*t*_{dep∕sbl}.
Contrary to the condensation and evaporation equation (Eq. 8),
the time evolution equation of mass (Eq. 11) is not stiff
because the curvature term is ignored and the solute effect does not exist.
Let us convert the equation to the time evolution equation of *m*^{2∕3}.
Then, in a situation when the ice particle is spherical and, at the same time, so small that the ventilation effect can be ignored,
then the equation reduces to $\text{d}{m}^{\mathrm{2}/\mathrm{3}}/\text{d}t=\mathrm{const}.$; i.e., the right-hand side (r.h.s.) does not depend on *m*.
Inspired by this fact, we adopt the forward Euler scheme to solve the time evolution equation of *m*^{2∕3}
even when the ice particle is not spherical or small.

The exchange of vapor and latent heat with moist air is calculated using Eqs. (71), (72), (74), (76), (78), and (79). Feedback from each ice super-particle is imposed only on the grid cell where the ice super-particle is located.

Δ*t*_{dep∕sbl} is restricted by the timescale of individual ice particle growth through deposition and sublimation,
and the phase relaxation time of supersaturation regarding deposition and sublimation.

### 5.5.5 Coalescence, riming, and aggregation

The stochastic process of coalescence, riming, and aggregation detailed in Sects. 4.1.8–4.1.11
is solved using the Monte Carlo algorithm of SDM (Shima et al., 2009).
The computational cost of this algorithm is proportional to the number of super-particles *O*(*N*_{s}),
which is achieved by an efficient collision candidate pair number reduction technique.
An additional advantage of this technique is the parallelizability of computation;
each super-particle belongs to only one candidate pair, and hence dependencies are eliminated.

Δ*t*_{collis} can be determined using the argument presented in the last paragraph of Sect. 5.1.3 in Shima et al. (2009);
however, here we repeat it in a slightly different way to provide a precise physical interpretation.
In short, the time step Δ*t*_{collis} is restricted by the mean free time of a particle,
i.e., the average waiting time for a particle between two successive coalescence/riming/aggregation events.
Let $\stackrel{\mathrm{\u203e}}{P}$ be the typical probability that a particle coalescence/rime/aggregate
with another particle within a small time interval Δ*t*_{collis}.
From Eq. (31), $\stackrel{\mathrm{\u203e}}{P}$ can be evaluated as

where ${N}_{\mathrm{r}}^{\prime}$ is the number of real particles in a volume Δ*V*, $\stackrel{\mathrm{\u203e}}{K}$ is
the typical value of the coalescence/riming/aggregation kernel *K*,
and *n*_{r} is the number concentration of real particles.
Requiring that $\stackrel{\mathrm{\u203e}}{P}<\mathrm{1}$ has to be satisfied,
we obtain

Here, we relate the above argument to that of Shima et al. (2009). Let $\stackrel{\mathrm{\u203e}}{{P}_{\mathrm{s}}}$ be the typical probability that a collision candidate super-particle pair coalescence/rime/aggregate after the pair number reduction technique is applied. Note that $\stackrel{\mathrm{\u203e}}{{P}_{\mathrm{s}}}$ is what Shima et al. (2009) evaluated in the last paragraph of Sect. 5.1.3. We can derive $\stackrel{\mathrm{\u203e}}{{P}_{\mathrm{s}}}\approx \stackrel{\mathrm{\u203e}}{P}$ as follows:

where ${N}_{\mathrm{s}}^{\prime}$ is the number of super-particles in the volume Δ*V*,
the first term *{*…*}* represents the scale-up factor due to the candidate pair number reduction,
and $\stackrel{\mathrm{\u203e}}{\mathit{\xi}}\approx {N}_{\mathrm{r}}^{\prime}/{N}_{\mathrm{s}}^{\prime}$ is the typical multiplicity.

In SDM, the multiple coalescence technique is used to make the algorithm robust to larger Δ*t*_{collis}.
Here, we clarify how we adapt it to riming and aggregation.
If it is a coalescence between a droplet *j* and $\stackrel{\mathrm{\u0303}}{\mathit{\gamma}}$ number of droplets *k*
(see Sect. 5.1.3 of Shima et al., 2009, for the definition of $\stackrel{\mathrm{\u0303}}{\mathit{\gamma}}$),
we modify Eqs. (33)–(35) by applying

If it is a coalescence between $\stackrel{\mathrm{\u0303}}{\mathit{\gamma}}$ number of droplets *j* and a droplet *k*,
we apply

Similarly, if it is a riming/aggregation between a particle *j* and $\stackrel{\mathrm{\u0303}}{\mathit{\gamma}}$ number of particles *k*,
we apply the following replacement to Eqs. (42)–(54)
and (60)–(70):

If it is a riming/aggregation between $\stackrel{\mathrm{\u0303}}{\mathit{\gamma}}$ number of particles *j* and a particle *k*,

What is not straightforward is the calculation of *V*_{max} used in the aggregation outcome formula.
For planar collector *j*, we consider that *V*_{max} is given by

For columnar collector *j*,

The exchange of the latent heat of fusion due to riming is calculated using Eqs. (74), (79), and (80). Feedback from each super-particle is imposed only on the grid cell where the super-particle is located.

## 5.6 Time integration of moist air fluid dynamics

Moist air fluid dynamics is governed by the compressible Navier–Stokes equation (Eqs. 71–81). In this study, as explained in the previous section, the four coupling terms from cloud microphysics denoted by $\partial \cdot /\partial t{|}_{\mathrm{cm}}$ are evaluated when calculating cloud microphysics.

We solve the compressible Navier–Stokes equation without the coupling terms using a finite volume method with an Arakawa-C staggered grid. For spatial discretization, the fourth-order central difference scheme is used for advection terms and the second-order central difference scheme is used for other spatial derivatives. To preserve the monotonicity, we apply the flux-corrected transport scheme of Zalesak (1979) to water vapor advection. For time integration, we use the three-step Runge–Kutta scheme of Wicker and Skamarock (2002). An artificial, fourth-order hyper-diffusion term is added to stabilize the calculation. For this study, we set the non-dimensional diffusion coefficient defined in Eq. (A132) of Nishizawa et al. (2015) as $\mathit{\gamma}={\mathrm{10}}^{-\mathrm{3}}$. For more details of the numerical schemes used for fluid dynamics, see Nishizawa et al. (2015) and Sato et al. (2015).

The time step Δ*t*_{dyn} must satisfy the CFL condition of acoustic waves.

The preceding sections described the basic equations and numerical implementation of SCALE-SDM. To evaluate our numerical model's performance, we conduct a 2-D simulation of an isolated cumulonimbus following the setup of Khain et al. (2004). In this section, we first describe the atmospheric conditions and numerical parameters used for the control case denoted by CTRL. To evaluate fluctuation, we conduct a 10-member ensemble of simulations by changing the pseudo-random number sequence. To investigate the simulation's numerical convergence, we will change the super-particle number concentration, grid sizes, and time steps of CTRL. Those ensembles are denoted by NSP, DX, and DT, respectively. Our choice of parameters is specified in the subsequent sections. Table 2 summarizes the model setup for all cases.

## 6.1 Control ensemble (CTRL)

In this section, we specify the atmospheric conditions and numerical parameters used for the CTRL ensemble.

### 6.1.1 Initial moist air conditions

The domain is 2-D (*x*–*z*),
60 km in the horizontal direction and
16 km in the vertical direction.

The initial atmospheric profile is horizontally uniform, and the
vertical moist air profile is given by sounding data from Midland, Texas, on 13 August 1999, as shown in Fig. 4 of Khain et al. (2004).
The cloud base and freezing level are at about 2.2 km
(14 ^{∘}C) and 4.1 km, respectively.
We consider that the wind is initially horizontal
and wind velocity increases from 4 m s^{−1} near the surface
to 7 m s^{−1} at 400 hPa,
and remains unchanged at higher levels.

### 6.1.2 Moist air boundary conditions

For the lateral boundaries, we impose periodic boundary conditions.
For the upper and lower boundaries, we set the vertical wind velocity *W* to zero,
i.e., a zero-fixed boundary condition for vertical momentum *ρ**W*,
and no flux boundary conditions for other prognostic variables.

### 6.1.3 Initial conditions of particles

Initially, the particles are distributed uniformly in space at random, and consist of pure ammonium bisulfate aerosol particles and mineral dust internally mixed with ammonium bisulfate.

The initial number-size distribution of the population of pure ammonium bisulfate particles is given by a bimodal log-normal distribution,

where ${r}_{\mathrm{dry}}^{\mathrm{sulf}}$ is the dry radius of the ammonium bisulfate component
and *N*^{sulf} is the accumulated number of particles smaller than ${r}_{\mathrm{dry}}^{\mathrm{sulf}}$ per unit volume of air.
The particle number concentrations are ${c}_{\mathrm{1}}^{\mathrm{sulf}}=\mathrm{270}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$ and ${c}_{\mathrm{2}}^{\mathrm{sulf}}=\mathrm{45}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$;
thus, the total particle number concentration is ${c}^{\mathrm{sulf}}={c}_{\mathrm{1}}^{\mathrm{sulf}}+{c}_{\mathrm{2}}^{\mathrm{sulf}}=\mathrm{315}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$.
The geometric mean radii are *r*_{1}=0.03 µm and *r*_{2}=0.14 µm,
with geometric standard deviations of *σ*_{1}=1.28 and *σ*_{2}=1.75, respectively.
This distribution is based on in situ maritime aerosol data as detailed in Sect. 2.2.3 of VanZanten et al. (2011),
but the number concentration is multiplied by 3.
As discussed in Sect. 2.4, we consider that
a droplet containing only soluble substances freezes only through a homogeneous freezing mechanism;
therefore, the freezing temperature of these particles is ${T}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C.
Therefore, pure ammonium bisulfate's initial distribution function can be calculated as

The other aerosol population consists of mineral dust internally mixed with ammonium bisulfate. We set the number concentration to ${c}^{\mathrm{dust}}=\mathrm{1}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$,
and for simplicity, set the mineral dust particle diameter to *d*^{dust}=1 µm initially (see, e.g., Fig. 3 of Hoose et al., 2010).
We assume that the size distribution of internally mixed ammonium bisulfate is the same as that of the pure ammonium bisulfate given by Eq. (95).
The probability density function of the freezing temperature *p*(*T*^{fz}) is given by Eq. (1).
Here, we use the INAS density formula from Niemand et al. (2012), but based on the discussion in Niedermeier et al. (2015),
we do not extrapolate the formula to lower or higher temperatures:

where ${T}_{\mathrm{max}}^{\mathrm{fz}}=-\mathrm{12}$ ^{∘}C and ${T}_{\mathrm{min}}^{\mathrm{fz}}=-\mathrm{36}$ ^{∘}C.
The mineral dust surface area is given by *A*^{insol}=*π*(*d*^{dust})^{2}.
As discussed in Sect. 2.4, we set ${T}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C
if the mineral dust is IN inactive and no INAS appears until ${T}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C.
Altogether, the mineral dust distribution function is given by

where *H*(*T*) is the Heaviside step function and ${P}_{\mathrm{INia}}:=P({T}^{\mathrm{fz}}\le -\mathrm{38}$ ^{∘}C) is the probability that
a single INAS does not appear until ${T}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C.
For *d*^{dust}=1 µm, *P*_{INia}≈0.056.

### 6.1.4 Boundary conditions for particles

We also impose periodic boundary conditions on particles for the lateral boundaries. If a particle crosses the upper or lower boundary, we remove that particle from the system.

### 6.1.5 Near-surface heating

Convective cloud development is triggered by a 20 min heating started from the beginning within a 10 km wide region
centered at *x*=5 km, and is expressed as

where $H=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{K}\phantom{\rule{0.125em}{0ex}}{\mathrm{h}}^{-\mathrm{1}}$, *x*_{0}=5 km, *w*=10 km, and *z*_{0}=0.5 km.
The heating has a parabolic shape in the horizontal direction
and decays exponentially in the vertical direction.

### 6.1.6 Grid size and time steps

We use a uniform grid throughout this study, with a grid size of $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z=\mathrm{62.5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$ in the CTRL case.
The time steps in the CTRL case are Δ*t*=0.4 s,
$\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{0.4}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$,
Δ*t*_{collis}=0.2 s,
$\mathrm{\Delta}{t}_{\mathrm{cnd}/\mathrm{evp}}=\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}}=\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$,
and Δ*t*_{dyn}=0.05 s.

### 6.1.7 Initialization of super-particles

Initially, the super-particles are distributed uniformly throughout the domain at random
with a number concentration of *c*^{SP}=128 per cell.
We consider half of them as pure ammonium bisulfate aerosol particles,
a few of them as IN inactive mineral dust particles internally mixed with ammonium bisulfate,
and the remainder to be IN active mineral dust particles internally mixed with ammonium bisulfate.

The multiplicity, ammonium bisulfate mass, and freezing temperature of each pure ammonium bisulfate super-particle is assigned as follows. For each pure ammonium bisulfate super-particle, we draw a random number uniformly in log-space from the interval $[{r}_{\mathrm{dry},\mathrm{min}}^{\mathrm{sulf}},{r}_{\mathrm{dry},\mathrm{max}}^{\mathrm{sulf}}]$ and determine the dry radius ${r}_{\mathrm{dry},i}^{\mathrm{sulf}}$. To accurately represent the size distribution given in Eq. (95), we set ${r}_{\mathrm{dry},\mathrm{min}}^{\mathrm{sulf}}=\mathrm{10.0}\phantom{\rule{0.125em}{0ex}}\mathrm{nm}$ and ${r}_{\mathrm{dry},\mathrm{max}}^{\mathrm{sulf}}=\mathrm{5.0}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m}$. From Eqs. (85) and (96), the super-particle's multiplicity is then given by

where *V*_{domain} is the total volume of the domain, *n* and *p* in Eq. (85)
in this case are given by

and *N*_{s}(0) in Eq. (85) is replaced by *N*_{s}(0)∕2
because we use half of the super-particles for pure ammonium bisulfate aerosol particles.
The ammonium bisulfate mass is calculated from the dry radius ${r}_{\mathrm{dry},i}^{\mathrm{sulf}}$
as ${m}_{\mathrm{1}i}^{\mathrm{sol}}=(\mathrm{4}\mathit{\pi}/\mathrm{3}){\mathit{\rho}}_{(\mathrm{NH}{)}_{\mathrm{4}}{\mathrm{HSO}}_{\mathrm{4}}}({r}_{\mathrm{dry},i}^{\mathrm{sulf}}{)}^{\mathrm{3}}$,
where ${\mathit{\rho}}_{(\mathrm{NH}{)}_{\mathrm{4}}{\mathrm{HSO}}_{\mathrm{4}}}=\mathrm{1.78}\phantom{\rule{0.125em}{0ex}}\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$.
The soluble aerosol particle freezing temperature is ${T}_{i}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C.

For IN inactive mineral dust super-particles, we use ${P}_{\mathrm{INia}}^{\mathrm{SP}}=\mathrm{0.05}$.
The mineral dust initially has the same size *d*^{dust}=1 µm. The dry radius ${r}_{\mathrm{dry},i}^{\mathrm{sulf}}$ is calculated
using the same procedure as the pure ammonium bisulfate aerosol particles; i.e.,
for each super-particle,
we draw a random number uniformly in log-space from the interval
$[{r}_{\mathrm{dry},\mathrm{min}}^{\mathrm{sulf}},{r}_{\mathrm{dry},\mathrm{max}}^{\mathrm{sulf}}]$.
The IN inactive mineral dust freezing temperature is ${T}_{i}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C.
From Eqs. (85) and (98),
an IN inactive mineral dust super-particle's multiplicity is then given by

Finally, we consider IN active mineral dust internally mixed with ammonium bisulfate.
The remaining super-particles, i.e., $(\mathrm{1}-{P}_{\mathrm{INia}}^{\mathrm{SP}})/\mathrm{2}$,
are used for this population.
The initial diameter of the mineral dust is *d*^{dust}=1 µm,
and the dry radius ${r}_{\mathrm{dry},i}^{\mathrm{sulf}}$ is determined as in the other populations.
We draw another random number uniformly from the interval
$[{T}_{\mathrm{min}}^{\mathrm{fz}},{T}_{\mathrm{max}}^{\mathrm{fz}}]$
and determine the freezing temperature ${T}_{i}^{\mathrm{fz}}$.
From Eqs. (85) and (98),
an IN active mineral dust super-particle's multiplicity is then given by

Note that multiplicity *ξ*_{i} is an integer variable.
We round the r.h.s. of Eqs. (101)–(105) to the nearest integer, and
if the r.h.s. is <1,
we draw a random number to decide whether to choose *ξ*_{i}=1 or *ξ*_{i}=0 to avoid sampling error.
If *ξ*_{i}=0, the super-particle will be removed from the system.

Assuming that all the particles are deliquescent,
we consider that the initial droplet radius *r*_{i} is equal to the equilibrium radius of condensation/evaporation growth equation (Eq. 8).
As the vapor profile is initially subsaturated relative to liquid water
and all particles contain soluble substances,
the growth equation (Eq. 8) has a unique, stable equilibrium solution.

### 6.1.8 Pseudo-random numbers

To evaluate the fluctuation, we conduct a 10-member ensemble of simulations by changing the pseudo-random number sequence.

Now, the atmospheric conditions and numerical parameters used for the CTRL ensemble have all been specified.

## 6.2 Other ensembles for investigating numerical convergence

We also try various other test cases by changing the CTRL ensemble's numerical parameters, and assess the sensitivity of results to numerical parameters. Our parameter selections are specified in the following sections and a summary is provided in Table 2.

### 6.2.1 NSP ensembles for super-particle number convergence

To investigate numerical convergence with respect to initial the super-particle number concentration *c*^{SP},
we vary *c*^{SP} as follows: 2, 4, …, 512 per cell.
Grid size and time steps are not changed.
These cases are respectively denoted by NSP002, NSP004, …, NSP512.
Note that NSP128 and CTRL are the same.

### 6.2.2 DX ensembles for grid convergence

To investigate numerical convergence with respect to the grid size, we run ensembles using different grid sizes.

The grid size of the DXx4 ensemble is 4 times that of CTRL:
$\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z=\mathrm{250}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$,
Δ*t*_{dyn}=0.2 s,
and $\mathrm{\Delta}t=\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{1.6}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$,
but other time steps, i.e., $\mathit{\{}\mathrm{\Delta}{t}_{\mathrm{collis}},\mathrm{\Delta}{t}_{\mathrm{cnd}/\mathrm{evp}},\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}}\mathit{\}}$,
are not changed.

The DXx2 ensemble's grid size is twice that of CTRL:
$\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z=\mathrm{125}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$,
Δ*t*_{dyn}=0.1 s,
and $\mathrm{\Delta}t=\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{0.8}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$.

Note that DXx1 and CTRL are the same.

The DX/2 ensemble has a grid size that is half that of CTRL:
$\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z=\mathrm{31.25}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$.
Δ*t*_{dyn}=0.025 s,
and $\mathrm{\Delta}t=\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{0.2}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$.

### 6.2.3 DT ensembles for time step convergence

To investigate numerical convergence with respect to the cloud microphysics time steps, we change the cloud microphysics time steps for CTRL without changing the time step for fluid dynamics.

The time steps for the DTx10 ensemble's cloud microphysics are 10 times that of CTRL:
$\mathrm{\Delta}t=\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{4.0}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$,
Δ*t*_{collis}=2.0 s, and
$\mathrm{\Delta}{t}_{\mathrm{cnd}/\mathrm{evp}}=\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}}=\mathrm{1.0}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$.
Δ*t*_{dyn} is not changed.

The time steps of the DTx5 ensemble are 5 times that of CTRL:
$\mathrm{\Delta}t=\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{2.0}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$,
Δ*t*_{collis}=1.0 s, and
$\mathrm{\Delta}{t}_{\mathrm{cnd}/\mathrm{evp}}=\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}}=\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$.

The time steps of the DTx2 ensemble are twice that of CTRL:
$\mathrm{\Delta}t=\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{0.8}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$,
Δ*t*_{collis}=0.4 s, and
$\mathrm{\Delta}{t}_{\mathrm{cnd}/\mathrm{evp}}=\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}}=\mathrm{0.2}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$.

Note that DTx1 and CTRL are the same.

The time steps of the DT/2 ensemble are half that of CTRL:
$\mathrm{\Delta}t=\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{0.2}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$,
Δ*t*_{collis}=0.1 s, and
$\mathrm{\Delta}{t}_{\mathrm{cnd}/\mathrm{evp}}=\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}}=\mathrm{0.05}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$.

The time steps of the DT/4 ensemble are one-quarter that of CTRL:
$\mathrm{\Delta}t=\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}{t}_{\mathrm{fz}/\mathrm{mlt}}=\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$,
Δ*t*_{collis}=0.05 s, and
$\mathrm{\Delta}{t}_{\mathrm{cnd}/\mathrm{evp}}=\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}}=\mathrm{0.025}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$.

From the 10 CTRL ensemble members, we selected the one that produced accumulated precipitation amounts closest to the mean value as the representative, hereafter referred to as the typical realization of CTRL. In this section, we analyze these results in detail.

## 7.1 Hydrometeor categorization

We do not categorize hydrometeors during the simulation, which is one of the salient features of our model because the artificial partitioning of hydrometeors could cause various artifacts. In contrast, when analyzing results, dividing hydrometeors into categories is useful to precisely understand the results.

In this study, we assume that hydrometeors completely freeze or melt instantaneously (see Sect. 4.1.4 and 4.1.5). Further, we assume that all particles contain soluble components and are hygroscopic. If not frozen, we assume that the particles are deliquescent even when humidity is low (see Sect. 4.1.6). We also introduced a limiter (Eq. 14) to prevent complete sublimation. Hence, all particles can be categorized as either droplets or ice particles with no ambiguity.

If a particle is a droplet and its radius *r* is <40 µm,
we consider it a cloud droplet. Otherwise,
the particle is considered a rain droplet.

If a particle is an ice particle with a rimed mass fraction satisfying ${m}^{\mathrm{rime}}/m>\mathrm{0.3}$,
we consider it a graupel particle. This criterion is
based on the riming categories in Fig. 5 of Mosimann et al. (1994), in which
0.3 corresponds to a densely rimed ice crystal.
If the maximum dimension of a graupel particle is >5 mm, we consider it a hailstone. However,
for the sake of simplicity, we consider hailstones as a subset of graupel
and they will not be distinguished in the figures.
If the ice particle is not a graupel particle
but is rather composed of >10 monomers, i.e., *n*^{mono}>10,
we consider the ice particle a snow aggregate.
Otherwise, we categorize the ice particle as a cloud ice particle.

## 7.2 Spatial structure of the cloud, water path, and precipitation amount

We first analyze the cloud's overall properties, and then, in the next section, we analyze the properties of individual ice particles.

Figure 1 shows how the cloud's spatial structure in the typical realization of CTRL evolved over time. The mixing ratio of cloud water, rainwater, cloud ice, graupel, and snow aggregates are plotted in fading white, yellow, blue, red, and green, respectively. See also Movie 1 in the video supplement.

Figure 2 shows how the amounts of hydrometeors in the atmosphere evolved over time. The domain-averaged cloud water, rainwater, cloud ice water, graupel water, and snow aggregate water paths are plotted in gray, yellow, blue, red, and green, respectively. Figure 3 shows the time evolution of domain-averaged accumulated precipitation amounts. The solid lines represent the typical realization of CTRL in both figures. The dark shades indicate the mean ± standard deviation that were calculated using the 10 CTRL ensemble members. The unbiased estimator was used to calculate the standard deviation. Pale shades indicate the maximum and minimum values of the 10 ensemble members.

The cloud started to form at approximately *t*=1200 s, and
at approximately *t*=1900 s, rain droplets started to be created through warm rain microphysics processes.
Soon after that, supercooled droplets near the cloud top started to freeze and
the number of supercooled cloud droplets quickly decreased because of
the Wegener–Bergeron–Findeisen process and riming.
At the same time, we also observed that convective cores
near the homogeneous freezing level (*z*≈9.3 km) containing high liquid water content were sustained until around *t*=5000 s.
For example, at *t*=3300 s, we observed a liquid water content of 2.1 g m^{−3}
at $(x,z)=(\mathrm{21.8}\phantom{\rule{0.125em}{0ex}}\mathrm{km},\mathrm{9.8}\phantom{\rule{0.125em}{0ex}}\mathrm{km})$, where $T=-\mathrm{37.5}$ ^{∘}C.
The existence of such high supercooled liquid water content down to the homogeneous freezing limit −38 ^{∘}C
are frequently observed in deep convective clouds (Rosenfeld and Woodley, 2000).
The ice particles quickly evolved into graupel particles through riming,
and then fell toward the surface. When crossing the freezing level at approximately *z*=4.1 km,
the graupel instantaneously melted into rain droplets, based on our model.
The peak of the rainwater path at *t*=2800 s was created by graupel melting.
The first rain droplet reached the surface at about *t*=2800 s, and
heavy precipitation was sustained for 1200 s, followed by weak precipitation.
At the end of the simulation (*t*=5400 s),
the domain-averaged accumulated precipitation amount was 1.2 mm.
An anvil cloud was created between *z*=10 km and *z*=12 km.
The anvil cloud was mostly composed of cloud ice particles,
with a small amount of snow aggregates that increased slowly over time
through the aggregation of cloud ice particles.
The maximum updraft and downdraft speeds were 39.0 and 21.9 m s^{−1},
which were observed at
$(t,x,z)=(\mathrm{2340}\phantom{\rule{0.125em}{0ex}}\mathrm{s},\mathrm{12.8}\phantom{\rule{0.125em}{0ex}}\mathrm{km},\mathrm{11.1}\phantom{\rule{0.125em}{0ex}}\mathrm{km})$
and
$(\mathrm{1620}\phantom{\rule{0.125em}{0ex}}\mathrm{s},\mathrm{9.5}\phantom{\rule{0.125em}{0ex}}\mathrm{km},\mathrm{4.1}\phantom{\rule{0.125em}{0ex}}\mathrm{km})$,
respectively.

Our model successfully simulated the life cycle of a cumulonimbus typically observed in nature (see, e.g., Chap. 8 of Cotton et al., 2010). At the same time, our results are limited because the simulation was conducted in 2-D; the turbulence characteristics are different in 2-D and 3-D. Furthermore, the convection was initiated from a stratified, non-turbulent atmosphere; however, this is unrealistic. Following Lasher-Trapp et al. (2005), imposing a spin-up period to develop turbulence in the boundary layer before initiating the deep convection would be desirable.

## 7.3 Ice particle morphology and fall speeds

Now, we analyze the properties of individual ice particles in the typical realization of CTRL.

Figure 4 shows the mass–dimension relationship of ice particles at *t*=2040 s (towering stage), 3000 s (mature stage), and 5400 s (dissipating stage).
The 2-D mass densities of cloud ice particles, graupel particles, and snow aggregates are plotted in fading blue, red, and green, respectively.
The horizontal axis represents the maximum ice particle dimension *D*.
The vertical axis represents the normalized ice particle mass *m*^{*},
which is defined by the ratio of ice particle mass *m* to
the mass of a spherical ice particle with the same maximum dimension *D*
and the true density of ice ${\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$:

Note that ${m}^{*}\le \mathrm{1}$ always holds.
To calculate 2-D mass densities, we divided the 2-D *D*–*m*^{*} space into
100×100 bins, accumulated the masses of ice particles in each bin,
and divided the total masses by the area of each bin measured in log _{10}(*D*) and log _{10}(*m*^{*}).
The colored slopes in Fig. 4 represent mass–dimension relationship formulas from various studies, and
M96, HK87, K89, M90, and LH74 indicate
Mitchell (1996), Heymsfield and Kajikawa (1987), Kajikawa (1989), Mitchell et al. (1990), and Locatelli and Hobbs (1974), respectively.
Note that “crystals with sector like branches (M96)” and “stellar crystals with broad arms (M96)”
consists of two slopes, respectively, but both are not continuous.
See also Movie 2 in the video supplement.

Figure 5 shows the relationship between ice particle aspect ratios and dimensions
at *t*=2040, 3000, and 5400 s.
The horizontal axis represents the maximum ice particle dimension *D*, and
the vertical axis represents the ice particle aspect ratio *ϕ*.
The 2-D mass densities of cloud ice particles, graupel particles, and snow aggregates are plotted
in the same manner as Fig. 4, except for differences in the vertical axis.
Note that if *ϕ*>1 (*ϕ*<1), ice particles are columnar (planar).
See also Movie 3 in the video supplement.

Figure 6 shows the relationship between ice particle apparent densities and dimensions
at *t*=2040, 3000, and 5400 s.
The horizontal axis represents the maximum ice particle dimension *D*, and
the vertical axis represents ice particle apparent density *ρ*^{i}.
The 2-D mass densities of cloud ice particles, graupel particles, and snow aggregates are plotted
in the same manner as Fig. 4, except for differences in the vertical axis.
See also Movie 4 in the video supplement.

Figure 7 shows the relationship between ice particle velocities and dimensions
at *t*=2040, 3000, and 5400 s.
The horizontal axis represents the maximum ice particle dimension *D*, and
the vertical axis represents ice particle terminal velocity *v*^{∞}.
The 2-D mass densities of cloud ice particles, graupel particles, and snow aggregates are plotted in the same manner as Fig. 4, except for differences in the vertical axis. The colored slopes in Fig. 7 represent the velocity–dimension relationship formulas from various studies, and SC85, W08, H72, KH83, A72, and H02 indicate
Starr and Cox (1985), Westbrook et al. (2008), Heymsfield (1972), Knight and Heymsfield (1983), Auer (1972), and Heymsfield et al. (2002), respectively. “Stokes' law for ice spheres” is based on the Stokes' terminal velocity for spherical ice particles with
the true ice density. We use the dynamic viscosity at a temperature of −20 ^{∘}C,
i.e., $\mathit{\mu}=\mathrm{1.630}\times {\mathrm{10}}^{-\mathrm{5}}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$.
The two slopes of W08 are based on the analytical formula of Westbrook et al. (2008) for <100 µm ice particles.
For “hexagonal plates”, $L/\mathrm{2}a=\mathrm{0.05}$ is assumed, with *L* being the height of the hexagonal prism
and $a=D/\mathrm{2}$ being the hexagon's maximal radius. The effective radius is calculated using the
horizontal orientation model from Roscoe (1949).
For “hexagonal columns”, $L/\mathrm{2}a=\mathrm{20}$ is assumed, and the effective radius is calculated using the
random orientation model of Hubbard and Douglas (1993). In both cases,
we use the dynamic viscosity at a temperature of −20 ^{∘}C.
See also Movie 5 in the video supplement.

At *t*=2040 s (towering stage), cloud glaciation had just started,
and a small amount of planar and columnar cloud ice particles and graupel particles can be observed.
The two horizontal red bands at *ϕ*=0.8 and $\mathit{\varphi}=\mathrm{1}/\mathrm{0.8}$ in Fig. 5 were created
because of our assumption that riming growth eventually makes ice particles quasi-spherical.

At *t*=3000 s (mature stage), many hailstones (graupel particles >5 mm) can be observed
in the cloud's middle layer. We also have many columnar cloud ice particles and a small number of snow aggregates
in the upper part of the cloud.
Those cloud ice particles were columnar because our model's inherent growth ratio Γ(*T*)
is >1 at this height.
Many of the snow aggregates were spherical because our model assumed that
the aspect ratio *ϕ* approaches 1 as aggregation occurs.

At *t*=5400 s (dissipating stage), most of the graupel particles
had fallen away and only a small amount remained.
More columnar cloud ice particles and snow aggregates can be observed in the anvil.

The mass–dimension relationship shown in Fig. 4, and
the velocity–dimension relationship shown in Fig. 7 show
a reasonable agreement between our model's predicted results
and existing formulas based on laboratory measurements and observations.
In both figures, cloud ice particles, graupel particles, and snow aggregates
are distributed near the blue, red, and green slopes, respectively.
In Fig. 7, one might note that snow aggregates in our
model fall faster than those in the formulas; however, this bias can be explained by
the air density dependence of the fall speed. The green slopes in
Fig. 7 represent the formulas of LH74 and H02. LH74's
formulas are constructed from data measured between altitudes of 750
and 1500 m above sea level; hence, the density is approximately
1.1 kg m^{−3}. H02's formula is for temperature
and pressure of −10 ^{∘}C and 500 hPa;
hence, the density is approximately 0.66 kg m^{−3}.
In our simulation, most of the snow aggregates exist in the anvil cloud,
wherein the density is approximately 0.38 kg m^{−3}.
Khvorostyanov and Curry (2002) estimated that the terminal velocities of large ice
particles scale with the ambient density to the power of $-\mathrm{1}/\mathrm{2}$.
Therefore, we can incorporate the density dependence by multiplying
the LH74 formulas for aggregates by a factor of
$(\mathrm{0.38}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}/\mathrm{1.1}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}{)}^{\mathrm{1}/\mathrm{2}}\approx \mathrm{1.70}$
and that of H02 for
aggregates by a factor of
$(\mathrm{0.38}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}/\mathrm{0.66}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}{)}^{\mathrm{1}/\mathrm{2}}\approx \mathrm{1.32}$.
We confirmed that these corrections improve the agreement
between our model results and the formulas (see Fig. R2-1 in
the authors' response to anonymous referee no. 2).

However, at the same time, we also see several types of seemingly unrealistic ice particles,
representative examples of which are indicated by symbols in Figs. 1, 4–7:
The ice particle denoted by the circle at *t*=3000 s is a long, slowly falling hailstone.
The square at *t*=3000 s is a columnar cloud ice particle that is inconsistent with known mass–dimension relationships.
The cross at *t*=3000 s is a hailstone with an extremely low apparent density.
The triangle at *t*=5400 s is an extremely long graupel particle with a low apparent density.
In Sect. 9.1, we will investigate the causes of these odd behaviors in more detail,
but those issues could be attributed to uncertainties in ice microphysics process formulations.

Our numerical model uses three types of numerical parameters, namely the super-particle number concentration, grid size, and time steps. These parameters correspond to the resolution of aerosol/cloud/precipitation particle distribution in real space and attribute space, the spatial resolution of moist air, and temporal resolution. The numerical solution from our model approaches the true solution of time evolution equations (Eqs. 2–81) as the super-particle number approaches the number of real particles, and the grid size and time steps approach zero.

To confirm the numerical convergence of the cumulonimbus case, we conducted a series of simulations changing the numerical parameters of CTRL. These ensembles are referred to as NSP, DX, and DT (see Table 2). Our results suggest that the numerical parameters used for the CTRL case could produce accurate numerical results. In what follows, the detail of this analysis is presented. Then, we conduct a general discussion of our model's numerical convergence characteristics and computational cost.

## 8.1 NSP ensembles and super-particle number convergence

Numerical convergence regarding the initial super-particle number concentration *c*^{SP}
was investigated by varying the *c*^{SP} value of CTRL as follows: 2, 4, …, 512 per cell (see Table 2).
These cases are referred to as NSP002, NSP004, …, and NSP512, respectively. Note that NSP128 and CTRL are the same.

Figure 8 shows the accumulated
precipitation amount statistics at the end of the simulation (*t*=5400 s)
versus the initial super-particle number concentration *c*^{SP}.
The error bars indicate the mean and standard deviation calculated
from the 10 members of each NSP ensemble. The unbiased estimator was
used to calculate standard deviations. The crosses denote
the maximum and minimum values of the 10 ensemble members.

Figure 9 shows the statistics of the maximum
water path of each hydrometeor type during the simulation (i.e., the
maximum of each line in Fig. 2) versus the initial
super-particle number concentration *c*^{SP}. The error bars
indicate the mean and standard deviation from the 10
members of each NSP ensemble. The unbiased estimator was used for
calculating the standard deviations. The symbols indicate the maximum
and minimum values of each hydrometeor type in the 10 ensemble members.

Our model has two sources of fluctuation, namely atmospheric turbulence and SDM randomness.
Pseudo-random numbers are used for the Monte Carlo calculation of coalescence, riming, and aggregation,
and to initialize the super-particles. The standard deviation (i.e., fluctuation) caused by SDM randomness
decreases proportionally to the inverse of the square root of the super-particle number.
However, Figs. 8 and 9 show that the fluctuation is not sensitive to
the initial super-particle number concentration *c*^{SP}.
This indicates that fluctuations in all simulations are mostly dominated by atmospheric turbulence. One might note that the fluctuations are slightly increasing as *c*^{SP} increases.
This suggests that the super-particle number affects the turbulence characteristics; however, we leave that for further investigation in future work.

Figure 8 indicates that the accumulated precipitation amount is less sensitive to the super-particle number. However, Fig. 9 reveals
that the initial super-particle number concentration *c*^{SP} affects the maximum water path statistics.
The numerical convergence of the maximum cloud water path is noticeably slow.
This is closely related to the onset of warm rain through coalescence.
From Fig. 2, we determine that the maximum of the cloud water path coincides with
the emergence of rainwater. Therefore, a small shift of the warm rain onset time changes the maximum cloud water path;
however, it does not have a considerable impact on the overall properties of the simulated cloud.
The maximum water paths of all the other hydrometeor types do not show a significant difference
if *c*^{SP} is larger than 64 or 128 per cell (see also Table R2-1 of authors' response to anonymous referee no. 2).
When the number of super-particles was too low,
more rain droplets were produced because of an erroneous enhancement of coalescence
that reduced the amount of cloud droplets, cloud ice particles, and graupel particles.

To summarize, we conclude
that the numerical convergence regarding the super-particle number
is fairly well achieved at NSP128 (CTRL), i.e., *c*^{SP}=128 per cell.

## 8.2 DX ensembles and grid convergence

We investigated the numerical convergence with respect to the grid size by varying $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z$ of CTRL as follows: 31.25, 62.5, 125, and 250 m (see Table 2). These cases are referred to as DX/2, DXx1, DXx2, and DXx4, respectively. Note that DXx1 and CTRL are the same.

Figure 10 shows the accumulated precipitation amount statistics at the end of the simulation versus the grid size, plotted in the same way as Fig. 8, except for the difference in the horizontal axis.

Figure 11 shows the maximum water path statistics for each hydrometeor type during the simulation versus the grid size, plotted in the same way as Fig. 9, except for the difference in the horizontal axis.

The DX/2 ensemble is the highest grid resolution tested in this study, and a snapshot of the cloud from the DX/2 ensemble is shown in Fig. 12. The mixing ratios are plotted in the same manner as Fig. 1. See also Movie 6 in the video supplement.

Figure 10 shows that the accumulated precipitation amount increased from a grid size of 125 m to a grid size of 62.5 m, but no significant difference exists between the 62.5 and 31.25 m grids. Similar behavior can be observed in the maximum rainwater path in Fig. 11; however, no significant difference exists for other hydrometeor types. Comparing Fig. 12 (31.25 m) and Fig. 1 (62.5 m), the spatial structures of the clouds also look similar.

Therefore, we conclude that the numerical convergence with respect to the grid size is achieved at DXx1 (CTRL), i.e., $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z=\mathrm{62.5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$.

## 8.3 DT ensembles and time step convergence

We investigated the numerical convergence with respect to the cloud microphysics time steps by varying CTRL's cloud microphysics time steps by factors of 1∕4, 1∕2, 1, 2, 5, and 10 (see Table 2). These cases are referred to as DT/4, DT/2, DTx1, DTx2, DTx5, and DTx10, respectively. Note that DTx1 and CTRL are the same.

We found that DTx10 diverges at around *t*=1200 s
because of a numerical instability. Let us compare the results of the other five ensembles.

Figure 13 shows the statistics of the accumulated precipitation amounts at the end of the simulation versus the ratio of cloud microphysics time steps to CTRL, plotted in the same manner as Fig. 8, except for the difference in the horizontal axis.

Figure 14 shows the statistics of the maximum water path of each hydrometeor type during the simulation versus the ratio of cloud microphysics time steps to CTRL, plotted in the same manner as Fig. 9, except for the difference in the horizontal axis.

Both figures show no significant difference among the five ensembles;
therefore, we conclude that the numerical convergence with respect to the time steps
is already attained at DTx1 (CTRL), i.e.,
(Δ*t*, Δ*t*_{adv}, Δ*t*_{fz∕mlt}, Δ*t*_{collis}, Δ*t*_{cnd∕evp}, $\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}})=$ (0.4, 0.4, 0.4, 0.2, 0.1, 0.1 s).
Because DTx5 does not show any difference, time steps of 5–10 times as large could suffice.

Further discussion of numerical convergence characteristics is provided in Sect. 8.4.

## 8.4 Interpretation and computational cost

As confirmed in the preceding sections, the numerical parameters used for the CTRL ensemble (see Table 2) could produce an accurate numerical solution of the cumulonimbus case.

The CTRL ensemble's super-particle number concentration is
*c*^{SP}=128 per cell, which is comparable to various
previous studies (e.g., Andrejczuk et al., 2010; Sölch and Kärcher, 2010; Riechelmann et al., 2012; Arabas and Shima, 2013; Unterstrasser and Sölch, 2014; Unterstrasser et al., 2017; Grabowski et al., 2018; Jaruga and Pawlowska, 2018; Dziekan et al., 2019; Hoffmann et al., 2019).
Those studies reported that approximately 50–200
super-particles per cell are needed to accurately simulate clouds in two or three dimensions.
If the number of attributes is increased, we generally need more
super-particles to cover the higher-dimensional attribute space. In
this study, we used 5 attributes to represent ice particles,
which is relatively large compared to previous studies. Therefore,
achieving numerical convergence with a super-particle
number concentration as low as 128 per cell is a remarkable
result, revealing the efficacy of a particle-based cloud modeling approach. Another example of studies using many attributes is
Jaruga and Pawlowska (2018), which included 8 attributes to study
aqueous-phase oxidation of sulfur to sulfate and confirmed that
the results do not change significantly if the number concentration of
super-droplets is larger than 64 per cell.
However, the readers must be warned that the performance is sensitive to how
super-particles are initialized (Unterstrasser et al., 2017), as we discussed in Sect. 5.3.

The CTRL ensemble's grid size is $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z=\mathrm{62.5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$, which is highly dependent on the simulated cloud's energy injection scale. As we will discuss in Sect. 9.3.10, introducing SGS turbulence models should improve the grid convergence characteristics.

The time steps for cloud microphysical processes used in the CTRL ensemble
are (Δ*t*_{adv}, Δ*t*_{fz∕mlt}, Δ*t*_{collis}, Δ*t*_{cnd∕evp}, $\mathrm{\Delta}{t}_{\mathrm{dep}/\mathrm{sbl}})=$ (0.4, 0.4, 0.2, 0.1, 0.1 s).
As shown in Sect. 8.3,
time steps as large as 5–10 times could suffice.
In the following section, we discuss how those time steps are determined and whether the constraints could be relaxed.

To accurately trace the flow of moist air, Δ*t*_{adv} should be limited by the CFL condition of wind velocity.

To avoid a sudden release of latent heat, Δ*t*_{fz∕mlt} must also be restricted through the CFL condition.

Δ*t*_{collis} is the time step of coalescence, riming, and aggregation. As shown in Shima et al. (2009) and clarified in Sect. 5.5.5, Δ*t*_{collis} can be estimated
from the number concentration and size of real particles,
and that Δ*t*_{collis} does not depend on the numerical parameters such as super-particle number concentration or grid size.
To make the calculation robust to larger time steps,
a technique to allow multiple coalescence, riming, and aggregation occurrences is implemented in the SDM
(see Sect. 5.5.5); however,
this does not work properly if the collected super-particle's multiplicity is not sufficiently large.
Multiple coalescence/riming/aggregation of collector particles would not be an accurate approximation either.
These issues could be improved by introducing a recursive algorithm (Okawa, 2015),
which could allow us to use larger Δ*t*_{collis} values.

Δ*t*_{cnd∕evp} and Δ*t*_{dep∕sbl}
are determined by the phase relaxation time of supersaturation, ${\mathit{\tau}}_{\mathrm{phase}}\propto \mathrm{1}/\sum {\mathit{\xi}}_{i}{r}_{i}$ (e.g., Squires, 1952).
The timescale of CCN activation/deactivation is normally much smaller than the phase relaxation time; however,
our model is not constrained by the activation/deactivation timescale because
the condensation and evaporation equation (Eq. 8) is solved implicitly (see Sect. 5.5.3).
However, we explicitly calculate the exchange of vapor and latent heat with moist air (see Sect. 5.5.3);
thus, Δ*t*_{cnd∕evp} and Δ*t*_{dep∕sbl} must be smaller than the phase relaxation time.
Otherwise, numerical instability occurs (Árnason and Brown, 1971).
This restriction could be relaxed if we fully implicitly solve this coupled process of droplets and moist air.
Perhaps the approach described in Sect. 2.6 of Grabowski et al. (2018)
for mitigating spurious cloud-edge supersaturations could also be used for this purpose.

We used the first-order operator splitting scheme to separate the calculation (Table 1). Employing higher-order operator splitting and/or tendencies would also improve numerical convergence characteristics.

Lastly, we discuss SCALE-SDM's actual computational cost.
Calculating one realization of the CTRL case required approximately 10 h using 160 Intel Xeon E5-2650v3 CPU cores.
To compare computational cost, we also tried the two-moment bulk scheme of Seiki and Nakajima (2014) implemented on SCALE.
This took approximately 20 min, which is about 30 times faster than SDM.
As SDM's computational cost scales linearly with the number of super-particles and the number concentration of super-particles for the CTRL case was 128 per cell, it is a plausible result.
We can solve the same mathematical model using a multi-dimensional bin scheme.
Let us also estimate the bin model's computational cost.
The effective number of attributes we used for ice particles is 5,
implying that the bin space is five-dimensional. If we assume that 10–100 bins are needed for each axis,
the total number of bins becomes 10^{5}–100^{5}.
For the binary collision calculation, most bin models assess all the combinations of the bins.
In this case, the computational cost scales
with the square of the number of bins, i.e., 10^{10}–100^{10}.
However, we can reduce the cost of bin models by introducing a collision pair number reduction technique
similar to that of SDM (Sato et al., 2009).
Therefore, if we enhance the efficiency by using this algorithm,
the computational cost of bin models scales linearly with the number of bins, i.e., 10^{5}–100^{5}.
However, this is still much larger than 100, i.e., the computational cost of SDM.

In SCALE-SDM, super-particles are distributed all over the simulated domain. If we use super-particles only inside the clouds by employing, e.g., the Twomey super-droplet methodology (Grabowski et al., 2018), computational costs could be considerably reduced.

Results of the typical realization of CTRL presented in Sect. 7 show that the life cycle of a cumulonimbus was successfully simulated and the predicted mass–dimension and velocity–dimension relationships agree fairly well with the existing formulas based on laboratory measurements and observations. At the same time, as indicated by the symbols in Figs. 1 and 4–7, our model produces several types of seemingly unrealistic ice particles. Another issue is the underestimation of columnar ice particle terminal velocities. As stated in Sect. 4.1.3, we did not properly implement the formula of Böhm (1989, 1992c, 1999) in our model. Further, not all of the elementary cloud microphysics processes critical for mixed-phase clouds are incorporated in our model yet. In this section, we explore the possible improvements and further sophistication of the model.

## 9.1 Origins of odd particles

Let us determine the origins of the four types of odd particles denoted by the symbols in Figs. 1 and 4–7. Once determined, we then modify the time evolution equations to resolve three of the four issues in effect.

### 9.1.1 Long, slowly falling hailstones

The ice particle denoted by the circle at *t*=3000 s is an example of hailstones that are too long and slow-falling.
The attributes related to this ice particle's morphology are
$\mathit{\{}a,c,{\mathit{\rho}}^{\mathrm{i}},{m}^{\mathrm{rime}}/m,{n}^{\mathrm{mono}}\mathit{\}}=\mathit{\{}\mathrm{2.58}\phantom{\rule{0.125em}{0ex}}\mathrm{mm},\mathrm{36.1}\phantom{\rule{0.125em}{0ex}}\mathrm{mm},\mathrm{1.12}\times {\mathrm{10}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}},\mathrm{0.98},\mathrm{213}\mathit{\}}$.
Therefore, this ice particle is categorized as a hailstone.
However, the aspect ratio is >10, which is unrealistically long for a hailstone.
Because of the odd shape, its terminal velocity ${v}^{\mathrm{\infty}}=\mathrm{4.25}\times {\mathrm{10}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$
is also much smaller than that of typical hailstones (Fig. 7).
It is located at $(x,z)=(\mathrm{20.0}\phantom{\rule{0.125em}{0ex}}\mathrm{km},\mathrm{4.11}\phantom{\rule{0.125em}{0ex}}\mathrm{km})$, which is near the freezing level (see Fig. 1).

This odd particle was caused by a problem with the riming density formula (Eqs. 55–57).
By analyzing this particle's history, we found that it was created by only a single riming event
between a graupel particle and a similarly sized rain droplet.
We can explain the mechanism as follows: consider the riming of a quasi-spherical columnar graupel particle with a radius of 1 mm
and a rain droplet with a radius slightly smaller than 1 mm.
Assume also that the ambient temperature is slightly lower than
0 ^{∘}C. Then, from Eqs. (55)–(57),
${\mathit{\rho}}_{\mathrm{rime}}=\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$.
In other words, the apparent volume of the rimed rain droplet expands 10-fold.
Because of the filling-in model we employed for riming outcome (see Sect. 4.1.10),
the resultant ice particle became a long columnar hailstone: $(a,c)=(\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{mm},\mathrm{11}\phantom{\rule{0.125em}{0ex}}\mathrm{mm})$.

However, ${\mathit{\rho}}_{\mathrm{rime}}=\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$ must be reconsidered.
Equation (56) has a global maximum of approximately 0.95 g cm^{−3} at around *Y*=3.7 and then quickly decreases,
becoming $<\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$ at around *Y*=5.5.
Then, from Eq. (55), ${\mathit{\rho}}_{\mathrm{rime}}=\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$ for *Y*>5.5.
Consequently, considering the definition $Y:=(-{r}_{k}{v}_{\mathrm{imp}}/{T}_{j}^{\mathrm{sfc}})/\left(\mathrm{\mu}\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}/{}^{\circ}\mathrm{C}\right)$,
${\mathit{\rho}}_{\mathrm{rime}}=\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$ frequently happens near the freezing level. For example,
*Y*=1000 for *r*_{k}=1 mm, ${v}_{\mathrm{imp}}=\mathrm{1}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, and ${T}_{j}^{\mathrm{sfc}}=-\mathrm{1}$ ^{∘}C.
However, *ρ*_{rime} would be much larger and even closer to ${\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}$ in such a situation in reality
because the rimed droplet freezes slowly. Therefore, we argue that
Eq. (56) is valid only up to *Y*=3.5 and levels off after that.
This correction can be made by replacing the *Y* value in
Eq. (56) (but not in Eq. 57) with

In Sect. 9.1.5, we will confirm that this correction eliminates those long hailstones (Figs. 15–18).

Additionally, the same problem occurs if a quasi-spherical planar graupel particle and a slightly smaller rain droplet collide and rime near the freezing level. However, it is less evident than with the previous case because the equatorial radius grows as the square root of the volume (Eq. 53). Regardless, this problem can also be addressed using the above correction.

### 9.1.2 Columns with steep mass–dimension relationship

The square at *t*=3000 s indicates another odd particle.
If we look around the square in Figs. 4 and 5,
we see that this particle belongs to a population of columnar cloud ice particles
that have a steeper mass–dimension relationship than observed.
The attributes related to this cloud ice particle's morphology are
$\mathit{\{}a,c,{\mathit{\rho}}^{\mathrm{i}},{m}^{\mathrm{rime}}/m,{n}^{\mathrm{mono}}\mathit{\}}$
$=\mathit{\{}\mathrm{24.9}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m},$ 138.8 µm, $\mathrm{269.7}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}},$
0.29, 1*}*. Its terminal velocity is
${v}^{\mathrm{\infty}}=\mathrm{1.46}\times {\mathrm{10}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and it is located at $(x,z)=(\mathrm{20.3}\phantom{\rule{0.125em}{0ex}}\mathrm{km},\mathrm{11.0}\phantom{\rule{0.125em}{0ex}}\mathrm{km})$ (see Fig. 1).

As in the previous case, we found that a single riming event
between a cloud ice particle and a cloud droplet followed by depositional growth created this columnar ice particle type.
We can explain the mechanism as follows: consider a quasi-spherical columnar ice particle with a radius of 10 µm
and a supercooled droplet with a radius slightly smaller than 10 µm.
Assuming an impact velocity of ${\mathrm{10}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$
and ambient temperature of −10 ^{∘}C,
then, from Eqs. (55)–(57),
$Y={\mathrm{10}}^{-\mathrm{2}}$ and ${\mathit{\rho}}_{\mathrm{rime}}=\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}\mathrm{g}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{-\mathrm{3}}$.
In other words, the apparent volume of the rimed droplet expands 10-fold
and creates a columnar graupel particle: $(a,c)=(\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m},\mathrm{110}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m})$
because of our riming outcome model's filling-in assumption.
Then, through subsequent depositional growth, this columnar graupel particle turns back into a columnar cloud ice particle.

Contrary to the previous case, the low riming density is reasonable. Instead, we must reconsider the filling-in model. We assumed that the ice particle's maximum dimension is preserved. However, this is not realistic for riming between an ice particle and a similarly sized droplet, as our thought experiment revealed. Generalizing the idea, we consider that the frozen droplet's diameter is preserved if the diameter is larger than the ice particle's maximum dimension. That is, we propose to replace Eq. (50) with

and Eq. (54) with

In Sect. 9.1.5, we will confirm that those columns that follow an extremely steep mass–dimension relationship can be eliminated using this correction (Figs. 15–18).

### 9.1.3 Low-density hailstones

The cross at *t*=3000 s represents a hailstone with an extremely low apparent density.
The attributes related to this hailstone's morphology are
$\mathit{\{}a,c,{\mathit{\rho}}^{\mathrm{i}},{m}^{\mathrm{rime}}/m,{n}^{\mathrm{mono}}\mathit{\}}$
$=\mathit{\{}\mathrm{12.6}\phantom{\rule{0.125em}{0ex}}\mathrm{mm},$ 15.0 mm, $\mathrm{10.7}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}},$
0.85, 1 585 116*}*. Its terminal velocity is
${v}^{\mathrm{\infty}}=\mathrm{4.31}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and it is located at $(x,z)=(\mathrm{10.6}\phantom{\rule{0.125em}{0ex}}\mathrm{km},\mathrm{11.5}\phantom{\rule{0.125em}{0ex}}\mathrm{km})$ (see Fig. 1).
What is unusual here is the very low apparent density ${\mathit{\rho}}^{\mathrm{i}}=\mathrm{10.7}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}$.
This particle is composed of many monomers *n*^{mono}=1 585 116,
and we set the limiting value of aggregate density in Eq. (62) to
${\mathit{\rho}}_{\mathrm{crt}}^{\mathrm{i}}=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}}$.
Thus, we can conclude that this hailstone is created by the repeated aggregation between graupel particles.

Lump graupel particles with apparent densities as low as 50 kg m^{−3} were reported in Locatelli and Hobbs (1974).
Therefore, a hailstone with an apparent density of
10 kg m^{−3} is not extremely unrealistic.
However, our aggregation model is crude.
Following Morrison and Grabowski (2010), we assumed that
collision–aggregation collection efficiency is a fixed constant of *E*_{agg}=0.1
regardless of morphology or temperature.
Therefore, this could cause the accumulation of graupel particles
near the limiting value ${\mathit{\rho}}_{\mathrm{crt}}^{\mathrm{i}}$ in Fig. 6.
There should be further detailed investigation to assess
our aggregation model's applicability to graupel particles.
See also Sect. 9.3.8, which provides a discussion
to refine our aggregation model.

### 9.1.4 Long graupel particles

The triangle at *t*=5400 s is an extremely long graupel particle with a low apparent density.
The attributes related to this cloud ice particle's morphology are
$\mathit{\{}a,c,{\mathit{\rho}}^{\mathrm{i}},{m}^{\mathrm{rime}}/m,{n}^{\mathrm{mono}}\mathit{\}}$
$=\mathit{\{}\mathrm{31.9}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m},$ 438.2 µm, $\mathrm{19.4}\phantom{\rule{0.125em}{0ex}}\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{3}},$
0.53, 26 679*}*. Its terminal velocity is
${v}^{\mathrm{\infty}}=\mathrm{1.02}\times {\mathrm{10}}^{-\mathrm{3}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, it is located at $(x,z)=(\mathrm{29.7}\phantom{\rule{0.125em}{0ex}}\mathrm{km},\mathrm{6.4}\phantom{\rule{0.125em}{0ex}}\mathrm{km})$ (see Fig. 1), and
the ambient temperature is $T=-\mathrm{14.4}$ ^{∘}C.

The particle is created by a sublimation of a graupel particle.
The inherent growth ratio Γ(*T*) proposed by Chen and Lamb (1994a)
was used to calculate the deposition and sublimation process as described in Sect. 4.1.7.
Γ(*T*)<1 if *T* is in the range of approximately
[−20 ^{∘}C,−10 ^{∘}C] and
[−5 ^{∘}C,0 ^{∘}C];
therefore, in this temperature range, ice particles grow to become planar through deposition
and shrink to become columnar by sublimation.

However, Γ(*T*) was derived from measurements of depositional growth; hence, it is questionable whether it is applicable for sublimation.
According to Harrington et al. (2019) and references therein,
Γ(*T*) should be considered as unity for sublimation,

thus, the aspect ratios of ice particles are preserved during sublimation.

In Sect. 9.1.5, we will confirm that those long graupel particles can be eliminated using this correction (Figs. 15–18).

### 9.1.5 Results after corrections

In the preceding sections, we proposed three corrections to the time evolution equations (Eqs. 107–110) to avoid the creation of ice particles with unrealistic morphologies.

We incorporated the proposed corrections into our model to create a new revision, SCALE-SDM 0.2.5-2.2.1. To assess the validity of these corrections, we conducted the same simulations as the typical realization of CTRL using the new model. By comparing these results (Figs. 15–18) to the original results (Figs. 4–7), we confirm that the three types of odd ice particles no longer exist, as we intended. See also Movies 7–11 in the video supplement. Note that we left the issue of low-density hailstones for future studies. These corrections have little effect on the overall cloud properties, i.e., spatial structure (Movie 7 in the video supplement), the time evolution of the water path (Fig. 19), and the accumulated precipitation amount (Fig. 20).

## 9.2 Fix of underestimated terminal velocities of columnar ice particles

As explained in Sect. 4.1.3, our model did not properly implement the ice particle terminal velocity formula of Böhm (1989, 1992c, 1999). In this section, we fix the problem and assess its impact on this study.

Noting that the area ratio *q*_{i}≤1 always holds in our model, Böhm's formula
${v}_{\text{B\xf6hm}}^{\mathrm{\infty}}({m}_{i},{\mathit{\varphi}}_{i},{d}_{i},{q}_{i};{\mathit{\rho}}_{i},{T}_{i})$
can be summarized as follows:

In our model (SCALE-SDM 0.2.5-2.2.0/2.2.1), we assumed that the characteristic length *d*_{i} is given by the maximum dimension ${D}_{i}=\mathrm{2}max({a}_{i},{c}_{i})$,
and the area ratio *q*_{i} is given by the area ratio regarding the circumcircle ${q}_{i}^{\mathrm{cc}}={A}_{i}/{A}_{i}^{\mathrm{cc}}$ (Eq. 4).
However, in Böhm's theory, they are defined by

i.e., for columnar particles, minor axis is used for *d*_{i},
and the area ratio regarding circumscribed ellipse is used for *q*_{i}.
Figure 1 in Böhm (1989) suggests ${q}_{i}={q}_{i}^{\mathrm{ce}}$.
It is not clearly specified, but from the second equality of Eq. 17 in Böhm (1992c), we can confirm that *d*_{i}=2*a*_{i}.

For planar ice particles (*ϕ*_{i}<1), ${v}_{\text{B\xf6hm}}^{\mathrm{\infty}}({d}_{i}=\mathrm{2}{a}_{i},{q}_{i}={q}_{i}^{\mathrm{ce}})$ and
${v}_{\text{B\xf6hm}}^{\mathrm{\infty}}({d}_{i}={D}_{i},{q}_{i}={q}_{i}^{\mathrm{cc}})$ yield the same results, because 2*a*_{i}=*D*_{i} and ${q}_{i}^{\mathrm{ce}}={q}_{i}^{\mathrm{cc}}$ hold.
However, for columnar ice particles (*ϕ*_{i}>1), ${v}_{\text{B\xf6hm}}^{\mathrm{\infty}}({D}_{i},{q}_{i}^{\mathrm{cc}})$ always underestimates the fall velocity.
From Eqs. (111)–(121), we can derive
${v}_{\text{B\xf6hm}}^{\mathrm{\infty}}(\mathrm{2}{a}_{i},{q}_{i}^{\mathrm{ce}})/{v}_{\text{B\xf6hm}}^{\mathrm{\infty}}({D}_{i},{q}_{i}^{\mathrm{cc}})={\mathit{\varphi}}_{i}^{\mathrm{3}/\mathrm{4}}$ for *X*≪1,
and ${v}_{\text{B\xf6hm}}^{\mathrm{\infty}}(\mathrm{2}{a}_{i},{q}_{i}^{\mathrm{ce}})/{v}_{\text{B\xf6hm}}^{\mathrm{\infty}}({D}_{i},{q}_{i}^{\mathrm{cc}})={\mathit{\varphi}}_{i}^{\mathrm{7}/\mathrm{8}}$ for *X*≫1.
Therefore, if *ϕ*_{i}=2, the ratio ${v}_{\text{B\xf6hm}}^{\mathrm{\infty}}(\mathrm{2}{a}_{i},{q}_{i}^{\mathrm{ce}})/{v}_{\text{B\xf6hm}}^{\mathrm{\infty}}({D}_{i},{q}_{i}^{\mathrm{cc}})$
is in the rage of 1.68–1.83. If *ϕ*_{i}=10, the range is 5.62–7.50, and if *ϕ*_{i}=20, it is 9.46–13.75.
We also confirmed that Böhm's original definition ${v}_{\text{B\xf6hm}}^{\mathrm{\infty}}(\mathrm{2}{a}_{i},{q}_{i}^{\mathrm{ce}})$
agrees well with the formulas of Westbrook et al. (2008), and Heymsfield and Westbrook (2010)
(see also Fig. R2-2 in the authors' response to anonymous referee no. 2).

Therefore, the correction (Eq. 122) generally increases
the fall speed of columnar ice particles, and the increase factor is
larger for longer particles. Then, through the ventilation effects
(Eqs. 11 and 15), the diffusional growth of columnar ice particles is
enhanced. Due to this mechanism, we observed the creation of extremely long
ice particles with aspect ratio *ϕ*_{i}>100 if we incorporate the
correction (Eq. 122) to SCALE-SDM 0.2.5-2.2.1. However, this is
unrealistic. The maximum aspect ratio reported is approximately 30 in
Auer and Veal (1970) (Fig. 12 therein) and 15.77 in Um et al. (2015).
In nature, such an extreme-shaped ice particle would be
shattered spontaneously or by collision. However, for the moment, we fix
this issue in an ad hoc way. We do not allow an ice particle to grow
by diffusion slenderer than *ϕ*_{i}=40 by imposing a limiter to the
effective inherent growth ratio Γ^{*} as follows:

We incorporated the corrections (Eqs. 122 and 123) into SCALE-SDM 0.2.5-2.2.1 to create a revision, SCALE-SDM 0.2.5-2.2.2. To assess the impact of these corrections, we conducted the same simulation as the typical realization of CTRL using the new model. We observed that the precipitation was developed a few minutes faster, but the total precipitation amount was almost the same as the previous versions (Fig. 20). Figure 19 compares the time evolution of water paths. Here, a noticeable decrease in the graupel water path can be observed, which is attributed to the increased fall speed of columnar graupel particles (i.e., densely rimed columns). This, in turn, increased the rainwater path. The time evolution of other hydrometeor water paths (cloud, cloud ice, and snow) was almost unchanged. Ice particle morphology distributions resemble closely to the previous results, except for the vanishing of cloud ice particles with relatively slow terminal velocities (Figs. R2-5–R2-8 in the authors' response to anonymous referee no. 2. See also Movies 13–16 in the video supplement). The corrections also do not alter the spatial structure of the cloud (Movie 12 in the video supplement).

## 9.3 Further sophistication of the model

Our model is based on a kinetic description, i.e., individual dynamics of particles and their stochastic collisions. However, a quantitative understanding of mixed-phase cloud microphysics is a long-standing meteorological issue, and a kinetic description of mixed-phase cloud microphysics has not been established. Further, our model does not incorporate several elementary processes that are critical for mixed-phase clouds. In this section, we explore the possibilities of further refining and sophisticating our model. Readers can also refer to Chen and Lamb (1994a, b), Misumi et al. (2010), Hashino and Tripoli (2007, 2008, 2011a, b), Jensen and Harrington (2015), Sölch and Kärcher (2010), Brdar and Seifert (2018), and Seifert et al. (2019), as these are modeling studies closely relevant to our study.

### 9.3.1 Ice nucleation pathways

There are various ice nucleation pathways (e.g., Kanji et al., 2017); however, in this study, we only considered condensation/immersion freezing and homogeneous freezing, as these are the dominant mechanisms in mixed-phase clouds.

Based on the singular hypothesis (Levine, 1950), we considered
that each insoluble particle has its own freezing temperature *T*^{fz}
that can be determined by INAS formulas.
In the model evaluation experiments, we assumed that ice nuclei consist of mineral dust and used the INAS formula of Niemand et al. (2012).
Formulas from Wex et al. (2015) and
Ullrich et al. (2017) can be used for biogenic substances and soot, respectively.

The singular hypothesis ignores the time dependence of ice nucleation; thus,
we assumed that particles initiate freezing immediately after
the temperature drops below *T*^{fz} and the ambient air becomes saturated over liquid water.
However, the time dependence of ice nucleation could be critical for clouds with long lifetimes,
known as the “stochastic hypothesis”. The soccer ball model of Niedermeier et al. (2011, 2014, 2015),
which is based on classical nucleation theory, could be used
to incorporate time dependence. Then, instead of the freezing temperature *T*^{fz},
the contact angle of the surface site *θ* must be treated as an attribute.

Note that our requirement that the ambient water vapor must be supersaturated over liquid water would be too restrictive for immersion freezing. Even under an unsaturated condition, it is reasonable to allow immersion freezing if the droplet is sufficiently large, for instance, larger than 1 µm in radius.

To express homogeneous freezing, we assigned a fixed freezing temperature of ${T}^{\mathrm{fz}}=-\mathrm{38}$ ^{∘}C
to all the IN inactive particles and ignored the time dependence of ice nucleation.
However, this is not appropriate for the homogeneous freezing of deliquescent aerosol particles
because homogeneous ice nucleation is suppressed when solute concentration increases.
Additionally, the time dependence of ice nucleation could be also critical
because the probability that a droplet freezing homogeneously is proportional to the liquid water volume.
These effects can all be incorporated using the model of Koop et al. (2000).

Condensation/immersion freezing of deliquescent IN particles can also be incorporated
by considering the depression of the freezing temperature *T*^{fz} by the solute
(see Wex et al., 2014, and references therein). Alternatively, a model based on
classical nucleation theory proposed by Khvorostyanov and Curry (2004, 2005)
can be used to incorporate time dependence.

The formation of ice directly from the vapor phase onto an IN particle is known as
deposition freezing. This can be observed at $<-\mathrm{25}$ ^{∘}C
and in air that is below water saturation. Marcolli (2014) suggested that the phenomena conventionally known as deposition freezing could be reinterpreted as pore condensation and freezing.
We can use the temperature-dependent and saturation-ratio-dependent INAS formula
proposed by Steinke et al. (2015) to incorporate this process.
Here, INAS density *n*_{S} is a function of *x*_{therm}, and
*x*_{therm} is a function of temperature *T* and saturation ratio over ice *S*^{i}.
We can assign *x*_{therm,i} to each particle as an attribute.
We consider that freezing occurs when ${x}_{\mathrm{therm}}(T,{S}^{\mathrm{i}})>{x}_{\mathrm{therm},i}$.

A crude model of pre-activation is incorporated in our model by inhibiting complete sublimation (see Eq. 14 and the explanation that follows). Pre-activation denotes “the capability of particles or materials to nucleate ice at lower relative humidities or higher temperatures compared to their intrinsic ice nucleation efficiency after having experienced an ice nucleation event or low temperature before” (Marcolli, 2017). Intensive sophistication based on laboratory studies is required; however, particle-based models are suitable for exploring the atmospheric relevance of pre-activation. Conversely, one might want to switch off pre-activation in our model, which is possible by resetting the particles as deliquescent aerosol particles when complete sublimation occurs.

Contact freezing is another ice nucleation mechanism in which
solid particles can initiate freezing
upon contacting the surface of a supercooled droplet.
Contact freezing occurs at temperatures greater than that of
the same particle immersed in a droplet (e.g., Shaw et al., 2005); therefore,
it might also be relevant to mixed-phase clouds.
To explain the scavenging of aerosol particles by droplets,
we must consider Brownian diffusion and phoretic forces.
This process can be incorporated into our model
by introducing the collision–coalescence kernels
detailed in Sect. 17.4.2 of Pruppacher and Klett (1997).
Then, based on the results of Shaw et al. (2005),
as suggested by Will H. Cantrell (personal communication, 2017),
contact freezing could be expressed by
increasing the particle's *T*^{fz}
by 4.5 ^{∘}C in each single particle–droplet collision event.
Another possibility is using laboratory data from Niehaus et al. (2014), who
measured the freezing efficiency of various insoluble particles.
This can be interpreted as the probability
that each particle–droplet collision results in a freezing event.

It is also known that the evaporation of a droplet could lead to inside-out contact freezing (e.g., Durant and Shaw, 2005); however, there are still substantial uncertainties.

### 9.3.2 Onset of melting

We assumed that ice particles start melting immediately after
the ambient temperature reaches > 0 ^{∘}C.
However, evaporative cooling delays the melting onset.
For example, at a relative humidity of 50 %,
melting starts at +4 ^{∘}C.
We can incorporate this effect by considering ice particle surface temperatures,
as discussed in Rasmussen and Pruppacher (1982).

### 9.3.3 Partially frozen/melted particles

After the onset of freezing or melting, we assumed that complete freezing/melting occurs instantaneously.

However, as shown in Murray and List (1972), the freezing time of millimeter-sized droplets could be of the order of 100 s. We can explicitly incorporate this process using the time evolution equation summarized in Sect. 16.1.4 of Pruppacher and Klett (1997), which is derived from a quasi-steady assumption of vapor and thermal diffusion around a partially frozen droplet.

We also assumed that rimed supercooled droplets freeze instantaneously; however, wet growth of graupel particles is critical to accurately predict hailstone formation. We can use the model from Rasmussen and Heymsfield (1987) to incorporate the wet growth process.

Depending on the relative humidity and warming rate, the melting time of spherical ice particles with radii of approximately 300–400 µm ranges between 20 and 70 s (Rasmussen and Pruppacher, 1982). A large hailstone could escape complete melting and reach the ground. The shedding of droplets could also occur if a partially melted hailstone contains excess meltwater, which could affect the raindrop size distribution below the cloud. Partially melted snow aggregates could create a layer of stronger radar reflectivity below the melting level, known as the “bright band”. We can explicitly incorporate these processes using the model summarized in Phillips et al. (2007).

Additionally, to complete the model, all other time evolution equations must be extended to make them compatible with partially frozen/melted particles, which would require some effort.

### 9.3.4 Condensation and evaporation

In SCALE-SDM, we assumed that water vapor's diffusivity
in air and moist air's thermal conductivity in Eq. (9) are fixed constants,
${D}_{\mathrm{v}}=\mathrm{2.52}\times {\mathrm{10}}^{-\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and
$k=\mathrm{2.55}\times {\mathrm{10}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{J}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{K}}^{-\mathrm{1}}$,
which are the values for *T*=20 ^{∘}C and *p*=1000 hPa.
However, this approximation is erroneous, particularly because diffusivity *D*_{v} is
inversely proportional to pressure.
In the case of the initial profile we used for model evaluation,
$T=-\mathrm{44}$ ^{∘}C and *p*=250 hPa at *z*=10 km.
Thus, ${D}_{\mathrm{v}}=\mathrm{6.08}\times {\mathrm{10}}^{-\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$,
which is about 2.4 times larger than we assumed.
The temperature and pressure dependence of water vapor's diffusivity in air *D*_{v},
and the temperature dependence of moist air's thermal conductivity *k*
must be considered. The formulas summarized in Sect. 13.1 of
Pruppacher and Klett (1997) can be used.

We considered the ventilation effect for deposition and sublimation but not for condensation and evaporation, even though it also enhances the growth and evaporation of larger droplets. We can include this effect by using the model described in Sect. 13.2.3 of Pruppacher and Klett (1997).

For cloud droplets, we must also consider kinetic correction to
*D*_{v} and *k*. See, e.g.,
Sect. 13.1 of Pruppacher and Klett (1997) and Kogan (1991).

In our model, aerosol particle hygroscopicity is expressed by
Raoult's law with the van 't Hoff factor *i* (Low, 1969); however,
using the kappa parameterization of Petters and Kreidenweis (2007) would be more convenient.

### 9.3.5 Deposition and sublimation

There are many issues around Γ(*T*),
which represents the primary growth habit of ice crystals.
Considering the amount of data used for the fitting,
the proposed shape of Γ(*T*) is subject to large uncertainties (see Fig. 3 of Chen and Lamb, 1994a).
The applicable range is also unclear. We set Γ(*T*)=1 for small ice crystals
*D*<10 µm. As discussed in Sect. 9.1.4,
Γ(*T*)=1 should be used for sublimation (Harrington et al., 2019, and references therein).
We might need to use some other form of Γ(*T*) for graupel particles and snow aggregates.
Connolly et al. (2012) had to adjust Γ(*T*) somewhat arbitrarily to obtain a better agreement.

Further, as shown by Kumai (1982) and Bailey and Hallett (2004),
at $T<-\mathrm{20}$ ^{∘}C, both plates and columns can be
created at the same temperature depending on the saturation ratio over ice *S*^{i},
and polycrystals can also be created. Therefore, for $T<-\mathrm{20}$ ^{∘}C,
Γ might better be considered a function of both *T* and *S*^{i},
and formation of polycrystals must be somehow incorporated into our model.
We can employ the mathematical model from Hashino and Tripoli (2008),
which extends Chen and Lamb (1994a)'s model
to describe these behaviors.

Harrington et al. (2019) reformulated the model from Chen and Lamb (1994a), and
their model does not rely on Γ(*T*), predicting the
aspect ratio evolution using the “facet-based hypothesis”.
The model is as good as Chen and Lamb's original
model at liquid saturation, and further, it can be applied to a wider range of
environmental conditions, such as low supersaturation and low
pressure. However, it is still unclear how well the model would
work for polycrystals or irregular ice particles.

We used Chen and Lamb (1994a)'s deposition density formula;
however, as discussed in Jensen and Harrington (2015),
their formula does not capture the wind tunnel data of Takahashi et al. (1991) very well.
Instead, Jensen and Harrington (2015) proposed a simple formula:
${\mathit{\rho}}_{\mathrm{dep}}={\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}\mathrm{\Gamma}\left(T\right)$ for Γ<1;
${\mathit{\rho}}_{\mathrm{dep}}={\mathit{\rho}}_{\mathrm{true}}^{\mathrm{i}}/\mathrm{\Gamma}\left(T\right)$ for Γ>1.
Their idea to relate deposition density *ρ*_{dep} to axis growth ratio
is plausible, but its dependence on *S*^{i} is lost.
Because *ρ*_{dep} accounts for the secondary growth habit,
dependence on *S*^{i} must be reconsidered.

In our model, each ice particle is approximated by a porous spheroid $(a,c,{\mathit{\rho}}^{\mathrm{i}})$.
We used spheroid capacitance *C*(*a*,*c*) to evaluate *C* in Eq. (11).
However, the spheroid (*a*,*c*) represents the ice particle's spatial extent,
and it might have a more detailed internal structure,
which is represented by the apparent density *ρ*^{i}.
The actual ice particle capacitance also depends on the internal structure.
Westbrook et al. (2008) accurately calculated the capacitance of realistic ice particles
by directly simulating the trajectories of diffusing water molecules.
Thus, we can use their formulas to refine our model's accuracy.
For example, they showed that the capacitance of snow aggregates can be
approximated by $C=D/\mathrm{4}$, which is half that of a sphere.

As with condensation and evaporation, we assumed that water vapor's diffusivity
in air *D*_{v} and moist air's thermal conductivity *k* in Eq. (12) are fixed constants,
but this must be revised.

Demange et al. (2017) constructed a sophisticated phase field model for ice crystal growth that successfully reproduced the formation of diverse ice crystal shapes. This model could help us construct a more accurate kinetic description of the deposition and sublimation processes.

### 9.3.6 Coalescence

For the collision efficiency of collision–coalescence ${E}_{\mathrm{coal}}^{\mathrm{collis}}$,
we used a modified table of Hall (1980) proposed in Seeßelberg et al. (1996) and Bott (1998).
However, the table of Pinsky et al. (2001) is more comprehensive and reliable.
It is based on numerical results but supported by the laboratory experiments of Vohl et al. (2007).
Another option is to use the formula of Böhm (1992b, 1999, 2004).
It is interesting to note that Böhm's formula (1992b, 1999) predicts that the collision–coalescence kernel *K*_{coal}
does not vanish for equal size droplets due to wake capture effect,
but caution must be taken because his theory has an error (Böhm, 2004).

We assumed that the coalescence efficiency is unity, ${E}_{\mathrm{coal}}^{\mathrm{coal}}=\mathrm{1}$, for simplicity; however, it can be much smaller than 1 for large droplets. Straub et al. (2010) proposed a simple formula based on their numerical results. We can also use the formula of Seifert et al. (2005), which compiles the formulas of Low and List (1982) and Beard and Ochs (1995).

### 9.3.7 Riming

For the collection efficiency of collision–riming *E*_{rime},
when a large spherical ice particle collects a supercooled droplet,
we used the formula from Beard and Grover (1974) with a mixed Froude number
(Eqs. 40 and 41).
von Blohn et al. (2009) demonstrated that the formula underestimates the efficiency if the spherical ice particle is large,
but Eq. (11) in their paper seems to be incorrect, and thus we did not consider this.

When a large droplet collects ice particles, we used the original formula from Beard and Grover (1974), approximating the ice particle as spherical. To consider the ice particle shape, we can use the formulas from Lew and Pruppacher (1983) for a large droplet collecting small columns, and Lew et al. (1985) for a large droplet collecting small planar crystals.

Beard and Grover (1974)'s formula is valid only for *p*<0.1,
where *p* is the size ratio of the collector ice/droplet and collected droplet/ice.
We forcibly applied the formula beyond this range, which
increases the collection efficiency of riming between small similar size droplets and ice particles,
as ${E}_{\mathrm{BG}\mathrm{74}}(p,{N}_{\mathrm{Re}},{N}_{\mathrm{St}})\approx {p}^{\mathrm{2}}/(\mathrm{1}+{p}^{\mathrm{2}})$ for *N*_{St}≪1.
This must be corrected.

When an ice particle collects a droplet, we employed the filling-in model and preserved the ice particle's maximum dimension. However, if the collector is a snow aggregate, we should use the similarity model proposed by Seifert et al. (2019). Unrimed/rimed snow aggregates have fractal structures, and Seifert et al. (2019) found a universal self-similar relation in snow aggregate growth through riming. The similarity model considers the maximum dimension's increase during the early stages of riming, which could lead to more rapid ice particle growth due to riming.

### 9.3.8 Aggregation

We assumed that collision–aggregation's collection efficiency is given by a constant *E*_{agg}=0.1,
following Morrison and Grabowski (2010), but this is a simplification.
*E*_{agg} should be larger for large particles because of the interlocking mechanism, and near water-saturated conditions.
*E*_{agg} can be decomposed into ${E}_{\mathrm{agg}}={E}_{\mathrm{agg}}^{\mathrm{collis}}{E}_{\mathrm{agg}}^{\mathrm{stick}}$, where ${E}_{\mathrm{agg}}^{\mathrm{collis}}$ is collision efficiency and
${E}_{\mathrm{agg}}^{\mathrm{stick}}$ is sticking efficiency.
For ${E}_{\mathrm{agg}}^{\mathrm{collis}}$, we can use the formula of Böhm (1989, 1992a, b, c, 1994, 1999, 2004).
For ${E}_{\mathrm{agg}}^{\mathrm{stick}}$