Articles | Volume 13, issue 9
https://doi.org/10.5194/gmd-13-4107-2020
https://doi.org/10.5194/gmd-13-4107-2020
Development and technical paper
 | Highlight paper
 | 
08 Sep 2020
Development and technical paper | Highlight paper |  | 08 Sep 2020

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

Shin-ichiro Shima, Yousuke Sato, Akihiro Hashimoto, and Ryohei Misumi
Abstract

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 Pinsky2018; 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 Milbrandt2015; Khain and Pinsky2018; 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., Beheng2010). 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 Klett1997; Hashino and Tripoli2007, 2008, 2011a, b; Khvorostyanov and Curry2014; Khain and Pinsky2018). 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 Harrington2015). 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 Ramanadham1954; Hardy1963; Srivastava1967). 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 Pfister2004; Shirgaonkar and Lele2006; Andrejczuk et al.2008, 2010; Shima et al.2009; Sölch and Kärcher2010; Riechelmann et al.2012; Brdar and Seifert2018; Seifert et al.2019; Jaruga and Pawlowska2018; Grabowski and Abade2017; 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 Pawlowska2017). 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 Rutland2000) 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 (Smoluchowski1916; Alfonso and Raga2017; Dziekan and Pawlowska2017). 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 −20C. 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 Milbrandt2015; Milbrandt and Morrison2016) 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. 24, 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 a={r,{mαsol},{mβinsol},Tfz,a,c,ρi,mrime,nmono,v}, 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 {a,c,ρi} 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 {a,c,ρi} 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)πr3.

2.3 Masses of soluble and insoluble substances

Let mαsol, α=1,2,,Nsol be the masses of soluble substances contained in the particle, and let mβinsol, β=1,2,,Ninsol 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” (Levine1950), we consider that each insoluble particle has its own freezing temperature Tfz, and that a supercooled droplet freezes as soon as the ambient temperature T decreases below Tfz. The freezing process is described in detail in Sect. 4.1.4.

Each particle's Tfz is directly connected to the ice nucleation active surface site (INAS) density concept (e.g., Fletcher1969; Connolly et al.2009; Niemand et al.2012; Hoose and Möhler2012).

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 nS(T) gives the accumulated number of INAS per unit surface area of the insoluble substance. Therefore, nS(T) is a function that increases as T decreases. The freezing temperature Tfz corresponds to the highest temperature at which the first INAS appears on the insoluble substance's surface. Let Ainsol be the insoluble substance's surface area. Then, the probability that Tfz is larger than T can be calculated as P(Tfz>T)=1-exp[-AinsolnS(T)]. The probability density function of Tfz then becomes

(1) p ( T ) = - d P ( T fz > T ) d T = - A insol d n S d T e - A insol n S .

We can determine Tfz 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 Tfz is the highest of all.

It is possible that a single INAS does not appear until −38C, 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 Tfz=-38C. If a particle contains only soluble substances, we also set Tfz=-38C.

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 Lee1966; Kikuchi et al.2013). Following the strategies of Chen and Lamb (1994a, b), Misumi et al. (2010), and Jensen and Harrington (2015), let us approximate each ice particle as a porous spheroid, which is characterized by three variables, namely equatorial radius a, polar radius c, and apparent density ρi. That is, the ice particle's apparent volume is V=(4π/3)a2c, and its mass can be evaluated as m=ρiV. The two radii a and c represent the ice particle's spatial extent and ρi represents its internal structure. Let us define the aspect ratio as ϕ:=c/a. A spheroid is considered a prolate spheroid if ϕ>1, and columns could be approximated by prolate spheroids. In contrast, plates and dendrites are approximated by oblate spheroids, i.e., ϕ<1. If an ice particle is hollowed out or intricately branched, ρi becomes smaller than the ice crystal's true density ρtruei916.8 kg m−3.

2.6 Rime mass and number of monomers

Following Brdar and Seifert (2018) we introduce two additional ice particle attributes, namely rime mass mrime and number of monomers nmono. Rime mass mrime records the mass of ice a particle has obtained through the riming process. The number of monomers nmono is an integer representing the number of primary ice crystals in the particle. In this study, mrime and nmono 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 mrime or nmono, 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 a={r,{mαsol},{mβinsol},Tfz,a,c,ρi,mrime,nmono,v}. We need the mass of insoluble substances {mβinsol,β=1,2,,Ninsol} (and corresponding INAS densities) to specify freezing temperature Tfz. However, as described in Sect. 4.1, time evolution equations do not depend on {mβinsol}. Rime mass mrime and the number of monomers nmono do not affect time evolution either. Particle velocity v is a diagnostic attribute. Therefore, the attributes directly relevant to time evolution are reduced to {r,{mαsol},Tfz,a,c,ρi}. 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 U=(U,V,W), density of dry air ρd, density of water vapor ρv, density of moist air ρ:=ρd+ρv, specific humidity qv:=ρv/ρ, mass of dry air per unit mass of moist air qd:=ρd/ρ, temperature T, pressure P, and potential temperature of moist air θ:=T/Π:=T/(P/P0)R/cp. Here, P0=1000 hPa is a reference pressure; Rd, Rv, and R:=qdRd+qvRv are the gas constants of dry air, water vapor, and moist air, respectively; and cpd, cpv, and cp:=qdcpd+qvcpv 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: G:={U,ρ,qv,θ,P,T}.

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 {{xi(t),ai(t)},i=1,2,,Nrwp}. Here, Nrwp 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 Ir(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

(2) d d t ( m i v i ) = F i drg - m i g z ^ , d x i d t = v i ,

where mi is the particle's mass, Fidrg is the force of drag from moist air, g is Earth's gravity, and z^ is the unit vector in the z-axis direction. Note that -Fidrg 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

(3) v i = U i - z ^ v i , d x i d t = v i ,

where Ui:=U(xi) is the ith particle's ambient wind velocity, and vi is the terminal velocity, which is a function of attributes ai and the state of the ambient air Gi.

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 Pruppacher1977) though that of micrometer-sized droplets is approximately 10−5 s (see, e.g., Eq. 1 of Chen et al.2018, and the discussion that follows). The acceleration of particles can be considered by explicitly solving the motion equation (see, e.g., Naumann and Seifert2015), 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): vi=vBeard(min(ri,3.5mm);ρi,Pi,Ti), where ρi:=ρ(xi) and Pi:=P(xi) are the density and pressure of ambient moist air, respectively. This formula applies to droplets with radii smaller than 3.5mm. 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.5mm 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): vi=vBöhm(mi,ϕi,di,qi;ρi,Ti), where di is the characteristic length, and qi is the area ratio.

In Böhm's theory, di is defined by 2ai, and qi is defined by the area ratio regarding circumscribed ellipse qice:=Ai/Aice, where Ai is the projected area perpendicular to the flow direction, and Aice is the area of the circumscribed ellipse of Ai, i.e., the area of the smallest ellipse that completely contains Ai.

However, in this study, we start from a slightly different definition of di and qi, which we adopted mistakenly:

(4) d i = D i := 2 max ( a i , c i ) , q i = q i cc := A i / A i cc ,

where Di is the maximum dimension, qicc is the area ratio regarding circumcircle, and Aicc is the area of the circumcircle of Ai, i.e., the area of the smallest circle that completely contains Ai.

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 Aicc=πmax(ai,ci)2. The projected area Ai can be roughly evaluated by the area of the circumscribed ellipse Aice=πaimax(ai,ci); however, we must subtract pores and indentations at boundaries from Aice. We assume that the ratio Ai/Aice is a power of the volume fraction ρii/ρtruei, and that the exponent κ is a function of the aspect ratio ϕi:

(5) A i = A i ce ρ i i ρ true i κ ( ϕ i ) .

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

(6) κ ( ϕ i ) = exp ( - ϕ i ) .

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 Ai/Aice is equal to the volume fraction ρii/ρtruei. Thus, κ(ϕi=0)=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 ρii/ρtruei will not affect the ratio Ai/Aice. Thus, κ(ϕi)=0.

For ϕi≈1, Jensen and Harrington (2015) argued that (ρii/ρtruei)κ=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 miDiβ and AiDiβ/s hold. Thus, by the definition of apparent density, ρii=mi/((4/3)πai2ci)Diβ-3. From Eq. (5), Diβ/s=Di2Di(β-3)κ. Hence, κ=(2s-β)/{s(3-β)} holds. Schmitt and Heymsfield (2010) estimated that (β,s)=(2.22,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 (β,s)=(2.20,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., ri>0; (2) the ambient water vapor is supersaturated over liquid water, i.e., ei>esw(Ti); and (3) the ambient temperature is lower than the particle's freezing temperature, i.e., Ti<Tifz. Here, ei:=e(xi) and Ti:=T(xi) are the ambient vapor pressure and temperature of the ith particle, respectively, and esw(T) 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 ρtruei. Therefore, attributes are initiated as follows: ri=0, ai=ci=ri(ρw/ρtruei)1/3, ρii=ρtruei, nimono=1, and mirime=0. The primed variables here denote values after the update, and ρw is the density of liquid water. {mαisol}, {mβiinsol}, and Tifz 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 0C, we consider that melting occurs immediately. Thus, the attributes are updated as follows: ri=(ai2ciρii/ρw)1/3 and ai=ci=ρii=nimono=mirime=0. {mαisol}, {mβiinsol}, and Tifz 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

(7) d m i d t = 4 π r i D v ( ρ v i - ρ v i sfc ) .

Here, Dv is water vapor's diffusivity in air, ρvi:=ρv(xi) is the ambient moist air's water vapor density, and ρvisfc 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 Tisfc and ambient temperature Ti are close to each other, i.e., (Tisfc-Ti)/Ti1, Eq. (7) can be reduced to

(8) r i d r i d t = 1 ρ w ( F k w + F d w ) S i w - e s i w , eff e s w ( T i ) ,

where Siw:=ei/esw(Ti) is the ambient saturation ratio over liquid water, and

(9) F k w = L v R v T i - 1 L v k T i , F d w = R v T i D v e s w ( T i ) ,

where Lv is the latent heat of vaporization, k is the thermal conductivity of moist air, and esiw,eff is the effective saturation vapor pressure regarding the ith droplet's surface. Following Köhler's theory (Köhler1936), an approximate formula of esiw,eff can be derived as

(10) e s i w , eff e s w ( T i ) = 1 + a ( T i ) r i - b m α i sol r i 3 ,

where a3.3×10-5cmK/Ti, b4.3cm3αIαmαisol/Mαsol, Iα is the van 't Hoff factor, which represents the degree of ionic dissociation, and Mα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 Klett1997). 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 Shima2017; Hoffmann2017; 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., Nakaya1954; Hallett and Mason1958; Kobayashi1961). 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):

(11) d m i d t = 4 π C D v ( ρ v i - ρ v i sfc ) f vnt = 4 π C S i i - 1 F k i + F d i f vnt ,

where Sii:=ei/esi(Ti) is the ambient saturation ratio over ice, and esi(T) is the saturation vapor pressure over ice at temperature T,

(12) F k i = L s R v T i - 1 L s k T i , F d i = R v T i D v e s i ( T i ) ,

where Ls is the latent heat of sublimation, C=C(ai,ci) is the electric capacitance of the spheroid, and fvnt is the particle-averaged ventilation coefficient.

The exact form of capacitance C(ai,ci) is given by Chen and Lamb (1994a). C(2ai+ci)/3 gives a good approximation for ϕi≈1.

The coefficient fvnt accounts for the ventilation effect, i.e., the enhancement of diffusional growth by air flow. Hall and Pruppacher (1976) suggested that fvnt could be described by

(13) f vnt = b 1 + b 2 X γ ,

where (b1,b2,γ)=(1.0,0.14,2) for X≤1, (b1,b2,γ)=(0.86,0.28,1) for X>1, X=NSc1/3(NReii)1/2, NSc=μ/(ρDv) is the Schmidt number, NReii=ρviDi/μ is the Reynolds number of ice particle i, and μ is the dynamic viscosity of moist air.

Note that mi 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 dmi as follows:

(14) d m i = max ( d m i , m min i - m i ) ,

where mmini 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 ρtruei. This is a crude representation of pre-activation (see, e.g., Marcolli2017, for a review). Each particle keeps the memory of ice activation until the ambient temperature rises above 0C. A particle with mmini ice grows immediately after the ambient air is supersaturated over ice, irrespective of its freezing temperature Tifz.

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:

(15) d c i d a i = Γ ( T i ) f vnt c i a i = : Γ * c i a i ,

where fvnt is the primary growth habit's ventilation coefficient, and Γ* is the effective inherent growth ratio, including the ventilation effect.

For purely diffusional growth, dci/dai=ci/ai 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 (Baran2012; Korolev and Isaac2003; 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

(16) Γ ( T ) = Γ ( - 30 C ) 1.28 , for T < - 30 C .

The ventilation coefficient fvnt 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 fvnt of the form

(17) f vnt = b 1 + b 2 X γ c i / C 1 / 2 b 1 + b 2 X γ a i / C 1 / 2 .

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

(18) d V i = d m i ρ dep , for d m i 0 (deposition) .

Deposition density ρdep can be expressed as

(19) ρ dep = ρ true i , for Γ ( T i ) < 1 a i < 100 µ m ; ρ dep CL 94 , otherwise .

Here, following Jensen and Harrington (2015), we assume that planar crystal branching does not occur if the equatorial radius ai is smaller than 100 µm. ρdepCL94 is an empirical formula of deposition density proposed by Chen and Lamb (1994a),

(20) ρ dep CL 94 = ρ true i exp - 3 max ( Δ ρ i - 0.05 g m - 3 , 0 ) Γ ( T i ) g m - 3 ,

where Δρi:=ρvi-ρvisfc. From Eq. (11), Δρi becomes

(21) Δ ρ i = S i i - 1 D v ( F k i + F d i ) .

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

(22) ( Δ ρ i ) = min ( S i i , e s w ( T i ) / e s i ( T i ) ) - 1 D v ( F k i + F d i ) .

For sublimation, the particle volume change dVi is given by

(23) d V i = d m i ρ sbl , for d m i < 0 (sublimation) ,

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

(24) ρ sbl = ρ i i .

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

(25) ρ i i ( t + d t ) = m i + d m i V i + d V i ,

where dmi is given in Eqs. (11) and (14), and dVi is given in Eqs. (18) and (23).

From Eq. (15) and the definition of volume Vi=(4π/3)ai2ci, after dt, the two radii become

(26)ai(t+dt)=aiexpdlogVi2+Γ*,(27)ci(t+dt)=ciexpΓ*dlogVi2+Γ*.

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 (Baran2012; Korolev and Isaac2003; 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{ai(t+dt),ci(t+dt)}<1µm,

(28)ρii(t+dt)=ρtruei,(29)ai(t+dt)=ci(t+dt)=mi+dmi(4π/3)ρii(t+dt)13,

where primed variables indicate values after correction.

For simplicity, we assume that the rime mass fraction mirime/mi does not change through sublimation, following Brdar and Seifert (2018):

(30) m i rime ( t + d t ) = m i rime , for d m i 0 ; m i rime m i + d m i m i , for d m i < 0 .

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 Pawlowska2017). 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+dt) is given by

(31) P j k = K ( a j , a k ; G ) d t Δ V ,

where the function K(aj,ak;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

(32) K coal = E coal ( r j , r k ) π ( r j + r k ) 2 | v j - v k | ,

where Ecoal(rj,rk) is the collection efficiency of collision–coalescence, which can be decomposed into Ecoal=EcoalcollisEcoalcoal. Here, collision efficiency Ecoalcollis 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 Ecoalcoal represents the fraction of collisions that result in permanent coalescence. In this study, we assume Ecoalcoal=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:

(33)rj=(rj3+rk3)13,(34)mαjsol=mαjsol+mαksol,α=1,2,,Nsol(35)mβjinsol=mβjinsol+mβkinsol,β=1,2,,Ninsol(36)Tjfz=max(Tjfz,Tkfz),

where primed values indicate the resultant droplet. Here, we assumed that the resultant particle's Tjfz is given by max(Tjfz,Tkfz), 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

(37) K rime = E rime A g | v j - v k | ,

where Erime is the collision–riming collection efficiency and Ag is the geometric cross-sectional area of j and k.

Figure 1 of Wang and Ji (2000) defines Ag for riming, but calculating it rigorously for porous spheroid models is impossible. Thus, we approximate Ag by

(38) A g = π ( a j + r k ) { max ( a j , c j ) + r k } - ( A j ce - A j ) ;

i.e., the indentation of the ice particle (Ajce-Aj) is subtracted from the area of an ellipse with semi-axes (aj+rk) and {max(aj,cj)+rk}. Therefore, if rkaj,cj, then AgAj. At the other extreme, if rkaj,cj, then Agπ(aj+rk){max(aj,cj)+rk}.

To evaluate collision–riming collection efficiency Erime, we combine formulas proposed by Beard and Grover (1974) and Erfani and Mitchell (2017).

If vj<vk, we consider droplet k as the collector and adopt the formula of Beard and Grover (1974):

(39) E rime = E BG 74 ( p i / w , N Re k w , N St i / w ) ,

where pi/w:=rji/rk, rji:=(aj2cj)1/3, NRekw=ρvk2rk/μ is the Reynolds number of droplet k, NSti/w=(pi/w)2ρjiNRekwCSC/(9ρ) is the Stokes impaction parameter when droplet k is collecting an ice particle, and CSC is the Cunningham slip correction factor.

If vjvk, 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 NSti/w with the mixed Froude number NmFr following Hall (1980), Rasmussen and Heymsfield (1985), and Heymsfield and Pflaum (1985). For columnar and planar ice particles, we use formulas EEM17clm and EEM17pln 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),

(40) E rime = ϕ j E BG 74 ( p w / i , N Re j i , N mFr ) + ( 1 - ϕ j ) E EM 17 pln ( N Re j i , N mFr ) .

For ϕj>1 (columnar),

(41) E rime = 1 ϕ j E BG 74 ( p w / i , N Re j i , N mFr ) + 1 - 1 ϕ j E EM 17 clm ( N Re j clm , N mFr ) .

Here, pw/i:=1/pi/w=rk/rji, NmFr=(vj-vk)vk/(gDj/2), and NRejclm=ρvj2aj/μ is the Reynolds number based on the width of column 2aj. 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(aj,cj)<rk, we assume that the resultant ice particle is spherical with the true ice density:

(42)ρji=ρtruei,(43)aj=cj=mj+mk(4π/3)ρtruei13,(44)mjrime=mjrime+mk,(45)njmono=njmono,(46)mαjsol=mαjsol+mαksol,α=1,2,,Nsol,(47)mβjinsol=mβjinsol+mβkinsol,β=1,2,,Ninsol,(48)Tjfz=max(Tjfz,Tkfz),

where primed values indicate the resultant ice particle.

If max(aj,cj)rk, we preserve the ice particle's maximum dimension, i.e., Dj=Dj, 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 Lamb1994b; Morrison and Grabowski2008, 2010; Jensen and Harrington2015; Morrison and Milbrandt2015). As graupels have an aspect ratio of approximately 0.8 (Heymsfield1978), we preserve the minor dimension if 0.8<ϕj1/0.8=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 1.0<ϕj1.25 (columnar but quasi-spherical),

(49)ρji=mj+mkVj+mk/ρrime,(50)aj=aj,(51)cj=Vj+mk/ρrime(4π/3)aj2.

Other attributes are updated using Eqs. (44)–(48). For ϕj>1.25 (columnar) and 0.8<ϕj1.0 (planar but quasi-spherical),

(52)ρji=mj+mkVj+mk/ρrime,(53)aj=Vj+mk/ρrime(4π/3)cj12,(54)cj=cj.

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:

(55) ρ rime = max { min { ρ rime HP 85 ( Y ) , 0.91 g cm - 3 } , 0.1 g cm - 3 } ,

where Y:=(-rkvimp/Tjsfc)/(µmms-1/C), vimp is impact velocity, and Tjsfc is the surface temperature of ice particle j.
https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-g01 where A=0.30gcm-3, B1=0.44, B2=-0.03115, B3=-1.7030, B4=0.9116, and B5=-0.1224.

Impact velocity can be calculated using the formula of Rasmussen and Heymsfield (1985): vimp=|vj-vk|max{fRH85(NReji, NStw/i),0}, where NStw/i=(pw/i)2ρwNReji/(9ρ) is the Stokes impaction parameter when an ice particle collects a droplet. Because the fRH85 given in Rasmussen and Heymsfield (1985) becomes slightly negative around 0.1<NStw/i<1.0, we impose a limiter to ensure it is positive. Surface temperature Tjsfc can be evaluated as

(58) T j sfc = T j + L s D v k Δ ρ j ,

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=min(Y,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

(59) K agg = E agg A j 1 2 + A k 1 2 2 | v j - v k | ,

where Eagg is the collision–aggregation collection efficiency. Following Morrison and Grabowski (2010), we assume that the efficiency is given by a constant, Eagg=0.1, in this study. Field et al. (2006) confirmed that Eagg=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 ρjki=(mj+mk)/(Vj+Vk) is closer to the true density of ice ρtruei, 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 ρjki,min, which can be evaluated using Eq. (61), which we will derive shortly.

In contrast, if ρjki 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 ρcrti=10kgm-3. This choice of value is roughly consistent with observations by Magono and Nakamura (1965). If ρjki is closer to ρcrti, we consider that the apparent density of the resultant aggregate is closer to the maximum possible density ρjki,max. Let us assume ρjki,max=ρjki.

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

Without loss of generality, assume that DjDk. For ϕj≤1 (planar),

(60) a j = a j ,

because we assumed the maximum dimension is preserved. The longest possible minor axis length is cj+min(ak,ck); hence, the largest possible volume becomes Vmax=(4π/3)aj2{cj+min(ak,ck)}. The minimum possible apparent density ρjki,min then becomes

(61) ρ j k i , min = m j + m k V max .

The resultant particle's apparent density is given by a weighted average of ρjki,max=ρjki and ρjki,min:

(62) ρ j i = ( ρ true i - ρ j k i ) ρ j k i , max + ( ρ j k i - ρ crt i ) ρ j k i , min ρ true i - ρ crt i ,

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

(63)cj=mj+mkρji(4π/3)aj2,(64)mjrime=mjrime+mkrime,(65)njmono=njmono+nkmono,(66)mαjsol=mαjsol+mαksol,α=1,2,,Nsol,(67)mβjinsol=mβjinsol+mβkinsol,β=1,2,,Ninsol,(68)Tjfz=max(Tjfz,Tkfz).

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

(69) c j = c j .

If approximating the largest possible particle using an ellipsoid, the largest possible volume becomes Vmax=(4π/3)cj{aj+min(ak,ck)}max(aj,ak,ck). Then, the resultant ice particle's apparent density ρji can be calculated using Eqs. (61) and (62). Then, the minor axis is updated by

(70) a j = m j + m k ρ j i ( 4 π / 3 ) c j 1 2 ,

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

Note that our aggregation outcome model does not produce particles lighter than ρcrti=10kgm-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:

(71)ρt+(ρU)=ρtcm,(72)ρqvt+(ρqvU)=ρqvtcm+Dv2(ρqv),(73)ρUt+(ρUU)=-P-ρgz^+ρUtcm+μ2U,(74)ρθt+(ρθU)=ρθtcm+kcp2θ,(75)P=ρRT=P0ρθRP0cp/cp-R,

where the four terms with the form /t|cm represent cloud microphysics coupling terms. ρ/t|cm=ρqv/t|cm is the source of vapor:

(76) ρ t cm = ρ q v t cm = s v + s s .

Here, sv and ss are sources of vapor through condensation/evaporation and deposition/sublimation, respectively:

(77)sv(x,t)=-iIr(t)δ3(x-xi(t))dmidtcnd/evp,(78)ss(x,t)=-iIr(t)δ3(x-xi(t))dmidtdep/sbl,

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

ρθ/t|cm represents heating due to the phase transition of water:

(79) ρ θ t cm = - L v s v + L s s s + L f s f c p Π ,

where Lf is the latent heat of fusion, and sf is the production rate of liquid water through freezing, melting, or riming. Let tnfz be the time of the nth freezing event and infz be the index of the frozen droplet. Similarly, let tnmlt and inmlt be the time and melted ice particle of the nth melting event, respectively. Let tnrime and inrime be the time and rimed droplet of the nth riming event, respectively. Then, sf is given by

(80) s f ( x , t ) = - freezing event n δ 3 ( x - x i n fz ( t ) ) δ ( t - t n fz ) m i n fz ( t ) + melting event n δ 3 ( x - x i n mlt ( t ) ) δ ( t - t n mlt ) m i n mlt ( t ) - riming event n δ 3 ( x - x i n rime ( t ) ) δ ( t - t n rime ) m i n rime ( t ) .

ρU/t|cm is the drag force from the particles. From Eq. (2), we can derive Fidrg=migz^+d(mivi)/dt. The terminal velocity assumption does not mean that the second term vanishes because mi and vi are still time dependent. However, even if a droplet accelerated from 0 to 10 m s−1 in 100 s through rapid precipitation development, the contribution of the second term is much smaller than that of the first term: 10ms-1/100sg. Thus, we finally obtain

(81) ρ U t cm = - i I r ( t ) δ 3 ( x - x i ( t ) ) F i drg - i I r ( t ) δ 3 ( x - x i ( t ) ) m i ( t ) g z ^ .

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. 7181) 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 ρqv, 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: ρ, ρqv, 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 Glmn 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 {{xi(t),ai(t)},i=1,2,,Nrwp} by a population of super-particles: {{ξi(t),xi(t),ai(t)},i=1,2,,Nswp} (see, e.g., Fig. 4 of Grabowski et al.2019). A super-particle is characterized by multiplicity ξi, position xi, and attributes ai. We consider that the ith super-particle represents ξi real particles {xi,ai}. Note that multiplicity ξi is an integer and is time dependent. Nswp 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(a,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:

(82) n ( a , x , t ) = i I r ( t ) δ d ( a - a i ( t ) ) δ 3 ( x - x i ( t ) ) ,

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

(83) n ( a , x , t ) = i I s ( t ) ξ i ( t ) δ d ( a - a i ( t ) ) δ 3 ( x - x i ( t ) ) = N s ( t ) ξ = 1 ξ p ( ξ , a , x , t ) ,

where p(ξ,a,x,t) is the probability density that a super-particle has multiplicity ξ, attributes a, and position x at time t; Is(t) is the set of super-particle indices existing in the domain at time t; and Ns(t):=#Is(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(ξ,a,x,t=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., ξ=ξ(a,x). Then, Eq. (83) at t=0 reduces to

(84) n ( a , x , 0 ) = N s ( 0 ) ξ ( a , x ) p ( a , x ) .

If we set ξ(a,x) as a constant, the probability density function must be proportional to the initial distribution function of real particles: p(a,x)n(a,x,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:

(85) ξ ( a , x ) = n ( a , x , 0 ) N s ( 0 ) p , p ( a , x ) = p = const .

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 Shima2013; 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., Niederreiter1978). 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 {{ξi,xi,ai}} and Glmn are updated from time t to tt.

Let Δtadv, Δtfz∕mlt, Δtcnd∕evp, Δtdep∕sbl, and Δtcollis 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 Δtdyn 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. 7681), and update moist air from Glmn(t) to Glmn(t). Then, we update super-particles {{ξi,xi,ai}} from t to tt. 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. 7681), and update the moist air from Glmn(t) to Glmn(tt). Table 1 shows an example of the calculation order.

Table 1An example of the calculation order when updating the system state from t to tt. 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.

Download Print Version

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 Gi:=G(xi) 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 Δtadv. We normally select a short enough Δtadv 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 Ulmn on the face grid. We then interpolate Ulmn 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 Δtfz∕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 Δtcnd∕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 r2 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 Δtcnd∕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 Δtdep∕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 m2∕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 dm2/3/dt=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 m2∕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.

Δtdep∕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.84.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(Ns), 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.

Δtcollis 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 Δtcollis 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 P be the typical probability that a particle coalescence/rime/aggregate with another particle within a small time interval Δtcollis. From Eq. (31), P can be evaluated as

(86) P N r K Δ t collis Δ V n r K Δ t collis ,

where Nr is the number of real particles in a volume ΔV, K is the typical value of the coalescence/riming/aggregation kernel K, and nr is the number concentration of real particles. Requiring that P<1 has to be satisfied, we obtain

(87) Δ t collis < 1 / ( n r K ) .

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

(88) P s N s ( N s - 1 ) 2 / N s 2 ξ K Δ t collis Δ V N s N r N s K Δ t collis Δ V P ,

where Ns 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 ξNr/Ns is the typical multiplicity.

In SDM, the multiple coalescence technique is used to make the algorithm robust to larger Δtcollis. Here, we clarify how we adapt it to riming and aggregation. If it is a coalescence between a droplet j and γ̃ number of droplets k (see Sect. 5.1.3 of Shima et al.2009, for the definition of γ̃), we modify Eqs. (33)–(35) by applying

(89) ( r k 3 , m α k sol , m β k insol ) γ ̃ ( r k 3 , m α k sol , m β k insol ) .

If it is a coalescence between γ̃ number of droplets j and a droplet k, we apply

(90) ( r j 3 , m α j sol , m β j insol ) γ ̃ ( r j 3 , m α j sol , m β j insol ) .

Similarly, if it is a riming/aggregation between a particle j and γ̃ number of particles k, we apply the following replacement to Eqs. (42)–(54) and (60)–(70):

(91) ( m k , V k , m k rime , n k mono , m α k sol , m β k insol ) γ ̃ ( m k , V k , m k rime , n k mono , m α k sol , m β k insol ) .

If it is a riming/aggregation between γ̃ number of particles j and a particle k,

(92) ( m j , V j , m j rime , n j mono , m α j sol , m β j insol ) γ ̃ ( m j , V j , m j rime , n j mono , m α j sol , m β j insol ) .

What is not straightforward is the calculation of Vmax used in the aggregation outcome formula. For planar collector j, we consider that Vmax is given by

(93) V max = ( 4 π / 3 ) a j 2 { c j + γ ̃ min ( a k , c k ) } , ( 4 π / 3 ) a j 2 { γ ̃ c j + min ( a k , c k ) } .

For columnar collector j,

(94) V max = ( 4 π / 3 ) c j { a j + γ ̃ min ( a k , c k ) } max ( a j , a k , c k ) , ( 4 π / 3 ) c j { γ ̃ a j + min ( a k , c k ) } max ( a j , a k , c k ) .

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. 7181). In this study, as explained in the previous section, the four coupling terms from cloud microphysics denoted by /t|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 γ=10-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 Δtdyn 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 2Summary of numerical experiments for model evaluation. The domain is two-dimensional (xz): 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.

Download Print Version | Download XLSX

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 (xz), 60 km in the horizontal direction and 16 km in the vertical direction.

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

6.1.2 Moist air boundary conditions

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

6.1.3 Initial conditions of particles

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

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

(95) d N sulf d log r dry sulf = a = 1 2 1 2 π c a sulf log σ a exp - log r dry sulf - log r a 2 2 log 2 σ a ,

where rdrysulf is the dry radius of the ammonium bisulfate component and Nsulf is the accumulated number of particles smaller than rdrysulf per unit volume of air. The particle number concentrations are c1sulf=270cm-3 and c2sulf=45cm-3; thus, the total particle number concentration is csulf=c1sulf+c2sulf=315cm-3. The geometric mean radii are r1=0.03 µm and r2=0.14 µm, with geometric standard deviations of σ1=1.28 and σ2=1.75, respectively. This distribution is based on in situ maritime aerosol data as detailed in Sect. 2.2.3 of VanZanten et al. (2011), but the number concentration is multiplied by 3. As discussed in Sect. 2.4, we consider that a droplet containing only soluble substances freezes only through a homogeneous freezing mechanism; therefore, the freezing temperature of these particles is Tfz=-38C. Therefore, pure ammonium bisulfate's initial distribution function can be calculated as

(96) n sulf ( log r dry sulf , T fz ) = d N sulf d log r dry sulf δ ( T fz - ( - 38 C ) ) .

The other aerosol population consists of mineral dust internally mixed with ammonium bisulfate. We set the number concentration to cdust=1cm-3, and for simplicity, set the mineral dust particle diameter to ddust=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(Tfz) 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:

(97) n S ( T ) = 0 , for T > T max fz ; n S Niemand ( T ) , for T max fz T > T min fz ; n S Niemand ( T min fz ) , for T min fz T ;

where Tmaxfz=-12C and Tminfz=-36C. The mineral dust surface area is given by Ainsol=π(ddust)2. As discussed in Sect. 2.4, we set Tfz=-38C if the mineral dust is IN inactive and no INAS appears until Tfz=-38C. Altogether, the mineral dust distribution function is given by

(98) n dust ( d dust , log r dry sulf , T fz ) = δ ( d dust - 1 µ m ) c dust c sulf d N sulf d log r dry sulf [ p ( T fz ) H ( T fz + 38 C ) + P INia δ ( T fz + 38 C ) ] ,

where H(T) is the Heaviside step function and PINia:=P(Tfz-38C) is the probability that a single INAS does not appear until Tfz=-38C. For ddust=1 µm, PINia≈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

(99)ρθtsfc=ρHmax(W,0),(100)W=-4w2(x-x0)2-w22exp-z-z0z0,

where H=10Kh-1, x0=5 km, w=10 km, and z0=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 Δx=Δy=Δz=62.5m in the CTRL case. The time steps in the CTRL case are Δt=0.4 s, Δtadv=Δtfz/mlt=0.4s, Δtcollis=0.2 s, Δtcnd/evp=Δtdep/sbl=0.1s, and Δtdyn=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 cSP=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 [rdry,minsulf,rdry,maxsulf] and determine the dry radius rdry,isulf. To accurately represent the size distribution given in Eq. (95), we set rdry,minsulf=10.0nm and rdry,maxsulf=5.0µm. From Eqs. (85) and (96), the super-particle's multiplicity is then given by

(101) ξ i = n sulf ( log r dry , i sulf , T i fz ) N s ( 0 ) / 2 V domain log ( r dry , max sulf / r dry , min sulf ) δ ( T fz - ( - 38 C ) ) = d N sulf d log r dry sulf log r dry , i sulf log ( r dry , max sulf / r dry , min sulf ) c SP / 2 ,

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

(102)n=nsulf(logrdry,isulf,Tifz),(103)p=δ(Tfz-(-38C))Vdomainlog(rdry,maxsulf/rdry,minsulf),

and Ns(0) in Eq. (85) is replaced by Ns(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 rdry,isulf as m1isol=(4π/3)ρ(NH)4HSO4(rdry,isulf)3, where ρ(NH)4HSO4=1.78gcm-3. The soluble aerosol particle freezing temperature is Tifz=-38C.

For IN inactive mineral dust super-particles, we use PINiaSP=0.05. The mineral dust initially has the same size ddust=1 µm. The dry radius rdry,isulf 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 [rdry,minsulf,rdry,maxsulf]. The IN inactive mineral dust freezing temperature is Tifz=-38C. From Eqs. (85) and (98), an IN inactive mineral dust super-particle's multiplicity is then given by

(104) ξ i = c dust c sulf d N sulf d log r dry sulf log r dry , i sulf log ( r dry , max sulf / r dry , min sulf ) c SP / 2 P INia P INia SP .

Finally, we consider IN active mineral dust internally mixed with ammonium bisulfate. The remaining super-particles, i.e., (1-PINiaSP)/2, are used for this population. The initial diameter of the mineral dust is ddust=1 µm, and the dry radius rdry,isulf is determined as in the other populations. We draw another random number uniformly from the interval [Tminfz,Tmaxfz] and determine the freezing temperature Tifz. From Eqs. (85) and (98), an IN active mineral dust super-particle's multiplicity is then given by

(105) ξ i = c dust c sulf d N sulf d log r dry sulf log r dry , i sulf p ( T i fz ) log ( r dry , max sulf / r dry , min sulf ) ( T max fz - T min fz ) ( c SP / 2 ) ( 1 - P INia SP ) .

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 ri 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 cSP, we vary cSP 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: Δx=Δy=Δz=250m, Δtdyn=0.2 s, and Δt=Δtadv=Δtfz/mlt=1.6s, but other time steps, i.e., {Δtcollis,Δtcnd/evp,Δtdep/sbl}, are not changed.

The DXx2 ensemble's grid size is twice that of CTRL: Δx=Δy=Δz=125m, Δtdyn=0.1 s, and Δt=Δtadv=Δtfz/mlt=0.8s.

Note that DXx1 and CTRL are the same.

The DX/2 ensemble has a grid size that is half that of CTRL: Δx=Δy=Δz=31.25m. Δtdyn=0.025 s, and Δt=Δtadv=Δtfz/mlt=0.2s.

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: Δt=Δtadv=Δtfz/mlt=4.0s, Δtcollis=2.0 s, and Δtcnd/evp=Δtdep/sbl=1.0s. Δtdyn is not changed.

The time steps of the DTx5 ensemble are 5 times that of CTRL: Δt=Δtadv=Δtfz/mlt=2.0s, Δtcollis=1.0 s, and Δtcnd/evp=Δtdep/sbl=0.5s.

The time steps of the DTx2 ensemble are twice that of CTRL: Δt=Δtadv=Δtfz/mlt=0.8s, Δtcollis=0.4 s, and Δtcnd/evp=Δtdep/sbl=0.2s.

Note that DTx1 and CTRL are the same.

The time steps of the DT/2 ensemble are half that of CTRL: Δt=Δtadv=Δtfz/mlt=0.2s, Δtcollis=0.1 s, and Δtcnd/evp=Δtdep/sbl=0.05s.

The time steps of the DT/4 ensemble are one-quarter that of CTRL: Δt=Δtadv=Δtfz/mlt=0.1s, Δtcollis=0.05 s, and Δtcnd/evp=Δtdep/sbl=0.025s.

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 mrime/m>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., nmono>10, we consider the ice particle a snow aggregate. Otherwise, we categorize the ice particle as a cloud ice particle.

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f01

Figure 1Typical 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.

Download

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

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

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

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

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

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f02

Figure 2Time 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.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f03

Figure 3Time 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.

Download

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 ρtruei:

(106) m * := ρ i a 2 c ρ true i ( D / 2 ) 3 , D = 2 max ( a , c ) .

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

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

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

Figure 7 shows the relationship between ice particle velocities and dimensions at t=2040, 3000, and 5400 s. The horizontal axis represents the maximum ice particle dimension D, and the vertical axis represents ice particle terminal velocity v. The 2-D mass densities of cloud ice particles, graupel particles, and snow aggregates are plotted in the same manner as Fig. 4, except for differences in the vertical axis. The colored slopes in Fig. 7 represent the velocity–dimension relationship formulas from various studies, and SC85, W08, H72, KH83, A72, and H02 indicate Starr and Cox (1985), Westbrook et al. (2008), Heymsfield (1972), Knight and Heymsfield (1983), Auer (1972), and Heymsfield et al. (2002), respectively. “Stokes' law for ice spheres” is based on the Stokes' terminal velocity for spherical ice particles with the true ice density. We use the dynamic viscosity at a temperature of −20C, i.e., μ=1.630×10-5kgm-1s-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/2a=0.05 is assumed, with L being the height of the hexagonal prism and a=D/2 being the hexagon's maximal radius. The effective radius is calculated using the horizontal orientation model from Roscoe (1949). For “hexagonal columns”, L/2a=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 −20C. See also Movie 5 in the video supplement.

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f04

Figure 4Mass–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.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f05

Figure 5Aspect 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.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f06

Figure 6Apparent 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 ρi. 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.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f07

Figure 7Velocity–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.

Download

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 ϕ=1/0.8 in Fig. 5 were created because of our assumption that riming growth eventually makes ice particles quasi-spherical.

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

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

The mass–dimension relationship shown in Fig. 4, and the velocity–dimension relationship shown in Fig. 7 show a reasonable agreement between our model's predicted results and existing formulas based on laboratory measurements and observations. In both figures, cloud ice particles, graupel particles, and snow aggregates are distributed near the blue, red, and green slopes, respectively. In Fig. 7, one might note that snow aggregates in our model fall faster than those in the formulas; however, this bias can be explained by the air density dependence of the fall speed. The green slopes in Fig. 7 represent the formulas of LH74 and H02. LH74's formulas are constructed from data measured between altitudes of 750 and 1500 m above sea level; hence, the density is approximately 1.1 kg m−3. H02's formula is for temperature and pressure of −10C and 500 hPa; hence, the density is approximately 0.66 kg m−3. In our simulation, most of the snow aggregates exist in the anvil cloud, wherein the density is approximately 0.38 kg m−3. Khvorostyanov and Curry (2002) estimated that the terminal velocities of large ice particles scale with the ambient density to the power of -1/2. Therefore, we can incorporate the density dependence by multiplying the LH74 formulas for aggregates by a factor of (0.38kgm-3/1.1kgm-3)1/21.70 and that of H02 for aggregates by a factor of (0.38kgm-3/0.66kgm-3)1/21.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, 47: 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. 281) 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 cSP was investigated by varying the cSP 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 cSP. 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 cSP. 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 cSP. This indicates that fluctuations in all simulations are mostly dominated by atmospheric turbulence. One might note that the fluctuations are slightly increasing as cSP 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 cSP 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 cSP 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., cSP=128 per cell.

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f08

Figure 8Statistics 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 cSP. 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.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f09

Figure 9Statistics 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 cSP. 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.

Download

8.2 DX ensembles and grid convergence

We investigated the numerical convergence with respect to the grid size by varying Δx=Δy=Δ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., Δx=Δy=Δz=62.5m.

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f10

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

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f11

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

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f12

Figure 12Spatial 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.

Download

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, Δtadv, Δtfz∕mlt, Δtcollis, Δtcnd∕evp, Δtdep/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.

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f13

Figure 13Statistics 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.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f14

Figure 14Statistics 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.

Download

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 cSP=128 per cell, which is comparable to various previous studies (e.g., Andrejczuk et al.2010; Sölch and Kärcher2010; Riechelmann et al.2012; Arabas and Shima2013; Unterstrasser and Sölch2014; Unterstrasser et al.2017; Grabowski et al.2018; Jaruga and Pawlowska2018; Dziekan et al.2019; Hoffmann et al.2019). Those studies reported that approximately 50200 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 Δx=Δy=Δz=62.5m, 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 tadv, Δtfz∕mlt, Δtcollis, Δtcnd∕evp, Δtdep/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, Δtadv should be limited by the CFL condition of wind velocity.

To avoid a sudden release of latent heat, Δtfz∕mlt must also be restricted through the CFL condition.

Δtcollis is the time step of coalescence, riming, and aggregation. As shown in Shima et al. (2009) and clarified in Sect. 5.5.5, Δtcollis can be estimated from the number concentration and size of real particles, and that Δtcollis 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 (Okawa2015), which could allow us to use larger Δtcollis values.

Δtcnd∕evp and Δtdep∕sbl are determined by the phase relaxation time of supersaturation, τphase1/ξiri (e.g., Squires1952). 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, Δtcnd∕evp and Δtdep∕sbl must be smaller than the phase relaxation time. Otherwise, numerical instability occurs (Árnason and Brown1971). 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 10100 bins are needed for each axis, the total number of bins becomes 1051005. 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., 101010010. 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., 1051005. 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 47, 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 47. 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 {a,c,ρi,mrime/m,nmono}={2.58mm,36.1mm,1.12×102kgm-3,0.98,213}. 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=4.25×10-1ms-1 is also much smaller than that of typical hailstones (Fig. 7). It is located at (x,z)=(20.0km,4.11km), 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 0C. Then, from Eqs. (55)–(57), ρrime=0.1gcm-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)=(1mm,11mm).

However, ρrime=0.1gcm-3 must be reconsidered. Equation (56) has a global maximum of approximately 0.95 g cm−3 at around Y=3.7 and then quickly decreases, becoming <0.1gcm-3 at around Y=5.5. Then, from Eq. (55), ρrime=0.1gcm-3 for Y>5.5. Consequently, considering the definition Y:=(-rkvimp/Tjsfc)/(µmms-1/C), ρrime=0.1gcm-3 frequently happens near the freezing level. For example, Y=1000 for rk=1 mm, vimp=1ms-1, and Tjsfc=-1C. However, ρrime would be much larger and even closer to ρtruei 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

(107) Y = min ( Y , 3.5 ) .

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

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 {a,c,ρi,mrime/m,nmono} ={24.9µm, 138.8 µm, 269.7kgm-3, 0.29, 1}. Its terminal velocity is v=1.46×10-2ms-1 and it is located at (x,z)=(20.3km,11.0km) (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 10-2ms-1 and ambient temperature of −10C, then, from Eqs. (55)–(57), Y=10-2 and ρrime=0.1gcm-3. In other words, the apparent volume of the rimed droplet expands 10-fold and creates a columnar graupel particle: (a,c)=(10µm,110µ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

(108) a j = max ( a j , r k ( ρ w / ρ rime ) 1 / 3 ) ,

and Eq. (54) with

(109) c j = max ( c j , r k ( ρ w / ρ rime ) 1 / 3 ) .

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

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 {a,c,ρi,mrime/m,nmono} ={12.6mm, 15.0 mm, 10.7kgm-3, 0.85, 1 585 116}. Its terminal velocity is v=4.31ms-1 and it is located at (x,z)=(10.6km,11.5km) (see Fig. 1). What is unusual here is the very low apparent density ρi=10.7kgm-3. This particle is composed of many monomers nmono=1 585 116, and we set the limiting value of aggregate density in Eq. (62) to ρcrti=10kgm-3. Thus, we can conclude that this hailstone is created by the repeated aggregation between graupel particles.

Lump graupel particles with apparent densities as low as 50 kg m−3 were reported in Locatelli and Hobbs (1974). Therefore, a hailstone with an apparent density of 10 kg m−3 is not extremely unrealistic. However, our aggregation model is crude. Following Morrison and Grabowski (2010), we assumed that collision–aggregation collection efficiency is a fixed constant of Eagg=0.1 regardless of morphology or temperature. Therefore, this could cause the accumulation of graupel particles near the limiting value ρcrti 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 {a,c,ρi,mrime/m,nmono} ={31.9µm, 438.2 µm, 19.4kgm-3, 0.53, 26 679}. Its terminal velocity is v=1.02×10-3ms-1, it is located at (x,z)=(29.7km,6.4km) (see Fig. 1), and the ambient temperature is T=-14.4C.

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 [−20C,−10C] and [−5C,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,

(110) Γ ( T i ) = 1 , for d m i < 0 (sublimation) ;

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

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

9.1.5 Results after corrections

In the preceding sections, we proposed three corrections to the time evolution equations (Eqs. 107110) 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. 1518) to the original results (Figs. 47), 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).

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f15

Figure 15This 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. 107110) proposed to avoid the creation of ice particles with unrealistic morphologies. See also Movie 8 in the video supplement.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f16

Figure 16This 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. 107110) proposed to avoid the creation of ice particles with unrealistic morphologies. See also Movie 9 in the video supplement.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f17

Figure 17This 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. 107110) proposed to avoid the creation of ice particles with unrealistic morphologies. See also Movie 10 in the video supplement.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f18

Figure 18This 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. 107110) proposed to avoid the creation of ice particles with unrealistic morphologies. See also Movie 11 in the video supplement.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f19

Figure 19Changes 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.

Download

https://gmd.copernicus.org/articles/13/4107/2020/gmd-13-4107-2020-f20

Figure 20Changes 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.

Download

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 qi≤1 always holds in our model, Böhm's formula vBöhm(mi,ϕi,di,qi;ρi,Ti) can be summarized as follows:

(111)X=8migρπμ2max(ϕi,1)qi1/4,(112)X=X1+(X/X0)21+1.6(X/X0)2,(113)X0=2.8×106,for ice particles,k=min{max(0.82+0.18ϕi,0.85),0.37+0.63ϕi,(114)1.33max(logϕi,0)+1.19},(115)Γ=max{1,min(1.98,3.76-8.41ϕi+9.18ϕi2-3.53ϕi3)},(116)CDP=max(0.292kΓ,0.492-0.200/ϕi),(117)CDO=4.5k2max(ϕi,1),(118)β=1+CDP6kXCDP1/21/2-1,(119)γ=CDO-CDP4CDP,(120)NRe=6kCDPβ21+2βe-βγ(2+β)(1+β),(121)vBöhm=μNReρdi.

In our model (SCALE-SDM 0.2.5-2.2.0/2.2.1), we assumed that the characteristic length di is given by the maximum dimension Di=2max(ai,ci), and the area ratio qi is given by the area ratio regarding the circumcircle qicc=Ai/Aicc (Eq. 4). However, in Böhm's theory, they are defined by

(122) d i = 2 a i , q i = q i ce = A i / A i ce ;

i.e., for columnar particles, minor axis is used for di, and the area ratio regarding circumscribed ellipse is used for qi. Figure 1 in Böhm (1989) suggests qi=qice. It is not clearly specified, but from the second equality of Eq. 17 in Böhm (1992c), we can confirm that di=2ai.

For planar ice particles (ϕi<1), vBöhm(di=2ai,qi=qice) and vBöhm(di=Di,qi=qicc) yield the same results, because 2ai=Di and qice=qicc hold. However, for columnar ice particles (ϕi>1), vBöhm(Di,qicc) always underestimates the fall velocity. From Eqs. (111)–(121), we can derive vBöhm(2ai,qice)/vBöhm(Di,qicc)=ϕi3/4 for X≪1, and vBöhm(2ai,qice)/vBöhm(Di,qicc)=ϕi7/8 for X≫1. Therefore, if ϕi=2, the ratio vBöhm(2ai,qice)/vBöhm(Di,qicc) is in the rage of 1.681.83. If ϕi=10, the range is 5.627.50, and if ϕi=20, it is 9.4613.75. We also confirmed that Böhm's original definition vBöhm(2ai,qice) 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:

(123) Γ * = 1 for d m i 0 ϕ i > 40 .

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 (Levine1950), we considered that each insoluble particle has its own freezing temperature Tfz 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 Tfz 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 Tfz, 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 Tfz=-38C 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 Tfz 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 <-25C 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 nS is a function of xtherm, and xtherm is a function of temperature T and saturation ratio over ice Si. We can assign xtherm,i to each particle as an attribute. We consider that freezing occurs when xtherm(T,Si)>xtherm,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” (Marcolli2017). 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 Tfz by 4.5C 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 Shaw2005); 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 +4C. 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 300400 µm ranges between 20 and 70 s (Rasmussen and Pruppacher1982). 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, Dv=2.52×10-5m2s-1 and k=2.55×10-2Jm-1s-1K-1, which are the values for T=20C and p=1000 hPa. However, this approximation is erroneous, particularly because diffusivity Dv is inversely proportional to pressure. In the case of the initial profile we used for model evaluation, T=-44C and p=250 hPa at z=10 km. Thus, Dv=6.08×10-5m2s-1, which is about 2.4 times larger than we assumed. The temperature and pressure dependence of water vapor's diffusivity in air Dv, 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 Dv 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 (Low1969); 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 Lamb1994a). 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<-20C, both plates and columns can be created at the same temperature depending on the saturation ratio over ice Si, and polycrystals can also be created. Therefore, for T<-20C, Γ might better be considered a function of both T and Si, 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: