A new method for large-eddy simulations of clouds with Lagrangian droplets including the effects of turbulent collision

In this paper, a new Lagrangian cloud model (LCM) is introduced in which the flow field is simulated by large-eddy simulation, and the droplets are treated as Lagrangian particles responding to the simulated flow field. In order to handle the extremely large number of droplets within a cloud, the concept of a super-droplet, which represents a large number of real droplets of the same size, is introduced, and the number of contributing real droplets is called the weighting factor. A novel method is developed to realize the collision/coalescence of droplets, in which the consequent variation of the droplet spectrum is represented in terms of the modification of the radius and weighting factor of super-droplets, while keeping the number of super-droplets unchanged. Using an idealized single cloud and trade wind cumuli, the LCM is shown to reproduce the general features of shallow cumulus clouds in agreement with traditional bulk models. The droplet spectrum simulated by the LCM, using collision kernels with and without the effects of turbulence, also shows a pattern consistent with the spectral bin model. Furthermore, the sensitivity of the LCM to two model parameters, the time step and the number of super-droplets, is examined.


Introduction
To simulate particle-laden turbulent flows, the Lagrangian model, which calculates the trajectories of individual particles under a turbulent velocity field simulated by large-eddy simulation (LES) or direct numerical simulation (DNS), is more natural than the Eulerian model, which calculates the variation of the particle concentration field (see, e.g., [1,2]). It is also important to note that the dynamics and microphysics of a cloud are directly related to the physical process of individual droplets, such as condensation, collision, transport and settling in the Lagrangian cloud model (LCM).
An LCM can realize the initiation of rain naturally by the motion of droplets settling gravitationally in turbulent flows, since the settling velocity is calculated directly. The transport of droplets, responding to the velocity field in a cloud, can reflect the different response of individual droplets depending on different inertial and gravitational forces, without resorting to the assumption of eddy diffusivity. Meanwhile, as pointed out by Andrejczuk et al [3], a Lagrangian model can avoid difficulties associated with the conversion from cloud condensation nuclei (CCN) into cloud droplets and the artificial broadening of the droplet spectrum in the spectral bin model. It should also be mentioned that particle-turbulence interaction, such as preferential concentration, droplet acceleration associated with the intermittency of turbulence and interaction of droplets with supersaturation fluctuation, which may exert an influence on cloud microphysics, can be understood only from Lagrangian perspectives (e.g. [4][5][6]).
Nonetheless, most of the cloud models developed so far followed the Eulerian model, using either the bulk model or the spectral bin method for microphysics. The bulk method calculates a 3 small number of bulk variables, such as vapor mass, cloud water mass and rain mass [7][8][9][10]. It does not require much computational cost, but depends heavily on the assumptions of the basic microphysical processes, such as the assumed size distribution of droplets. The spectral bin method solves the integral-differential equation that describes the evolution of droplet spectra, called the stochastic collection equation (SCE), assuming that the distribution of cloud droplets within a grid box is uniform [11][12][13].
In the spectral bin model the evolution of the droplet spectrum is mainly controlled by the collection kernel K (R, r ), which represents the probability that a droplet of radius R collects a droplet of radius r per unit time, given that both are present with unit concentration. To solve the SCE, the droplet spectrum must be discretized into bins. If the number of bins is not sufficient, the spectral bin method suffers from the numerical diffusion of the droplet spectrum, which leads to an earlier onset of rainfall [3,14]. On the other hand, the computational cost for calculating the SCE at every grid point with a sufficiently large number of bins becomes insurmountably high.
The progress in LCMs could not be realized until recently because of two challenging problems. One is how to handle the extremely large number of droplets within a cloud, which makes it practically impossible to realize the motion of all droplets. However, the total number of droplets is important in the simulation of a cloud, because it determines the total liquid water content (LWC). The other problem is how to realize the collisions of Lagrangian droplets, followed by coalescence, whose process is highly complicated and sensitive to the fine-scale turbulence structure.
Only two groups are known to the authors so far who attempted to develop an LCM recently: Andrejczuk et al [3,15] and Shima et al [16]. In both the cases, each simulated droplet is assumed to represent a large number of droplets; it is referred to as a 'super-droplet' by Shima et al [16]. As to the collision process, Shima et al [16] used a Monte Carlo scheme for the collision of simulated droplets, in which collisions between randomly sampled pairs of droplets are calculated, while assuming that the resulting droplet statistics are equivalent to those from the collision between all super-droplets. On the other hand, Andrejczuk et al [3] used a scheme in which collisions are assumed to occur between super-droplets within a subdivided grid box, and new super-droplets are created as the outcome of the collisions.
In this work, we have developed a new LCM. Similarly to Andrejczuk et al [3,15] and Shima et al [16], in our model the flow field is simulated by LES and the concept of 'superdroplets' is applied. Thereby a single super-droplet represents a large number of real droplets of the same size, and the number of contributing droplets to a super-droplet is called a 'weighting factor'. However, in contrast to the existing LCMs, the droplet collision in our LCM is realized in a radically different way. Instead of considering individual collision events between superdroplets explicitly, we choose a statistical approach in which the growth of the droplet size after a large number of collision events is determined by the background droplet spectrum and turbulence, and the outcome is described in terms of the modification of the radius and weighting factor. In this approach the creation/annihilation of droplets after collision/coalescence is not necessary. A detailed description is presented in section 2.
The collection kernel plays a critical role in the calculation of droplet growth by collision/coalescence in the present scheme, as in the case of the spectral bin model. The gravitational kernel with collision efficiencies from Hall [11] (hereafter called the Hall kernel) has been widely used, which considers only gravitational collision independent of turbulence. However, recent studies of the collection kernel through either DNS [17][18][19] or a kinematic/stochastic representation of turbulence [20] reveal that turbulence enhances the collision frequency significantly. The main mechanisms of enhancing the collision frequency are the enhanced relative velocity between colliding droplets, the preferential concentration and the enhanced collision efficiency. Ayala et al [17] developed a new collision kernel including the effects of turbulence based on the analysis of DNS data (hereafter called AW kernel), and Wang and Grabowski [21] showed using a spectral bin model that the AW kernel can explain the observed droplet broadening between 10 and 50 µm, which could not be predicted from the Hall kernel.
The purpose of this paper is to propose a new LCM, which has the potential to provide new insight for the understanding of cloud physics and to validate it. Section 2 discusses the fluid and particle models and their numerical implementation. Following that, the developed LCM is applied to the simulation of shallow cumulus convection of warm clouds: for an idealized single cloud and trade wind cumuli corresponding to the Barbados Oceanographic and Meteorological Experiment (BOMEX; [22]) shallow cumulus case. The case of a single cloud, idealized as a rising air bubble, is presented in section 3. Here the characteristics of cloud dynamics and microphysics are analyzed and compared with the traditional bulk model results. The simulation of the LCM is carried out with different collection kernels, the Hall and the AW kernel, in order to clarify the effects of turbulence on droplet growth. The idealized single cloud is also used for sensitivity tests of model parameters such as the time step and the number of super-droplets. In section 4, the validity of the model is further examined for the BOMEX data in comparison with other model results. Concluding remarks are given in section 5.

Fluid model
The LES model used in this work is PALM (PArallelized Large-Eddy Simulation Model; [23]), which is designed for optimal performance and high scalability on massively parallel computer architectures [23]. PALM has been successfully applied to various geophysical applications including the convective boundary layer (see, e.g., [24][25][26]). It has also been applied to simulate the motion of Lagrangian particles in turbulent flows, such as particle settling, dispersion and footprint problems [27,28].
PALM is based on the non-hydrostatic incompressible Boussinesq equation and the equations for the conservation of mass, energy and moisture, which are filtered over the grid size, as Here u i are the velocity components (u, v, w); x i are the Cartesian coordinates (x, y, z), f i are the Coriolis parameters, u g is the geostrophic wind, ρ 0 is the density of dry air, p * is the socalled perturbation pressure, g is the gravitational acceleration and θ v = θ (1 + 0.608q − q l ) is the virtual potential temperature with potential temperature θ , specific humidity q and liquid water mixing ratio q l . Furthermore, L e is the latent heat of evaporation, c p is the specific heat of dry air at constant pressure, = ( p/ p 0 ) R d /c p is the Exner function with the reference pressure p 0 = 1000 hPa and the gas constant R d , and is the source/sink term of q l through the condensation/evaporation process. In the present model q l and are calculated by the particle model described in section 2.2. All the variables in the equations are filtered ones, but the overbar indicating filtered quantities is omitted except for the subgrid-scale (SGS) flux terms for readability. The variables with the double prime represent the SGS components. Molecular diffusion and radiation processes are generally neglected, but radiative cooling is included in one of our simulation setups (see section 4.1). The filtering is realized by averaging over discrete grid volumes after Schumann [29]. The SGS turbulence is parameterized according to Deardorff [30] that includes a prognostic equation for the SGS turbulent kinetic energy (SGS-TKE) e = 1 2 u i u i : Here ε is the dissipation rate. The SGS fluxes in the filtered equations (1) and (3)-(5) are parameterized using the SGS eddy viscosity and diffusivity K m and K h as Here K m and K h are calculated by where c m = 0.1 and s = 3 √ x y z. The SGS mixing length l is calculated by where z is the distance from the bottom. The dissipation rate ε is given by The model equations are discretized by finite differences. In the present study, a fifthorder advection scheme of Wicker and Skamarock [31] and a third-order Runge-Kutta time step scheme [32] are used. To ensure a divergence-free wind field, a predictor-corrector method is used where the Poisson equation for the perturbation pressure is solved after every time step. For a more detailed description of PALM, see [23,33].
In addition to LCM, PALM also has the option for a traditional one moment bulk cloud model (see, e.g., [23,34]). In this model prognostic equations for the liquid water potential temperature and total water content are solved following Ogura [35] and Orville [36]. The condensation scheme follows Cuijpers and Duynkerke [37] and is an all-or-nothing method, where a grid volume is regarded as unsaturated when the total water content is below the saturation value, and saturated otherwise. In the saturated case q l is calculated by the difference between the total water content and the saturation value.

Particle model
Within the model each particle is defined by a set of features (position, velocity components, radius, weighting factor, etc), which are stored as components of a FORTRAN90-derived data type variable. One element of this data type defines a complete particle with its entire features and the whole set of particles is stored in a one-dimensional (1D) array of this derived data type. The particle model of PALM is parallelized using the domain-decomposition concept, as in the fluid model. Therefore the particles are handled by and stored in the memory of the processors assigned to those subdomains where the particles are located. The particle model includes a sorting algorithm which speeds up the performance by avoiding cache misses. For this purpose the particles are sorted at each time step. The sorting determines all particles belonging to the same grid box and also sorts them according to their size which accelerates the calculation of the collision/coalescence process. For more details on the parallelization and sorting technique, see the appendix.

Definition and characteristics of a super-droplet.
In the present model, in order to handle the extremely large number of droplets, the concept of a super-droplet is introduced, similarly to Andrejczuk et al [3,15] and Shima et al [16]. A single super-droplet represents a large number of real droplets of the same size, and the number of droplets contributing to a super-droplet is called the 'weighting factor' A n (t). A n (t) can be defined equivalently as the ratio of the total volume of a large number of real droplets represented by a super-droplet to that of a super-droplet.
The super-droplets experience all cloud microphysical processes in the model, and represent the LWC L by where ρ l is the density of liquid water, V is the volume of a grid box, N p is the number of super-droplets in a grid box and r n is the droplet radius. Note that q l obtained from (15) is used to calculate θ v (=θ (1 + 0.608q − q l )).

Advection of droplets.
Since ρ 0 /ρ l 1, the motion of each droplet can be simulated by where X i is the droplet position, V i is the droplet velocity, u i (X i ) is the fluid velocity at the droplet position and τ p is the droplet relaxation time. The inertial response time τ p is calculated using the nonlinear drag law [38]: which yields a good prediction of the terminal velocity [39]. The nonlinear drag law is used because the Stokes law is only valid for droplets with a radius up to 30 µm [40]. Here Re p is the droplet Reynolds number and ν is the molecular viscosity of air. For these calculations, the velocity components of the flow field are interpolated to the current droplet position from the surrounding eight grid points of the respective PALM grid box by bilinear interpolation. SGS velocity fluctuations for the droplet advection are not considered.

2.2.3.
Droplet growth due to condensation/evaporation. The droplet growth due to the condensation/evaporation process is calculated using the equation suggested by Mason [41] as where S is supersaturation, e v the vapor pressure and e s the saturation vapor pressure over a plane surface, F k is the term associated with heat conduction, R v is the individual gas constant for water, T is the temperature, K is the coefficient of thermal conductivity of air, F d is the term associated with vapor diffusion and D is the coefficient of diffusion of water vapor in air. This approach neglects the temperature difference between the droplet and its environment as well as the solution and curvature effects on the droplet's equilibrium vapor pressure; thus the activation process of aerosol particles is not implemented. The radius of all super-droplets is initially given by r = 1 µm, which corresponds to the typical size of the smallest droplet/activated CCN. It also corresponds to the typical size of the smallest bin in the spectral bin model [8][9][10]. We take this simplified condition as an initial approach, because the activation process is not important in the evolution of shallow convection and the droplet growth by collision/coalescence, which are the main concern in this work. The activation process, which can actually be more naturally represented in the LCM, can be included in the next stage, as was done by Andrejczuk et al [15].
To evaluate (18), e v and e s are calculated from q and θ, which are computed by the fluid model (3) and (4), and interpolated to the droplet position. The changes of q and θ are induced by the change in LWC from the condensation/evaporation of all droplets within the respective 8 grid box. The change in LWC by condensation/evaporation per unit time is calculated by the particle model as where r * is the droplet radius modified by condensation/evaporation after one time step.

Droplet growth due to collision and coalescence.
Firstly, we reinterpret the collision/ coalescence of two droplets as the growth of the larger droplet and the disappearance of the smaller droplet. According to this interpretation, after each collision/coalescence between two super-droplets, the larger super-droplet increases its radius while keeping the weighting factor invariant and the smaller super-droplet decreases its weighting factor while keeping its radius invariant. For example, if a collision occurs between the super-droplets with the radius R and r (R > r ) and the weighting factors In this way, we can avoid the creation/annihilation of super-droplets after collision/coalescence. In the present model the evolution of droplet spectra is accomplished by calculating dR/dt and A, which are described in the following. Secondly, we calculate the increase in the radius of the collector droplet (dR/dt) using a statistical approach, instead of considering individual collision events between super-droplets, as was done by Andrejczuk et al [3] and Shima et al [16]. That is, we consider the growth by collision/coalescence given the background droplet spectrum and turbulence. It also implies the growth of the mean radius of real droplets belonging to a super-droplet, resulting from a large number of collisions with real droplets within a grid box during t. The present approach is conceptually similar to the mean field theory in which the n-body system is replaced by a onebody system with a chosen good external field that replaces the interaction of all the particles with an arbitrary particle [42].
Following this approach, the increase of the radius of a collector droplet can be calculated using the continuous growth equation as (e.g. [43]) where K (R, r ) is the collection kernel and n(r ) the number density of droplets. In the present model we assume that the droplet distribution is uniform within a grid box and thus neglect the effect of the distance between droplets, as in the spectral bin model. Accordingly, the radius of a collector droplet after one time step by collision/coalescence R * is calculated in the model by where N R is the total number of super-droplets with r n < R within a grid box. In the present model, the generation of the variance of the droplet size within a super-droplet after collision/coalescence is ignored. According to the stochastic collection model (e.g. [43]), the total number of coalescence to this particular droplet with R per unit time experienced by droplets with r (r < R) is proportional to (e.g. [43]) 9 Therefore we assume that the decrease of the weighting factor of each collected super-droplet with r n < R is given by with a constant α which results from the requirement of mass conservation. The conservation of LWC during the collision/coalescence process requires This leads to the evaluation of α at the collision/coalescence of a collector droplet with the size R as This collision/coalescence process is repeated for all super-droplets, except the smallest one, starting from the largest one during one time step. For example, after the super-droplets are sorted according to their size as {(r 1 , A 1 ), . . . , (r N p , A N p )}, the changes in the radius of the largest super-droplet r N p and A n (n = 1, . . . , N p − 1) are calculated by (21), (23) and (25) first.
The resulting values of A n (n = 1, . . . , N p − 1) are then used to calculate the changes of r N p −1 and A n (n = 1, . . . , N p − 2). The process is repeated until A n (n = 1, 2) is used to calculate r 2 and A 1 . As a result, for each super-droplet, except the smallest one, the increase of the radius by (21) is calculated only once during one time step. On the other hand the decrease of the weighting factor by (23) and (25) is repeated as many times as the number of super-droplets larger than itself. In this way the growth of the droplet radius by collision according to (21) and the decrease of the total number of droplets by coalescence to other droplets according to (23) are satisfied simultaneously for a given super-droplet, while LWC is conserved. The numerical technique for droplet sorting is described in the appendix.

Collection kernel
In this work, two different collection kernels are applied: the traditional Hall kernel, which considers only gravitational collision, and the AW kernel, which includes the effects of turbulence on droplet collision parameterized based on the analysis of DNS results. The Hall kernel is expressed in the following form (e.g. [43]): where V t is the terminal velocity of a droplet, and E g is the collision efficiency. The coalescence efficiency is assumed to be unity. The terminal velocities used in (26) are calculated by semiempirical formulae of Beard [44]. The AW kernel incorporates the effects of turbulence by extending the Hall kernel as where | W Rr | is the average radial relative velocity between two droplets with radius R and r , g Rr is the radial distribution function and η E is the turbulent enhancement factor for the collision efficiency E g . It includes three turbulence effects: | W Rr | describes how the relative velocity between two colliding droplets is modified in a turbulent flow field, g Rr represents the effect of preferential concentration on the droplet distribution due to the turbulent flow and η E describes how turbulence affects the flow field around droplets and thus influences the collision efficiency. For | W Rr | and g Rr , Ayala et al [17] developed a parameterization based on theoretical formulations and compared it with DNS results. These parameterizations are provided as functions of the dissipation rate ε and the Taylor-microscale Reynolds number Re λ (= u λ/ν = u 2 √ 15/(εν); u is the flow rms fluctuation velocity and calculated from ε in the AW kernel). The enhancement factor η E is obtained in terms of ε by analyzing the hybrid DNS results of Wang and Grabowski (see table 1 in [21]). The information on ε and Re λ is provided by LES of the flow field at each grid box in the present model.
Note that, in the absence of turbulence effects, | W Rr | = | V t (R) − V t (r )|/2 and g Rr = η E = 1, and the AW kernel is reduced to the Hall kernel. It should also be mentioned that the AW kernel has been developed and validated using (hybrid) DNS at relatively low values of Re λ , and its validity needs to be confirmed for higher values of Re λ , corresponding to observed clouds (e.g. [45]). Figure 1 shows the ratio of the AW kernel to the Hall kernel for ε = 1 × 10 −2 and 4 × 10 −2 m 2 s −3 . A noticeable enhancement occurs for droplets less than 100 µm, and the overall enhancement is moderate with a value ranging from 1.0 to 5.0. Here the effect of preferential concentration g Rr has the largest influence, followed by the turbulence enhancement on the collision efficiency η E and the radial relative velocity | W Rr | [46].

Particle-turbulence interaction
It is important to remark that the particle-turbulence interaction associated with the fine-scale turbulence structure (or the SGS turbulence), such as the droplet acceleration associated with the intermittency of turbulence, preferential concentration and supersaturation fluctuation (see, e.g., [4]), is not treated explicitly in the present LES model, and its effect is supposed to be included only through the SGS parameterization, such as the collection kernel. Only the effect of preferential concentration is parameterized in the collection kernel so far in the present model.
Based on the assumption of uniform droplet distribution within a grid box, only the droplet size distribution within a grid box, not the detailed information on droplet location within a grid, is relevant for determining the droplet growth from both condensation and collision, in the present model, as in the spectral bin model. The flow field generated by LES mainly contributes to the transport of droplets, and the motion of droplets affects the fluid motion only thermodynamically; that is, modifications of the grid-averaged values θ , q and q l by condensation/evaporation, represented by in (3) and (4), lead to the change of θ v in (1). Here is calculated by the change of total volume of droplets from condensation/evaporation in a grid box in the particle model (19). This implies that both the evolution of droplet spectra and the particle-fluid coupling are not sensitive to the exact location of droplets within the grid in the present model. Therefore, the SGS velocity fluctuation of droplets (e.g. [48,49]) and the higher-order interpolation schemes for the velocity field (e.g. [50,51]) are not likely to play a significant role in the present model, unlike the case of DNS simulating fine-scale particle-turbulence interaction.

Simulations of an idealized single cloud
In the first step, the simulation of an idealized single cloud in the form of a rising warm air bubble is used to examine the new LCM scheme. The bulk model is also carried out for comparison and hereafter called BULK. Furthermore the LCM is simulated with two different collection kernels, the Hall kernel and the AW kernel, with the aim to investigate the effects of turbulence on the droplet growth. These two simulations are called LCM HALL and LCM AW, respectively.

Simulation setup
A 2D rising warm air bubble is triggered by an initial potential temperature difference θ * given by  [22]. These are the same profiles as those used in the LES intercomparison experiment [52]. Note that there is no background wind for this idealized case. The BOMEX case will be discussed further in section 4. The LCM simulations use a constant time step of 0.1 s, while the BULK simulation uses a variable time step with an average value of around  3 s. The very small time step is required in the LCM simulation for the accurate calculation of condensation given by (18). Note that the typical length and time scales of large eddies in moderate cumulus convection are estimated as 10 2 m and 10 2 s [4]. The super-droplets are released at the beginning of the simulation and uniformly distributed all over the model domain, up to a height of 2800 m 4 . The average distance of the super-droplets is initially 4.5 m, yielding a total number of roughly 1.5 × 10 8 and about 87 super-droplets per grid box. Using a weighting factor of 9 × 10 9 a droplet concentration of approximately 100 cm −3 is represented. All initial droplets have the same radius r = 1 µm. For an overview of the main simulation parameters see table 1. The simulated time for all simulations is 30 min, during which the cloud reaches the inversion layer and fully develops.   averages are compared with it too. For the vertical profiles, w 2 and θ 2 are plotted instead of w and θ . Here the cloud reaches the maximum rising velocity at t ≈ 1200 s before impinging on the inversion layer at z = 1500 m.

Cloud evolution
Both models reproduce the typical pattern of an idealized single convective cloud (e.g. [11,12,14,[53][54][55]). At t = 1200 s the cross section of the vertical velocity shows the dipole pattern associated with a rising air bubble. The maximum values of θ and L appear at the center of the cloud, since the entrainment and mixing of environmental dry air decrease the values near the boundaries. The negative θ near the top of the cloud is due to evaporation of droplets through the mixing of the cloudy air with the entrained dry air. The maximum value of e appears near the top of the cloud, across which strong entrainment occurs, and the maximum value of RH appears in the lower part of the cloud, which is affected by the vertical advection of humid air with lower θ from below. The cloud spreads horizontally after impinging on the inversion layer at t = 1600 s. The upward motion comes to an end, and the lower portion of the cloud is occupied by negative values of θ because of intensive evaporation. The results indicate that the LCM reproduces the evolution of the idealized cloud consistent with the bulk model.
Nevertheless, some differences in cloud evolution appear between the two models owing to the difference in cloud microphysics schemes. Unlike LCM, BULK cannot get RH larger than 100% due to the all-or-nothing condensation scheme in which the difference between the total water content and the saturation value is instantly converted into q l . Therefore at t = 1200 s the BULK simulation shows larger values of L, thus inducing larger θ due to the release of latent heat and consequently larger w and e. On the other hand, at t = 1600 s, L in BULK becomes smaller than in that LCM, because droplet growth by collision at a later stage makes the evaporation in LCM smaller than that in BULK. A similar contrast is also found in the comparison of cumulus cloud models with the one-moment bulk model and spectral bin models [14,53].
The distributions of corresponding variables in the horizontal cross section at z = 1400 m and t = 1200 s, shown in figure 5, reveal the turbulent nature of cumulus convection. The length scale of fluctuations is comparable to the cloud itself. The location of the maximum w, L, θ and RH, which appear in the most interior region of the cloud, also fluctuates accordingly. However, the values of SGS-TKE are lower in this most interior region. Meanwhile, the fluctuation of LWC in LCM is stronger in magnitude and smaller in the horizontal scale than the fluctuation of LWC in BULK.
The LCM with the AW kernel produces almost the same results as the results for LCM HALL shown in figures 3 and 4, implying that the effect of microphysics on the dynamical process of shallow convection is insignificant in this case. Figure 6 shows the distribution of the mean radius of the droplets at each grid point. The maximum values, which appear in the center of the cloud, are more than two times larger in the AW kernel panel (d) than in the Hall kernel panel (c) at t = 1600 s, although no noticeable difference is observed at t = 1200 s panels (a) and (b). One can also notice that in both cases the mean radius in the center of this idealized cloud is several times larger than that near the edge of the cloud, which can be understood from the higher level of supersaturation and dissipation rate at the center, as can be expected from figures 3 and 4.  The difference in droplet growth is more clearly identified in the mass density distribution g(ln r ), shown in figure 7. Here g(ln r ) is obtained using the mass of the super-droplets inside the cloud, multiplied by the respective weighting factor. The cloud is defined as the region where L > 1 × 10 −2 g m −3 . Both distributions show almost identical shapes at t = 800 s, but a difference starts appearing for larger radius in the spectrum at t = 1200 s. From t = 800 to 1200 s the radius with the maximum g(ln r ) increases to 20 µm, while the spectral band becomes narrower, suggesting the dominance of the condensational process. From t = 1200 to 1600 s, however, the spectral transfer of droplets occurs from r ≈ 20 µm to r > 50 µm, and leads to the bimodal distribution, which characterizes the typical collision/coalescence process (e.g. [56]). This transition to the bimodal distribution is clearly identified in the case of the AW kernel with the second peak near 100 µm.

Effects of turbulence on microphysics
Using a spectral bin model, Xue et al [57] showed that the AW kernel enhances the droplet growth significantly in comparison with the Hall kernel, especially in the region between 20 and 100 µm, and thus shortens the rain initiation time by 20-40%. For example, when ε = 300 cm 2 s −3 and u = 202 cm s −1 , the percentage of mass for droplets larger than 100 µm at t = 1800 s is 76 and 0.32%, respectively, for the AW and Hall kernels. Although a direct comparison is difficult because of different conditions, the present result, shown in figure 7, appears to be consistent with these results. Using parameterizations of the autoconversion and accretion processes based on the AW kernel for a bulk microphysical model, Seifert et al [46] also showed that rain formation is significantly enhanced when including the turbulence effects on collisions in LES of trade wind cumuli.

Sensitivity to model parameters
In order to verify the validity of the new model it is important to examine the sensitivity of the model solution to the numerical parameters of the model. The parameters that are most likely to affect the model results are the time step size and the number of super-droplets. It is important to find optimum values, which produce accurate results with minimum computational demands (e.g. computing time and memory space), considering the high computational demand of the present model. More accurate model results are expected with a smaller time step and a larger number of super-droplets. Meanwhile, the computational demands increase drastically with decreasing time step size and increasing the number of super-droplets.
The time step size is expected to influence the condensation/evaporation process by affecting the level of supersaturation. On the other hand, the number of super-droplets is important in the calculation of the collision/coalescence process, because the number of superdroplets per grid box has to be large enough to properly represent the droplet distribution for the collision/coalescence process. For the sensitivity tests we repeated LCM HALL with different time step sizes and super-droplet numbers, as shown in table 2. The weighting factor is modified according to the change in the number of super-droplets in order to keep the number of real droplets invariant. The number of super-droplets within a grid box is changed accordingly; for example, from ≈186 super-droplets for the largest super-droplet number (HALL LN) to ≈26 droplets for the smallest one (HALL SN2) (see table 2).   ST) are very close, although the first peak near 20 µm is slightly smaller in HALL ST. Small differences for larger droplets (>30 µm) imply that the effect on the collision/coalescence process is only small. On the other hand, the mass density distribution with the longer time step (HALL LT) shows not only a decrease of the first peak near 20 µm, but also a significant modification for larger droplets, indicating that large time steps affect the droplet growth due to collision/coalescence as well.
The mass density distributions after 1600 s in figure 8(b) show that the results with different numbers of super-droplets are quite similar, except for the case with the smallest number of super-droplets (HALL SN2). In particular, the largest droplet radii in HALL SN2 are much smaller than those in the other simulations. This suggests that an initial number of 26 superdroplets per grid box is too small to provide a sufficient droplet statistic.
Overall, the sensitivity studies show that a time step of 0.1 s and an initial number of ≈80 super-droplets per grid box are approximately the limiting values for robust model results. The information on the optimum values of these parameters is important, because the computational demands increase drastically with a decrease in time step size and an increase in the number of super-droplets in terms of computing time and memory requirement.
To further validate the presented LCM scheme, we carry out a simulation for the shallow cumulus case BOMEX [22]. It is a trade wind cumulus case whose behavior is remarkably steady, and for which there were no apparent complications from precipitation or mesoscale advection. Furthermore, the BOMEX observation was used for an intercomparison of LES (Siebesma et al [52]), so the simulation of this case allows us to examine whether the new model produces results that are consistent with those from existing models. All LES models participating in the intercomparison study used a simple one-moment bulk model similar to our bulk model [52].

Simulation setup
The initial and boundary conditions of all our BOMEX simulations closely follow those used by the LES intercomparison study (Siebesma et al [52]), which simulated the BOMEX observations (see [22]). The initial profiles of the potential temperature, the specific humidity and the horizontal velocity components have already been presented in figure 2. Furthermore, an initial profile for the TKE is prescribed as 1 − z/3000 m 2 s −2 . The surface boundary conditions include a constant surface sensible heat flux of 8 × 10 −3 K m s −1 , a constant surface latent heat flux of 5.2 × 10 −5 ms −1 and a total momentum flux determined by u * = 0.28 ms −1 . Large-scale forcing terms such as subsidence, radiative cooling and low-level drying are also prescribed following Siebesma et al [52]. To trigger convection, small random perturbations of the potential temperature and specific humidity fields are initially added to the lowest 1600 m. All simulations use a model domain of 6400 m × 6400 m × 3000 m. Siebesma et al [58] used a grid spacing of 100 m in the horizontal and 40 m in the vertical direction. However, a much higher resolution with a grid size as small as 20 m is inevitable in the LCM in order to perform the collision/coalescence process properly under the assumption of homogeneity of droplet distribution within a grid box. On the other hand, it is known that some quantities from the LES of trade wind cumuli are sensitive to the grid resolution, especially the size of the clouds, the cloud cover and the average LWC [59,60]. Therefore, we perform the simulations of LCM and BULK with a grid spacing of 20 m, hereafter referred to as HR LCM and HR BULK, and perform an additional BULK with standard resolution equivalent to that used in the intercomparison [52], hereafter referred to as SR BULK, that is, 100 m in the horizontal and 40 m in the vertical direction. Matheou et al [60] suggested that grid convergence is achieved for a grid spacing of 20 m for the case of nonprecipitating shallow cumulus convection. Only the Hall kernel is used as a collection kernel in the LCM in this case. The average time step is about 4 s for the SR BULK simulation and about 1.6 s for the HR BULK simulation. The HR LCM case uses a fixed time step of 0.1 s. An overview of the main simulation parameters is given in table 3.
The settings of the super-droplets are close to those used for the idealized case in section 3.1: all super-droplets are released at the beginning of the simulation, they are uniformly distributed all over the model domain up to a height of 2800 m and have the same initial radius of r = 1 µm. For this case the initial super-droplet distance is 4 m. With a weighting factor of 1 × 10 10 this results in a represented droplet concentration of approximately 150 cm −3 and a total number of super-droplets of roughly 1.7 × 10 9 (about 125 per grid box).   [52]. The gray bands around this mean have a width that is twice the standard deviation of the model results of Siebesma et al [52]. Figure 9 shows the temporal evolution of the total cloud cover, the liquid water path (LWP) and the vertically integrated TKE for our simulations in comparison with the results of Siebesma et al [52]. The total cloud cover is defined by Siebesma et al [52] as the fraction of vertical columns that contain liquid water, and the same definition is applied to the cases of SR BULK and HR BULK. In the HR LCM simulation, however, a small amount of water exists in every grid box because of the initially uniform distribution of super-droplets. Thus, in this case, only model columns that contain at least one grid cell with L > 1 × 10 −2 g m −3 are defined as cloudy  [52]. The gray bands around this mean have a width that is twice the standard deviation of the model results of Siebesma et al [52].

Comparison with the intercomparison results
columns. As reported by Siebesma et al [52], all models are in a spin-up phase during the first two hours of the simulation. Good agreement is found between the present simulations and Siebesma et al [52] in the cases of LWP and TKE. The cloud cover is still in a good agreement with that in Siebesma et al [52], while the resolution is the same as that in Siebesma et al [52] (SR BULK). However, the cloud cover from the high-resolution models, HR LCM and HR BULK, is much higher than that in [52]. According to the recent sensitivity experiments of Matheou et al [60], the cloud cover increases substantially with increasing model resolution in LES of shallow convection. They performed LES of trade wind cumuli from the RICO (Rain in Cumulus over the Ocean, [61]) field experiment and observed an increase of the cloud cover with grid spacing from about 7% for x = 80 m to about 25% for x = 20 m. Even though these simulations are not directly comparable to ours, a total cloud cover of 20-25% in our high-resolution simulations is consistent with that found by Matheou et al [60] for a grid spacing of 20 m.
A selection of mean profiles averaged over 1 h (hours 2-3) is shown in figure 10. Siebesma et al [52] found excellent model-to-model agreement in the vertical profiles of potential temperature and specific humidity, but a rather large variation in LWC. They attributed this to the low cloud cover in the simulations, which yields a strong dependence of LWC on a relatively few saturated grid points. The vertical profiles of potential temperature and specific humidity from the present simulations also show good agreement with Siebesma et al [52], but LWC shows large variability.
Matheou et al [60] showed that the LWC in the upper parts of clouds increases with grid resolution, which is related to the increase of cloud cover. In agreement with this, the present high-resolution simulation results, HR BULK and HR LCM, also show the larger LWC in the upper parts of clouds. It is found in particular that the LWC from HR LCM is substantially larger in the upper part than other results. It should be mentioned that none of the models participating in the intercomparison used a sophisticated microphysics scheme such as the spectral bin model and therefore the growth of the droplet size was not taken into account. As shown in the idealized single-cloud experiment (figures 2 and 3), the growth of the droplet size suppresses evaporation, and this may lead to higher LWC in the upper parts of clouds.
In conclusion, the present simulation results, including LCM, show good agreement with the previous intercomparison results of the LES of shallow convection [52]. The disagreement seen in the cloud cover and in the vertical profiles of LWC may be attributed to the insufficient resolution and the lack of microphysics scheme in the models participating in the intercomparison.

Conclusion
In this paper, we introduce a new LCM in which the flow field is calculated by LES and the droplets are treated as Lagrangian particles, and we show that the new model produces results that are consistent with those of existing cloud models. LCM simulations of an idealized single cloud and trade wind cumuli (the BOMEX case) reproduce general features of shallow cumulus clouds that are in agreement with traditional bulk models. The droplet spectrum simulated by the LCM, using the collection kernels with and without the effect of turbulence (the Hall kernel and the Ayala-Wang kernel), is also consistent with those predicted by the spectral bin model. The sensitivity tests of the LCM for two model parameters, the time step and the number of super-droplets, show that a maximum time step size of 0.1 s and a minimum number of about 80 initial super-droplets per grid box are necessary for the convergence of solutions.
In the new LCM, in order to handle the extremely large number of droplets within a cloud, the concept of a super-droplet, which represents a large number of real droplets of the same size, is introduced, and the number of contributing real droplets is called a weighting factor. A novel method is developed to realize the collision/coalescence of droplets, in which the consequent variation of the droplet spectrum is represented in terms of the modification of the radius and weighting factor of super-droplets while keeping the number of superdroplets invariant. Furthermore, we use a statistical approach to calculate the droplet growth by collision/coalescence given the background droplet spectrum and turbulence, in contrast to Andrejczuk et al [3] and Shima et al [16] who considered individual collision events between super-droplets.
This work illustrates that the new LCM can be a valid and promising approach for the simulation of clouds. Once it is developed successfully, the LCM, which is a more natural approach to simulate particle-laden turbulent flows such as clouds, has many advantages over traditional cloud models. In the LCM, the dynamics and microphysics of a cloud are directly related to the physical processes of individual droplets; for example, precipitation and condensation are naturally realized by individual droplets, and the evolution of the droplet spectrum is determined directly by the growth of individual droplets.
Meanwhile, for more realistic simulation of clouds, further improvement of the model in various aspects needs to be made in the future. The collision/coalescence scheme and the collection kernel should be verified definitely and further elaborated, possibly by using information from DNS with high Reynolds numbers. It is also necessary to include the activation process of CCN in the LCM. The efficiency of the numerical code must be improved to overcome the extremely large number of super-droplets. Finally, for a more rigorous verification of the model, the LCM results must be directly compared with observation data: especially for observations including precipitation, in which the advantage of the Lagrangian model can be more clearly identified (e.g. the Rain in Cumulus over the Ocean (RICO) field study [61,62]), or for observations with high-resolution turbulence measurements (e.g., ACTOS (Airborne Cloud Turbulence Observations System) [63]). dimensioned as P j (N ), where N is the number of particles in use. At the initial stage of a model run, the particle data are sorted following the sequence of grid point data stored in the memory. For example, P j (1 : 10) contains data of those particles located in grid box (1,1,1), P j (11 : 20) those of particles in (2,1,1) and so on. If the number of particles in each box is known, the particles of this box (i.e. their respective array index) can be easily accessed. The second advantage of storing particles in this order is that it is computationally very efficient because it takes advantage of the cache mechanism avoiding time-consuming cache misses. If, for example, the new particle coordinate x for the next time step is calculated following (16) in a loop like DO n = 1, N P(n)%x = P(n)%x + P(n)%speed_x * EXP(-delta_t/tau_p) & + u_int * (1.0 -EXP(-delta_t/tau_p)) ENDDO where speed_x is the particle velocity in the x-direction, delta_t is the time step, tau_p is the particle relaxation time τ p and u_int is the u-component of the velocity field of the fluid as interpolated from the neighboring grid points, then the u-velocity grid point data are accessed in this loop in the order in which they are stored in memory.
During the simulation, due to turbulent mixing, the correlation between the particles and the grid boxes gets lost. Although the data of the first particle are still stored in P j (1), the particle now probably stays in a box far away from that where it was initially started. Carrying out the above loop now requires more or less random access to the memory to obtain the velocity data needed for interpolation. This causes a large number of cache misses which significantly slows down the computational speed, especially for large grids whose data do not fit into the cache. Therefore a re-sorting of the particles is performed after each time step to avoid the cache miss. The sorting is carried out in such a way that the original correlation between the sequence of particle data and the sequence of grid point data stored in memory is re-established. This method yields a great advantage, especially for the calculation of the collision/coalescence process. To determine the particles in the vicinity of a particle under consideration would need a nested loop search like DO n = 1, N DO k = 1, N IF ( k /= n ) THEN IF ( ABS( P(k)%x -P(n)%x ) < threshold ) THEN ...

ENDIF ENDDO ENDDO
if the particles are not stored in a consecutive order and would require an order of N 2 operations. If particles are sorted according to the grid boxes, however, the particle collisions within these boxes can be calculated very fast using a loop like DO n = n_start(i,j,k), n_end(i,j,k) ... ENDDO where n_start and n_end are the lower and upper indices of those particles within grid box (i, j, k). This method only requires a number of operations of the order of N .
To further simplify and speed up the calculation of the collision/coalescence process, the particles are additionally sorted according to their size within each grid box. In our model the collision/coalescence process is calculated using the continuous growth equation (see (20)). Thus the growth of each particle R depends on the number and size of all particles with r < R within the same grid box. Following that, each calculation would need a nested loop search, similar to the one shown above, to determine all smaller particles. If the particles are already sorted according to their size, however, the necessary information can be easily accessed without further calculations. Therefore the additional particle sorting is done at every time step in each grid box.