Continuous growth of cloud droplets in cumulus cloud

Turbulence and cloud are key to understanding atmospheric dynamics, forecasting weather, and predicting the future climate. Because the minimum grid size in numerical weather (climate) predictions is of the order of a few kilometers under optimal conditions when using high performance computers, the degrees of freedom of turbulent flow and cloud microphysical processes, such as cloud condensation nucleation, condensation/evaporation, and collision/coalescence/breakup of water droplets below the grid size, are adequately modeled or parametrized with certain hypotheses. For example, both the bin and super droplet methods use the sub-grid model for the turbulent motion below the grid scale and the prescribed functional form for the collision kernel [1–6]. Although many efforts have been made to understand the dynamics of turbulence and cloud formation, and to develop more accurate models [1–9], our knowledge of these degrees of freedom is limited because turbulence is still not fully understood [10, 11], and the direct observation and measurement of cloud microphysical processes is difficult [12–15].

Then it is quite natural and reasonable to study the microphysical processes and turbulent transport from more fundamental view point according to first principle, which is to numerically compute all the degrees of motion as accurately with respect to the fundamental physics as possible. This approach requires both Eulerian and Lagrangian computations for the turbulence and the cloud droplets, respectively, and needs large computational resources. Recent progress in high performance computing has driven this trend. Vaillancourt, Yau, and Grabowski [16, 17] were the first to undertake this research for studying the broadening of the droplet spectrum (size distribution). They introduced a small parcel, which had a cubic domain of (15 cm)³ inside the quasi-adiabatic cloud core and derived two sets of equations for calculating the variables of the cloud and turbulent flow. A macroscopic set of equations was used to describe the evolution of the various means within the parcel, and a microscopic set of equations was used to describe the dynamics of each droplet and the turbulent fluctuations of velocity, temperature, and the water vapor mixing ratio. It was found that the increase in the variance of the droplet radius was too small to explain the observed broadening of the droplet spectrum by the sedimentation, randomness of the particle position and/or the initial variation of the droplet radius and turbulence.

Then droplet collision is considered to be a key process for the broadening of the droplet spectrum. The collision process of particles in turbulent flow, and its statistical properties, have been studied extensively since Saffman and Turner [18]. Many studies have analytically and numerically examined particle density [19, 20], collision frequency [21–24], and various influences on collision efficiency such as sedimentation [25–27], the clustering and spatial distribution of the inertial particles [28–36], and the hydrodynamic interaction [37–39]. It was found that the collision of cloud droplets with a similar size depends on the Stokes number ${St}={\tau }_{p}/{\tau }_{K}$, where ${\tau }_{p}$ and ${\tau }_{K}$ are the characteristic time of the droplet and the Kolmogorov time of turbulent flow, respectively, and becomes frequent for ${St}=O(1)$, at which the droplet radius is about 30 μm. Sedimentation due to gravity and the clustering of droplets with ${St}=O(1)$ tends to enhance collisions.

However, these findings raise further questions, for example, as to how fast the cloud droplets grow from the water drops around the condensation nuclei to about 30 μm beyond which the collision/coalescence becomes important, how the size distribution of the aerosols affects the evolution of the droplet spectrum, and the condensation process in the supersaturated air dominates this early stage of growth but how can the supersaturated state be maintained? For the last question, a plausible way is to consider the adiabatic cooling by the updraft inside a cloud. However, at the same time the updraft engulfs the unsaturated air in the cloud and mixes it with saturated air and cloud droplets, the entrainment, which tends to slow down the droplet growth. Therefore, the updraft effects are not onefold but bring other factors into the droplet growth [40–50].

There are many accumulated knowledges for individual processes, but the whole is not necessarily a simple sum of elementary processes. What is necessary in the research for cloud evolution is a study that explores the continuous growth of cloud droplets from sub-micron to hundreds micron meter radius, a whole story of the rain drop formation. The goal of our research is to reconstruct the whole process, to elucidate the real feature of the phenomena and to extract the useful knowledge for better numerical modeling.

From the above arguments it is reasonable to infer a scenario about the rain drop formation that cloud droplets grow by condensation without an appreciable broadening of the droplet spectrum up to about 30 μm radius and then evolve into rain drops by the collision-coalescence. If this scenario is applied to the core region of cumulus cloud where the updraft exists, it suggests that the domain of the droplet growth dominated by the collision-coalescence is not necessarily the same as that dominated by the condensation. Therefore the computation of the cloud microphysics within a small domain fixed inside cloud as used in the previous studies is not appropriate for seamless computation of the continuous growth of the cloud droplets to the rain drops.

Our purpose in this paper is to present a new method that is able to seamlessly simulate the continuous growth of cloud droplets from about 10 μm to a few hundred μm with a minimum set of processes which are essential for the rain drop formation. And we examine its consistency from the view points of the numerical and physical properties. In order to study the evolution of the system in more natural and closer circumstances to the realistic conditions, we consider a scenario that a parcel which contains a large number of droplets and turbulence rises inside a typical maritime cumulative cloud as studied by Takahashi [51]. In the time domain, the evolution process to be studied is about 20 min. However, it should be stressed that in this paper we do not intend to accurately simulate all the processes in the cloud microphysics. Instead, we focus on to show effectiveness and wider capability of the new method for incorporating various detailed dynamics in the long time integration. As a first step, we studied the continuous growth, without collisions, of cloud droplets by condensation within a small parcel advected by the updraft. The paper is organized as follows. Section 2 describes the mathematical formulation of the problem, the equations of the fluid, temperature, the water vapor mixing ratio, and the droplets. Section 3 describes the direct numerical simulation (DNS) and the various numerical parameters. The results are presented in section 4, and the summary and discussions are given in section 5.

Here, we consider a small cubic box of length ${L}_{{\rm{box}}}$ which is moving with an updraft velocity $\bar{W}(t)$ inside a cumulus cloud above the tropical sea (see figure 1) [51]. The box contains many cloud droplets that undergo advection by turbulent flow and growth or decay due to condensation, evaporation or collision processes. A coordinate system ${\boldsymbol{x}}=(x,y,z)$ whose origin is fixed on the sea surface is introduced in such a way that the z axis is set along the vertical upward direction and the x and y axes are in the horizontal directions. The environmental conditions for the cumulus cloud in this study are specified by the environmental temperature T_e(z) and the environmental water vapor mixing ratio Q_ve(z), respectively. We assume that ${L}_{{\rm{box}}}$ is much smaller than the typical length ${L}_{{\rm{cloud}}}$ of the cumulus cloud. The turbulent velocity ${\boldsymbol{u}}$, temperature T, and water vapor mixing ratio Q_v are well mixed on a scale larger than ${L}_{{\rm{box}}}$, but their fluctuations appear to be observable and statistically homogeneous (but axially symmetric) when seen on a scale below ${L}_{{\rm{box}}}$ [42–47]. Therefore, it is reasonable to impose periodic boundary conditions in three directions for the three fields. We also assume that the fluid is incompressible and the buoyancy force is expressed in terms of the Boussinesq approximation. A local coordinate system ${{\boldsymbol{x}}}^{\prime }=({x}^{\prime },{y}^{\prime },{z}^{\prime })$ is introduced in such a way that the system is attached to the cubic box and moves with the updraft $\bar{W}(t)$. Two coordinate systems are related to each other by the time-dependent Galilean transform:

$$\begin{eqnarray}&&{\boldsymbol{x}}={{\boldsymbol{x}}}^{\prime }+{\int }_{0}^{t}\bar{W}({t}^{\prime }){\rm{d}}{t}^{\prime }\ {{\boldsymbol{e}}}_{z}={{\boldsymbol{x}}}^{\prime }+\bar{H}(t){{\boldsymbol{e}}}_{z},\qquad {\boldsymbol{u}}={{\boldsymbol{u}}}^{\prime }+\bar{W}(t){{\boldsymbol{e}}}_{z},\end{eqnarray} \tag{ 1 }$$

where $\bar{H}(t)$ and ${{\boldsymbol{e}}}_{z}$ are the altitude of the box and the unit vector in the vertical direction, respectively. Then the equation of motion for the velocity field with respect to this coordinate is:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}^{\prime} }{\partial t}+{\boldsymbol{u}}^{\prime} \cdot {\rm{\nabla }}^{\prime} {\boldsymbol{u}}^{\prime} =-\displaystyle \frac{1}{{\rho }_{a}}{\rm{\nabla }}^{\prime} p^{\prime} +{\nu }_{a}{\rm{\nabla }}{^{\prime} }^{2}{\boldsymbol{u}}^{\prime} +\left(B-\displaystyle \frac{{\rm{d}}\bar{W}}{{\rm{d}}t}\right){{\boldsymbol{e}}}_{z},\quad \end{eqnarray} \tag{ 2 }$$

where ${p}^{\prime }$ is the pressure, ${\rho }_{a}$ and ${\nu }_{a}$ are the density and kinematic viscosity of dry air, respectively, and B is the buoyancy force which is given below. The updraft velocity $\bar{W}$ is determined by the requirement that the volume average of the fluid velocity over the cubic box ${V}_{{\rm{box}}}$ vanishes. Taking the volume average of the Navier–Stokes equation over ${V}_{{\rm{box}}}$ and using the periodic boundary condition and the assumption $\bar{{\boldsymbol{u}}^{\prime} }(0)=0$, we obtain:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\bar{{\boldsymbol{u}}^{\prime} }}{{\rm{d}}t}=\left(\bar{B}-\displaystyle \frac{{\rm{d}}\bar{W}}{{\rm{d}}t}\right){{\boldsymbol{e}}}_{z}=0\ \ \longrightarrow \ \ \bar{W}(t)=\bar{W}(0)+{\int }_{0}^{t}\bar{B}({t}^{\prime })\ {\rm{d}}{t}^{\prime },\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&\bar{H}(t)=\bar{H}(0)+{\int }_{0}^{t}\bar{W}({t}^{\prime })\ {\rm{d}}{t}^{\prime },\end{eqnarray} \tag{ 4 }$$

where the volume average is defined by:

$$\begin{eqnarray}&&\bar{A}\equiv \displaystyle \frac{1}{{V}_{{\rm{box}}}}{\int }_{{V}_{{\rm{box}}}}A({\boldsymbol{x}}^{\prime} )\ {\rm{d}}{\boldsymbol{x}}^{\prime} .\end{eqnarray} \tag{ 5 }$$

With respect to this moving frame, the mean fluid velocity is zero so that the homogeneity of turbulent fields holds, but it should be noted that the mean velocity of the cloud droplets inside the box is not necessarily zero. The environmental temperature and water vapor mixing ratio at the altitude of the parcel are given by:

$$\begin{eqnarray}&&{T}_{{\rm{e}}}(\bar{H}(t))={T}_{{\rm{e}}0}-{\int }_{{H}_{0}}^{\bar{H}(t)}{{\rm{\Gamma }}}_{{\rm{e}}}(z){\rm{d}}z,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{Q}_{\mathrm{ve}}(\bar{H}(t))={Q}_{\mathrm{ve}0}-{\int }_{{H}_{0}}^{\bar{H}(t)}{{\rm{\Gamma }}}_{\mathrm{ve}}(z){\rm{d}}z,\end{eqnarray} \tag{ 7 }$$

respectively. Note that the integration is made with respect to the z coordinate fixed to the sea surface. Here T_e0 and Q_ve0 are the temperature and the water vapor mixing ratio of the environment at ${H}_{0}\equiv \bar{H}(t=0)$, respectively, and ${{\rm{\Gamma }}}_{{\rm{e}}}(z)$ and ${{\rm{\Gamma }}}_{\mathrm{ve}}(z)$ are the lapse rate of the air and the rate of change of the water vapor mixing ratio $| {\rm{d}}{Q}_{\mathrm{ve}}/{\rm{d}}z|$, respectively.

Zoom In Zoom Out Reset image size

The incompressibility condition, the equations of the fluctuating turbulent velocity, the temperature T, and the water vapor mixing ratio Q_v are given by:

$$\begin{eqnarray}&&{\rm{\nabla }}^{\prime} \cdot {\boldsymbol{u}}^{\prime} =0,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}^{\prime} }{\partial t}+{\boldsymbol{u}}^{\prime} \cdot {\rm{\nabla }}^{\prime} {\boldsymbol{u}}^{\prime} =-\displaystyle \frac{1}{{\rho }_{a}}{\rm{\nabla }}^{\prime} p^{\prime} +{\nu }_{a}{\rm{\nabla }}{^{\prime} }^{2}{\boldsymbol{u}}^{\prime} +(B-\bar{B}){{\boldsymbol{e}}}_{z}+{{\boldsymbol{f}}}^{\prime },\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}+{\boldsymbol{u}}^{\prime} \cdot {\rm{\nabla }}^{\prime} T={\kappa }_{T}{\rm{\nabla }}{^{\prime} }^{2}T+\displaystyle \frac{{L}_{{\rm{v}}}}{{c}_{p}}{C}_{{\rm{d}}}-{\rm{\Gamma }}(t)\bar{W}(t),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {Q}_{{\rm{v}}}}{\partial t}+{\boldsymbol{u}}^{\prime} \cdot {\rm{\nabla }}^{\prime} {Q}_{{\rm{v}}}={\kappa }_{{\rm{v}}}{\rm{\nabla }}{^{\prime} }^{2}{Q}_{{\rm{v}}}-{C}_{{\rm{d}}},\end{eqnarray} \tag{ 11 }$$

respectively, where ${{\boldsymbol{f}}}^{\prime }$ is the external force that is introduced to maintain the turbulence in a steady state, the ${\kappa }_{T}$ and ${\kappa }_{{\rm{v}}}$ are the molecular diffusivity coefficients of the temperature and water vapor mixing ratio, respectively, L_v is the latent heat of vaporization, and c_p is the specific heat of the air [1, 10]. The C_d is the condensation rate, which is discussed below, and the last term of the right-hand side of equation (10) denotes the adiabatic cooling due to the updraft and ${\rm{\Gamma }}(t)={\rm{\Gamma }}(\bar{H}(t))$ is the rate of change of the mean temperature inside the parcel at $\bar{H}(t)$ which is described below. Because the box size ${L}_{{\rm{box}}}$ is much smaller than the cloud size, the contributions from the vertical fluctuating fluid velocity ${\rm{\Gamma }}{u}_{3}^{\prime }$ to the adiabatic cooling is neglected. The ${\rm{\Gamma }}$ is assumed to be given by the sum of two contributions as

$$\begin{eqnarray}&&{\rm{\Gamma }}={{\rm{\Gamma }}}_{s}+{{\rm{\Gamma }}}_{{\rm{ent}}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}{{\rm{\Gamma }}}_{s}(t)=\left\{\begin{array}{ll}{{\rm{\Gamma }}}_{{\rm{d}}}\left(1+\displaystyle \frac{{L}_{{\rm{v}}}{\bar{Q}}_{{\rm{v}}}(t)}{{R}_{{\rm{d}}}\bar{T}(t)}\right){\left(1+\displaystyle \frac{\varepsilon {L}_{{\rm{v}}}^{2}{\bar{Q}}_{{\rm{v}}}(t)}{{c}_{p}{R}_{{\rm{d}}}{\bar{T}}^{2}(t)}\right)}^{-1} & (S\gt 0),\\ {{\rm{\Gamma }}}_{{\rm{d}}}{(1+0.846{\bar{Q}}_{{\rm{v}}}(t))}^{-1} & (S\lt 0),\end{array}\right.\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{ent}}}(t)=\mu \left(\bar{T}(t)-{T}_{{\rm{e}}}(\bar{H}(t))+\displaystyle \frac{{L}_{{\rm{v}}}}{{c}_{p}}[{\bar{Q}}_{{\rm{v}}}(t)-{Q}_{\mathrm{ve}}(\bar{H}(t))]\right){\left(1+\displaystyle \frac{{L}_{{\rm{v}}}^{2}{\bar{Q}}_{{\rm{v}}}(t)}{{c}_{p}{R}_{{\rm{v}}}{\bar{T}}^{2}(t)}\right)}^{-1},\end{eqnarray} \tag{ 14 }$$

where ${{\rm{\Gamma }}}_{s}$ is the effect of the moist air and ${{\rm{\Gamma }}}_{{\rm{d}}}$ is the dry adiabatic lapse rate. When the parcel is ascending, the moist air is very close to being supersaturated ($S\gt 0$) and when it is descending the parcel of air is undersaturated ($S\lt 0$). The ${{\rm{\Gamma }}}_{{\rm{ent}}}$ is the effect due to the entrainment [1, 2] and μ is the entrainment strength parameter.

For the cloud droplets, we followed the formalism developed by Vaillancourt, Yau and Grabowski [16]. The microscopic dynamics of the cloud droplets are assumed to be specified by the position ${{\boldsymbol{X}}}_{j}^{\prime }(t)$, the velocity ${{\boldsymbol{V}}}_{j}^{\prime }(t)$, and the radius R_j(t) of the jth droplet, with respect to the comoving frame. Then, the buoyancy force $B({{\boldsymbol{x}}}^{\prime },t)$, the condensation rate ${C}_{{\rm{d}}}({{\boldsymbol{x}}}^{\prime },t)$, and the liquid water mixing ratio ${q}_{{\rm{l}}}({{\boldsymbol{x}}}^{\prime },t)$, are given by

$$\begin{eqnarray}&&B({{\boldsymbol{x}}}^{\prime },t)=g\left(\displaystyle \frac{T({\boldsymbol{x}}^{\prime} ,t)-{T}_{{\rm{e}}}(\bar{H}(t))}{{T}_{{\rm{e}}}(\bar{H}(t))}+\varepsilon ({Q}_{{\rm{v}}}({\boldsymbol{x}}^{\prime} ,t)-{Q}_{\mathrm{ve}}(\bar{H}(t)))-{q}_{{\rm{l}}}({\boldsymbol{x}}^{\prime} ,t)\right),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{C}_{{\rm{d}}}({\boldsymbol{x}}^{\prime} ,t)=\displaystyle \frac{4\pi {\rho }_{{\rm{l}}}{K}_{s}(t)}{{\rho }_{a}{({\rm{\Delta }}x)}^{3}}\displaystyle \sum _{j=1}^{{N}_{{\rm{\Delta }}}({\boldsymbol{x}}^{\prime} ,t)}{R}_{j}(t)S({\boldsymbol{X}}{^{\prime} }_{j}+\bar{H}(t){{\boldsymbol{e}}}_{z},t),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{q}_{{\rm{l}}}({\boldsymbol{x}}^{\prime} ,t)=\displaystyle \frac{4\pi {\rho }_{{\rm{l}}}}{3{\rho }_{a}{\rm{\Delta }}{x}^{3}}\displaystyle \sum _{j=1}^{{N}_{{\rm{\Delta }}}({\boldsymbol{x}}^{\prime} ,t)}{R}_{j}^{3}(t),\end{eqnarray} \tag{ 17 }$$

where g is the gravitational acceleration, and $\varepsilon ={R}_{{\rm{d}}}/{R}_{{\rm{v}}}=0.61$, where R_d and R_v are the gas constants for dry air and water vapor, respectively. ${N}_{{\rm{\Delta }}}({\boldsymbol{x}}^{\prime} ,t)$ is the number of droplets within the grid cell ${({\rm{\Delta }}x)}^{3}$, ${\rho }_{{\rm{l}}}$ is the liquid density, respectively, S is the supersaturation, and K_s is a temperature-dependent diffusion coefficient specified below.

The cloud droplets are assumed to move around inside the box under the turbulence, the gravity force and the Stokes drag. The use of the Stokes drag for the cloud droplet is justified for droplets with a radius smaller than 30 μm in this study because the droplet Reynolds number is smaller than unity. Then, the equations of motion of the cloud droplets with respect to the comoving frame are:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\boldsymbol{X}}{^{\prime} }_{j}(t)}{{\rm{d}}t}={\boldsymbol{V}}{^{\prime} }_{j}(t),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\boldsymbol{V}}{^{\prime} }_{j}(t)}{{\rm{d}}t}=\displaystyle \frac{1}{{\tau }_{j}(t)}({\boldsymbol{u}}^{\prime} ({\boldsymbol{X}}{^{\prime} }_{j}(t),t)-{\boldsymbol{V}}{^{\prime} }_{j}(t))-\left(g+\displaystyle \frac{{\rm{d}}\bar{W}(t)}{{\rm{d}}t}\right){{\boldsymbol{e}}}_{z},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{R}_{j}(t)\displaystyle \frac{{\rm{d}}{R}_{j}(t)}{{\rm{d}}t}={K}_{s}(t)S({{\boldsymbol{X}}}_{j}(t),t),\qquad {{\boldsymbol{X}}}_{j}(t)={\boldsymbol{X}}{^{\prime} }_{j}(t)+\bar{H}(t){{\boldsymbol{e}}}_{z},\end{eqnarray} \tag{ 20 }$$

where ${\tau }_{j}(t)$ and $S({{\boldsymbol{X}}}_{j}(t),t)$ are the relaxation time and the supersaturation at the particle position defined by

$$\begin{eqnarray}&&{\tau }_{j}(t)=\displaystyle \frac{2}{9}\displaystyle \frac{{\rho }_{{\rm{l}}}}{{\rho }_{a}}\displaystyle \frac{{R}_{j}^{2}(t)}{{\nu }_{a}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&S({{\boldsymbol{X}}}_{j}(t),t)=\displaystyle \frac{{Q}_{{\rm{v}}}({{\boldsymbol{X}}}_{j}(t),t)}{{Q}_{\mathrm{vs}}({{\boldsymbol{X}}}_{j}(t))}-1,\end{eqnarray} \tag{ 22 }$$

respectively, and:

$$\begin{eqnarray}&&{K}_{s}(t)={\kappa }_{{\rm{v}}}\displaystyle \frac{{\rho }_{\mathrm{vs}}(\bar{T}(t))}{{\rho }_{{\rm{l}}}}{\left[1+\displaystyle \frac{{\kappa }_{{\rm{v}}}{L}_{{\rm{v}}}{\rho }_{\mathrm{vs}}(\bar{T}(t))}{{\kappa }_{T}\bar{T}(t)}\left(\displaystyle \frac{{L}_{{\rm{v}}}}{{R}_{{\rm{v}}}\bar{T}(t)}-1\right)\right]}^{-1},\end{eqnarray} \tag{ 23 }$$

where ${\rho }_{\mathrm{vs}}(\bar{T}(t))$ is the saturated vapor density. We assume that the supersaturated water vapor mixing ratio Q_vs at the droplet position is determined by Tetens' formula:

$$\begin{eqnarray}&&{Q}_{\mathrm{vs}}({{\boldsymbol{X}}}_{j}(t))=\varepsilon \displaystyle \frac{610.78}{P(\bar{T}(t),\bar{H}(t))}\mathrm{exp}\left(17.269\displaystyle \frac{T({{\boldsymbol{X}}}_{j}(t))-273.16}{T({{\boldsymbol{X}}}_{j}(t))-35.86}\right).\end{eqnarray} \tag{ 24 }$$

The pressure and the altitude are related by:

$$\begin{eqnarray}&&P(\bar{T}(t),\bar{H}(t))={P}_{{\rm{sea}}}{\left(1-\displaystyle \frac{0.0065\bar{H}(t)}{\bar{T}(t)+0.0065\bar{H}(t)}\right)}^{5.257},\end{eqnarray} \tag{ 25 }$$

where ${P}_{{\rm{sea}}}$ is the pressure at the sea level and set to be 1013.25 (hPa).

The set of equations (8)–(25) is numerically integrated over about 20 min. For cloud droplets with a radius of the order of 10 μm the representative time increment ${\rm{\Delta }}t$ is about 1 ms, and therefore the total number of time steps for one run is 1.2 million, which requires long time integration and large computational resources. Before proceeding to heavy computations it is important to examine the consistency of the above microscopic equations with the macroscopic behavior that has been found by many previous studies. Also knowing in advance the average behavior of the above microscopic equations is helpful and useful to choose proper parameters for DNSs of the full set of equations such as the environmental temperature, entrainment strength, and so on. For this purpose, we studied a parcel model which is obtained by taking the volume average of the equations derived in the previous section and described in detail in the appendix. It was found that the time evolution of the altitude and updraft velocity of the parcel and of the mean droplet radius was consistent with the results by the bin method [51] and with those in the literature [57–60]. The results of the above preliminary study were encouraging in that the set of the microscopic equations was expected to adequately simulate the detailed evolution of the cloud droplets and turbulent field inside the parcel as well as the mean behavior.

The equations under consideration are equations (8)–(17) for the turbulent air flow, equations (18)–(25) for the droplets. The numerical methods used were as follows. For the spatial derivatives, we used a hybrid code, which employs the spectral method for the velocity field and the combined compact finite difference method (CCFD) for the temperature and the water vapor mixing ratio. We developed an efficient three-dimensional fast Fourier transform (3DFFT) using the two-dimensional domain decomposition, and we also developed an efficient solver for the 8th order CCFD using a three-dimensional domain decomposition. The numerical accuracy and performance of the hybrid code have been carefully examined in [52] and were found to be very satisfactory. The velocity and scalar fields at the droplet position were linearly interpolated, and the contributions from the condensation term and the buoyancy force were distributed on the surrounding grid points using the same weighting as that of the linear interpolation, i.e., the particle in the cell method (PIC). The time integration was made in terms of the predictor corrector method.

The environmental temperature and water vapor mixing ratio in the air are taken from observation data at Hawaii [51] and approximately given by

$$\begin{eqnarray}{T}_{{\rm{e}}}(P)={T}_{0}+{\int }_{{P}_{0}}^{P}{{\rm{\Gamma }}}_{{\rm{e}}}({p}^{\prime }){\rm{d}}{p}^{\prime },\quad {{\rm{\Gamma }}}_{{\rm{e}}}=\left\{\begin{array}{ll}0.0588\ {\rm{K}}\ {(\mathrm{hPa})}^{-1} & (P\gt 780\ {\rm{hPa}})\\ -0.0333\ {\rm{K}}\ {(\mathrm{hPa})}^{-1} & (P\lt 780\ {\rm{hPa}}),\end{array}\right.\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}{Q}_{\mathrm{ve}}(P)={Q}_{\mathrm{ve}0}+{\int }_{{P}_{0}}^{P}{{\rm{\Gamma }}}_{\mathrm{ve}}({p}^{\prime }){\rm{d}}{p}^{\prime },\quad {{\rm{\Gamma }}}_{\mathrm{ve}}=\left\{\begin{array}{ll}3.92\times {10}^{-5}\text{}{(\mathrm{hPa})}^{-1} & (P\gt 780\ {\rm{hPa}})\\ 1.33\times {10}^{-4}\text{}{(\mathrm{hPa})}^{-1} & (P\lt 780\ {\rm{hPa}}),\end{array}\right.\end{eqnarray} \tag{ 27 }$$

where ${{\rm{\Gamma }}}_{\mathrm{ve}}(P)=| {\rm{d}}{Q}_{\mathrm{ve}}(P)/{\rm{d}}P|$ and the argument z of ${{\rm{\Gamma }}}_{{\rm{e}}}(z)$ of equation (6) and ${{\rm{\Gamma }}}_{\mathrm{ve}}(z)$ of equation (7) is now switched to the pressure P as ${{\rm{\Gamma }}}_{{\rm{e}}}(P)$ and ${{\rm{\Gamma }}}_{\mathrm{ve}}(P)$ using equation (25) as is often used. Note that $P\lt {P}_{0}$ and P₀ is the pressure at $z=H(0)$, the initial position of the box. The distributions of T_e(P) (red, color online) and Q_ve(P) (blue) are plotted in figure 2 as functions of the pressure (left axis) as well as functions of the altitude (right axis). It is seen that the inversion occurs at 780 hPa (about 2200 m above the sea surface).

Zoom In Zoom Out Reset image size

A cubic box with size ${L}_{{\rm{box}}}=25.6\;{\rm{cm}}$ is initially at rest $\bar{W}(0)=0$ ms⁻¹ at $\bar{H}(0)=500\;{\rm{m}}$, at which the mean temperature and mean water vapor mixing ratio in the environment are:

$$\begin{eqnarray}&&{T}_{{\rm{e}}0}=293\ {\rm{K}},\qquad {Q}_{\mathrm{ve}}=15.6\ {\rm{g}}\ {(\mathrm{Kg})}^{-1},\end{eqnarray} \tag{ 28 }$$

respectively. Inside the box, the mean temperature is set to be slightly higher than the environmental temperature as:

$$\begin{eqnarray}&&{T}_{0}=293.25\ {\rm{K}}.\end{eqnarray} \tag{ 29 }$$

The mean water vapor mixing ratio is assumed to be saturated and is determined by Tetens' formula for the saturated vapor pressure as:

$$\begin{eqnarray}&&{e}_{{\rm{sat}}}=610.78\times \mathrm{exp}\left(17.269\displaystyle \frac{{T}_{0}-273.16}{{T}_{0}-35.86}\right),\end{eqnarray} \tag{ 30 }$$

and then the mean water vapor mixing ratio is:

$$\begin{eqnarray}&&{Q}_{{\rm{v}}0}=\varepsilon \displaystyle \frac{{e}_{{\rm{sat}}}}{{P}_{0}-{e}_{{\rm{sat}}}}\approx \varepsilon \displaystyle \frac{{e}_{{\rm{sat}}}}{{P}_{0}}.\end{eqnarray} \tag{ 31 }$$

The K_s, which is defined by equation (23) and used in equation (16), was treated as a constant determined by T₀ and P₀ for most runs because the temperature dependence was found to be very weak. The entrainment strength μ in (14) is chosen to be $7.5\times {10}^{-5}$ cm⁻¹.

The fluctuating turbulent velocity ${{\boldsymbol{u}}}^{\prime }$ and the fluctuations of temperature θ and water vapor mixing ratio q_v are defined as follows:

$$\begin{eqnarray}&&{{\boldsymbol{u}}}^{\prime }={{\boldsymbol{u}}}^{\prime }({{\boldsymbol{x}}}^{\prime },t),\qquad T={T}_{0}(t)+\theta ({{\boldsymbol{x}}}^{\prime },t),\qquad {Q}_{{\rm{v}}}={Q}_{{\rm{v}}0}(t)+{q}_{{\rm{v}}}({{\boldsymbol{x}}}^{\prime },t).\end{eqnarray} \tag{ 32 }$$

The fluctuations are assumed to obey multivariate Gaussian statistics with zero means. The kinetic energy spectrum and the variance spectra of the temperature and water vapor mixing ratio are defined by:

$$\begin{eqnarray}&&E(t)=\displaystyle \frac{\langle {{{\boldsymbol{u}}}^{\prime }}^{2}\rangle }{2}=\displaystyle \frac{3}{2}{u}_{0}^{2}={\int }_{0}^{\infty }E(k,t){\rm{d}}k,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{E}_{\theta }(t)=\langle {\theta }^{2}\rangle ={\theta }_{0}^{2}={\int }_{0}^{\infty }{E}_{\theta }(k,t){\rm{d}}k,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{E}_{q}(t)=\langle {q}_{{\rm{v}}}^{2}\rangle ={q}_{0}^{2}={\int }_{0}^{\infty }{E}_{q}(k,t){\rm{d}}k,\end{eqnarray} \tag{ 35 }$$

respectively, where $\langle \ \rangle$ is the ensemble average and ${u}_{0},{\theta }_{0}$ and q₀ are the root mean square velocity, temperature, and water vapor mixing ratio, respectively. The initial kinetic energy spectrum is:

$$\begin{eqnarray}&&E(k,0)=16\sqrt{\displaystyle \frac{2}{\pi }}{u}_{0}^{2}{k}_{0}^{5}\mathrm{exp}\left(-\displaystyle \frac{{k}^{2}}{2{k}_{0}^{2}}\right),\end{eqnarray} \tag{ 36 }$$

where k₀ is the peak wavenumber and is set to be ${k}_{0}{L}_{\mathrm{box}}=4$ and ${u}_{0}=10$ cms⁻¹. The spectra ${E}_{\theta }(k,0)$ and ${E}_{q}(k,0)$ are of the same functional form as equation (36), with amplitudes of ${\theta }_{0}/{T}_{0}={10}^{-2}$ and ${q}_{0}/{Q}_{{\rm{v}}0}={10}^{-4}$, respectively.

The cloud droplets are assumed to be monodisperse, with an identical radius of 10 μm and are uniformly distributed inside the box, with zero initial velocity. The mean number density of the droplets n_p is 125 cm⁻³, which is a typical value in a cloud. Various parameters are summarized in table 1.

Table 1.  Various parameters. N_p and n_p are the total number and the mean number density of droplets, respectively.

Box length	${L}_{{\rm{box}}}$	25.6 (cm)
Number of droplets	Np	221
Droplet number density	np	125 (cm−3)
Initial temperature at 500 m	Te0	293 (K)
Pressure at sea surface	${P}_{{\rm{sea}}}$	1013.25 (hPa)
Gas constant of dry air	Rd	$287\times {10}^{4}$ (cm2 s−2 K−1)
Gas constant of water vapor	Rv	$462\times {10}^{4}$ (cm2 s−2 K−1)
Specific heat for constant pressure	cp	$1006\times {10}^{4}$ (cm2 s−2 K−1)
Dry adiabatic lapse rate	${{\rm{\Gamma }}}_{{\rm{d}}}$	$9.76\times {10}^{-5}$ (K cm−1)
Latent heat of vaporization	Lv	$2.46\times {10}^{10}$ (cm2 s−2)
Initial saturated vapor pressure	${e}_{{\rm{sat}}}$	$2.34\times {10}^{4}$ (g cm−1 s−2)
Density of saturated water vapor	${\rho }_{{\rm{vs}}}$	$17.30\times {10}^{-6}$ (g cm−3)
Density of dry air	${\rho }_{a}$	$1.21\times {10}^{-3}$ (g cm−3)
Kinematic viscosity of dry air	${\nu }_{a}$	$1.5\times {10}^{-1}$ (cm2 s−1)
Thermal diffusivity	${\kappa }_{T}$	$2.22\times {10}^{-1}$ (cm2 s−1)
Diffusivity of water vapor	${\kappa }_{{\rm{v}}}$	$2.16\times {10}^{-1}$ (cm2 s−1)

In this computation, the buoyancy force is not strong enough to maintain the turbulence in a statistically steady state and random force is applied to the turbulence. The solenoidal random force is assumed to be a Gaussian random with zero mean and to be white noise [55]. It is applied to the low wavenumber band as:

$$\begin{eqnarray}&&\langle {\boldsymbol{f}}({\boldsymbol{k}},t)\rangle =0,\quad \langle {f}_{i}({\boldsymbol{k}},t){f}_{j}(-{\boldsymbol{k}},s)\rangle ={P}_{{ij}}({\boldsymbol{k}})\displaystyle \frac{F(k)}{4\pi {k}^{2}}\delta (t-s),\quad \end{eqnarray} \tag{ 37 }$$

where ${P}_{{ij}}({\boldsymbol{k}})={\delta }_{{ij}}-\tfrac{{k}_{i}{k}_{j}}{{k}^{2}}$ and F(k) is the spectrum of the random force and is given by:

$$\begin{eqnarray}F(k)=\left\{\begin{array}{ll}{c}_{f} & (1\leqslant {{kL}}_{{\rm{box}}}\leqslant 4)\\ 0 & ({\rm{otherwise}}),\end{array}\right.\end{eqnarray} \tag{ 38 }$$

where c_f is a constant and is set to be 70 cm²s⁻³.

The turbulence in the box is statistically homogeneous and axially symmetric due to the buoyancy force in the vertical direction. Therefore, it is reasonable to expand the statistical quantities in terms of the Legendre polynomials [56]. However, we consider only the zeroth order coefficients (i.e., the isotropic parts) for simplicity of analysis. The integral scale, Taylor microscale and the Taylor microscale Reynolds number are defined by:

$$\begin{eqnarray}&&{ \mathcal L }(t)=\displaystyle \frac{3\pi }{4}\displaystyle \frac{{\int }_{0}^{\infty }{k}^{-1}E(k,t){\rm{d}}k}{{\int }_{0}^{\infty }E(k,t){\rm{d}}k},\qquad \lambda (t)=\sqrt{\displaystyle \frac{5\langle {u}_{x}^{2}\rangle }{\langle {({\partial }_{x}{u}_{x})}^{2}\rangle }},\qquad {R}_{\lambda }(t)=\displaystyle \frac{{u}_{{\rm{rms}}}(t)\lambda (t)}{{\nu }_{a}},\end{eqnarray} \tag{ 39 }$$

respectively, and the rate of kinetic energy dissipation per unit mass and the mean Kolmogorov length of the turbulent flow are defined by:

$$\begin{eqnarray}&&\epsilon (t)=2{\nu }_{a}{\int }_{0}^{\infty }{k}^{2}E(k,t){\rm{d}}k,\qquad \eta (t)={\left(\displaystyle \frac{{\nu }_{a}^{3}}{\epsilon (t)}\right)}^{1/4},\end{eqnarray} \tag{ 40 }$$

respectively.

The target range of time evolution for the present DNS is about 20 min and the total number of time steps with the time increment ${\rm{\Delta }}t\quad =$ 1 ms is 1.2 million. Therefore, we examine the small system with the grid number $N={128}^{3}$ as the first step in future large-scale computations. Five computations were made and the numerical and turbulence parameters are summarized in table 2. Runs 1, 2, 3a, and 4 used four different values of μ to examine the entrainment effects. Run 3b was made using $\bar{W}=0$ to compare with the case of zero updraft. The computers used were the K-computer at the Research Organization for Information Science and Technology in Kobe and the FX10 at Nagoya University and the University of Tokyo.

Table 2.  Numerical and mean turbulence parameters. Updraft is On for $\bar{W}\ne 0$, Off for $\bar{W}=0$. μ is the entrainment strength.

	Updraft	μ	${R}_{\lambda }$	E		${ \mathcal L }$	λ	η
		×10−5 (cm−1)		(cm2 s−2)	(cm2 s−3)	(cm)	(cm)	(cm)
Run 1	On	2.5	86	73.3	32.5	5.58	1.84	0.101
Run 2	On	5.0	86	73.8	33.0	5.58	1.84	0.101
Run 3a	On	7.5	86	73.6	32.8	5.58	1.84	0.101
Run 3b	Off	7.5	86	73.6	32.8	5.58	1.84	0.101
Run 4	On	10.0	86	73.3	32.8	5.58	1.84	0.101

4.1. Mean values

Figures 3 and 4 show the time evolution of the mean updraft velocity $\bar{W}(t)$ and the altitude $\bar{H}(t)$ of the parcel for Run 3a (red, color online), respectively. Also the curves of the parcel model (green) are plotted to see the turbulence effects. The $\bar{W}(t)$ attains a maximum of about 8 ms⁻¹, and $\bar{H}(t)$ quickly reaches 2500 m at about 600 s and then slowly increases in time with the oscillation, which are compared to the results of 2.5 ms⁻¹ at 5 min and 7.7 ms⁻¹ at 25 min that were obtained by the bin method in Takahashi [51]. The period of oscillation is about 330 s which is longer than 310 s for the parcel model and about 50% longer than the Brunt Väisälä period that is about 220 s. This deviation is reasonable when considering the difference between the linear stability analysis and the fully nonlinear computation. The difference between the DNS and the parcel model in figures 3 and 4 is due to the fluctuations of the temperature, the water vapor mixing ratio and the droplets radius generated by the turbulent motion. The turbulence tends to slow down the evolution of $\bar{W}(t)$ and $\bar{H}(t)$ and to slightly decrease the oscillation amplitudes. Figure 5 shows the time evolution of the mean temperature $\bar{T}(t)$ (red, color online) inside the parcel and that of the environment (blue, color online) at the same altitude as that of the parcel for Run 3a. When the parcel ascends, $\bar{T}(t)$ decreases below the environmental temperature and increases according to the parcel'(s) descent, and then begins to oscillate. The time evolution of the mean water vapor mixing ratio is similar to that of the mean temperature, but is always larger than the environmental value (figure not shown). Figure 6 indicates the time evolution of each component (temperature B_T, water vapor mixing ratio B_q, and the liquid water mixing ratio B_ql) and the total buoyancy force B for Run 3a. It can be seen that the temperature term dominates the buoyancy force and the liquid water mixing ratio becomes comparable to the temperature variation at later times due to the growth of the cloud droplets. The water vapor mixing ratio makes a small contribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Turbulent flow

Now we consider the evolution of the turbulent flow inside the parcel. To examine the updraft effects, we compare Run 3a for $\bar{W}\ne 0$ with Run 3b for $\bar{W}(t)=0$. However, the differences in the one point statistics of the fluctuating velocity and temperature, $E,\epsilon$ and E_θ are very small for the two cases. Therefore, we describe the updraft effects for the cases when they are appreciable.

Figures 7 and 8 show the time evolution of ${R}_{\lambda }(t)$ and $\epsilon (t)$ for Run 3a, respectively. The turbulent flow is in a statistically steady state and ${R}_{\lambda }(t)$ fluctuates around the mean value of 86 during the computation and the time averaged kinetic energy dissipation rate per unit mass $\epsilon (t)$ is 33 cm² s⁻³. These values mean that the turbulence in this study is very weak when compared to the values in the actual cloud because of the computational limitation. The updraft effects are seen in E_q(t) as shown in figures 9 and 10 for Runs 3a and 3b, respectively. The fluctuation amplitude of E_q(t) for Run 3a is slightly suppressed when compared to the case in the absence of the updraft. This is because the water vapor is more systematically condensed on the droplets when the supersaturated state is maintained by the adiabatic cooling due to the updraft. At times after t = 600 s E_q(t) begins to oscillate and its fluctuation amplitude increases when the parcel is close to the local minimum or maximum altitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although the updraft effects on the one point statistics are appreciable only in the water vapor, the spectrum is more sensitive to the updraft. Figures 11–13 show the time variation of three spectra $E(k,t),{E}_{\theta }(k,t)$, and ${E}_{q}(k,t)$ for Runs 3a and 3b at 760 s (red), 770 s (green), 840 s (blue), 890 s (magenta), 920 s (cyan), 930 s (black). As for the kinetic energy spectrum, there is almost no difference in $E(k,t)$. This is because (1) the mean buoyancy force (which corresponds to the Fourier mode at the origin, $B({\boldsymbol{k}}=0,t)$) generating the uniform updraft is subtracted from the buoyancy force as ${B}^{\prime }=B-\bar{B}$ so that no direct effects of the updraft can be observed in the turbulence and (2) the amplitude of the buoyancy force is smaller by about two order in magnitude than the random force and masked by the random force. On the other hand, ${E}_{\theta }(k,t)$ and ${E}_{q}(k,t)$ are affected by the updraft. They are approximately constant at ${{kL}}_{\mathrm{box}}\lt 10$ and oscillate in time. By comparing the time variation of these spectra with $\bar{H}(t)$ in figure 4, we see that the one period of the oscillation of the spectra roughly corresponds to the half the oscillation period of the parcel altitude. However, this correspondence becomes obscure at latter times. The condensation term ${C}_{{\rm{d}}}({\boldsymbol{x}},t)$ acts as the source of the scalar fluctuations at scales smaller than the Kolmogorov scale (note that Sc = 0.72 so that $\bar{\eta }\approx {\bar{\eta }}_{{\rm{B}}}$, where ${\eta }_{{\rm{B}}}$ is the Batchelor length). The fluctuations of the water vapor and temperature generated by ${C}_{{\rm{d}}}({\boldsymbol{x}},t)$ are advected by the turbulent flow and spread over the entire domain of the box. In the wavenumber space, the excitation of ${C}_{{\rm{d}}}({\boldsymbol{k}},t)$ at high wavenumbers (figure not shown) is quickly transferred toward the low wavenumbers and begins to oscillate in time. The same mechanism can be applied to ${E}_{\theta }(k,t)$ and ${E}_{q}(k,t)$ but without oscillation when the updraft is absent as seen in figures 12(b) and 13(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The spectra ${E}_{\theta }(k,t)$ and ${E}_{q}(k,t)$ have quite large cusps near the high wavenumber cut off, which is in sharp contrast to $E(k,t)$. The cusp has two origins: the first is due to the fact that the fluctuations are injected at small scales below the computational grids, while the second is due to the fact that the contributions of C_d are treated as the point force. This behavior has also been observed and is unavoidable—in the DNS, which used the point particle approximation to the droplet [22]. It is difficult to make a quantitative assessment of the two contributions, and therefore we do not consider this problem further.

4.3. Cloud droplets

Figures 14 and 15 show a comparison of the time evolution of the mean droplet radius and the variance of the droplet radius around the mean for Run 3a (red: color online) and Run 3b (blue), respectively. Also the result of the parcel model is plotted in figure 14 for comparison. The effects of the updraft are very dramatic when compared to the case with no updraft. Starting from 10 μm the $\bar{R}(t)$ grows to about 24 μm by 600 s at which $\bar{H}=2500\;{\rm{m}}$ and begins to oscillate according to the parcel oscillation at the inversion layer. This result is not far from the data of 30 μm at 2000 m above the sea as reported in [51]. Although the oscillation amplitude gradually decreases, the maximum of $\bar{R}(t)$ remains almost unchanged. In contrast, the $\bar{R}(t)$ for the zero updraft case grows very slowly. When compared to the parcel model, the turbulence tends to make the evolution slower and to decrease the oscillation amplitudes.

Zoom In Zoom Out Reset image size

The evolution of the distribution of the cloud droplet radius is a main target of the present simulation. Figure 15 shows the comparison of the evolution of the variance of the droplet radius fluctuations for Runs 3a and 3b. It is interesting to observe that the variance for Run 3a becomes smaller than that for Run 3b. This is because all the droplets grow uniformly in the box due to the decrease of the saturated vapor pressure when the parcel ascends. Figures 16(a) and (b) show the time evolution of the droplet spectrum for Runs 3a and 3b, respectively. The peak of the droplet spectrum for Run 3a moves rapidly to the right with increase of time and then oscillates at around 20 μm, while the droplet spectrum for Run 3b moves slightly to the right and the spectrum width is a little wider than that for Run 3a, which is consistent with figure 15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 17 indicates comparison of the time evolution of the liquid water content (LWC), which is defined by:

$$\begin{eqnarray}&&{m}_{{\rm{l}}}(t)=\displaystyle \frac{4\pi }{3}\displaystyle \frac{{\rho }_{{\rm{l}}}}{{L}_{{\rm{box}}}^{3}}\displaystyle \sum _{j=1}^{{N}_{p}}{R}_{j}^{3}(t).\end{eqnarray} \tag{ 41 }$$

When compared to the parcel model (green) the LWC by the DNS (red) increases slowly and begins to oscillate around 5 gm⁻³, which is consistent with the value observed in the cloud. Figure 18 shows the variation of the mean vertical velocity of the cloud droplets with respect to the moving parcel frame. When the updraft exists, the mean vertical velocity is negative at all times and the amplitude slowly increases, and then begins to oscillate between about −8 and −4 cm s⁻¹, which means that the cloud droplets are slowly falling inside the parcel. However, because the mean updraft velocity is oscillating between about 8 ms and −6 ms⁻¹, the ensemble of cloud droplets is, as a whole, comoving with the parcel. The droplet Reynolds number, based on the mean droplet vertical velocity with respect to the air inside the parcel, is about 0.2, and therefore the use of Stokes' law for the droplet drag in this study is justified.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is well known and discussed in the literature that the entrainment affects the evolution of the updraft velocity and the parcel altitude, which in turn leads to different evolution history of droplets inside the parcel [3]. Since the parcel is assumed to be small and advected inside the whole cloud, unlike the LES model, the entrainment strength can not be determined within the present microscopic approach but is given as the parameter [40–50]. Then it is useful and necessary for the future development of the present method and for large scale computations at very high Reynolds numbers with larger box size to know the variation of the above results to change of the entrainment strength. Different entrainment strength may, roughly speaking, be considered to be the difference in the initial position of the parcel, core region for low μ and near interface region between cloud and clear air for high μ [3].

We performed four computations Runs 1, 2, 3a, and 4 by changing μ while keeping other parameters the same, and these are summarized in table 2. Although these values of μ are larger than the values, for example, those between $2.0\times {10}^{-6}$ and $2.0\times {10}^{-5}\quad {\mathrm{cm}}^{-1}$ used in [53, 54], they were chosen partially by the reason for the numerical stability. Since our main purpose in the present study is to examine the consistency of the new method, we do not enter the detailed arguments for the parameter. The effects of the entrainment on the parcel altitude and the mean droplet radius are compared in figures 19 and 20, respectively. The general trend is that when the entrainment μ decreases the oscillation amplitudes become larger and the oscillation period becomes gradually shorter and the migration of the mean parcel altitude tends to be higher. This trend is observed also in $\bar{W}(t)$ and the LWC (figure not shown) and is consistent with the previous studies [57–60].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we examined the effects of the temperature dependence of the diffusion constant K_s on cloud droplet growth. However, the effects on the mean radius of the cloud droplets were found to be very weak in this study.

We have developed a new numerical method that simulates the detailed dynamics of cloud droplets from the microscopic view point. Novel feature of the present method is the introduction of the parcel advected by the updraft that is self-consistently determined within the theory, so that the seamless computation for the continuous growth of the droplets along the parcel trajectory inside the cloud becomes possible. We call the set of equations of motion and growth of the radius of a large number of cloud droplets, and the equations of turbulent flow, temperature, and water vapor mixing ratio as the cloud microphysics simulator.

Continuous evolution of the droplets and turbulence during 20 min has successfully been simulated by the cloud microphysics simulator which has previously not been possible in simulations resolving all the degrees of freedom of cloud droplets. By the seamless simulation we found that (1) the mean radius of the cloud droplets grows from 10 μm to about 24 μm, (2) the turbulence slows down the evolution of the system when compared to the prediction by the parcel model, and (3) the width of the droplet spectrum becomes narrower than that in the absence of the updraft.

It is worthwhile and important to examine the effects of fluctuations due to the turbulence. In the parcel model all the droplets are assumed to be exposed to the same supersaturation and have the same radius and the same condensation rate, which leads to a systematic build up of the buoyancy force. On the other hand, the turbulent flow advects the droplets, temperature and water vapor so that the particle positions and radii change randomly in space and time. These fluctuations prevent the buoyancy force that is uniform inside the box from building up, as a result the evolution of the mean $\bar{B}$ becomes slower than in the case of the parcel model, and so do for the means of $T({\boldsymbol{x}},t),{Q}_{{\rm{v}}}({\boldsymbol{x}},t),{C}_{{\rm{d}}}({\boldsymbol{x}},t)$ and ${q}_{{\rm{l}}}({\boldsymbol{x}},t)$. This is the reason for the slow down of the mean evolution of the cloud droplets by turbulence although the effects are small for the case of the condensation alone. When the collision-coalescence is included, the effects of turbulence would become stronger because the variation of the mean and fluctuations of $\bar{{q}_{{\rm{l}}}}$ become larger. It should be noted that the total water mass ($M(t)=\int [{M}_{{\rm{d}}}{Q}_{{\rm{v}}}({\boldsymbol{x}},t)+{q}_{{\rm{l}}}({\boldsymbol{x}},t)]{\rm{d}}{\boldsymbol{x}}$ where M_d is the mass of dry air inside the parcel) is conserved (indeed, it was numerically conserved within 0.2% accuracy) in the present set up even when the collision coalescence is included.

In the present study, the updraft generated by the mean buoyancy force is treated as a time dependent uniform flow in the vertical direction which can only indirectly drive the turbulence through the variation of the parcel height. When a whole cloud is considered, there exists non-uniformity of the updraft and the large scale shear appears as driving force to the turbulence. Since the parcel size is so small when compared to the whole cloud, such a mechanism is not incorporated within the present framework and must be introduced as the external force. The random external force in this study was weak in the sense that ${R}_{\lambda }=86$ was very small when compared to the one observed in the cloud, which is of the order of 10³, but quite strong to mask the effects of the buoyancy force ${B}^{\prime }\equiv B-\bar{B}$. The dominance of the random force over the buoyancy force is due to the artificial nature of the white in time property. In the future study a random force with finite memory should be considered [3].

Since the detailed dynamics of the cloud droplets are naturally incorporated in the cloud microphysics simulator, the process of collision-coalescence of the droplets with effects such as the offset distance, the relative droplet sizes, the surface tension, the Reynolds number dependence of the drag, and the hydrodynamic interaction are to be examined. Also the DNS of turbulent flow with a large number of grid points and cloud droplets is required to study the early phase of formation of rain drops and cloud. To conclude the cloud microphysics simulator developed in this study is very successful and promising. Its application to the above issues and further developments are great challenges in turbulence and cloud physics research.

The authors thank Prof. Watanabe for his valuable discussions and comments. This research used the computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science, through the High Performance Computing Infrastructure (HPCI) System Research project (hp120147, hp140024, hp150088). The computational support provided by Japan High Performance Computing and Networking, Large-scale Data Analyzing and Information Systems (JHPCN) (jh130010, jh140025, jh150012), and High Performance Computing (HPC) at Nagoya University are also appreciated. TG was supported by Grants-in-Aid for Scientific Research Nos.15H02218 and 24360068, and IS by No. 15H02218, respectively, from the Ministry of Education, Culture, Sports, Science, and Technology of Japan.

Parcel model is a set of equations of means constructed by taking the volume average of the equations derived in section 2. Then all of the spatial derivative terms of equations (9)–(11) vanish by the periodic boundary condition. Replacing the average of the product of the fluctuating quantities by the product of their means, we obtain the parcel model as:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\bar{W}(t)}{{\rm{d}}t}=\bar{B}(t),\qquad \bar{H}(t)=\bar{H}(0)+{\int }_{0}^{t}\bar{W}({t}^{\prime })\ {\rm{d}}{t}^{\prime },\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\bar{T}(t)}{{\rm{d}}t}=\displaystyle \frac{{L}_{{\rm{v}}}}{{c}_{p}}{\bar{C}}_{{\rm{d}}}(t)-{\rm{\Gamma }}(t)\bar{W}(t),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\bar{Q}}_{{\rm{v}}}(t)}{{\rm{d}}t}=-{\bar{C}}_{{\rm{d}}}(t),\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\bar{R}}^{2}(t)}{{\rm{d}}t}=2{K}_{s}(t)\ \bar{S}(t),\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\bar{V}}_{z}(t)}{{\rm{d}}t}=-\displaystyle \frac{1}{\bar{\tau }(t)}{\bar{V}}_{z}(t)-(g+\bar{B}(t)),\end{eqnarray} \tag{ 46 }$$

where

$$\begin{eqnarray}&&\bar{B}(t)=g\left(\displaystyle \frac{\bar{T}(t)-{T}_{{\rm{e}}}(\bar{H}(t))}{{T}_{{\rm{e}}}(\bar{H}(t))}+\varepsilon [{\bar{Q}}_{{\rm{v}}}(t)-{Q}_{\mathrm{ve}}(\bar{H}(t))]-{\bar{q}}_{{\rm{l}}}(t)\right),\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}&&\bar{S}(t)=\displaystyle \frac{{\bar{Q}}_{{\rm{v}}}(t)}{{Q}_{\mathrm{vs}}(\bar{T}(t))}-1,\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{\bar{C}}_{{\rm{d}}}(t)=4\pi {n}_{p}\ \bar{R}(t)\ {K}_{s}(t)\ \bar{S}(t)\left(\displaystyle \frac{{\rho }_{{\rm{l}}}}{{\rho }_{a}}\right),\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{\bar{q}}_{{\rm{l}}}(t)=\displaystyle \frac{4\pi {\bar{R}}^{3}(t)}{3}{n}_{p}\left(\displaystyle \frac{{\rho }_{{\rm{l}}}}{{\rho }_{a}}\right),\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{\bar{Q}}_{\mathrm{vs}}(t)=\varepsilon \times \displaystyle \frac{610.78}{P(\bar{T}(t),\bar{H}(t))}\mathrm{exp}\left(17.269\ \displaystyle \frac{\bar{T}(t)-273.16}{\bar{T}(t)-35.86}\right),\end{eqnarray} \tag{ 51 }$$

where ${\rm{\Gamma }}$ is given by equations (12)–(14). The numerical computations were made for the same constants, and the same initial and environmental conditions as those in DNS with various entrainment strengths. The results for the entrainment strength $\mu =7.5\times {10}^{-5}\quad {\mathrm{cm}}^{-1}$ are plotted in the green lines in figures 3, 4, 14 and 17.

Equations (42)–(45) are the same as equations (11)–(15) of [16] for the macroscopic variables, although the precise meaning of 'macroscopic' is not given in their derivation. In our derivation, the volume average is defined in equation (5) and the mean values correspond to their macroscopic variables. The differences are: (1) the updraft velocity $\bar{W}$ in [16] is constant and given as parameters and the buoyancy force is dropped, while in the present study it is self-consistently determined through the buoyancy force term; (2) the saturation vapor pressure is given as a function of the altitude using Tetens' formula, but it is not certain whether or not it is constant in [16]; (3) T_e and Q_ve are taken from the observation data of [51], while they are treated as constants in [16]; (4) ${\rm{\Gamma }}(t)$ has two contributions from the effects of the moist air and the entrainment, while it is a constant in [16].