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

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.

1 Introduction

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 Attributes of atmospheric particles

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 a. 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.

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³.

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}-\exp [-A^{\mathrm{insol}}{n}_{\mathrm{S}}(T\left)\right].$ The probability density function of T^fz then becomes

$$\begin{array}{}\text{(1)}& p\left(T\right)=-{\displaystyle \frac{\text{d}P(T^{\mathrm{fz}}>T)}{\text{d}T}}=-A^{\mathrm{insol}}{\displaystyle \frac{\text{d}{n}_{\mathrm{S}}}{\text{d}T}}{e}^{-A^{\mathrm{insol}}{n}_{\mathrm{S}}}.\end{array}$$

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 ρⁱ. 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=ρⁱV. The two radii a and c represent the ice particle's spatial extent and ρⁱ 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, ρⁱ becomes smaller than the ice crystal's true density ${\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}}\approx \mathrm{916.8}$ kg m⁻³.

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.

3 Variables for moist air

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₀=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{\}}$.

4 Time evolution equations of mixed-phase clouds

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

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\text{d}}{\text{d}t}}\left(m_i{\mathit{v}}_i\right)={\mathit{F}}_i^{\mathrm{drg}}-m_{i}g\widehat{\mathit{z}},{\displaystyle \frac{\text{d}{\mathit{x}}_i}{\text{d}t}}={\mathit{v}}_i,\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\mathit{v}}_i={\mathit{U}}_i-\widehat{\mathit{z}}{v}_i^{\mathrm{\infty }},{\displaystyle \frac{\text{d}{\mathit{x}}_i}{\text{d}t}}={\mathit{v}}_i,\end{array}$$

where U_i:=U(x_i) is the ith 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⁻⁵ 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}\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öhm}}^{\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 2a_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:

$$\begin{array}{}\text{(4)}& d_i=D_i:=\mathrm{2}max(a_i,c_i),q_i=q_i^{\mathrm{cc}}:=A_i/A_i^{\mathrm{cc}},\end{array}$$

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:

$$\begin{array}{}\text{(5)}& A_i=A_i^{\mathrm{ce}}{\left({\displaystyle \frac{{\mathit{\rho }}_i^{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}}}}\right)}^{\mathit{\kappa }\left({\mathit{\varphi }}_i\right)}.\end{array}$$

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

$$\begin{array}{}\text{(6)}& \mathit{\kappa }\left({\mathit{\varphi }}_i\right)=\exp (-{\mathit{\varphi }}_i).\end{array}$$

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 ith 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 ith 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

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\text{d}{m}_i}{\text{d}t}}=\mathrm{4}\mathit{\pi }{r}_i{D}_{\mathrm{v}}({\mathit{\rho }}_{\mathrm{v}i}-{\mathit{\rho }}_{\mathrm{v}i}^{\mathrm{sfc}}).\end{array}$$

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

$$\begin{array}{}\text{(8)}& r_i{\displaystyle \frac{\text{d}{r}_i}{\text{d}t}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}^{\mathrm{w}}(F_{\mathrm{k}}^{\mathrm{w}}+F_{\mathrm{d}}^{\mathrm{w}})}}\left\{S_i^{\mathrm{w}}-{\displaystyle \frac{e_{\mathrm{s}i}^{\mathrm{w},\mathrm{eff}}}{e_{\mathrm{s}}^{\mathrm{w}}\left(T_i\right)}}\right\},\end{array}$$

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

$$\begin{array}{}\text{(9)}& F_{\mathrm{k}}^{\mathrm{w}}=\left({\displaystyle \frac{L_{\mathrm{v}}}{R_{\mathrm{v}}{T}_i}}-\mathrm{1}\right){\displaystyle \frac{L_{\mathrm{v}}}{k{T}_i}},F_{\mathrm{d}}^{\mathrm{w}}={\displaystyle \frac{R_{\mathrm{v}}{T}_i}{D_{\mathrm{v}}{e}_{\mathrm{s}}^{\mathrm{w}}\left(T_i\right)}},\end{array}$$

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 ith 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

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{e_{\mathrm{s}i}^{\mathrm{w},\mathrm{eff}}}{e_{\mathrm{s}}^{\mathrm{w}}\left(T_i\right)}}=\mathrm{1}+{\displaystyle \frac{a\left(T_i\right)}{r_i}}-{\displaystyle \frac{b\left(\left\{m_{\mathit{\alpha }i}^{\mathrm{sol}}\right\}\right)}{r_i^{\mathrm{3}}}},\end{array}$$

where $a\approx \mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{cm}\mathrm{K}/T_i$, $b\approx \mathrm{4.3}{\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):

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\text{d}{m}_i}{\text{d}t}}=\mathrm{4}\mathit{\pi }C{D}_{\mathrm{v}}({\mathit{\rho }}_{\mathrm{v}i}-{\mathit{\rho }}_{\mathrm{v}i}^{\mathrm{sfc}})\stackrel{\mathrm{‾}}{f_{\mathrm{vnt}}}=\mathrm{4}\mathit{\pi }C{\displaystyle \frac{S_i^{\mathrm{i}}-\mathrm{1}}{F_{\mathrm{k}}^{\mathrm{i}}+F_{\mathrm{d}}^{\mathrm{i}}}}\stackrel{\mathrm{‾}}{f_{\mathrm{vnt}}},\end{array}$$

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,

$$\begin{array}{}\text{(12)}& F_{\mathrm{k}}^{\mathrm{i}}=\left({\displaystyle \frac{L_{\mathrm{s}}}{R_{\mathrm{v}}{T}_i}}-\mathrm{1}\right){\displaystyle \frac{L_{\mathrm{s}}}{k{T}_i}},F_{\mathrm{d}}^{\mathrm{i}}={\displaystyle \frac{R_{\mathrm{v}}{T}_i}{D_{\mathrm{v}}{e}_{\mathrm{s}}^{\mathrm{i}}\left(T_i\right)}},\end{array}$$

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{‾}}{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{‾}}{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{‾}}{f_{\mathrm{vnt}}}$ could be described by

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{‾}}{f_{\mathrm{vnt}}}=b_{\mathrm{1}}+b_{\mathrm{2}}{X}^{\mathit{\gamma }},\end{array}$$

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 dm_i as follows:

$$\begin{array}{}\text{(14)}& \text{d}{m}_i=max(\text{d}{m}_i,m_{\min }^{\mathrm{i}}-m_i),\end{array}$$

where $m_{\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_{\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:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\text{d}{c}_i}{\text{d}{a}_i}}=\mathrm{\Gamma }\left(T_i\right)f_{\mathrm{vnt}}{\displaystyle \frac{c_i}{a_i}}=:{\mathrm{\Gamma }}^{*}{\displaystyle \frac{c_i}{a_i}},\end{array}$$

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

$$\begin{array}{}\text{(16)}& \mathrm{\Gamma }\left(T\right)=\mathrm{\Gamma }(-\mathrm{30}{}^{\circ }\mathrm{C})\approx \mathrm{1.28},\text{for}T<-\mathrm{30}{}^{\circ }\mathrm{C}.\end{array}$$

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

$$\begin{array}{}\text{(17)}& f_{\mathrm{vnt}}={\displaystyle \frac{b_{\mathrm{1}}+b_{\mathrm{2}}{X}^{\mathit{\gamma }}{\left(c_i/C\right)}^{\mathrm{1}/\mathrm{2}}}{b_{\mathrm{1}}+b_{\mathrm{2}}{X}^{\mathit{\gamma }}{\left(a_i/C\right)}^{\mathrm{1}/\mathrm{2}}}}.\end{array}$$

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 dV_i is given by

$$\begin{array}{}\text{(18)}& \text{d}{V}_i={\displaystyle \frac{\text{d}{m}_i}{{\mathit{\rho }}_{\mathrm{dep}}}},\text{for}\text{d}{m}_i\ge \mathrm{0}\text{(deposition)}.\end{array}$$

Deposition density ρ_dep can be expressed as

$$\begin{array}{}\text{(19)}& {\mathit{\rho }}_{\mathrm{dep}}=\left\{\begin{array}{ll}{\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}},& \text{for}\mathrm{\Gamma }\left(T_i\right)<\mathrm{1}\bigwedge a_i<\mathrm{100}\mathrm{\mu }\mathrm{m};\\ {\mathit{\rho }}_{\mathrm{dep}}^{\mathrm{CL}\mathrm{94}},& \text{otherwise}.\end{array}\right.\end{array}$$

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),

$$\begin{array}{}\text{(20)}& {\mathit{\rho }}_{\mathrm{dep}}^{\mathrm{CL}\mathrm{94}}={\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}}\exp \left[-{\displaystyle \frac{\mathrm{3}max(\mathrm{\Delta }{\mathit{\rho }}_i-\mathrm{0.05}\mathrm{g}{\mathrm{m}}^{-\mathrm{3}},\mathrm{0})}{\mathrm{\Gamma }\left(T_i\right)\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}}}\right],\end{array}$$

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

$$\begin{array}{}\text{(21)}& \mathrm{\Delta }{\mathit{\rho }}_i={\displaystyle \frac{S_i^{\mathrm{i}}-\mathrm{1}}{D_{\mathrm{v}}(F_{\mathrm{k}}^{\mathrm{i}}+F_{\mathrm{d}}^{\mathrm{i}})}}.\end{array}$$

Here, following Miller and Young (1979), we limit ρ_vi by water saturation and replace the Δρ_i in Eq. (20) with

$$\begin{array}{}\text{(22)}& (\mathrm{\Delta }{\mathit{\rho }}_i{)}^{\downarrow }={\displaystyle \frac{min(S_i^{\mathrm{i}},e_{\mathrm{s}}^{\mathrm{w}}(T_i)/e_{\mathrm{s}}^{\mathrm{i}}(T_i\left)\right)-\mathrm{1}}{D_{\mathrm{v}}(F_{\mathrm{k}}^{\mathrm{i}}+F_{\mathrm{d}}^{\mathrm{i}})}}.\end{array}$$

For sublimation, the particle volume change dV_i is given by

$$\begin{array}{}\text{(23)}& \text{d}{V}_i={\displaystyle \frac{\text{d}{m}_i}{{\mathit{\rho }}_{\mathrm{sbl}}}},\text{for}\text{d}{m}_i<\mathrm{0}\text{(sublimation)},\end{array}$$

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.,

$$\begin{array}{}\text{(24)}& {\mathit{\rho }}_{\mathrm{sbl}}={\mathit{\rho }}_i^{\mathrm{i}}.\end{array}$$

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

$$\begin{array}{}\text{(25)}& {\mathit{\rho }}_i^{\mathrm{i}}(t+\text{d}t)={\displaystyle \frac{m_i+\text{d}{m}_i}{V_i+\text{d}{V}_i}},\end{array}$$

where dm_i is given in Eqs. (11) and (14), and dV_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 dt, the two radii become

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }{a}_i(t+\text{d}t)=a_i\exp \left({\displaystyle \frac{\text{d}\log V_i}{\mathrm{2}+{\mathrm{\Gamma }}^{*}}}\right),\text{(27)}& {\displaystyle }& {\displaystyle }{c}_i(t+\text{d}t)=c_i\exp \left({\displaystyle \frac{{\mathrm{\Gamma }}^{*}d\log V_i}{\mathrm{2}+{\mathrm{\Gamma }}^{*}}}\right).\end{array}$$

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}\mathrm{\mu }\mathrm{m}$,

$$\begin{array}{}\text{(28)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_i^{\mathrm{i}\prime }(t+\text{d}t)={\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}},\text{(29)}& {\displaystyle }& {\displaystyle }{a}_i^{\prime }(t+\text{d}t)=c_i^{\prime }(t+\text{d}t)={\left({\displaystyle \frac{m_i+\text{d}{m}_i}{(\mathrm{4}\mathit{\pi }/\mathrm{3}){\mathit{\rho }}_i^{\mathrm{i}\prime }(t+\text{d}t)}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\end{array}$$

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):

$$\begin{array}{}\text{(30)}& m_i^{\mathrm{rime}}(t+\text{d}t)=\left\{\begin{array}{ll}{m}_i^{\mathrm{rime}},& \text{for}\text{d}{m}_i\ge \mathrm{0};\\ m_i^{\mathrm{rime}}{\displaystyle \frac{m_i+\text{d}{m}_i}{m_i}},& \text{for}\text{d}{m}_i<\mathrm{0}.\end{array}\right.\end{array}$$

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 dt. 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

$$\begin{array}{}\text{(31)}& P_{jk}=K({\mathit{a}}_j,{\mathit{a}}_k;\mathit{G}){\displaystyle \frac{\text{d}t}{\mathrm{\Delta }V}},\end{array}$$

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

$$\begin{array}{}\text{(32)}& K_{\mathrm{coal}}=E_{\mathrm{coal}}(r_j,r_k)\mathit{\pi }(r_j+r_k{)}^{\mathrm{2}}|v_j^{\mathrm{\infty }}-v_k^{\mathrm{\infty }}|,\end{array}$$

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:

$$\begin{array}{}\text{(33)}& {\displaystyle }& {\displaystyle }{r}_j^{\prime }=(r_j^{\mathrm{3}}+r_k^{\mathrm{3}}{)}^{\frac{\mathrm{1}}{\mathrm{3}}},\text{(34)}& {\displaystyle }& {\displaystyle }{m}_{\mathit{\alpha }j}^{\mathrm{sol}\prime }=m_{\mathit{\alpha }j}^{\mathrm{sol}}+m_{\mathit{\alpha }k}^{\mathrm{sol}},\mathit{\alpha }=\mathrm{1},\mathrm{2},\mathrm{\dots },N^{\mathrm{sol}}\text{(35)}& {\displaystyle }& {\displaystyle }{m}_{\mathit{\beta }j}^{\mathrm{insol}\prime }=m_{\mathit{\beta }j}^{\mathrm{insol}}+m_{\mathit{\beta }k}^{\mathrm{insol}},\mathit{\beta }=\mathrm{1},\mathrm{2},\mathrm{\dots },N^{\mathrm{insol}}\text{(36)}& {\displaystyle }& {\displaystyle }{T}_j^{\mathrm{fz}\prime }=max(T_j^{\mathrm{fz}},T_k^{\mathrm{fz}}),\end{array}$$

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

$$\begin{array}{}\text{(37)}& K_{\mathrm{rime}}=E_{\mathrm{rime}}{A}_{\mathrm{g}}|v_j^{\mathrm{\infty }}-v_k^{\mathrm{\infty }}|,\end{array}$$

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

$$\begin{array}{}\text{(38)}& A_{\mathrm{g}}=\mathit{\pi }(a_j+r_k)\mathit{\{}max(a_j,c_j)+r_k\mathit{\}}-(A_j^{\mathrm{ce}}-A_j);\end{array}$$

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):

$$\begin{array}{}\text{(39)}& E_{\mathrm{rime}}=E_{\mathrm{BG}\mathrm{74}}(p^{\mathrm{i}/\mathrm{w}},N_{\mathrm{Re}k}^{\mathrm{w}},N_{\mathrm{St}}^{\mathrm{i}/\mathrm{w}}),\end{array}$$

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),

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{E}_{\mathrm{rime}}=& {\mathit{\varphi }}_j{E}_{\mathrm{BG}\mathrm{74}}(p^{\mathrm{w}/\mathrm{i}},N_{\mathrm{Re}j}^{\mathrm{i}},N_{\mathrm{mFr}})\\ & +(\mathrm{1}-{\mathit{\varphi }}_j)E_{\mathrm{EM}\mathrm{17}}^{\mathrm{pln}}(N_{\mathrm{Re}j}^{\mathrm{i}},N_{\mathrm{mFr}}).\end{array}\end{array}$$

For ϕ_j>1 (columnar),

$$\begin{array}{}\text{(41)}& \begin{array}{rl}{E}_{\mathrm{rime}}=& {\displaystyle \frac{\mathrm{1}}{{\mathit{\varphi }}_j}}{E}_{\mathrm{BG}\mathrm{74}}(p^{\mathrm{w}/\mathrm{i}},N_{\mathrm{Re}j}^{\mathrm{i}},N_{\mathrm{mFr}})\\ & +\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\varphi }}_j}}\right)E_{\mathrm{EM}\mathrm{17}}^{\mathrm{clm}}(N_{\mathrm{Re}j}^{\mathrm{clm}},N_{\mathrm{mFr}}).\end{array}\end{array}$$

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 2a_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:

$$\begin{array}{}\text{(42)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_j^{\mathrm{i}\prime }={\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}},\text{(43)}& {\displaystyle }& {\displaystyle }{a}_j^{\prime }=c_j^{\prime }={\left({\displaystyle \frac{m_j+m_k}{(\mathrm{4}\mathit{\pi }/\mathrm{3}){\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\text{(44)}& {\displaystyle }& {\displaystyle }{m}_j^{\mathrm{rime}\prime }=m_j^{\mathrm{rime}}+m_k,\text{(45)}& {\displaystyle }& {\displaystyle }{n}_j^{\mathrm{mono}\prime }=n_j^{\mathrm{mono}},\text{(46)}& {\displaystyle }& {\displaystyle }{m}_{\mathit{\alpha }j}^{\mathrm{sol}\prime }=m_{\mathit{\alpha }j}^{\mathrm{sol}}+m_{\mathit{\alpha }k}^{\mathrm{sol}},\mathit{\alpha }=\mathrm{1},\mathrm{2},\mathrm{\dots },N^{\mathrm{sol}},\text{(47)}& {\displaystyle }& {\displaystyle }{m}_{\mathit{\beta }j}^{\mathrm{insol}\prime }=m_{\mathit{\beta }j}^{\mathrm{insol}}+m_{\mathit{\beta }k}^{\mathrm{insol}},\mathit{\beta }=\mathrm{1},\mathrm{2},\mathrm{\dots },N^{\mathrm{insol}},\text{(48)}& {\displaystyle }& {\displaystyle }{T}_j^{\mathrm{fz}\prime }=max(T_j^{\mathrm{fz}},T_k^{\mathrm{fz}}),\end{array}$$

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),

$$\begin{array}{}\text{(49)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_j^{\mathrm{i}\prime }={\displaystyle \frac{m_j+m_k}{V_j+m_k/{\mathit{\rho }}_{\mathrm{rime}}}},\text{(50)}& {\displaystyle }& {\displaystyle }{a}_j^{\prime }=a_j,\text{(51)}& {\displaystyle }& {\displaystyle }{c}_j^{\prime }={\displaystyle \frac{V_j+m_k/{\mathit{\rho }}_{\mathrm{rime}}}{(\mathrm{4}\mathit{\pi }/\mathrm{3})a_j^{\mathrm{2}}}}.\end{array}$$

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),

$$\begin{array}{}\text{(52)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_j^{\mathrm{i}\prime }={\displaystyle \frac{m_j+m_k}{V_j+m_k/{\mathit{\rho }}_{\mathrm{rime}}}},\text{(53)}& {\displaystyle }& {\displaystyle }{a}_j^{\prime }={\left\{{\displaystyle \frac{V_j+m_k/{\mathit{\rho }}_{\mathrm{rime}}}{(\mathrm{4}\mathit{\pi }/\mathrm{3})c_j}}\right\}}^{\frac{\mathrm{1}}{\mathrm{2}}},\text{(54)}& {\displaystyle }& {\displaystyle }{c}_j^{\prime }=c_j.\end{array}$$

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:

$$\begin{array}{}\text{(55)}& {\mathit{\rho }}_{\mathrm{rime}}=max\mathit{\{}min\mathit{\{}{\mathit{\rho }}_{\mathrm{rime}}^{\mathrm{HP}\mathrm{85}}\left(Y\right),\mathrm{0.91}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}\mathit{\}},\mathrm{0.1}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}\mathit{\}},\end{array}$$

where $Y:=(-r_k{v}_{\mathrm{imp}}/T_j^{\mathrm{sfc}})/\left(\mathrm{\mu }\mathrm{m}\mathrm{m}{\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}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$, B₁=0.44, $B_{\mathrm{2}}=-\mathrm{0.03115}$, $B_{\mathrm{3}}=-\mathrm{1.7030}$, B₄=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

$$\begin{array}{}\text{(58)}& T_j^{\mathrm{sfc}}=T_j+{\displaystyle \frac{L_{\mathrm{s}}{D}_{\mathrm{v}}}{k}}\mathrm{\Delta }{\mathit{\rho }}_j,\end{array}$$

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

$$\begin{array}{}\text{(59)}& K_{\mathrm{agg}}=E_{\mathrm{agg}}{\left(A_j^{\frac{\mathrm{1}}{\mathrm{2}}}+A_k^{\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}|v_j^{\mathrm{\infty }}-v_k^{\mathrm{\infty }}|,\end{array}$$

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{‾}}{{\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},\min }$, which can be evaluated using Eq. (61), which we will derive shortly.

In contrast, if $\stackrel{\mathrm{‾}}{{\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}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. This choice of value is roughly consistent with observations by Magono and Nakamura (1965). If $\stackrel{\mathrm{‾}}{{\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},\max }$. Let us assume ${\mathit{\rho }}_{jk}^{\mathrm{i},\max }=\stackrel{\mathrm{‾}}{{\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),

$$\begin{array}{}\text{(60)}& a_j^{\prime }=a_j,\end{array}$$

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_{\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},\min }$ then becomes

$$\begin{array}{}\text{(61)}& {\mathit{\rho }}_{jk}^{\mathrm{i},\min }={\displaystyle \frac{m_j+m_k}{V_{\max }}}.\end{array}$$

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

$$\begin{array}{}\text{(62)}& {\mathit{\rho }}_j^{\mathrm{i}\prime }={\displaystyle \frac{({\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}}-\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{jk}^{\mathrm{i}}}){\mathit{\rho }}_{jk}^{\mathrm{i},\max }+(\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{jk}^{\mathrm{i}}}-{\mathit{\rho }}_{\mathrm{crt}}^{\mathrm{i}}){\mathit{\rho }}_{jk}^{\mathrm{i},\min }}{{\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}}-{\mathit{\rho }}_{\mathrm{crt}}^{\mathrm{i}}}},\end{array}$$

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

$$\begin{array}{}\text{(63)}& {\displaystyle }& {\displaystyle }{c}_j^{\prime }={\displaystyle \frac{m_j+m_k}{{\mathit{\rho }}_j^{\mathrm{i}\prime }(\mathrm{4}\mathit{\pi }/\mathrm{3}){a_j^{\prime }}^{\mathrm{2}}}},\text{(64)}& {\displaystyle }& {\displaystyle }{m}_j^{\mathrm{rime}\prime }=m_j^{\mathrm{rime}}+m_k^{\mathrm{rime}},\text{(65)}& {\displaystyle }& {\displaystyle }{n}_j^{\mathrm{mono}\prime }=n_j^{\mathrm{mono}}+n_k^{\mathrm{mono}},\text{(66)}& {\displaystyle }& {\displaystyle }{m}_{\mathit{\alpha }j}^{\mathrm{sol}\prime }=m_{\mathit{\alpha }j}^{\mathrm{sol}}+m_{\mathit{\alpha }k}^{\mathrm{sol}},\mathit{\alpha }=\mathrm{1},\mathrm{2},\mathrm{\dots },N^{\mathrm{sol}},\text{(67)}& {\displaystyle }& {\displaystyle }{m}_{\mathit{\beta }j}^{\mathrm{insol}\prime }=m_{\mathit{\beta }j}^{\mathrm{insol}}+m_{\mathit{\beta }k}^{\mathrm{insol}},\mathit{\beta }=\mathrm{1},\mathrm{2},\mathrm{\dots },N^{\mathrm{insol}},\text{(68)}& {\displaystyle }& {\displaystyle }{T}_j^{\mathrm{fz}\prime }=max(T_j^{\mathrm{fz}},T_k^{\mathrm{fz}}).\end{array}$$

For ϕ_j>1 (columnar), the polar axis length is preserved

$$\begin{array}{}\text{(69)}& c_j^{\prime }=c_j.\end{array}$$

If approximating the largest possible particle using an ellipsoid, the largest possible volume becomes $V_{\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

$$\begin{array}{}\text{(70)}& a_j^{\prime }={\left\{{\displaystyle \frac{m_j+m_k}{{\mathit{\rho }}_j^{\mathrm{i}\prime }(\mathrm{4}\mathit{\pi }/\mathrm{3})c_j^{\prime }}}\right\}}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

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}\mathrm{kg}{\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:

$$\begin{array}{}\text{(71)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{U}\right)={\left.{\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}\right|}_{\mathrm{cm}},\text{(72)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\rho }{q}_{\mathrm{v}}}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }{q}_{\mathrm{v}}\mathit{U}\right)={\left.{\displaystyle \frac{\partial \mathit{\rho }{q}_{\mathrm{v}}}{\partial t}}\right|}_{\mathrm{cm}}+D_{\mathrm{v}}{\mathrm{\nabla }}^{\mathrm{2}}\left(\mathit{\rho }{q}_{\mathrm{v}}\right),\text{(73)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\rho }\mathit{U}}{\partial t}}+\mathrm{\nabla }\cdot (\mathit{\rho }\mathit{U}\otimes \mathit{U})=-\mathrm{\nabla }P-\mathit{\rho }g\widehat{\mathit{z}}+{\left.{\displaystyle \frac{\partial \mathit{\rho }\mathit{U}}{\partial t}}\right|}_{\mathrm{cm}}+\mathit{\mu }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{U},\text{(74)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\rho }\mathit{\theta }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{\theta }\mathit{U}\right)={\left.{\displaystyle \frac{\partial \mathit{\rho }\mathit{\theta }}{\partial t}}\right|}_{\mathrm{cm}}+{\displaystyle \frac{k}{c_{\mathrm{p}}}}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\theta },\text{(75)}& {\displaystyle }& {\displaystyle }P=\mathit{\rho }RT=P_{\mathrm{0}}{\left({\displaystyle \frac{\mathit{\rho }\mathit{\theta }R}{P_{\mathrm{0}}}}\right)}^{c_{\mathrm{p}}/\left(c_{\mathrm{p}}-R\right)},\end{array}$$

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:

$$\begin{array}{}\text{(76)}& {\left.{\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}\right|}_{\mathrm{cm}}={\left.{\displaystyle \frac{\partial \mathit{\rho }{q}_{\mathrm{v}}}{\partial t}}\right|}_{\mathrm{cm}}=s_{\mathrm{v}}+s_{\mathrm{s}}.\end{array}$$

Here, s_v and s_s are sources of vapor through condensation/evaporation and deposition/sublimation, respectively:

$$\begin{array}{}\text{(77)}& {\displaystyle }& {\displaystyle }{s}_{\mathrm{v}}(\mathit{x},t)=-\sum _{i\in I_{\mathrm{r}}\left(t\right)}{\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_i(t\left)\right){\left.{\displaystyle \frac{\text{d}{m}_i}{\text{d}t}}\right|}_{\mathrm{cnd}/\mathrm{evp}},\text{(78)}& {\displaystyle }& {\displaystyle }{s}_{\mathrm{s}}(\mathit{x},t)=-\sum _{i\in I_{\mathrm{r}}\left(t\right)}{\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_i(t\left)\right){\left.{\displaystyle \frac{\text{d}{m}_i}{\text{d}t}}\right|}_{\mathrm{dep}/\mathrm{sbl}},\end{array}$$

where δ³(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:

$$\begin{array}{}\text{(79)}& {\left.{\displaystyle \frac{\partial \mathit{\rho }\mathit{\theta }}{\partial t}}\right|}_{\mathrm{cm}}=-{\displaystyle \frac{L_{\mathrm{v}}{s}_{\mathrm{v}}+L_{\mathrm{s}}{s}_{\mathrm{s}}+L_{\mathrm{f}}{s}_{\mathrm{f}}}{c_p\mathrm{\Pi }}},\end{array}$$

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 nth 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 nth melting event, respectively. Let $t_n^{\mathrm{rime}}$ and $i_n^{\mathrm{rime}}$ be the time and rimed droplet of the nth riming event, respectively. Then, s_f is given by

$$\begin{array}{}\text{(80)}& \begin{array}{rl}{s}_{\mathrm{f}}(\mathit{x},t)=& -\sum _{\text{freezing event}n}{\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_{i_n^{\mathrm{fz}}}(t\left)\right)\mathit{\delta }(t-t_n^{\mathrm{fz}})m_{i_n^{\mathrm{fz}}}\left(t\right)\\ & +\sum _{\text{melting event}n}{\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_{i_n^{\mathrm{mlt}}}(t\left)\right)\mathit{\delta }(t-t_n^{\mathrm{mlt}})m_{i_n^{\mathrm{mlt}}}\left(t\right)\\ & -\sum _{\text{riming event}n}{\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_{i_n^{\mathrm{rime}}}(t\left)\right)\mathit{\delta }(t-t_n^{\mathrm{rime}})m_{i_n^{\mathrm{rime}}}\left(t\right).\end{array}\end{array}$$

$\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⁻¹ in 100 s through rapid precipitation development, the contribution of the second term is much smaller than that of the first term: $\mathrm{10}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}/\mathrm{100}\mathrm{s}\ll g$. Thus, we finally obtain

$$\begin{array}{}\text{(81)}& \begin{array}{rl}{\left.{\displaystyle \frac{\partial \mathit{\rho }\mathit{U}}{\partial t}}\right|}_{\mathrm{cm}}& =-\sum _{i\in I_{\mathrm{r}}\left(t\right)}{\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_i(t\left)\right){\mathit{F}}_i^{\mathrm{drg}}\\ & \approx -\left[\sum _{i\in I_{\mathrm{r}}\left(t\right)}{\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_i(t\left)\right)m_i\left(t\right)\right]g\widehat{\mathit{z}}.\end{array}\end{array}$$

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.

5 Numerical schemes and implementation

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 ρU are defined on the faces of each grid cell. To simplify the notation, we use 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 ith 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 x and time t. The following relation then holds:

$$\begin{array}{}\text{(82)}& n(\mathit{a},\mathit{x},t)=〈\sum _{i\in I_{\mathrm{r}}\left(t\right)}{\mathit{\delta }}^d(\mathit{a}-{\mathit{a}}_i(t\left)\right){\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_i(t\left)\right)〉,\end{array}$$

where 〈⋯〉 denotes the mean, and δ^d(a) is the d-dimensional Dirac delta function. Super-particles reproduce the behavior of particles in expectation:

$$\begin{array}{}\text{(83)}& \begin{array}{rl}n(\mathit{a},\mathit{x},t)& =〈\sum _{i\in I_{\mathrm{s}}\left(t\right)}{\mathit{\xi }}_i\left(t\right){\mathit{\delta }}^d(\mathit{a}-{\mathit{a}}_i(t\left)\right){\mathit{\delta }}^{\mathrm{3}}(\mathit{x}-{\mathit{x}}_i(t\left)\right)〉\\ & =N_{\mathrm{s}}\left(t\right)\sum _{\mathit{\xi }=\mathrm{1}}^{\mathrm{\infty }}\mathit{\xi }p(\mathit{\xi },\mathit{a},\mathit{x},t),\end{array}\end{array}$$

where $p(\mathit{\xi },\mathit{a},\mathit{x},t)$ is the probability density that a super-particle has multiplicity ξ, attributes a, and position x at time 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 x based on the probability density function p(a,x), and determine the super-particle's multiplicity ξ by using a deterministic function of a and x, i.e., $\mathit{\xi }=\mathit{\xi }(\mathit{a},\mathit{x})$. Then, Eq. (83) at t=0 reduces to

$$\begin{array}{}\text{(84)}& n(\mathit{a},\mathit{x},\mathrm{0})=N_{\mathrm{s}}\left(\mathrm{0}\right)\mathit{\xi }(\mathit{a},\mathit{x})p(\mathit{a},\mathit{x}).\end{array}$$

If we set ξ(a,x) 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 ν_const-init in Unterstrasser et al. (2017).

Instead, we can set p(a,x) as a constant (i.e., uniform sampling). Multiplicity then becomes proportional to the initial distribution function of real particles:

$$\begin{array}{}\text{(85)}& \mathit{\xi }(\mathit{a},\mathit{x})={\displaystyle \frac{n(\mathit{a},\mathit{x},\mathrm{0})}{N_{\mathrm{s}}\left(\mathrm{0}\right)p}},p(\mathit{a},\mathit{x})=p=\mathrm{const}.\end{array}$$

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.

Table 1 An example of the calculation order when updating the system state from t to t+Δt. We first calculate the fluid dynamics and then calculate cloud microphysics. Each process is integrated one time step forward at a time. Processes lagging in time are calculated preferentially.

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² 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{‾}}{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{‾}}{P}$ can be evaluated as

$$\begin{array}{}\text{(86)}& \stackrel{\mathrm{‾}}{P}\approx N_{\mathrm{r}}^{\prime }\stackrel{\mathrm{‾}}{K}{\displaystyle \frac{\mathrm{\Delta }{t}_{\mathrm{collis}}}{\mathrm{\Delta }V}}\approx n_{\mathrm{r}}\stackrel{\mathrm{‾}}{K}\mathrm{\Delta }{t}_{\mathrm{collis}},\end{array}$$

where $N_{\mathrm{r}}^{\prime }$ is the number of real particles in a volume ΔV, $\stackrel{\mathrm{‾}}{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{‾}}{P}<\mathrm{1}$ has to be satisfied, we obtain

$$\begin{array}{}\text{(87)}& \mathrm{\Delta }{t}_{\mathrm{collis}}<\mathrm{1}/\left(n_{\mathrm{r}}\stackrel{\mathrm{‾}}{K}\right).\end{array}$$

Here, we relate the above argument to that of Shima et al. (2009). Let $\stackrel{\mathrm{‾}}{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{‾}}{P_{\mathrm{s}}}$ is what Shima et al. (2009) evaluated in the last paragraph of Sect. 5.1.3. We can derive $\stackrel{\mathrm{‾}}{P_{\mathrm{s}}}\approx \stackrel{\mathrm{‾}}{P}$ as follows:

$$\begin{array}{}\text{(88)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{P_{\mathrm{s}}}& \approx \left\{{\displaystyle \frac{N_{\mathrm{s}}^{\prime }(N_{\mathrm{s}}^{\prime }-\mathrm{1})}{\mathrm{2}}}/\left[{\displaystyle \frac{N_{\mathrm{s}}^{\prime }}{\mathrm{2}}}\right]\right\}\stackrel{\mathrm{‾}}{\mathit{\xi }}\stackrel{\mathrm{‾}}{K}{\displaystyle \frac{\mathrm{\Delta }{t}_{\mathrm{collis}}}{\mathrm{\Delta }V}}\\ & \approx N_{\mathrm{s}}^{\prime }{\displaystyle \frac{N_{\mathrm{r}}^{\prime }}{N_{\mathrm{s}}^{\prime }}}\stackrel{\mathrm{‾}}{K}{\displaystyle \frac{\mathrm{\Delta }{t}_{\mathrm{collis}}}{\mathrm{\Delta }V}}\\ & \approx \stackrel{\mathrm{‾}}{P},\end{array}\end{array}$$

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{‾}}{\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{̃}}{\mathit{\gamma }}$ number of droplets k (see Sect. 5.1.3 of Shima et al., 2009, for the definition of $\stackrel{\mathrm{̃}}{\mathit{\gamma }}$), we modify Eqs. (33)–(35) by applying

$$\begin{array}{}\text{(89)}& (r_k^{\mathrm{3}},m_{\mathit{\alpha }k}^{\mathrm{sol}},m_{\mathit{\beta }k}^{\mathrm{insol}})\to \stackrel{\mathrm{̃}}{\mathit{\gamma }}(r_k^{\mathrm{3}},m_{\mathit{\alpha }k}^{\mathrm{sol}},m_{\mathit{\beta }k}^{\mathrm{insol}}).\end{array}$$

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

$$\begin{array}{}\text{(90)}& (r_j^{\mathrm{3}},m_{\mathit{\alpha }j}^{\mathrm{sol}},m_{\mathit{\beta }j}^{\mathrm{insol}})\to \stackrel{\mathrm{̃}}{\mathit{\gamma }}(r_j^{\mathrm{3}},m_{\mathit{\alpha }j}^{\mathrm{sol}},m_{\mathit{\beta }j}^{\mathrm{insol}}).\end{array}$$

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

$$\begin{array}{}\text{(91)}& \begin{array}{rl}& (m_k,V_k,m_k^{\mathrm{rime}},n_k^{\mathrm{mono}},m_{\mathit{\alpha }k}^{\mathrm{sol}},m_{\mathit{\beta }k}^{\mathrm{insol}})\\ \to \stackrel{\mathrm{̃}}{\mathit{\gamma }}& (m_k,V_k,m_k^{\mathrm{rime}},n_k^{\mathrm{mono}},m_{\mathit{\alpha }k}^{\mathrm{sol}},m_{\mathit{\beta }k}^{\mathrm{insol}}).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(92)}& \begin{array}{rl}& (m_j,V_j,m_j^{\mathrm{rime}},n_j^{\mathrm{mono}},m_{\mathit{\alpha }j}^{\mathrm{sol}},m_{\mathit{\beta }j}^{\mathrm{insol}})\\ \to \stackrel{\mathrm{̃}}{\mathit{\gamma }}& (m_j,V_j,m_j^{\mathrm{rime}},n_j^{\mathrm{mono}},m_{\mathit{\alpha }j}^{\mathrm{sol}},m_{\mathit{\beta }j}^{\mathrm{insol}}).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(93)}& V_{\max }=\left\{\begin{array}{l}(\mathrm{4}\mathit{\pi }/\mathrm{3})a_j^{\mathrm{2}}\mathit{\{}{c}_j+\stackrel{\mathrm{̃}}{\mathit{\gamma }}min(a_k,c_k\left)\mathit{\right\}},\\ (\mathrm{4}\mathit{\pi }/\mathrm{3})a_j^{\mathrm{2}}\mathit{\{}\stackrel{\mathrm{̃}}{\mathit{\gamma }}{c}_j+min(a_k,c_k\left)\mathit{\right\}}.\end{array}\right.\end{array}$$

For columnar collector j,

$$\begin{array}{}\text{(94)}& V_{\max }=\left\{\begin{array}{l}(\mathrm{4}\mathit{\pi }/\mathrm{3})c_j\mathit{\{}{a}_j+\stackrel{\mathrm{̃}}{\mathit{\gamma }}min(a_k,c_k\left)\mathit{\right\}}max(a_j,a_k,c_k),\\ (\mathrm{4}\mathit{\pi }/\mathrm{3})c_j\mathit{\{}\stackrel{\mathrm{̃}}{\mathit{\gamma }}{a}_j+min(a_k,c_k\left)\mathit{\right\}}max(a_j,a_k,c_k).\end{array}\right.\end{array}$$

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.

6 Design of numerical experiments for model evaluation: 2-D simulation of an isolated cumulonimbus

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.

Table 2 Summary of numerical experiments for model evaluation. The domain is two-dimensional (x–z): 60 km in the horizontal direction and 16 km in the vertical direction. The initial profile of moist air is given by sounding data from Midland, Texas, on 13 August 1999, as shown in Fig. 4 of Khain et al. (2004). The particles are initially distributed uniformly in space at random and consist of pure ammonium bisulfate aerosol particles and mineral dust internally mixed with ammonium bisulfate. The numerical parameters used in each case are listed in the table, and values changed from the CTRL case are in bold. We conducted a 10-member ensemble of simulations for each case by changing the pseudo-random number sequence to evaluate fluctuation. The vertical dots indicate NSP008, NSP016, NSP032, and NSP064.

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⁻¹ near the surface to 7 m s⁻¹ 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,

$$\begin{array}{}\text{(95)}& {\displaystyle \frac{\text{d}{N}^{\mathrm{sulf}}}{\text{d}\log r_{\mathrm{dry}}^{\mathrm{sulf}}}}=\sum _{a=\mathrm{1}}^{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}}}{\displaystyle \frac{c_a^{\mathrm{sulf}}}{\log {\mathit{\sigma }}_a}}\exp \left[{\displaystyle \frac{-{\left(\log r_{\mathrm{dry}}^{\mathrm{sulf}}-\log r_a\right)}^{\mathrm{2}}}{\mathrm{2}{\log }^{\mathrm{2}}{\mathit{\sigma }}_a}}\right],\end{array}$$

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}{\mathrm{cm}}^{-\mathrm{3}}$ and $c_{\mathrm{2}}^{\mathrm{sulf}}=\mathrm{45}{\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}{\mathrm{cm}}^{-\mathrm{3}}$. The geometric mean radii are r₁=0.03 µm and r₂=0.14 µm, with geometric standard deviations of σ₁=1.28 and σ₂=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

$$\begin{array}{}\text{(96)}& n^{\mathrm{sulf}}(\log r_{\mathrm{dry}}^{\mathrm{sulf}},T^{\mathrm{fz}})={\displaystyle \frac{\text{d}{N}^{\mathrm{sulf}}}{\text{d}\log r_{\mathrm{dry}}^{\mathrm{sulf}}}}\mathit{\delta }(T^{\mathrm{fz}}-(-\mathrm{38}{}^{\circ }\mathrm{C}\left)\right).\end{array}$$

The other aerosol population consists of mineral dust internally mixed with ammonium bisulfate. We set the number concentration to $c^{\mathrm{dust}}=\mathrm{1}{\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:

$$\begin{array}{}\text{(97)}& n_{\mathrm{S}}\left(T\right)=\left\{\begin{array}{ll}\mathrm{0},& \text{for}T>T_{\max }^{\mathrm{fz}};\\ n_{\mathrm{S}}^{\mathrm{Niemand}}\left(T\right),& \text{for}{T}_{\max }^{\mathrm{fz}}\ge T>T_{\min }^{\mathrm{fz}};\\ n_{\mathrm{S}}^{\mathrm{Niemand}}\left(T_{\min }^{\mathrm{fz}}\right),& \text{for}{T}_{\min }^{\mathrm{fz}}\ge T;\end{array}\right.\end{array}$$

where $T_{\max }^{\mathrm{fz}}=-\mathrm{12}$ ^∘C and $T_{\min }^{\mathrm{fz}}=-\mathrm{36}$ ^∘C. The mineral dust surface area is given by A^insol=π(d^dust)². 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

$$\begin{array}{}\text{(98)}& \begin{array}{rl}{n}^{\mathrm{dust}}& (d^{\mathrm{dust}},\log r_{\mathrm{dry}}^{\mathrm{sulf}},T^{\mathrm{fz}})\\ =& \mathit{\delta }(d^{\mathrm{dust}}-\mathrm{1}\mathrm{\mu }\mathrm{m}){\displaystyle \frac{c^{\mathrm{dust}}}{c^{\mathrm{sulf}}}}{\displaystyle \frac{\text{d}{N}^{\mathrm{sulf}}}{\text{d}\log r_{\mathrm{dry}}^{\mathrm{sulf}}}}\\ & \left[p\right(T^{\mathrm{fz}}\left)H\right(T^{\mathrm{fz}}+\mathrm{38}{}^{\circ }\mathrm{C})+P_{\mathrm{INia}}\mathit{\delta }(T^{\mathrm{fz}}+\mathrm{38}{}^{\circ }\mathrm{C}\left)\right],\end{array}\end{array}$$

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

$$\begin{array}{}\text{(99)}& {\displaystyle }& {\displaystyle }{\left.{\displaystyle \frac{\partial \mathit{\rho }\mathit{\theta }}{\partial t}}\right|}_{\mathrm{sfc}}=\mathit{\rho }Hmax(W,\mathrm{0}),\text{(100)}& {\displaystyle }& {\displaystyle }W=\left(-{\displaystyle \frac{\mathrm{4}}{w^{\mathrm{2}}}}\right)\left[(x-x_{\mathrm{0}}{)}^{\mathrm{2}}-{\left({\displaystyle \frac{w}{\mathrm{2}}}\right)}^{\mathrm{2}}\right]\exp \left[-{\displaystyle \frac{z-z_{\mathrm{0}}}{z_{\mathrm{0}}}}\right],\end{array}$$

where $H=\mathrm{10}\mathrm{K}{\mathrm{h}}^{-\mathrm{1}}$, x₀=5 km, w=10 km, and z₀=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}\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}\mathrm{s}$, Δt_collis=0.2 s, $\mathrm{\Delta }{t}_{\mathrm{cnd}/\mathrm{evp}}=\mathrm{\Delta }{t}_{\mathrm{dep}/\mathrm{sbl}}=\mathrm{0.1}\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},\min }^{\mathrm{sulf}},r_{\mathrm{dry},\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},\min }^{\mathrm{sulf}}=\mathrm{10.0}\mathrm{nm}$ and $r_{\mathrm{dry},\max }^{\mathrm{sulf}}=\mathrm{5.0}\mathrm{\mu }\mathrm{m}$. From Eqs. (85) and (96), the super-particle's multiplicity is then given by

$$\begin{array}{}\text{(101)}& \begin{array}{rl}{\mathit{\xi }}_i& ={\displaystyle \frac{n^{\mathrm{sulf}}(\log r_{\mathrm{dry},i}^{\mathrm{sulf}},T_i^{\mathrm{fz}})}{N_{\mathrm{s}}\left(\mathrm{0}\right)/\mathrm{2}}}{\displaystyle \frac{V_{\mathrm{domain}}\log (r_{\mathrm{dry},\max }^{\mathrm{sulf}}/r_{\mathrm{dry},\min }^{\mathrm{sulf}})}{\mathit{\delta }(T^{\mathrm{fz}}-(-\mathrm{38}{}^{\circ }\mathrm{C}\left)\right)}}\\ & ={\displaystyle \frac{\text{d}{N}^{\mathrm{sulf}}}{\text{d}\log r_{\mathrm{dry}}^{\mathrm{sulf}}}}\left(\log r_{\mathrm{dry},i}^{\mathrm{sulf}}\right){\displaystyle \frac{\log (r_{\mathrm{dry},\max }^{\mathrm{sulf}}/r_{\mathrm{dry},\min }^{\mathrm{sulf}})}{c^{\mathrm{SP}}/\mathrm{2}}},\end{array}\end{array}$$

where V_domain is the total volume of the domain, n and p in Eq. (85) in this case are given by

$$\begin{array}{}\text{(102)}& {\displaystyle }& {\displaystyle }n=n^{\mathrm{sulf}}(\log r_{\mathrm{dry},i}^{\mathrm{sulf}},T_i^{\mathrm{fz}}),\text{(103)}& {\displaystyle }& {\displaystyle }p={\displaystyle \frac{\mathit{\delta }(T^{\mathrm{fz}}-(-\mathrm{38}{}^{\circ }\mathrm{C}\left)\right)}{V_{\mathrm{domain}}\log (r_{\mathrm{dry},\max }^{\mathrm{sulf}}/r_{\mathrm{dry},\min }^{\mathrm{sulf}})}},\end{array}$$

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}\mathrm{g}{\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},\min }^{\mathrm{sulf}},r_{\mathrm{dry},\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

$$\begin{array}{}\text{(104)}& {\mathit{\xi }}_i={\displaystyle \frac{c^{\mathrm{dust}}}{c^{\mathrm{sulf}}}}{\displaystyle \frac{\text{d}{N}^{\mathrm{sulf}}}{\text{d}\log r_{\mathrm{dry}}^{\mathrm{sulf}}}}\left(\log r_{\mathrm{dry},i}^{\mathrm{sulf}}\right){\displaystyle \frac{\log (r_{\mathrm{dry},\max }^{\mathrm{sulf}}/r_{\mathrm{dry},\min }^{\mathrm{sulf}})}{c^{\mathrm{SP}}/\mathrm{2}}}{\displaystyle \frac{P_{\mathrm{INia}}}{P_{\mathrm{INia}}^{\mathrm{SP}}}}.\end{array}$$

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_{\min }^{\mathrm{fz}},T_{\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

$$\begin{array}{}\text{(105)}& \begin{array}{rl}{\mathit{\xi }}_i=& {\displaystyle \frac{c^{\mathrm{dust}}}{c^{\mathrm{sulf}}}}{\displaystyle \frac{d{N}^{\mathrm{sulf}}}{d\log r_{\mathrm{dry}}^{\mathrm{sulf}}}}\left(\log r_{\mathrm{dry},i}^{\mathrm{sulf}}\right)p\left(T_i^{\mathrm{fz}}\right)\\ & {\displaystyle \frac{\log (r_{\mathrm{dry},\max }^{\mathrm{sulf}}/r_{\mathrm{dry},\min }^{\mathrm{sulf}})(T_{\max }^{\mathrm{fz}}-T_{\min }^{\mathrm{fz}})}{(c^{\mathrm{SP}}/\mathrm{2})(\mathrm{1}-P_{\mathrm{INia}}^{\mathrm{SP}})}}.\end{array}\end{array}$$

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}\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}\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}\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}\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}\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}\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}\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}\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}\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}\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}\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}\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}\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}\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}\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}\mathrm{s}$.

7 Typical behavior of CTRL ensemble

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.

Figure 1 Typical realization of CTRL cloud spatial structures at t=2040, 2460, 3000, 4200, and 5400 s. The mixing ratio of cloud water, rainwater, cloud ice, graupel, and snow aggregates are plotted in fading white, yellow, blue, red, and green, respectively. The symbols indicate examples of unrealistic predicted ice particles (Sects. 7.3 and 9.1). See also Movie 1 in the video supplement.

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⁻³ at $(x,z)=(\mathrm{21.8}\mathrm{km},\mathrm{9.8}\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⁻¹, which were observed at $(t,x,z)=(\mathrm{2340}\mathrm{s},\mathrm{12.8}\mathrm{km},\mathrm{11.1}\mathrm{km})$ and $(\mathrm{1620}\mathrm{s},\mathrm{9.5}\mathrm{km},\mathrm{4.1}\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.

Figure 2 Time evolution of the domain-averaged water path in the CTRL ensemble. The cloud water, rainwater, cloud ice water, graupel water, and snow aggregate water paths are plotted in gray, yellow, blue, red, and green, respectively. The solid line represents the typical realization of CTRL. Dark shades indicate the mean ± standard deviation, and pale shades indicate the maximum and minimum values of the 10 ensemble members.

Figure 3 Time evolution of domain-averaged accumulated precipitation amounts in the CTRL ensemble. The solid line represents the typical realization of CTRL. The dark shading indicates the mean ± standard deviation, and the pale shading indicates the maximum and minimum values of the 10 ensemble members.

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}}$:

$$\begin{array}{}\text{(106)}& m^{*}:={\displaystyle \frac{{\mathit{\rho }}^{\mathrm{i}}{a}^{\mathrm{2}}c}{{\mathit{\rho }}_{\mathrm{true}}^{\mathrm{i}}(D/\mathrm{2}{)}^{\mathrm{3}}}},D=\mathrm{2}max(a,c).\end{array}$$

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 ₁₀(D) and log ₁₀(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 ρⁱ. 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}}\mathrm{kg}{\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.

Figure 4 Mass–dimension relationship of ice particles in the typical realization of CTRL at t=2040, 3000, and 5400 s. 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 and vertical axes represent the maximum ice particle dimension D and the normalized ice particle mass m^*, respectively. The colored slopes represent various mass–dimension relationship formulas. The symbols indicate examples of unrealistically predicted ice particles (Sects. 7.3 and 9.1). See also Movie 2 in the video supplement.

Figure 5 Aspect ratio–dimension relationship of ice particles in the typical realization of CTRL at t=2040, 3000, and 5400 s. The vertical axis represents the ice particle aspect ratio ϕ. This figure is the same as Fig. 4, except for the vertical axis. The symbols indicate examples of unrealistically predicted ice particles (Sects. 7.3 and 9.1). See also Movie 3 in the video supplement.

Figure 6 Apparent density–dimension relationship of ice particles in the typical realization of CTRL at t=2040, 3000, and 5400 s. The vertical axis represents the ice particle apparent density ρⁱ. This figure is the same as Fig. 4, except for the vertical axis. The symbols indicate examples of unrealistically predicted ice particles (Sects. 7.3 and 9.1). See also Movie 4 in the video supplement.

Figure 7 Velocity–dimension relationship of ice particles in the typical realization of CTRL at t=2040, 3000, and 5400 s. The vertical axis represents the ice particle terminal velocity v^∞. This figure is the same as Fig. 4, except for the vertical axis. The colored slopes represent various velocity–dimension relationship formulas. The symbols indicate examples of unrealistically predicted ice particles (Sects. 7.3 and 9.1). 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⁻³. H02's formula is for temperature and pressure of −10 ^∘C and 500 hPa; hence, the density is approximately 0.66 kg m⁻³. In our simulation, most of the snow aggregates exist in the anvil cloud, wherein the density is approximately 0.38 kg m⁻³. 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}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}/\mathrm{1.1}\mathrm{kg}{\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}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}/\mathrm{0.66}\mathrm{kg}{\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.

8 Numerical convergence characteristics

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.

Figure 8 Statistics of NSP ensemble accumulated precipitation amounts. The vertical axis represents the accumulated precipitation at the end of the simulation (t=5400 s), and the horizontal axis represents 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 crosses denote maximum and minimum values of the 10 ensemble members.

Figure 9 Statistics of NSP ensemble maximum water paths for each hydrometeor type. The vertical axis represents the maximum water path of each hydrometeor type during the simulation (i.e., the maximum of each line in Fig. 2), and the horizontal axis represents 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 symbols denote the maximum and minimum values of the 10 ensemble members.

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}\mathrm{m}$.

Figure 10 Statistics of DX ensemble accumulated precipitation amounts. The horizontal axis represents the grid size $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z$. This figure has the same form as Fig. 8, except for the horizontal axis.

Figure 11 Statistics of the DX ensemble maximum water paths for each hydrometeor type. The horizontal axis represents the grid size $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z$. This figure has the same form as Fig. 9, except for the horizontal axis.

Figure 12 Spatial structure of the cloud at t=3000 s from a DX/2 ensemble member. This figure is formatted the same as Fig. 1, except for the grid resolution. See also Movie 6 in the video supplement.

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.

Figure 13 Statistics of DT ensemble accumulated precipitation amounts. The horizontal axis represents the ratio of cloud microphysics time steps to CTRL. This figure has the same form as Fig. 8, except for the horizontal axis.

Figure 14 Statistics of the DT ensemble maximum water path for each hydrometeor type. The horizontal axis represents the ratio of cloud microphysics time steps to CTRL. This figure has the same form as Fig. 9, except for the horizontal axis.

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}\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⁵–100⁵. 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¹⁰–100¹⁰. 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⁵–100⁵. 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.

9 Improvement of the model

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}\mathrm{mm},\mathrm{36.1}\mathrm{mm},\mathrm{1.12}\times {\mathrm{10}}^{\mathrm{2}}\mathrm{kg}{\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}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ is also much smaller than that of typical hailstones (Fig. 7). It is located at $(x,z)=(\mathrm{20.0}\mathrm{km},\mathrm{4.11}\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}\mathrm{g}{\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}\mathrm{mm},\mathrm{11}\mathrm{mm})$.

However, ${\mathit{\rho }}_{\mathrm{rime}}=\mathrm{0.1}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$ must be reconsidered. Equation (56) has a global maximum of approximately 0.95 g cm⁻³ at around Y=3.7 and then quickly decreases, becoming $<\mathrm{0.1}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$ at around Y=5.5. Then, from Eq. (55), ${\mathit{\rho }}_{\mathrm{rime}}=\mathrm{0.1}\mathrm{g}{\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}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}/{}^{\circ }\mathrm{C}\right)$, ${\mathit{\rho }}_{\mathrm{rime}}=\mathrm{0.1}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$ frequently happens near the freezing level. For example, Y=1000 for r_k=1 mm, $v_{\mathrm{imp}}=\mathrm{1}\mathrm{m}{\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

$$\begin{array}{}\text{(107)}& Y^{\downarrow }=min(Y,\mathrm{3.5}).\end{array}$$

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}\mathrm{\mu }\mathrm{m},$ 138.8 µm, $\mathrm{269.7}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}},$ 0.29, 1}. Its terminal velocity is $v^{\mathrm{\infty }}=\mathrm{1.46}\times {\mathrm{10}}^{-\mathrm{2}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ and it is located at $(x,z)=(\mathrm{20.3}\mathrm{km},\mathrm{11.0}\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}}\mathrm{m}{\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}\mathrm{g}{\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}\mathrm{\mu }\mathrm{m},\mathrm{110}\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

$$\begin{array}{}\text{(108)}& a_j^{\prime }=max(a_j,r_k({\mathit{\rho }}^{\mathrm{w}}/{\mathit{\rho }}_{\mathrm{rime}}{)}^{\mathrm{1}/\mathrm{3}}),\end{array}$$

and Eq. (54) with

$$\begin{array}{}\text{(109)}& c_j^{\prime }=max(c_j,r_k({\mathit{\rho }}^{\mathrm{w}}/{\mathit{\rho }}_{\mathrm{rime}}{)}^{\mathrm{1}/\mathrm{3}}).\end{array}$$

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}\mathrm{mm},$ 15.0 mm, $\mathrm{10.7}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}},$ 0.85, 1 585 116}. Its terminal velocity is $v^{\mathrm{\infty }}=\mathrm{4.31}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ and it is located at $(x,z)=(\mathrm{10.6}\mathrm{km},\mathrm{11.5}\mathrm{km})$ (see Fig. 1). What is unusual here is the very low apparent density ${\mathit{\rho }}^{\mathrm{i}}=\mathrm{10.7}\mathrm{kg}{\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}\mathrm{kg}{\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⁻³ were reported in Locatelli and Hobbs (1974). Therefore, a hailstone with an apparent density of 10 kg m⁻³ 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}\mathrm{\mu }\mathrm{m},$ 438.2 µm, $\mathrm{19.4}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}},$ 0.53, 26 679}. Its terminal velocity is $v^{\mathrm{\infty }}=\mathrm{1.02}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, it is located at $(x,z)=(\mathrm{29.7}\mathrm{km},\mathrm{6.4}\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,

$$\begin{array}{}\text{(110)}& \mathrm{\Gamma }\left(T_i\right)=\mathrm{1},\text{for}\text{d}{m}_i<\mathrm{0}\text{(sublimation)};\end{array}$$

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).

Figure 15 This figure is the same as Fig. 4 but shows results from SCALE-SDM 0.2.5-2.2.1, which incorporates the three corrections (Eqs. 107–110) proposed to avoid the creation of ice particles with unrealistic morphologies. See also Movie 8 in the video supplement.

Figure 16 This figure is the same as Fig. 5 but shows results from SCALE-SDM 0.2.5-2.2.1, which incorporates the three corrections (Eqs. 107–110) proposed to avoid the creation of ice particles with unrealistic morphologies. See also Movie 9 in the video supplement.

Figure 17 This figure is the same as Fig. 6 but shows results from SCALE-SDM 0.2.5-2.2.1, which incorporates the three corrections (Eqs. 107–110) proposed to avoid the creation of ice particles with unrealistic morphologies. See also Movie 10 in the video supplement.

Figure 18 This figure is the same as Fig. 7 but shows results from SCALE-SDM 0.2.5-2.2.1, which incorporates the three corrections (Eqs. 107–110) proposed to avoid the creation of ice particles with unrealistic morphologies. See also Movie 11 in the video supplement.

Figure 19 Changes in the domain-averaged water path before and after corrections. The long dashed, solid, and short dashed lines represent the SCALE-SDM 0.2.5-2.2.0, -2.2.1, and -2.2.2, respectively.

Figure 20 Changes in accumulated precipitation amounts before and after corrections. The long dashed, solid, and short dashed lines represent the SCALE-SDM 0.2.5-2.2.0, -2.2.1, and -2.2.2, respectively.

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öhm}}^{\mathrm{\infty }}(m_i,{\mathit{\varphi }}_i,d_i,q_i;{\mathit{\rho }}_i,T_i)$ can be summarized as follows:

$$\begin{array}{ll}\text{(111)}& {\displaystyle }& {\displaystyle }X={\displaystyle \frac{\mathrm{8}{m}_{i}g\mathit{\rho }}{\mathit{\pi }{\mathit{\mu }}^{\mathrm{2}}max({\mathit{\varphi }}_i,\mathrm{1})q_i^{\mathrm{1}/\mathrm{4}}}},\text{(112)}& {\displaystyle }& {\displaystyle }{X}^{\prime }=X{\displaystyle \frac{\mathrm{1}+(X/X_{\mathrm{0}}{)}^{\mathrm{2}}}{\mathrm{1}+\mathrm{1.6}(X/X_{\mathrm{0}}{)}^{\mathrm{2}}}},\text{(113)}& {\displaystyle }& {\displaystyle }{X}_{\mathrm{0}}=\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{6}},\text{for ice particles},{\displaystyle }& {\displaystyle }k=min\mathit{\{}max(\mathrm{0.82}+\mathrm{0.18}{\mathit{\varphi }}_i,\mathrm{0.85}),\left(\mathrm{0.37}+{\displaystyle \frac{\mathrm{0.63}}{\sqrt{{\mathit{\varphi }}_i}}}\right),\\ \text{(114)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1.33}}{max(\log {\mathit{\varphi }}_i,\mathrm{0})+\mathrm{1.19}}}\mathit{\}},\text{(115)}& {\displaystyle }& {\displaystyle }\mathrm{\Gamma }=max\mathit{\{}\mathrm{1},min(\mathrm{1.98},\mathrm{3.76}-\mathrm{8.41}{\mathit{\varphi }}_i+\mathrm{9.18}{\mathit{\varphi }}_i^{\mathrm{2}}-\mathrm{3.53}{\mathit{\varphi }}_i^{\mathrm{3}}\left)\mathit{\right\}},\text{(116)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{DP}}=max(\mathrm{0.292}k\mathrm{\Gamma },\mathrm{0.492}-\mathrm{0.200}/\sqrt{{\mathit{\varphi }}_i}),\text{(117)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{DO}}=\mathrm{4.5}{k}^{\mathrm{2}}max({\mathit{\varphi }}_i,\mathrm{1}),\text{(118)}& {\displaystyle }& {\displaystyle }\mathit{\beta }={\left[\mathrm{1}+{\displaystyle \frac{C_{\mathrm{DP}}}{\mathrm{6}k}}{\left({\displaystyle \frac{X^{\prime }}{C_{\mathrm{DP}}}}\right)}^{\mathrm{1}/\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}-\mathrm{1},\text{(119)}& {\displaystyle }& {\displaystyle }\mathit{\gamma }={\displaystyle \frac{C_{\mathrm{DO}}-C_{\mathrm{DP}}}{\mathrm{4}{C}_{\mathrm{DP}}}},\text{(120)}& {\displaystyle }& {\displaystyle }{N}_{\mathrm{Re}}={\displaystyle \frac{\mathrm{6}k}{C_{\mathrm{DP}}}}{\mathit{\beta }}^{\mathrm{2}}\left[\mathrm{1}+{\displaystyle \frac{\mathrm{2}\mathit{\beta }{e}^{-\mathit{\beta }\mathit{\gamma }}}{(\mathrm{2}+\mathit{\beta })(\mathrm{1}+\mathit{\beta })}}\right],\text{(121)}& {\displaystyle }& {\displaystyle }{v}_{\text{Böhm}}^{\mathrm{\infty }}={\displaystyle \frac{\mathit{\mu }{N}_{\mathrm{Re}}}{\mathit{\rho }{d}_i}}.\end{array}$$

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

$$\begin{array}{}\text{(122)}& d_i=\mathrm{2}{a}_i,q_i=q_i^{\mathrm{ce}}=A_i/A_i^{\mathrm{ce}};\end{array}$$

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=2a_i.

For planar ice particles (ϕ_i<1), $v_{\text{Böhm}}^{\mathrm{\infty }}(d_i=\mathrm{2}{a}_i,q_i=q_i^{\mathrm{ce}})$ and $v_{\text{Böhm}}^{\mathrm{\infty }}(d_i=D_i,q_i=q_i^{\mathrm{cc}})$ yield the same results, because 2a_i=D_i and $q_i^{\mathrm{ce}}=q_i^{\mathrm{cc}}$ hold. However, for columnar ice particles (ϕ_i>1), $v_{\text{Böhm}}^{\mathrm{\infty }}(D_i,q_i^{\mathrm{cc}})$ always underestimates the fall velocity. From Eqs. (111)–(121), we can derive $v_{\text{Böhm}}^{\mathrm{\infty }}(\mathrm{2}{a}_i,q_i^{\mathrm{ce}})/v_{\text{Böhm}}^{\mathrm{\infty }}(D_i,q_i^{\mathrm{cc}})={\mathit{\varphi }}_i^{\mathrm{3}/\mathrm{4}}$ for X≪1, and $v_{\text{Böhm}}^{\mathrm{\infty }}(\mathrm{2}{a}_i,q_i^{\mathrm{ce}})/v_{\text{Böhm}}^{\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öhm}}^{\mathrm{\infty }}(\mathrm{2}{a}_i,q_i^{\mathrm{ce}})/v_{\text{Böhm}}^{\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öhm}}^{\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:

$$\begin{array}{}\text{(123)}& {\mathrm{\Gamma }}^{*}=\mathrm{1}\text{for}\text{d}{m}_i\ge \mathrm{0}\bigwedge {\mathit{\varphi }}_i>\mathrm{40}.\end{array}$$

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ⁱ. 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}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ and $k=\mathrm{2.55}\times {\mathrm{10}}^{-\mathrm{2}}\mathrm{J}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}{\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}}{\mathrm{m}}^{\mathrm{2}}{\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ⁱ, and polycrystals can also be created. Therefore, for $T<-\mathrm{20}$ ^∘C, Γ might better be considered a function of both T and Sⁱ, 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ⁱ is lost. Because ρ_dep accounts for the secondary growth habit, dependence on Sⁱ 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 ρⁱ. 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}}$, Pruppacher and Klett (1997, Sect. 16.2) provides a simple formula that depends solely on temperature. The $E_{\mathrm{agg}}^{\mathrm{stick}}$ formula provided by Phillips et al. (2015) is physically based and should thus be more reliable.

Calculating the attributes of the resultant ice particles is also not easy. Let $(a^{\prime },c^{\prime },{\mathit{\rho }}^{\mathrm{i}\prime },m^{\mathrm{rime}\prime },n^{\mathrm{mono}\prime })$ be the ice particle created by the aggregation of $(a_{\mathrm{1}},c_{\mathrm{1}},{\mathit{\rho }}_{\mathrm{1}}^{\mathrm{i}},m_{\mathrm{1}}^{\mathrm{rime}},n_{\mathrm{1}}^{\mathrm{mono}})$ and $(a_{\mathrm{2}},c_{\mathrm{2}},{\mathit{\rho }}_{\mathrm{2}}^{\mathrm{i}},m_{\mathrm{2}}^{\mathrm{rime}},n_{\mathrm{2}}^{\mathrm{mono}})$. For rime mass and number of monomers, $m^{\mathrm{rime}\prime }=m_{\mathrm{1}}^{\mathrm{rime}}+m_{\mathrm{2}}^{\mathrm{rime}}$ and $n^{\mathrm{mono}\prime }=n_{\mathrm{1}}^{\mathrm{mono}}+n_{\mathrm{2}}^{\mathrm{mono}}$ hold. To determine the remainder, $(a^{\prime },c^{\prime },{\mathit{\rho }}^{\mathrm{i}\prime })$, specifying two out of the three attributes is sufficient because of the conservation of the total mass. In this study, as in the case of riming, we assumed that the filling-in model can be applied to aggregation; i.e., the maximum dimension is conserved and only the minor axis grows. Therefore, $D^{\prime }=max(D_{\mathrm{1}},D_{\mathrm{2}})=max(a_{\mathrm{1}},c_{\mathrm{1}},a_{\mathrm{2}},c_{\mathrm{2}})$. However, one more attribute must be specified. In this study, instead of predicting minor axis growth, we predict the apparent density ρ^i′ by introducing an intuitive model that considers the compaction of fluffy snowflakes. Consequently, the fractal dimension of the mass–dimension relationship of snow aggregates predicted by our model is close to 2 (see the green shading in Figs. 4 and 15), which agrees well with various previous studies (e.g., Brown and Francis, 1995; Heymsfield et al., 2010; Mitchell, 1996; Schmitt and Heymsfield, 2010).

However, the filling-in assumption is not valid for aggregation. Higuchi (1960) introduced a parameter called the separation ratio: $s:=\mathrm{2}l/(D_{\mathrm{1}}+D_{\mathrm{2}})$, $s\in [\mathrm{0},\mathrm{1}]$, where l is the horizontal distance between the centers of the two particles. For an aggregation between two planar ice particles, the resultant ice particle's maximum dimension can be evaluated by $D^{\prime }=max\mathit{\{}{D}_{\mathrm{1}},D_{\mathrm{2}},(\mathrm{1}+s\left)\right(D_{\mathrm{1}}+D_{\mathrm{2}})/\mathrm{2}\mathit{\}}$. Our model corresponds to the special case of s=0, but it has been reported that s≈0.5–0.6 for two planar crystals and dendrites (Higuchi, 1960; Kajikawa and Heymsfield, 1989; Kajikawa et al., 2002), and s≈0.9 for spatial dendrites (Kajikawa et al., 2002). In contrast, s=0 for columnar ice crystals can be justified from Kajikawa (1995)'s observation that two needles of similar sizes tend to attach with their centers close (s≈0) and a right angle between their polar axes (crossed adhesion). Notably, the cross-adhesion displacement gives the largest possible volume V_max, which we used to calculate the apparent density ρ^i′ of the resultant ice particle by interpolation.

Another issue of the filling-in assumption is that it gradually makes snow aggregates quasi-spherical (see the green shading in Figs. 5 and 16). Measurements indicate that snow aggregates have an average aspect ratio of 0.6 (e.g., Korolev and Isaac, 2003) or smaller (Jiang et al., 2017).

Introducing the separation ratio s in our model is straightforward and could improve our model's accuracy. In general, this tends to reduce the mass–dimension relationship's fractal dimension, and their aspect ratio. Locatelli and Hobbs (1974) reported that aggregates of dendrites and aggregates of unrimed side planes had fractal dimensions of 1.4 (plotted in Figs. 4 and 15), which is smaller than 2.

In our model, the apparent density ρ^i′ after aggregation is predicted by the formula given in Eq. (62). It is natural to assume that there is a lower limit of apparent density; however, this is a crude expression of the idea and requires further validation and improvement. Also note that a contact angle model was used in Chen and Lamb (1994b) and Hashino and Tripoli (2011a) to determine the resultant ice particle.

Several numerical models can create detailed 3-D structures of snow aggregates consisting of primary ice crystals (e.g., Westbrook et al., 2004a, b; Maruyama and Fujiyoshi, 2005; Schmitt and Heymsfield, 2010). We can refine our aggregation outcome model by using the results of those more microscopic models that resolve snow aggregate structures. For example, Przybylo et al. (2019) and Dunnavan et al. (2019) intensively studied the geometry of aggregates using such numerical models.

9.3.9 Spontaneous/collisional breakup

Several mechanisms can induce the spontaneous/collisional breakup of hydrometeors. However, we did not consider any of them in the present study. In particular, rime splintering (Findeisen and Findeisen, 1943; Hallett and Mossop, 1974), and the collisional breakup of ice particles (Vardiman, 1978) are critical in mixed-phase clouds, as these processes are thought to be responsible for the large excess in the observed number concentration of ice particles to the number concentration of IN aerosol particles (e.g., Field et al., 2017).

First, a particle-based numerical algorithm for calculating spontaneous/collisional breakup processes has not yet been established. A simple strategy is to add more super-particles to the system when a breakup event occurs, but this could be computationally inefficient.

Mathematical models of spontaneous/collisional breakup processes are available from various studies. For the spontaneous breakup of rain droplets >6.5 mm, we can use the mathematical model from Kamra et al. (1991). For the collisional breakup of droplets, the models compiled and compared in Prat et al. (2012) can be used. For the shedding of excess meltwater, Phillips et al. (2007)'s model can be used. For rime splintering, the model summarized in Sect. 16.1.6 of Pruppacher and Klett (1997) can be used. Readers may also refer to Field et al. (2017) and the references cited therein. For the collisional breakup of ice particles, Phillips et al. (2017)'s model can be used.

9.3.10 Subgrid-scale turbulence

The grid size we tested for evaluating the model ranged from 31.25 to 250 m, and only flows that are larger than the chosen grid size can be resolved. A substantial portion of turbulence kinetic energy is accumulated in large scales, and small-scale turbulence is mostly driven by large-scale motions; therefore, SGS turbulence is of secondary importance to the phenomena. Nevertheless, SGS turbulence does affect moist air flow and atmospheric particle behavior. SGS turbulence should be appropriately incorporated to improve the model's grid convergence.

The Smagorinsky–Lilly model (Smagorinsky, 1963; Lilly, 1962; Brown et al., 1994; Scotti et al., 1993), which is already available in SCALE-SDM, can be used for the diffusion of moist air by SGS turbulence. However, we did not use it in this study because the model is designed for 3-D turbulence.

SGS turbulence can enhance particle collision, which can be incorporated by using the collision kernels proposed in Wang et al. (2008), Onishi and Seifert (2016), and Chen et al. (2018). Particle velocity fluctuations due to SGS turbulence can be modeled as an Ornstein–Uhlenbeck process (e.g., Pope, 1994; Schilling et al., 1996; Grabowski and Abade, 2017). The fluctuation of supersaturation through eddy-hopping and entrainment can be considered by introducing a new stochastic attribute (Grabowski and Abade, 2017; Abade et al., 2018) or by applying the linear eddy model to particles (Hoffmann et al., 2019).

10 Conclusions

Using SDM, we constructed a detailed numerical model of mixed-phase clouds based on a kinetic description, and subsequently demonstrated that a large-eddy simulation of a cumulonimbus that predicts ice particle morphology without assuming ice categories or mass–dimension relationships is possible. Our results strongly support the particle-based modeling methodology's efficacy for simulating mixed-phase clouds.

In our model, ice particles are approximated by porous spheroids. The elementary cloud microphysics processes that the model considers include advection and sedimentation; immersion/condensation and homogeneous freezing; melting; condensation and evaporation including the activation and deactivation of CCN; deposition and sublimation; and coalescence, riming, and aggregation. Moist air fluid dynamics is described using the compressible Navier–Stokes equation.

Our model successfully simulated the life cycle of a cumulonimbus, and the predicted mass–dimension and velocity–dimension relationships were comparable with existing formulas. Numerical convergence was achieved at a super-particle number concentration as low as 128 per cell, which consumed 30 times more computational time than a typical two-moment bulk model. We then fixed several issues of the original model and developed two updated versions: SCALE-SDM 0.2.5-2.2.1 (fix of the odd ice particle creation) and SCALE-SDM 0.2.5-2.2.2 (fix of the underestimated columnar ice terminal velocity).

A more detailed evaluation of the model to explore the applicability of the new approach is an essential step forward. Our results strongly indicate that ice particle morphology can be predicted more accurately by further developing particle-based models. However, from this study, we cannot quantify the extent to which the refined representation of mixed-phase cloud microphysics could improve the predictability of mixed-phase clouds' macroscopic properties. Such proficiency can be addressed by conducting a thorough comparison with observations and other models. In addition, further sophistication of the model is necessary. As discussed in Sect. 9.3, various elementary processes must be incorporated or refined in the model. In particular, rime splintering and the collisional breakup of ice particles are critical because these processes are thought to be responsible for secondary ice production. Therefore, establishing an accurate and efficient particle-based algorithm for spontaneous/collisional breakup is also crucial.

Particle-based model accuracy is more subject to cloud microphysics uncertainties than numerical errors. Therefore, a quantitative understanding of elementary cloud microphysics processes is becoming increasingly important. More laboratory, observational, and theoretical studies to advance our knowledge of cloud microphysics are desired in the future (Morrison et al., 2020). Additionally, we can go into a more microscopic description of cloud microphysics than kinetic description, i.e., to explicitly resolve droplet and ice particle shapes and deterministically consider their collisions (e.g., Demange et al., 2017; Wang and Ji, 2000; Westbrook et al., 2004a, b; Maruyama and Fujiyoshi, 2005; Schmitt and Heymsfield, 2010; Mazloomi Moqaddam et al., 2015). Such model studies would also be useful for refining kinetic descriptions.

Our model's computational cost is at least 1 or 2 orders of magnitude larger than that of bulk models. To further accelerate calculation, the use of SGS models discussed in Sect. 9.3.10 is crucial. Further reduction of the computational cost could also be achieved by using the Twomey super-droplet methodology described in Grabowski et al. (2018); however, it is vital to introduce dynamic load balancing. The acceleration achieved by those improvements might be insufficient to allow using particle-based cloud microphysics models in weather or climate models. Studies to construct a high-fidelity bulk model or another form of macroscopic cloud microphysics model must also be pursued (e.g., Noh et al., 2018; Morrison et al., 2020).

Appendix A: List of symbols

Table A1 summarizes important variables used in this study.

Table A1 List of symbols.

Appendix B: List of abbreviations

Table B1 summarizes important abbreviations that are used in this study.

Auer (1972)Beard and Grover (1974)Chen and Lamb (1994a)Erfani and Mitchell (2017)Heymsfield et al. (2002)Heymsfield (1972)Heymsfield and Kajikawa (1987)Heymsfield and Pflaum (1985)Kajikawa (1989)Knight and Heymsfield (1983)Locatelli and Hobbs (1974)Mitchell et al. (1990)Mitchell (1996)Starr and Cox (1985)Westbrook et al. (2008)

Table B1 List of abbreviations.