Turbulence and cloud droplets in cumulus clouds

In this paper, we report on the successful and seamless simulation of turbulence and the evolution of cloud droplets to raindrops over 10 minutes from microscopic viewpoints by using direct numerical simulation. Included processes are condensation–evaporation, collision–coalescence of droplets with hydrodynamic interaction, Reynolds number dependent drag, and turbulent flow within a parcel that is ascending within a self-consistently determined updraft inside a cumulus cloud. We found that the altitude and the updraft velocity of the parcel, the mean supersaturation, and the liquid water content are insensitive to the turbulence intensity, and that when the turbulence intensity increases, the droplet number density swiftly decreases while the spectral width of droplets rapidly increases. This study marks the first time the evolution of the mass density distribution function has been successfully calculated from microscopic computations. The turbulence accelerated to form a second peak in the mass density distribution function, leading to the raindrop formation, and the radius of the largest drop was over 300 μm at the end of the simulation. We also found that cloud droplets modify the turbulence in a way that is unlike the Kolmogorov–Obukhov–Corrsin theory. For example, the temperature and water vapor spectra at low wavenumbers become shallower than k − 5 / 3 in the inertial-convective range, and decrease slower than exponentially in the diffusive range. This spectra modification is explained by nonlinear interactions between turbulent mixing and the evaporation–condensation process associated with large numbers of droplets.


Introduction
Weather forecasting and accurate predictions of future climate conditions are becoming increasingly important for human activities, and numerical simulations with initial conditions generated from detailed data provided by state-of-the-art high-technology equipment play a central role in both weather and climate predictions. However, because the minimum grid size in global simulations is limited to a few kilometers, even when petascale high-performance supercomputers are used, it is necessary to use certain hypotheses when modeling or parameterizing turbulent motion and cloud microphysical processes smaller than the grid scale. Although there has been a great deal of research aimed at understanding the dynamics of cloud formation and turbulence, and to develop more accurate models [2,15,30,32,47,62,64,80], our knowledge of these processes is still limited because turbulence is still not fully understood [13,22,82], and it is difficult to conduct direct measurements of cloud microphysical properties [66,67,75]. Thus, turbulence and clouds remain as significant barriers to understanding atmospheric dynamics, forecasting weather, and predicting the future climate.
One reasonable way to overcome these barriers is to simulate turbulence and cloud microphysical processes from microscopic viewpoints, which means to numerically compute all of the motions of the degrees of freedom as faithfully as possible with respect to their fundamental physics. Although this approach requires large computational resources, research trends in this direction have been growing since [76,77] due to ongoing increases in computer power. A basic approach to such problems is to consider a small parcel inside a quasiadiabatic cloud core, and to examine the evolution of cloud droplets at ambient temperature and the water vapor mixing ratio under turbulent mixing conditions by using direct numerical simulation (DNS). The central issue in such research has been the growth of the droplet spectrum (size distribution), especially the broadening of the distribution over time. Turbulence fluctuations affect the spatial distribution of cloud condensation nuclei (CCN) and cloud droplets, thereby enhancing droplet clustering via intermittency [5,6,8,11,15,18,19,57,58,76], and supersaturation fluctuations broaden the size distribution of cloud droplets depending on the number density of the aerosols [10,55].
However, the studies conducted to date have only examined a few of the relevant microphysical processes, and almost no studies that include all of the key microphysical processes have been published. For example, collision-coalescence is not included in [10,55,76,77], and while the effects of turbulence on the collision kernel have been studied theoretically and numerically, the effects of coalescence are quite often neglected [3,4,51,72,78,79]. In addition, while hydrodynamic interaction effects on collision-coalescence have been examined in [46], condensation-evaporation was not included.
Since it is natural and important to fully examine and understand the evolution of droplets and turbulence inside clouds when all (or most) of the key microphysical processes are included under conditions that closely match actual circumstances, this study seeks the answers to a number of questions including, but not limited to: (1) How can the supersaturated state be maintained during the evolution of small droplets from a few μm to around 30 μm? (2) When does the collision-coalescence process becomes dominant over condensation growth? (3) How does strong turbulence affect the evolution of low order moments and the size distribution of the cloud droplets? and, (4) How is turbulence modified by the cloud droplets and the associated heat release? To answer these questions, it will be necessary to conduct an extremely large numerical simulation that examines essential key processes, from cloud droplets to raindrops, from first principles viewpoints.
It is also necessary to estimate the amount of numerical computations that will be required for the above research. To some extent, the space-time range in the microscopic simulation considered in the present context needs to be overlapped with that of the macroscopic simulation, so that the statistical data obtained by the microscopic computation can be connected to the macroscopic parameters, coefficients, or functions.
However, in space and time scales, there are large computational gaps between the microscopic and (semi) macroscopic approaches. To overcome these gaps, a very long integration time that allows sufficient spatial scale range for both turbulence and droplets is necessary. For example, in order to connect with the (semi) macroscopic simulation like the super droplet method which uses a 4 m cube grid for 30 minutes [64], it is roughly estimated that about 4000 3 grid points and more than 16 million time steps (about 10 000 large eddy turnover times) will be required. Although such a large computation is beyond the scope of the present study, we are optimistic from our experiences in large-scale DNS of turbulence [26,27,33,34] that it will be possible in the future. Accordingly, this paper describes some of the preliminary results from steps taken in that direction.
We will begin by discussing a newly developed model, named as a cloud microphysics simulator, that can seamlessly compute the evolution of turbulence and cloud droplets up to raindrop formation within a cumulus cloud from microscopic viewpoints [25] (hereafter referred to as GSS). The method shares with [29,42] in the introduction of a small cubic box ascending within a cumulus cloud, but differs in the following points: (1) the updraft velocity of the box is self-consistently determined by the latent heat release, (2) the supersaturated state is sustained by the ascending motion through the self-consistently determined mean buoyancy, (3) the dynamics of the cloud droplets is directly integrated instead of using the bin method, and (4) the turbulent air flow, temperature and water vapor mixing ratio within the box are computed by DNS. In addition to the above points, in the present study we newly include the following items: (5) the collision-coalescence process, (6) the hydrodynamic interaction effects in terms of Hall's table [31], and (7) the Reynolds number dependent drag for the droplets. On the other hand, the activation of CCN is not included to avoid an excessive amount of computational work, and the entrainment effects are introduced through one parameter.
In our simulation, we numerically and seamlessly integrate the dynamics of the turbulent flow, temperature, water vapor, and droplets with those microphysical processes over long time periods, and then examine the evolution and interaction of turbulence and droplets in order to obtain answers to the questions raised above. To the authors' best knowledge, this type of numerical study, which covers the entire evolutionary process, has never been conducted previously.
The remainder of this paper is organized as follows. Section 2 describes the numerical model, and section 3 describes the numerical setup. Section 4 presents the results of the simulation, and section 5 provides the summary and discussion.

Model description
While the system under consideration is basically the same as that used in GSS, there are some modifications. Accordingly, for the benefit of our readers and later arguments, we will provide a concise description of those modifications below.
We begin by considering a small air parcel ascending inside the core region of maritime cumulus cloud. The parcel is a cubic box with lengths of L box per side, and is assumed to be much smaller than the size of the entire cloud so that statistical properties of fluctuating quantities inside the box can be regarded as homogeneous and periodic boundary conditions in three directions are imposed on the flow field. The governing equations are those for the parcel, turbulent flow, and droplets.

Parcel
The mean altitude ( ) H t and the mean updraft velocity ( ) W t of the box are governed by the following equations: where the overbar denotes the volume average over the box, and B is the mean buoyancy force acting on the box. The parcel is assumed to vertically move without drag in the prescribed environmental temperature and water vapor mixing ratio (ˆ) T z e and (ˆ) Q z ve , respectively, where = (ˆˆˆ) x x y z , , represents a coordinate system with its origin fixed on the sea surface [49,50,63]. A local coordinate, = ( ) x x y z , , , moving in tandem with the box is introduced, and is related to the coordinate system fixed on the sea surface by the time-dependent Galilean transform as 3 z respectively, where e z is the unit vector in the vertical direction.

Flow field
The fluid inside the box is assumed to be incompressible under the Boussinesq approximation. The temperature T and the water vapor mixing ratio Q v are expressed as the sum of the mean and fluctuation as respectively. The evolution of the means is determined by where L v and c p are the latent heat and specific heat of air at constant pressure, respectively, and C d is the condensation rate defined below. The fluctuating turbulent velocity = ( ) u u u u , , 3 , temperature θ, and water vapor mixing ratio q are assumed to obey in the local coordinate system of the box, respectively, where p is the pressure, and f is the external force whose properties are described in section 3.2. The kinematic viscosity of dry air n a , as well as the diffusivity of temperature and water vapor k T and k v , are the same as those in GSS. The buoyancy force B is given by , and R d and R v are the gas constants for dry air and water vapor, respectively. The term with G( ) t in (9) denotes the cooling effect imposed by the ascending motion and is assumed to be the sum of two contributions as Here, G d is the dry adiabatic lapse rate, G ( ) t ent is the effect due to entrainment, and μ is the entrainment strength parameter [1,28,[38][39][40].
Droplets affect the evolution of fluid through the term C d in (9) and (10) (representing latent heat release and mass exchange through condensation, respectively), and the term q l in (11) (representing drag force due to the weight of the condensed water). These terms are assumed to be given by are the radius and the position for the jth droplet, respectively, D ( ) x N t , is the number of droplets in the grid cell D ( ) x 3 , S is the supersaturation described below in (23), and K s is a temperature-dependent diffusion coefficient that is assumed to be constant because its dependence on temperature is very weak (GSS).

Droplet
The cloud droplets are treated as point particles since the droplets studied in the present computation are assumed to be smaller than the grid size. The evolution equations for the jth cloud droplet are given by is the velocity of the jth droplet, t ( ) t j is the relaxation time for Stokes drag: , 2 0 j j a l 2 a and ξ is the nonlinear drag coefficient [61] given by and Re p is the particle Reynolds number defined as (15) and (19)) is the supersaturation at the droplet position defined by The saturation mixing ratio for water vapor Q vs at the droplet position is determined by Tetens' formula:

Collision-coalescence
Since the droplet number density is assumed to be small, only binary collisions are considered. Upon collision and coalescence, the mass and momentum of droplets are conserved. For collision detection between time ( ) t n and + D ( ) t t n , the trajectory of the droplet is calculated by the following linear interpolation: is defined as the separation distance between the centers of the jth and kth droplets at time t. According to [79], two droplets j and k are assumed to collide if both of the following inequalities are satisfied: In the present simulation, the right-hand sides of the above two inequalities are approximated by because the change in droplet radius within one time step is very small. The effects of the hydrodynamic interaction on the collision are counted by using the Hall table [31]. When geometric collision is detected for a certain particle pair according to the above detection scheme, collision efficiency of the pair due to hydrodynamic interaction is taken from the Hall table (table1 of [31]) based on the radii of two particles. If the efficiency is larger than where c is the offset and R R , j k are the radii of jth and kth particles,respectively, the particles are assumed to collide. Additionally, upon collision, coalescence is assumed to occur with a probability of unity [3,20,46,78].
The Hall table for the collision efficiency is computed under the following assumptions: (1) the droplet is a rigid sphere, (2) the flow field around the sphere is laminar, (3) the droplet Reynolds number is small, (4) the computation is based on the terminal velocity due to gravity alone, and (5) the coalescence efficiency is unity [31]. Our computation shares (1) and (5) with Hall's assumptions, but not the others; the flow is turbulent, the droplet Reynolds number, which is initially small, becomes larger as evolution progresses, and collisions are determined directly based on the actual trajectories of droplets through the turbulent flow [3,4,56]. The assumption of the coalescence efficiency being unity is acceptable when we consider that the most of droplets in the present study are smaller than 100 μm, and that the surface tension is high enough to prevent the droplets from fragmenting [54,60,85]. The efficiency in the Hall table is quite small, and below 0.5 for droplets with a radius of about 30 μm, which are computed for the laminar flows. However, since other studies have reported that the turbulence enhances droplet collisions, it is expected that the collision efficiency of droplets in turbulent flows would be between those for the Hall table and unity [5,6,8,9,11,14,18,19,[57][58][59].
In the above model, the activation of CCN is not included, because it would add an excessive amount of computational work. The collision-coalescence computation for cloud droplets in excess of 16 million droplets over a million time steps, as described later, consumes considerable computational time and it is not certain whether such a computation would be feasible, even if a high-performance supercomputer were to be used. Therefore, in this study, we only include collision-coalescence with the condensation-evaporation of cloud droplets (GSS), and focus on the development of efficient computational methods that are capable of accomplishing long-period integration. Without the activation of CCN, the number density of cloud droplets decreases due to collision-coalescence, and as the parcel continues to ascend and the temperature decreases, supersaturation is expected to increase [29]. Simulations are stopped before supersaturation becomes too large, as is done in bin-model simulations by [29].

Environmental and initial conditions
The same environmental profiles for the temperature and water vapor mixing ratio as those in GSS are used. These profiles are modeled from observations taken in Hawaii in [74]. The environmental temperature T e and water vapor mixing ratio Q ve are given as functions of pressure P by  respectively, where P 0 is the pressure at the initial height of the box ( = ( ) z H 0 ). An inversion layer exists above the 780 hPa level (corresponding to about 2200 m). The relationship between the pressure P and the vertical coordinateẑ is given by respectively. In the box, the initial value of the average temperature is set to be slightly higher than the environmental temperature as The mean water vapor mixing ratio is assumed to be saturated and is determined by Tetens' formula as follows: The initial conditions for the fluctuations of velocity, temperature, and water vapor mixing ratio are multivariate Gaussian random fields with a zero mean. Their spectra are for the kinetic energy spectrum, and follow same forms as shown above for the variance spectra of temperature q ( ) E k, 0 and water vapor mixing ratio ( ) E k, 0 q , except that u 0 is replaced by q 0 and q 0 , respectively. Here, k 0 is the peak wavenumber, while u 0 , q 0 and q 0 are the initial values of the root mean square of u, θ, and q, respectively. We set where ( ) n t is the mean number density, m(R) is the droplet mass with radius R and liquid water density r l , and . Similarly, the mass density distribution function of the droplet radius is given by is the mean cloud water mixing ratio. The initial value of the mean droplet number density is = = -( ) n N L 0 125 cm p box 3 3 , and the droplets with total number of = N 2 p 24 are randomly distributed in the box. The initial distribution of ( ) D m t , is the same as in [7], and is given by where m f is the mean mass of the number density distribution function and n ( ) G is defined as and Λ is the gamma function. We set n = 0 and

Random forcing
Turbulence is maintained in a statistically steady state by random forcing applied in the velocity field, and all wavenumber components within a shell   kL 1 4 box are subjected to the forcing. The forcing is solenoidal and produced from an Ornstein-Uhlenbeck stochastic process as where the correlation time ζ is set to be z = -( ) u k 0 1 and u 0 is the root-mean-square of the fluctuating velocity field at the initial time [16,17,84]. The term ( ) w k t , is a Gaussian random process [24] given by where F(k) is the forcing spectrum as follows: The factor z ( ) 2 in (39) is introduced so that the mean energy injection rate due to random forcing is approximately independent of the correlation time ζ and is controlled solely by the parameter c f . The values of c f used for experiments are summarized in table 2.

Collision detection
In order to efficiently accelerate the computation of collision detection, the so-called cell index method is used [72,79].
For a selected droplet j, the primary cell is defined as the cell that includes the droplet's position X j at each time step. Collision detection involving this selected droplet is restricted to the primary and adjacent cells. The size of the cell is chosen such that L cell is greater than the distance traveled by a droplet within a single time step.
To check the validity of the above collision detection scheme, the theoretical value of the collision rate is compared to the results of numerical simulations for certain droplet size distributions. These comparisons show that the collision detection scheme in the present simulation can reproduce the theoretical value with sufficient accuracy. The details of the validation are presented in appendix A.

DNS of flow and droplets
The equations for the turbulence are numerically computed using the pseudo spectral method, while the velocity and scalar fields at the droplet position are obtained by linear interpolation. The weighting used in the linear interpolation is also applied to distribute the contribution from the droplet positions on the surrounding grid points. The second-order Runge-Kutta method is used for the time integration. To avoid the momentum equation stiffness for a droplet with a small radius, a method that can analytically solve linear terms, the details of which are described in appendix B, is used.
It should be noted that the evolution of the means is separately computed as (5) and (6). Computation using the total value, such as q = + T T , leads to information loss because the mean value is quite large compared to the fluctuation amplitudes, such as q =´-T 1.7 10 0 0 4 and =q Q 10 0 v0 4 , which in turn yields cusps in the spectrum of the temperature variance near the wavenumber cut off (GSS). The parameters values used are the same as those shown in GSS table 1, except for those indicated in table 1 of the present paper.
Four computations are performed and the numerical and turbulence parameters are summarized in table 2. Runs 1 and 2 are conducted to check the significance of collision-coalescence for droplet growth. Run 2 uses the same parameters as Run 1, except that Run 2 is conducted without collision-coalescence. Runs 3 and 4 are conducted to determine the effects of turbulence intensity on the evolution of cloud droplets. The kinetic energy dissipation rate ò is increased from about 27 to 400 cm s 2 3 by changing the random force intensity. These values are consistent with the observed data that ranges from 10 to 1000 cm s 2 3 [30]. The Taylor micro scale Reynolds numbers are between 104 and 167, and are smaller by one order than in the actual cloud due to the present simulation size. The elementary processes directly related to droplet growth are condensation-evaporation, collision-coalescence, and the updraft.
The supercomputers used in this study were the K-computer at the Research Organization for Information Science and Technology (RIST) in Kobe and the Fujitsu FX100 installed at Nagoya University. In order to Density of liquid water r l 1.0 (g cm −3 ) Diffusion coefficient in (15) and (19) K s´-1.20 10 6 (cm 2 s −1 ) Entrainment strength μ´-2.5 10 5 (cm −1 ) perform long-time integration, parallelization of the computer program was developed intensively. The cubic box domain was decomposed into two dimensions for the velocity and scalar computation using a fast Fourier transform and for the droplets. Both message passing interface and open multi-processing were used for parallel computation of most of parts, including collision-coalescence. The grid number is set to N=512 for each direction, and the grid spacing is D = = x L N 1 mm box for all runs. The time increment is chosen so that the CFL condition for the turbulence is satisfied. The length of time integration is 600 s for Runs 1-4, while the number of time steps is 1.2 million for Runs 1 and 2, 1.5 million for Run 3, and 2.4 million for Run 4. From table 2, it can be seen that the large eddy turnover time T L is about 1.1 s for Runs 1 and 2, 0.71 s for Run 3, and 0.45 s for Run 4. Therefore, the total time integration width corresponds to about T 545 L for Runs 1 and 2, 845T L for Run 3, and 1333T L for Run 4. This computation used 8192 cores and the longest computational time for Run 4 was about 100 h. In our opinion, this central processing unit (CPU) time is acceptable and feasible for our present purposes.

Mean value and size distribution
We first compare the results for Run 1 and 2 in order to check the significance of the collision-coalescence.  1 , and is consistent with the moist adiabatic lapse rate at these temperatures. In the present model, the condensation proceeds and the cloud water mixing ratio increases, but the total water content (the sum of water vapor and cloud water) does not change. For Run 1, supersaturation increases and finally reaches 1 %. The small kink of the green curve in figure 1(d) at about 570 s arises as the parcel reaches the inversion layer at about 2200 m, and such kinks are not seen for other curves due to the rapid growth of the supersaturation. Runs are stopped at 600 s so that supersaturation does not become unrealistically large in comparison to actual observations [68].
Figures 2(a) and (b) show the evolution of the droplet number density and spectral width (the standard deviation of the number density distribution function of droplet radius ( ) F R t , ), respectively. For Run 1, the number density and spectral width first change gradually, and then begin to rapidly decrease and increase, respectively, after about eight minutes. Since the changes until this transition are similar to those for Run 2, where collision-coalescence is not included and droplets grow only by condensation, the eight-minute transition indicates that the dominant process for droplet growth changes from condensation to collisioncoalescence.
The time evolution of droplet size distribution shows the detailed characteristics of condensation and collision-coalescence growth. Figure 3 shows the time evolution of the probability density function (PDF) for the droplet radius, ( ) h R t , , which is related to the number density distribution function ( ) where ( ) n t is the mean number density. For Run 1 ( figure 3(a)), the PDF first narrows with time as the peak shifts rightward. This corresponds to the decrease of the spectral width in figure 2(b), and is similar to the result for Run 2 ( figure 3(b)). Since distribution narrowing is a well-known characteristic of condensation growth for which the growth rate of the droplet radius is proportional to the inverse of the radius, as seen in (19), these results add further support to the conclusion that condensation is the dominant process for the first stage of Run 1 droplet growth. After about six minutes, the PDF for Run 1 gradually broadens, showing a nearly exponential right tail. This broadening is caused by droplet collision-coalescence. Table 2. Numerical and mean turbulence parameters. Collision-coalescence efficiency is shown in the column with the same name, 0: not included, Hall: Hall's table for the collision efficiency is used [31]. l R is the Taylor microscale Reynolds number, E is the kinetic energy, ò is the mean energy dissipation rate per unit mass,  is the integral scale, λ is the Taylor microscale, η is the Kolmogorov length, h k max is the cut off wavenumber normalized by the Kolmogorov length, Dt is the time increment, and c f is the forcing amplitude. The time evolution of the mass density distribution function clearly shows different phases of collisioncoalescence growth in the rain initiation process. Figure 4 shows the time evolution of the mass density distribution function (defined in (36)) for Run 1. The red curve is the mean volume radius and the blue curve is the mass-weighted mean radius, or 'predominant radius', as defined in (11b) in [7]. The mean volume radius and the predominant radius correspond approximately to the primary and second peak of the distribution, respectively.    As droplets begin growing (primarily due to condensation), the cloud water and mean volume radius of the droplets increases. Accordingly, the area of the distribution increases and the primary peak shifts rightward, as shown in figure 4. From about eight minutes, the right tail of the primary peak gradually broadens due to droplet collision-coalescence. The second peak emerges from the right tail after nine minutes and grows rapidly. Near the end of the simulation, the predominant radius reaches 100 μm, the distribution broadens significantly due to collision-coalescence, and the radius of the largest drop is over 300 μm (drizzle).
The droplet collision growth characteristics described above are well-known in bin-model simulations. As described by Berry and Reinhardt [7], the gradual broadening of the right tail of the PDF and the mass density distribution function is caused by auto-conversion, in which the mass of similar size droplets near the primary peak is transferred to larger droplets by a weak accretion mechanism. The auto-conversion phase is followed by the accretion phase, where larger drizzle drops efficiently collect smaller cloud droplets near the primary peak. This corresponds to the emergence and rapid growth of the second peak in the mass density distribution function after nine minutes, as seen in figure 4.

Effects of turbulence intensity
It is argued that turbulence accelerates rain initiation by making droplet collision-coalescence more frequent. For example, Grabowski and Wang [29] applied the turbulent collision kernel model by [4] to bin-model simulations and showed that the time for warm rain initiation could be reduced by 25%-40%.
To investigate the effects of turbulence intensity on rain initiation in the present simulation, we conducted Runs 3 and 4, in which the forcing amplitude was increased such that the mean energy dissipation rate ò was about 100 and 400 cm s 2 3 , respectively. The results of Runs 3 and 4 are shown in figures 1-3, 5 and 6. From these results, it can be seen that the parcel's height, updraft velocity, and cloud water are very insensitive to the turbulence intensity. The dependency of the supersaturation, number density, spectral width, and the PDFs of the radius on the turbulence intensity becomes visible later. On the other hand, the mass density distribution functions at the final state show a clear tendency that indicates stronger turbulence results in faster raindrop formation (figures 5 and 6). To facilitate comparisons, the mass density distribution function for Run 1 at t=600 s is overplotted in figures 5 and 6, where it can be seen that the second peak becomes greater and the right tail lengthens as the turbulence intensity becomes stronger. This tendency is consistent with previous binmodel simulations [29], but the present study is the first to confirm it by DNS of turbulence with the microscopic dynamics of cloud droplets.
At the final period of the present simulation, the large droplet radius exceeds 300 μm whose fall speed is about 2.  scales smaller than r s would return to its steady state during t s after droplet sweep, on the other hand, those motions with larger scales do not and the effects of droplet sweep remain by the time of the next sweep. In other words, those large droplets would experience almost the same (frozen) turbulence on the scale larger than r s and the spatial configuration of those droplets at large scales differs from that of large box size. This implies less chance for those large droplets to collide with the droplets of similar or larger size and might slow down the evolution of the second peak of the mass pdf. The condition that the large droplets experience the different turbulence in the next sweep is > t T s L , that is, < u u fall rms , which means that the Reynolds number must be large. As for the temperature and water vapor mixing ratio, the spatial configuration of those droplets at scales larger than r s would also differ from the case of large box, implying that the spectral tail at low wavenumber would also be affected, but the condensation process for droplets of this size is expected to be very small. Indeed we observed that the evolution of the scalar spectra for Run 2 (without collision coalescence) is almost identical with that of Run 1 (figure not shown). In the present computation, since  » r 2 s and the simulation has not been performed beyond 600 s, we expect that the effects due to the large fall speed on the results are negligible. Figures 7(a)-(c) show the evolution of the kinetic energy spectrum ( ) E k t , , the variance spectra for the temperature q ( ) E k t , , and the water vapor mixing ratio ( ) E k t , q , respectively, for Run 1, while figures 8 and 9 show the same information for Runs 3 and 4, respectively. Since the statistical nature of a turbulence field is axially symmetric, all the spectra in what follows are those of the isotropic parts alone.

Spectra
The evolution of the kinetic energy spectrum is qualitatively similar to that for GSS (figures 11-13 therein), but there are no cusps in the spectra q ( ) E k t , and ( ) E k t , q near the maximum wavenumber. We found that the cusps, as found in GSS, are due to the information loss that arises from the computation for the total field (large mean plus small fluctuation). The kinetic energy spectrum evolves quickly and has already attained a steady state at 10 s, where it remains unchanged except in the far dissipation range. The spectral slope at low wavenumbers (   kL 2 2 0 box ) is close to -5 3, which matches the Kolmogorov spectrum in the inertial range and is consistent with previous turbulence DNSs [24,26,33,34,84]. On the other hand, the spectral tail in the dissipation range gradually rises and tends to have a slope close to −2 for Run 1, which has quite often been observed in previous works that use point particle approximation for the droplets [73]. This trend arises from the q l term of (11). Indeed, a DNS conducted without this term yielded an evolution similar to ( ) E k t , , as shown in figure 7(a) in the inertial range, but the spectrum quickly rolled off as exponentially without such a rise in the dissipation range (figure not shown). However, when the Reynolds number is increased, as seen in figures 8(a) and 9(a), the slow decay of ( ) E k t , in the dissipation range becomes less appreciable, as it is apparently buried by the relatively high excitation of kinetic energy. The spectra q ( ) E k t , and ( ) E k t , q rapidly attain the nearlyk 5 3 spectrum in the inertial-convective range at 10 s, and as time progresses their amplitudes increase at all wavenumbers. However, as time goes on, the spectra become shallower, roughly ask 1  , the spectra roll off roughly ask 2 at later times, which is much slower than an exponential decay. The transition of the spectra fromk 1 3 tok 2 like behavior occurs at around   kL 10 3 0 box . These observations are common for three runs in the present study, but it is difficult to determine the precise spectral form at low and high wavenumbers because the wavenumber ranges of the spectra are not wide enough to accurately determine whether they obey the power law with definite exponents. We consider it likely that the deviation of the spectra from the Obukhov-Corrsin spectrum is due to the terms arising from the condensation fluctuation This is based on the following facts: (1) the temperature equation is linear, so that the temperature fluctuation is given by the sum of contributions arising from terms of Gu 3 andC d , (2) the equations governing the temperature and the water vapor mixing ratio fluctuations share the condensation fluctuationC d , and (3) the large scale DNS shows that a temperature spectrum that is only excited by a uniform gradient (like Gu 3 ) obeys the Obkhov-Corrsin spectrumk 5 3 in the inertial convective range, followed by the exponential decay in the far diffusive range [26,27,81,83,84].
It is interesting to note that the power spectrum shallowness of the liquid water content (LWC) variance was seen in the stratocumulus cloud observations reported in [23], and mountaintop clouds in [48,69]. The kinetic energy power spectrum is close tof 5 3 and its rapid roll off begins at about h = (¯) f 1 8 8 d (about 50 Hz) in figure 4 of [69]. On the other hand, the LWC power spectrum in figure 12 of [69] only indicatesf 5 3 power law at frequencies below f sb (the frequency of the phase relaxation time). The spectrum slope becomes between −2/3 and −1 at < < f f f 10 sb d (∼40 Hz), after which the LWC spectrum begins to decay faster thanf 1 for  Although the cloud observed at mountain top differs from the cumulus cloud studied in the present study, does not correspond to the LWC power spectrum, and the Reynolds number in the present DNS is one order smaller than those in their experiments, still it is very interesting to compare ( ) E k t , q in the present study with theirs. For example, in figure 9(c) with the measured LWC power spectrum suggests that it may be possible to find common physics for the spectrum of the scalar fluctuations in clouds. We will argue this problem after examining the spectra that are related to the cloud droplet distribution and condensation fluctuations.

Spectra of droplet number density, condensation rate, and liquid water mixing ratio
The isotropic part of the spectra of droplet number density ( ) , and liquid water mixing ratio ( ) , d,  Figure 10 shows the evolution of ( )  The two-point correlation function of the droplet number density with separation r consists of the sum of two terms: x r x x r x r r r r n t n t n t n t n n t w n t n n t g n t , , , where ( ) r w is the correlation function of the droplets and = -( ) ( ) r r g w 1 for > r 0 is the radial distribution function satisfying the condition  ( ) r g 0 as  ¥ r . The ( ) r g term in the third line is related to the probability of finding a droplet within a small volume V d at r conditioned on one droplet within V d at the origin, and the second term is related to the probability that two droplets will reside within the same volume [41]. Therefore, the first term (correlated term) describes the spatial coherency of two droplets at nonzero separation, while the second term arises from the discreteness of the droplets and corresponds to spatially uncorrelated distribution, which is why we refer to it as an uncorrelated term. This was anticipated in [35,44]. In [35] the contributions of the second term of (48) are subtracted by considering the deviation of the particle concentration ( ) ,0, in which case the total number of particles is conserved. In [44] the Fourier transform of the first term of (48) is directly computed by using the particle positions that are most accurate, but requires a large amount of computation. Since the total number of droplets changes in time, and since various spectra related to droplet distribution are to be computed in the present study, we will first examine the spectrum of the second term, and then subtract it from the total spectrum to see the first term.
The spectrum for the uncorrelated part (the second term of the right-hand side of (48)) is However, in our computation, the droplet number density ( ) x n t , , the liquid water mixing ratio ( ) x q t , l and the condensation rate ( ) for each droplet are distributed on the grid points by using the linear weight in the particle in cell (PIC) method, in which the fluid variables are computed and are then treated as the continuum. This means that droplet position spatial accuracy is lost for scales smaller than the grid spacing p D = ( ) x K L 2 max box . Figure 11(d) shows pk L 4 2 box 3 (black) and the ensemble-average of ( ) ( ) E k n , 0 0 n (cyan) computed from the continuum number density after distributing the droplets on the grid points using the PIC method, where ensembles are calculated by changing the random droplet positions. The ratio of the black to cyan curve in figure 11(d) provides a measure of the filter function ( ) F k filter due to the PIC method. The function  ( ) F k filter is shown by a red curve in figure 11(d). Here, it can be seen that ( ) F k filter is close to unity for < kL 70 box and decreases for high wavenumbers > kL 100 box . Figure 11(a) shows the normalized spectrum ( ) ( ) E k t n t , n . Here, it can be seen that the various curves collapse well onto a single curve for > kL 50 box as expected. From similar arguments, we can deduce the amplitudes are in proportion to pk 4 2 for the uncorrelated part spectra as box is excellent, meaning that an uncorrelated spectrum that is proportional to k 2 arises from the discrete property of the droplets. The correlated part spectra ( ) E k t , a c are obtained by subtracting the uncorrelated spectrum as a a a c u c and shown in figure 12.
First, we will consider ( ) E k t , n . When turbulence mixes the droplets, ( ) x n t , fluctuates and deviates from a spatially uniform distribution, which means that the width of the PDF ofñ becomes wider so that increases with time. At later times, the number of droplets decreases due to collision-coalescence so that the amplitude of ( ) E k t , n begins to decrease after 480 s. The curves of ( ) E k t , n c at high wavenumbers show strong oscillations because E n c at high wavenumbers is not well established and is very close to E n uc , which means that the two spectra nearly cancel each other out. The same trend applies also to E C  figure 12, respectively, where it can be seen that they share similar spectral behavior, such as a comparable k 1 power law increase at low wavenumbers, followed by a plateau at < < kL 20 50 box , and a slow decay as a < < a k , 1 2 a a for = ( ) a C q , d l . From the above arguments, the total spectra = ( ) ( ) E k t a n C q , , , , a d l at later times in figure 10 can be interpreted as follows. For low wavenumbers is dominated by the correlated part E a c and the uncorrelated part is small. However, as the wavenumber increases, the uncorrelated part E a uc becomes comparable with E a c for < < kL 20 50 box and larger than it for < < kL 50 120 box , where the spectra behave as k 2 , before finally becoming horizontal for < kL 120 box due to the filtering effect caused by information loss. When the turbulence intensity increases, the amplitude of the correlated part of the spectra becomes larger relative to the uncorrelated part and prevails over a wider wavenumber range. In contrast, the uncorrelated spectra are bounded by A n (t).

Spectra modification
Next, we will consider spectra modifications, beginning with ( ) E k t , q . The equation for ( ) E k t , q is given by represents the scalar variance transfer function due to the convective term and ( ) F k t , q provides the input rate due to the condensation. The input term can be written symbolically as , , is Lagrangian response function of the scalar, and the time integral is taken along with the Lagrangian trajectory of the fluid particle. Since rigorous analysis is far beyond the scope of the present paper, we will proceed with the dimensional arguments. Thus, ( ) , , is the characteristic time of the water vapor mixing ratio at the wavenumber k. Substituting ( ) F k t , q into (55), and then by repeating the same analysis as the passage from (57) to (58), we obtain , .
It is important to note that the total spectrum ( ) For the range < kL 10 box , for example, see figure 9(a), the turbulence spectrum is nearlyk 5 3 , meaning that t µ - , as seen from figure 10(b), we obtain µ -( ) E k t k , q 2 which is steeper thank 5 3 and slower than an exponential decay. The spectral slope of ( ) E k t , q for < < kL 50 120 box at later times of figure 7(c) is close to k 2 . When the turbulence intensity increases, the transition wavenumber k * , which is determined by , shifts to the high wavenumber range. The wavenumber range ofẼ k C d expands and the plateau range moves toward high wavenumbers, so that the steeper slope of ( ) E k t , q for high wavenumbers < < kL 50 120 box appears (for example, compare figure 9(c) with figure 7(c)). The evolution of the temperature spectrum q ( ) E k t , can be explained in the same way, as far as the contributions from the heat exchange through the condensation-evaporation dominate the contributions due to the G( ) t u 3 term. Therefore, the overall functional form of the temperature spectrum is determined by the relative strength of the excitations. When the G( ) t u 3 term is dominant over theC d term, most of q ( ) E k t , would obeyk 5 3 , and only the part at the high end of the inertial convective to the diffusive ranges would be modified. Indeed, even though more careful analysis will be necessary in the future, it is interesting to see that the curves q ( ) E k t , at t=240s (blue curves) in figures 7(b), 8(b) and 9(b) are in such a state, and resemble the LWC power spectrum seen in figure 12 of [69], which is observed at mountaintops. The temperature spectrum peels off the Obukhov-Corrsin spectrumk 5 3 at k c (t) which decreases in time. We infer that k c (t) corresponds to the wavenumber of the scale break argued in the literature [23,36,45] and determined by the balance between the k 5 3 spectrum and the modification. However, since the Reynolds number in the present computation is not high enough and the width of the inertial range is insufficient, we do not argue this problem.
The above results show that the temperature and the water vapor mixing ratio are excited by the condensation-evaporation associated with a large number of small particles that are advected by the turbulent flow. The particles are discrete objects. At large scales, they can be regarded as a continuum, but at small scales their discrete nature emerges with a delta correlated distribution in space. In the spectral space, this discrete property appears as the spectrum of k 2 , the equipartition spectrum. In other words, the discreteness of the particle appears as if it is a kind of thermal noise with particle temperature A a , and the continuum and discrete properties crossover at k * . Therefore, the source term in the scalar equation has spectral support over the entire wavenumber range and consists of coherent and random noise. The latter prevails at small scales and is independent of the turbulent flow, so that the spectra in the far diffusive range take the form ofk 2 , which is universal. At small wavenumbers, the scalar spectrum is modified according to the correlated spectrum of the source term. While it is obvious that we need to develop a theory that will explain the correlated spectrum ( ) E k c , that chore is beyond the scope of the present study. A similar mechanism applies to the slow decay of the kinetic energy spectrum in the far dissipation range through the q l term in the buoyancy force, an example of which can be seen in figure 7(a). Since the turbulence is excited by a large scale external force, the kinetic energy spectrum at low wavenumbers is not affected by the q l term. However, in the far dissipation range, µ -( ) E k t k , 2 is observed when the turbulence intensity is relatively low.
Although the above explanation needs further examination, the fundamental physics for the modification of the turbulence spectra is consistent with DNS observations. This modification of turbulence spectra through the condensation-evaporation process and the liquid water mixing ratio is a peculiar characteristic of cloud turbulence.

Summary and discussion
In this paper, we reported on the development of a cloud microphysics simulator, a DNS model based on the microscopic dynamics that resolves droplet dynamics, turbulence, and the microscale thermodynamic effects embedded within a one-dimensional parcel model. The droplet dynamics include condensation-evaporation, collision-coalescence with the hydrodynamic interaction in terms of the Hall table, but not CCN activation. We also report on the successful integration of a set of equations that covers about 1300 large eddy turnover times, which has never been achieved before, from the viewpoint of the microscopic dynamics.
The results show that the evolution of the parcel altitude, updraft velocity and mean water vapor mixing ratio is insensitive to the inclusion of collisions and the turbulence intensity studied in the present DNSs. On the other hand, the evolution of the size distribution is affected by the collision process. Droplets first grow, primarily through condensation, during which the left tail of the peak of the PDF for droplet radius steepens and the peak moves toward larger radius. Then, the collision-coalescence processes gradually become dominant and the right tail of the mass density distribution function broadens due to autoconversion. The accretion phase finally sets in and droplets grow rapidly due to collision-coalescence. A second peak emerges and grows rapidly in the mass density distribution function and the largest droplet radius reaches about 300 μm after 10 minutes. The abovementioned features of condensation and collision-coalescence growth are well known from previous bin-model simulations [7,29,52,80]. However, to the best of our knowledge, this is the first study to reproduce these features via a DNS of turbulence over long time by the spectral method that includes the microscopic dynamics of cloud droplets such as Reynolds number dependent drag, processes of evaporation-condensation, collisioncoalescence, and hydrodynamic interaction.
The cloud microphysics simulator also reveals notable turbulence-droplet interaction features. We found that when the turbulence intensity becomes stronger, collision-coalescence is enhanced and the broadening of the droplet spectrum becomes faster, as seen in the previous studies. Furthermore, we found a new feature in the variance spectra of the temperature and the water vapor mixing ratio. Specifically, we found that, in contrast to the traditional Kolmogorov-Obukhov-Corrsin theory, the spectra are shallower than thek 5 3 spectrum at low wavenumbers and attain a plateau followed by the algebraic decay in the diffusive range. The physical mechanism for this modification was studied by examining the droplet related spectra, such as the droplet number density, the condensation rate, and the liquid water mixing ratio. We found that those spectra consist of correlated and uncorrelated parts, and that the contributions arising from the condensation process are responsible for the turbulence spectra modifications. The condensation spectrum is determined both by turbulent mixing at large scales and by a delta correlated distribution at small scales due to the discrete nature. The same mechanism applies to the spectrum of the liquid water mixing ratio in the buoyancy force. In this way, the droplets affect the dynamics of turbulence and two scalars with the turbulence relaxation time, which is a function of the wavenumber. The direct interaction between turbulence and droplets through the condensation process and buoyancy force are, to the best of the authors' knowledge, new findings obtained via the present DNS study.
In the direction of studying the whole process of cloud droplets to raindrops by synthesizing the key cloud microphysical processes, the present results of long-time simulations from the microscopic viewpoint are promising, and can be expected to provide motivation for conducting larger size simulations in order to examine the effects of higher Reynolds numbers. It is also possible and vital to extend the space and time domain of the microscopic simulation to larger and longer scales in order to connect to the latest advances in (semi) macroscopic modeling and computing.
However, at the same time, we need to anticipate the limitations of the present simulation. The cloud droplets are treated as the point particles, so that the time integration can not proceed beyond the time at which the droplet radius becomes comparable to the grid size. To extend the computation longer, the simulation code needs to resolve the cloud droplet on the grid points, which would consume a considerable computer resources. The entrainment effects are introduced as single macro scale parameter μ in the present study, and plays the important role to determine the mean evolution of the parcel, and affects indirectly the turbulent flow. However, the entrainment is one of the major mechanisms to directly drive small scale turbulence. In the present computation, the turbulence is maintained in a steady state by the random forcing, which is artificial. One method to drive the turbulence in a more realistic way is to introduce a forcing mechanism related to the shearing motion acting on the parcel at scales larger than it, implying that the turbulence forcing is to be introduced in conjunction with the entrainment. More reasonable methods to achieve the entrainment effects are to be developed. The hydrodynamic interaction is partially included by using Hall's table, but a more comprehensive table that includes the effects of high droplet Reynolds numbers and surface tension will need to be introduced. The CCN process should also be included in future studies. The preliminary results for the CCN effects on the evolution of the droplet spectrum by the microscopic simulation are quite satisfactory and will be reported in future works.
On the other hand, a new class of fundamental problems regarding turbulence, such as how the correlated spectra are related to the inertial particles, how we theoretically compute, to what extent they are universal, and what modification of the turbulence and scalars by the droplets arises, have also arisen. Studies of the above fundamental problems by using seamless simulations, using the above essential key processes at higher Reynolds numbers, are certain to be challenging targets for future research into cloud turbulence.
To validate the collision detection scheme used in the present model, we follow a way similar to that described in section 3 in [21], and compare two theoretical collision rates with the results of numerical simulations. In this validation, droplets do not coalesce after collision and are allowed to move independently of other droplets [79].
The first theoretical collision rate used for validation is that for droplets falling only under the force of gravity in air without turbulence. For bi-dispersed droplets with number density n i , radius R i and terminal velocity V i = ( ) i 1, 2 , the theoretical collision rate is given by Assuming that the drag force is given by Stokes drag, the terminal velocities for droplets with = . The radius of half of the droplets is set to 10 μm while that of the other half is set to 20 μm. The time step is D = t 1 ms, and the fluid is kept stationary during time-integration. After droplets reach their terminal velocities, the time average of the collision rate over 10 4 time steps is calculated. Twenty ensemble runs are conducted using different random seeds for the initial droplet positions. The ensemble average of the collision rate is 0.45425 cm −3 s −1 , which is different from the theoretical prediction by 0.1%. Therefore, the simulation results give a sufficiently accurate value for the collision rate.
The second theoretical collision rate used for validation is that for monodisperse droplets without inertia advected by homogeneous isotropic turbulence. In this case, the simulation results include numerical errors associated with interpolation of the velocity field to the position of the droplet. For monodisperse droplets with the number density n and radius R, the theoretical collision rate is given by [53,72]: where ò is the mean energy dissipation rate. Numerical simulations are conducted using the grid number N=128 for each direction, grid spacing D = x 1 mm, and the number of droplets = N 2 p 18 (corresponding to the number density of 125 cm −3 ). Droplets without inertia and condensation growth are advected by the turbulent flow that is maintained by the same forcing as discussed in (39). Droplets are monodisperse with R = 10 μm, the time step is set to D = t 1 ms, and = c 80 f cm 2 s −3 . After a statistically steady state is attained, the time averages of the collision rate and energy dissipation rate ò over 10 6 time steps are calculated. The average of ò is 101.3 cm 2 s −3 , and the collision rate is´-2.03 10 3 cm −3 s −1 . When compared to the theoretical value of =´-N 2.10 10 3 cm −3 s −1 given by (A.2), the result of the numerical simulation is different by 3%. The deviation from the theory is greater than that for droplets only under gravity, presumably because the simulation results include errors associated with interpolation of the velocity field. However, the deviation is small and the accuracy is sufficient for the purposes of the present study.
. It is important to examine the meaning of the ensemble. When turbulence exists, the distribution function for the droplet configuration (position and radius) and the turbulence fields at the droplet position are given by considering the joint ensemble of realizations of turbulence and droplets. Thus, the joint distribution function is written as , , , , , , , , which is the probability of finding droplets within f