**Research article**
23 Nov 2018

**Research article** | 23 Nov 2018

# Direct Lagrangian tracking simulation of droplet growth in vertically developing cloud

Yuichi Kunishima and Ryo Onishi

**Yuichi Kunishima and Ryo Onishi**Yuichi Kunishima and Ryo Onishi

- Center for Earth Information Science and Technology, Japan Agency for Marine-Earth Science and Technology 3173-25 Showa-machi, Kanazawa-ku Yokohama Kanagawa 236-0001, Japan

- Center for Earth Information Science and Technology, Japan Agency for Marine-Earth Science and Technology 3173-25 Showa-machi, Kanazawa-ku Yokohama Kanagawa 236-0001, Japan

**Correspondence**: Ryo Onishi (onishi.ryo@jamstec.go.jp)

**Correspondence**: Ryo Onishi (onishi.ryo@jamstec.go.jp)

Received: 30 Mar 2018 – Discussion started: 27 Apr 2018 – Revised: 04 Sep 2018 – Accepted: 15 Nov 2018 – Published: 23 Nov 2018

We present a direct Lagrangian simulation that computes key warm-rain processes in a vertically developing cloud, including cloud condensation nuclei (CCN) activation, condensational growth, collisional growth, and droplet gravitational settling. This simulation, which tracks the motion and growth of individual particles, is applied to a kinematic simulation of an extremely vertically elongated quasi-one-dimensional domain, after which the results are compared with those obtained from a spectral-bin model, which adopts the conventional Eulerian framework. The comparison results, which confirm good bulk statistical agreement between the Lagrangian and conventional spectral-bin simulations, also show that the Lagrangian simulation is free from the numerical diffusion found in the spectral-bin simulation. After analyzing the Lagrangian statistics of the surface raindrops that reach the ground surface, back-trajectory scrutiny reveals that the Lagrangian statistics of surface raindrops contains the information about the sky where the raindrops grow like the shape does for snow crystals.

- Article
(6684 KB) -
Supplement
(34852 KB) - BibTeX
- EndNote

Because clouds play central roles in weather and climate, numerous cloud microphysics models have been developed to investigate their physics and predict their development. Conventional models are divided into either bulk or spectral-bin types based on their microphysical representations. In bulk models, all of the microphysical processes, such as mixing ratios or number concentrations of cloud hydrometers (cloud, rain, cloud ice, snow, graupel, etc.) are described in terms of grid-averaged parameters. In contrast, in spectral-bin models, hydrometer size or mass distributions are modeled directly. Although spectral-bin models are more complex and prognose more variables than bulk models, both types are the same in the sense that they prognose grid-averaged values. In other words, all conventional models are based on the Eulerian framework.

However, in recent years, a number of new cloud models that are based on the Lagrangian framework, in which the motion and growth of individual droplets are tracked, have been developed. These Lagrangian cloud models are also divided into two groups: super-droplet models (Dziekan and Pawlowska, 2017; Riechelmann et al., 2012; Shima et al., 2009) and direct tracking models (Chen et al., 2018; Grabowski and Wang, 2013; Kumar et al., 2014; Onishi et al., 2015; Saito and Gotoh, 2018). The former introduces the multiplicity concept, in which each droplet represents multiple droplets that share approximately the same attributes and positions, while the latter directly computes the motion and growth of each droplet and can consider the influence of microscale flow and scalar field on the droplet motion and growth, even though their computational costs are extremely high. For example, Onishi et al. (2015) developed the Lagrangian cloud simulator (LCS), which adopts the Eulerian–Lagrangian framework and can provide reference data for cloud microphysical models. Their study revealed the turbulence enhancement of collisional growth on the temporal evolution of droplet size distributions for a periodic box (cubic) domain. Saito and Gotoh (2018) reported on a direct tracking simulation for a periodic box domain that additionally considers condensational growth while ignoring the cloud condensation nuclei (CCN) activation. Here, it should be noted that other existing direct tracking models still lack the CCN activation process and also rely on periodic box domains, which are useful, but unphysical in the sense that large drops repeatedly reenter the same domain. However, no direct tracking simulations have ever succeeded in representing all of the warm-rain processes, which include CCN activation, condensation–evaporation, collision, and gravitational settling.

With the above background in mind, this study aims to establish a direct Lagrangian model that can compute all of the abovementioned warm-rain processes and obtain Lagrangian statistics on droplet growth in a vertically developing cloud. First, a simple stochastic CCN activation model for direct Lagrangian simulations was developed and implemented in the LCS. The integrated LCS was then applied to a kinematic simulation for a vertically developing warm cloud. The resulting computational model adopts an extremely vertically elongated quasi-one-dimensional domain that can allow gravitational settling of large raindrops in physical (rather than periodic) space at a feasible computational cost.

The rest of this paper is organized as follows. In the next section, we describe the numerical methods, including the newly developed CCN activation model, used for direct Lagrangian simulations. Computational settings, including the domain setting, are described in Sect. 3. Then, results and discussion for both Eulerian and Lagrangian simulations and statistics are presented in Sect. 4, and our conclusions are given in Sect. 5.

## 2.1 Lagrangian cloud simulator (LCS) overview

The LCS (Onishi et al., 2015) adopts a hybrid Eulerian–Lagrangian framework in which flow motion and scalar transport are computed using a Eulerian method, and the particle motion and growth are obtained via a Lagrangian tracking method. Within that framework, cloud microphysics are computed by direct Lagrangian tracking, which follows the motion and growth of individual particles. Previously, the LCS was used primarily to focus on the collisional growth of droplets in a turbulent medium. In this study, in order to include the water droplet phase change, the humidity and temperature fields are calculated using the same Eulerian method that was used for the flow. In addition, to include the droplet activation (nucleation) process, this study presents a newly developed CCN activation model based on Twomey's modeling (Twomey, 1959). The resulting integrated LCS can track the key warm-rain processes of water droplets (CCN activation, condensational and collisional growth, and gravitational settling).

The following subsection describes the present kinematic simulation while Sect. 2.3 describes the present Eulerian methods used for the calculation of flow and scalar fields. Then, Sect. 2.4 describes the present Lagrangian method used for particle phase, which includes the newly developed CCN activation model designed for use with a Lagrangian framework.

## 2.2 Overview of quasi-one-dimensional kinematic simulation

### 2.2.1 Kinematic model

This study adopts the computational settings of the kinematic driver
(KiD) model (Shipway and Hill, 2012), which was originally designed
to provide an intercomparison framework for bulk and spectral-bin
cloud microphysics models. The flow and temperature profiles are prescribed
in a way that minimizes feedback between dynamics and microphysics
and thus facilitate straightforward comparisons. Here, we adopt the
simplest test, which is a shallow convection case (*warm-1*
case) in which a simple updraft is prescribed as

with *w*_{0}=2 m s^{−1}. This prescribed updraft lifts all the
air parcels by ($\int w\mathrm{d}t=$) 764 m. The duration and depth of the
simulation are 3600 s and 3000 m, respectively. Table 1
shows the initial temperature and moisture profiles, which are set to be
similar to those used in the Global Energy and Water Cycle Experiment (GEWEX)
Cloud System Study (GCSS) Rain In Cumulus over the Ocean (RICO) composite
intercomparison. Inside the surface layer below 764 m, the air is assumed to
be well mixed and neutrally stratified. More specifically, the potential
temperature Θ ($=T(p/{p}_{\mathrm{0}}{)}^{R/{C}_{\mathrm{p}}}$, where *T* is the
temperature, *p* is the pressure, *R* is the gas constant of dry air
(*R*=287 J kg^{−1} K^{−1}), *C*_{p} is the specific heat at
constant pressure (*C*_{p}=1004 J kg^{−1} K^{−1}) and
*p*_{0} ($=\mathrm{1.00}\times {\mathrm{10}}^{\mathrm{5}}$ Pa) is the reference pressure), and the initial
vapor mixing ratio *Q*_{υ} are set to the same values inside the
layer.

### 2.2.2 MSSG-Bin simulation (a spectral-bin simulation)

The Multi-Scale Simulator for the Geoenvironment (MSSG;
Matsuda et al., 2018; Takahashi et al., 2013) has a warm-bin–cold-bulk
hybrid cloud microphysics model named MSSG-Bin (Onishi and Takahashi, 2012) in which
the mass coordinate *m* was discretized as ${m}_{k}={\mathrm{2}}^{\mathrm{1}/s}{m}_{k-\mathrm{1}}$, and in
which
*s* was set to 16. The representative radius of the first bin was
2.7 µm and 528 classes were calculated, the largest class of which
had a representative radius of 5.4 mm. The vertical grid spacing was set to
Δ_{z}=25 m and the time interval was set to Δ*t*=1 s.
Following the KiD protocol, the initial CCN concentration was set to
5.0×10^{7} m^{−3} and kept constant for the entire simulation.

## 2.3 Flow and scalar phase in the Eulerian framework

In addition to solving the three-dimensional (3-D) continuity and
Navier–Stokes equations for incompressible flows based on a finite difference
method, the LCS solves the transport equations of scalars *ϕ* as

where *U* is the flow velocity, *κ*_{ϕ} is the diffusion
coefficient and *S*_{ϕ} is the source term. The LCS considers two scalars;
the *ϕ* scalar, which can be either the water vapor mixing ratio
*Q*_{υ} or the potential temperature Θ.

The spatial derivatives were calculated using fourth-order central
differences. The conservative scheme devised by Morinishi et al. (1998) was
used for the advection terms and the second-order Runge–Kutta scheme was used
for time integration. Since the temperature field for the KiD *warm-1*
simulation is fixed, only the transport equation for *Q*_{υ} was
calculated explicitly in this study. The diffusion coefficient for the water
vapor is ${\mathit{\kappa}}_{q}=\mathit{\nu}/\mathit{\text{Sc}}$, where *ν* is the kinematic viscosity
and *Sc* is the Schmidt number (*Sc* = 0.675). The vapor field was
coupled with the particle phase through the source term *S*_{q}, described in
Eq. (12) in Sect. 2.5.

## 2.4 Particle phase in the Lagrangian framework

### 2.4.1 Motion of particles

Under the limit of a large ratio of the density of the particle material to
that of the fluid (${\mathit{\rho}}_{\mathrm{p}}/\mathit{\rho}\gg \mathrm{1}$, where *ρ*_{p} and
*ρ* are the densities of particle and air, respectively), the governing
equation for the *i*th particle is given by

where *V*_{p} is the particle velocity, *U*(*x*_{p},*t*) is
the air velocity at the particle position, *u* is the disturbance flow
velocity caused by the surrounding particles, and *τ*_{p} is the
particle relaxation time defined as
${\mathit{\tau}}_{\mathrm{p}}=(\mathrm{2}/\mathrm{9})({\mathit{\rho}}_{\mathrm{p}}/\mathit{\mu}){r}^{\mathrm{2}},$ where *ρ*_{p} is
the particle density (*ρ*_{p}=1000 kg m^{−3}), *μ* the air
viscosity ($\mathit{\mu}=\mathrm{1.79}\times {\mathrm{10}}^{-\mathrm{5}}$ Pa⋅s at standard atmosphere), and *r*
the particle radius. The impulsive acceleration *f*_{coll}
represents the force resulting from collisions and *g* is the gravitational
acceleration vector ($g=(\mathrm{0},\mathrm{0},-g)$ with *g*=9.81 m s^{−2}). The
nonlinear drag coefficient *α* (Rowe and Henwood, 1961) is obtained by
$\mathit{\alpha}=\mathrm{1}+\mathrm{0.15}{\left({\mathit{\text{Re}}}_{\mathrm{p}}\right)}^{\mathrm{0.687}}$, where the
particle Reynolds number *Re*_{p} is defined by
${\mathit{\text{Re}}}_{\mathrm{p}}=\left(\left|U-{V}_{\mathrm{p}}\right|\phantom{\rule{0.125em}{0ex}}\mathrm{2}{r}_{\mathrm{p}}\right)/\mathit{\nu}$.

One important physical process that is often neglected because of the high
computational cost associated with its simulation is hydrodynamic interaction
(HI) among moving particles. While moving in a flow medium, a particle
induces a flow disturbance (denoted by *u* in Eq. 3) in its
vicinity that may reduce collisions between approaching particle pairs. The
significance of HI on the collisional growth of cloud droplets with
O(10 µm) sizes is well known and has actually been confirmed
recently by a Lagrangian tracking simulation (Onishi et al., 2015). This
study considers HI using the binary-based superposition method (BiSM)
(Onishi et al., 2013), which is one of the so-called superposition methods
(Ayala et al., 2007), and which considers HI based on a pair-wise technique
while assuming that interactions via three or more particles are negligible.
This assumption dramatically reduces computational costs while maintaining
reliability as long as the particle number concentration is small, as is
found in atmospheric clouds. It should be noted that the superposition
method, which assumes Stokes flows around moving particles, is valid for
small particles whose *Re*_{p} is smaller than unity. This
requirement typically corresponds to *r*<40 µm. For larger
droplets, associated errors will grow but the influence of HI will become
insignificant as the particle inertia increases. As a result, the
significance of HI relative to particle inertia would level off at a certain
particle size, the threshold of which would depend on the flow conditions.
While a more precise discussion on such issues could be done with
size-resolving simulations, they are not a central focus of this study.
Furthermore, the lubrication effect is not included in the superposition
method. This would cause a slight overestimation of the collisional growth
rate in the present simulations. Again, while this issue could be solved with
cutting-edge size-resolving simulations, such detailed discussions are out of
the scope of this study.

The second-order Runge–Kutta method was used for the time integration. The
flow velocity *U* at the droplet position was linearly interpolated from
the adjacent grid values. Since the Δ*t* time interval should be
smaller than the relaxation time, it should be set to a very small value if a
tiny particle is included in the domain. In order to avoid extremely small
time intervals, we designate particles whose relaxation times are smaller
than the time interval as tracer particles. In this study, we set the time
interval to $\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{3}}$ s and designate particles with a radii
smaller than 14 µm as tracers. This treatment does not
significantly alter the growth of particles with radii smaller than the
threshold because condensational growth is more dominant than collisional
growth. Note that all the CCN (dry aerosol) particles, whose sizes were not
treated explicitly, were also considered to be tracers.

The growth of the *i*th liquid droplet is calculated as

where the subscripts “act”, “cond”, and “coll” represent the growth due to activation, condensation, and collision, respectively. It should be noted that cond represents the change due to phase change and includes evaporation shrinkage.

### 2.4.2 Growth of particles

### CCN activation

Next, we explain our newly developed stochastic CCN activation method for direct Lagrangian tracking simulations. First, the number density of newly activated liquid droplets at the corresponding time step is diagnosed in each Euler grid cell (i.e., inside a ${\mathrm{\Delta}}_{x}\times {\mathrm{\Delta}}_{y}\times {\mathrm{\Delta}}_{z}$ volume, where Δ is the grid spacing and the subscript denotes the direction) using the bulk air saturation and bulk number liquid droplet density. Then, a stochastic judgment is made on whether or not to activate each CCN particle. If the decision is affirmative, another stochastic judgment is then made on which size is to be set.

The formulation of the liquid droplet activation process is based on the
relationship between the number of activated CCN *N*_{act} and the
supersaturation ratio
${S}_{\mathrm{w}}=({Q}_{\mathit{\upsilon}}-{{Q}_{\mathit{\upsilon}}}_{\mathrm{sw}})/{{Q}_{\mathit{\upsilon}}}_{\mathrm{sw}}$,
where ${{Q}_{\mathit{\upsilon}}}_{\mathrm{sw}}$ is the saturated mixing ratio with
respect to water. The saturated mixing ratio is obtained as
${{Q}_{\mathit{\upsilon}}}_{\mathrm{sw}}=\mathit{\epsilon}{e}_{\mathrm{s}}/\left(p-{e}_{\mathrm{s}}\right)$, where *ε* is the ratio
of the gas constants for dry air and water vapor (i.e.,
$\mathit{\epsilon}=R/{R}_{\mathit{\upsilon}}$=0.622) and *e*_{s} is the saturated vapor
pressure given by the Tetens formula ${e}_{\mathrm{s}}=\mathrm{611.2}\mathrm{exp}\mathit{\left\{}\mathrm{17.67}\right(T-$
$\mathrm{273.15})/(T-\mathrm{29.65}\left)\mathit{\right\}}$ (Murray, 1967). In Twomey (1959), the
relationship between *N*_{act} and *S*_{w} is written in the
form ${N}_{\mathrm{act}}=C{S}_{\mathrm{w}}^{k}$, where *C* and *k* depend on the
CCN type. If we define *S*_{max} as the supersaturation needed to
activate the total particle count *N*_{CCN}+*N*_{w}, where
*N*_{CCN} and *N*_{w} are the number concentrations of dry
CCN and liquid droplets, then *C* can be represented as
$C=({N}_{\mathrm{CCN}}+{N}_{\mathrm{w}}){S}_{\mathrm{max}}^{-k}$. Thus, the number of
activated CCN (tiny liquid droplets) can be expressed as

The number of newly activated droplets is calculated as

Assuming the maritime conditions given in Onishi and Takahashi (2012), the
parameters *k* and *S*_{max} were set to 0.6 and 0.008,
respectively. The conversion from a dry aerosol to a liquid droplet is
stochastically determined using the probability
*N*_{act}∕*N*_{CCN} for each dry aerosol particle. The
activated (nucleated) liquid droplet size (*r*_{act}) is then
determined stochastically so that its size distribution follows the
exponential probability distribution given in Soong (1974) as

where $\stackrel{\mathrm{\u203e}}{r}$ is the radius of an average mass droplet and it was set to 11 µm (Onishi and Takahashi, 2012). This stochastic procedure is introduced to avoid detail calculations of the CCN activation process, in which microscale supersaturation fluctuations and CCN density size fluctuations are to be properly represented. The present procedure reproduces realistic size distributions of activated droplets without those detail calculations.

The mass growth due to this activation process is then calculated as

### Condensational growth

The phase change for the *i*th particle is written as (Houze, 1993)

Here, *F*_{d} represents a thermodynamic term associated with vapor
diffusion and *F*_{k} is associated with heat conduction. The terms with
coefficients *α* and *β* represent the curvature effect and the
reduction in vapor pressure due to CCN hydrophilicity, respectively. This
study simply ignored the solute effect and adopted
$\mathit{\alpha}=\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{7}}$ m K and *β*=0 m^{3}. *R*_{υ} is the
gas constant for water vapor (*R*_{υ}=461.7 J kg^{−1} K^{−1}),
*κ*_{q} is the diffusion coefficient for water vapor, and *K* is the
coefficient of air thermal conductivity. It should be noted that
Eq. (9) also covers particle evaporation. When
undersaturated (i.e., *S*_{w}<0), a liquid droplet shrinks due to
evaporation. When a particle becomes smaller than the threshold, it reverts
to the CCN (dry aerosol) category. The threshold radius was set to
1 µm.

The mass growth rate is obtained from

### Collisional growth

It was assumed that colliding particles are immediately united by
conserving mass and momentum. In other words, the coalescence efficiency
was unity, and subsequent breakups (collisional breakups) were not
considered. This no-breakup assumption can be justified for shallow
clouds like the present KiD *warm-1* case, and
none of the other currently used spectral-bin models actually consider
breakups for such cases. A collision is judged from the trajectories
of a droplet pair by assuming linear particle movement for the time
interval Δ*t*.

## 2.5 Phase coupling

This study adopts one-way momentum coupling, thereby neglecting the momentum
feedback from the particle phase to the flow phase. Since the mass fraction
of the particles to flow is O(10^{−3}) in the present simulations, this is
easily justified by the dilute condition. The vapor and particle phases are
coupled in order to conserve the water content. The coupling is done via the
source term *S*_{q} in Eq. (2).

## 3.1 Quasi-one-dimensional computational domain

Table 2 shows the computational settings for
the present LCS simulations. Two differently sized domains were used. This
study focuses on the smaller size simulation (NORM), which has the
horizontal size of (${L}_{x}={L}_{y}=$) 0.01 m and a vertical size of 3000 m.
Roughly speaking, the influential length scale of HI relative to the particle
diameter is approximately 10 for small *Re*_{p} values
(Ayala et al., 2007; Onishi et al., 2013), although it would be smaller for larger
particles. The maximum raindrop diameter is around 10^{−3} m in the
present simulations. Since the horizontal domain size is 10 times larger than
the maximum diameter, the HI of a particle would not affect itself
unphysically through the horizontal periodic boundary condition even for the
largest drops.

The large size simulation (LARGE), the volume of which was 9 times larger than that of NORM, was performed in order to confirm the robustness of the NORM simulation as will be discussed in Sect. 4.1. The number of initial particles was set to 9 times larger than that for NORM in order to keep the particle number concentration the same.

The mean separation length ${l}_{\mathrm{sep}}\left(={n}_{p\mathrm{0}}^{-\mathrm{1}/\mathrm{3}}\right)$ was
$\mathrm{2.71}\times {\mathrm{10}}^{-\mathrm{3}}$ m. The horizontal grid spacing
Δ_{x} ($={L}_{x}/{n}_{x}$, where *n*_{x} is the number of grid points in
the *x* direction) was set similar to *l*_{sep}, while the vertical
grid spacing Δ_{z} ($={L}_{z}/{n}_{z}$) was set longer than the maximum
path of the largest drop within a single time interval (∼O(10) m s^{−1} × O(10^{−3}) s $\sim O\left({\mathrm{10}}^{-\mathrm{2}}\right)$ m).

## 3.2 Initial and boundary conditions

The air flow and the temperature field were prescribed as shown in Sect. 2.2.1. The pressure was determined by the hydrostatic equilibrium as

The obtained pressure and temperature are also listed in
Table 1. Air viscosity depends on temperature and decreases
by about 5 % from 297.9 K (domain bottom) to 280.4 K (domain top). This
study simply neglected this change and used a constant value:
$\mathit{\mu}=\mathrm{1.79}\times {\mathrm{10}}^{-\mathrm{5}}$ Pa s at 298 K and 1 atm. Initially, dry
aerosols (CCN) were randomly and uniformly distributed. The initial number
density was ${n}_{p\mathrm{0}}=\mathrm{5.00}\times {\mathrm{10}}^{\mathrm{7}}$ m^{−3}, which is the same value used
in the MSSG-Bin simulation.

Periodic boundary conditions were applied in all three directions for the
particle field, except for the particles that reached the ground surface
(*z*=0). Those particles (hereafter, surface raindrops) are removed from the
system immediately after reaching the ground. Due to the prescribed updraft
(see Eq. 1) some particles move out from the top boundary
and reenter the domain from the bottom (ground surface). Although this
sounds unphysical, it does not alter the physical result since all the
reentry particles are dry CCN (not activated liquid droplets). The
prescribed updraft lifts all the air parcels by ($\int w\mathrm{d}t=$)
764 m, while there is an over 1000 m vertical gap between the domain top
(*z*=3000 m as in Table 1) and cloud top (*z*∼2000 m)
as will be seen in Fig. 3. This simple periodic treatment for
the vertical direction precisely conserves the number of particles, thus
keeping the particle number concentration near the surface, until the first
raindrop reaches the surface and disappears.

Periodic boundary conditions were applied in horizontal directions
for the vapor field, while the Neumann boundary condition with zero
gradient was applied in the vertical direction. This restricts *Q*_{υ}
near the surface to its initial value.

## 3.3 Sensitivity of computational settings and computational performance

We performed 30 production runs for NORM and one for LARGE. Each run was intense and required massively parallel high-performance computing. One simulation for LARGE required 11 times longer to complete than was required for NORM, while the domain size and number of particles for LARGE were 9 times larger than those for NORM. In NORM, the ratio of the number of particles to that of grid points was 21.7, which means the particle calculation was much heavier. The particle calculation component actually occupied 98 % of the total computational cost.

Next we checked the effect of the computational settings used for grid
spacing and time interval on the results obtained by performing a finer-grid-resolution simulation for NORM with a half-size
grid, i.e., number of grids
8 (=2^{3}) times larger, and comparing the results to the
reference setting in Table 2. We also
performed an additional NORM simulation using time intervals that
were finer by a factor of 2, i.e., with $\mathrm{\Delta}t=\mathrm{6.25}\times {\mathrm{10}}^{-\mathrm{4}}$ s. The
results from those additional simulations were within the statistical
fluctuations summarized in Table 3, thus indicating the
validity of the present computational settings.

## 4.1 Bulk statistics

Figure 2 shows the temporal evolution of water paths obtained from the present LCS together with those from the MSSG-Bin. The threshold radius between cloud and rain was set to 40 µm. All the results from the 30 LCS runs are drawn. However, it should be noted that even though the LCS results look like a single thick line, each LCS run is slightly different from the others as denoted by standard deviations, as shown in Table 3.

The total liquid water path increases due to the supersaturation caused by
the prescribed updraft for *t*≤600 s, which remains constant after the
updraft stops until the first raindrop reaches the ground surface and is
removed from the system. Rainwater appears at around *t*∼600 s and
increases rapidly until *t*∼1350 s when the raindrops start to reach the
surface. These results are consistent between the LCS and MSSG-Bin results.
However, there is a difference in the maximum rainwater path, which can be
attributed to the difference in the CCN treatment. Following the KiD
protocol, the MSSG-Bin simulation kept *N*_{CCN} constant while the
LCS simulations did not. *N*_{CCN} in the LCS simulation decreased
with time since CCN particles were consumed for activation. This led to
larger raindrops, and consequently to larger rainwater amounts.

Table 3 shows a comparison of the bulk statistical values between the NORM and LARGE LCS runs. The values for LARGE range within the statistical variations in NORM. This suggests that the limited size of the NORM computational domain did not matter for bulk statistics.

Figure 3 shows the temporal evolution of the vertical profile of
liquid water content. The vertical spacing of the LCS result for the plot was
set to 25 m, which is the same vertical grid spacing used in the MSSG-Bin
simulation. The results show that the cloud bottom exists at around
*z*=600 m and that the cloud top is at *z*=2000 m. Raindrops appear from
the bottom at around *t*=1100 s. These results are consistent for both the
LCS and MSSG-Bin simulations. However, a remarkable difference is seen in the
raindrop settling path for $\mathrm{600}<t<\mathrm{1100}$ s. The LCS result shows a sharp
ridge, drawn by the contour line for $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{3}}$ kg m^{−3},
inside the cloud layer, while the MSSG-Bin result does not. The MSSG-Bin
somehow suffers from numerical diffusion because it adopts the Eulerian
approach, while the LCS that uses a Lagrangian tracking approach does not.
This clearly shows that the LCS can provide robust diffusionless numerical
results.

## 4.2 Lagrangian statistics

One of the strengths of direct Lagrangian simulations is that they
can provide Lagrangian statistics on histories of individual Lagrangian
particles. This study focuses on the Lagrangian statistics of surface
raindrops, which are the raindrops that reach the ground surface (*z*=0 m). The 30 NORM simulations obtained a total of 36 018 surface
raindrops, which means that each run obtained 1200 surface raindrops
on average.

Figure 4 shows the initial locations of particles that comprise the 36 018 surface raindrops used in the 30 NORM runs. The total number of constituent particles was 119 490 655, which means that each surface raindrop contained an average of 3320 particles. The prescribed updraft lifted the air by 764 m. The air parcel was lifted by 764 m, while the cloud top was about 2000 m as in Fig. 3. The CCN initially located below about 1200 m (2000–764 m) in height participated in the cloud layer formation but those between 900 and 1200 m did not participate in the formation of the surface raindrops. The particles that were initially below 0 m can be considered in the same manner as those transported into the domain along, or generated near, the surface. This figure shows that surface raindrops consist of the CCN particles initially found below 900 m in height. This kind of information, which cannot be obtained from the conventional Eulerian-based bulk or spectral-bin simulations, can be used in other studies such as investigations into the chemical compositions of surface raindrops.

Figure 5 plots the surface raindrop volume,
*V*_{sr}, against the number of constituent particles,
*N*_{memb}. The regression line shows
${V}_{\mathrm{sr}}=\mathrm{2.38}\times {\mathrm{10}}^{-\mathrm{14}}{N}_{\mathrm{memb}}$, which means that the
average diameter of each constituent is 36.1 µm. For example, a
surface raindrop with a 1 mm diameter, whose volume is
$\mathrm{5.24}\times {\mathrm{10}}^{-\mathrm{10}}$ m^{−3}, consists of 22 000 constituent particles.
Recalling that the mean diameter of the nucleated droplets was set to
22.0 µm (see Sect. 2.4.2), it can be
surmised that condensation caused each constituent to grow from 22.0 to
36.1 µm, and then to grow further due to collisions. This is
consistent with findings that show collisional growth is typically dominant
for droplets with diameters larger than 40 µm.

Figure 6 shows the probability density function (PDF) of the
surface raindrops when they reach 600 m height. Raindrops smaller than
200 µm in radius evaporated completely in the unsaturated air below
the cloud bottom (*z*∼600 m) before reaching the ground surface, and thus
could not become surface raindrops. That is, the unsaturated air near the
surface acts as a barrier that prevents small raindrops from reaching the
ground surface. This process cannot be represented in the parcel
(zero-dimensional) concept and it is thus often neglected in theoretical
works. For example, Kostinski and Shaw (2005) and Wilkinson (2016)
discussed, neglecting that process, the stochastic aspect of the onset of
surface precipitations and showed possible large deviations in the time
required for droplets to grow, via collisions, from 10 to 50 µm in
radius. The present result, however, suggests that the time required to grow
from 10 to 200 µm, rather than to 50 µm, in radius should
be measured for the discussion of the rapid onset of surface precipitations

Figure 7 shows the height against time collision history of the
first surface raindrop, which reached the ground surface first in each run.
The surface precipitation time *T*_{surf} was 1347 s. Here, we
identify the “mother” particle from the tens of thousands of constituent
particles making up the raindrop, each of which has its own ID number in the
LCS simulation. When two particles collide and unite, the collector particle
grows and the collected particle is removed, thereby conserving mass and
momentum. In this process, the particle with the highest attained height of
the raindrop constituents survives and its ID number is affixed to the
raindrop. That particle is then designated as the mother of the surface
raindrop. It should be noted that one particle attained its maximum height at
*t*∼600 s when the updraft halted. In Fig. 7a, the trajectory
of the mother particle is drawn as a line. The maximum height,
*H*_{max}, attained by the mother particle can be used to
characterize the surface raindrop. It should also be noted that the
trajectory of a particle with constant settling velocity is drawn in a
straight line, whereas the trajectory line of the mother particle appears to
be parabolic. This result is attributed to the increasing settling velocity
that results from collisional growth.

## 4.3 First 10 surface raindrops

Figure 8 shows the relationship between the maximum
attained height of the mother particle for each surface raindrop
*H*_{max} against its surface precipitation time *T*_{surf}
for all 30 NORM runs. In total, 36 018 dots are drawn in
the figure. The black dashed line is obtained by least-squares regression for
all the surface raindrops and is written as

The positive correlation between *H*_{max} and
*T*_{surf} means that larger *H*_{max} leads to larger
*T*_{surf} because it takes longer for a droplet to fall far to the
surface.

Interestingly enough, if restricted to the surface raindrops ranked as the
first 10 with the earliest *T*_{surf} for each run, the correlation
becomes negative as shown in Fig. 8b. Although the plots
are very scattered, a statistical test has proven that a positive correlation
is rejected with a 76 % confidence level, meaning that a negative
correlation is very likely. We have also verified that the confidence level
depends on computational conditions. If the evaporation rate is artificially
strengthened by setting *F*_{k}=0 in Eq. (9), the confidence
level becomes 99 % (virtually certain). This adds extra confidence to the
supposition that this correlation reversal is possible under specific
conditions.

In order to reveal the correlation reversal mechanism, we additionally
collected *T*_{600}, the time when the raindrop passes *z*=600 m (i.e., near
the cloud bottom). Figure 9 plots *H*_{max}
against *T*_{600}. The black dashed line is obtained by least-squares
regression for all the surface raindrops and is written as

The sensitivity of *H*_{max} to *T*_{600} is stronger, i.e., the
slope is steeper than that for *T*_{surf}. Raindrops with higher
*H*_{max} take longer to fall from *H*_{max} to the cloud
bottom (*z*∼600 m). However, they can also collect more droplets during
the longer fall and grow larger, and thus take less time to fall from 600 m
to the surface. This weakens the correlation between *T*_{600} and
*T*_{surf} and consequently results in the stronger sensitivity
(steeper slope in the linear form) of *H*_{max} to *T*_{600} than
to *T*_{surf}.

In contrast to Fig. 8, the correlation between
*H*_{max} and *T*_{600} for the first 10 is positive in
Fig. 9. The sensitivity of *H*_{max} to
*T*_{600} for the first 10 raindrops is somewhat stronger than that for all
the other surface raindrops. That is, the correlation reversal happens while
the first 10 surface raindrops fall from the cloud bottom to the surface. These
facts suggest the plausible scenario that larger drops among the first 10
surface raindrops, with higher *H*_{max}, pass the cloud bottom
later and then overtake other raindrops before they can reach the surface.
This kind of overtake would only occur if specific conditions are met.

This scenario can be confirmed in the Supplement (see Movie S1). The movie view consists of three windows; domain side view
(left window), liquid water path (right top window) and look-up view from the
ground (right bottom window). The look-up view shows that red spheres
(raindrops about 500 µm in diameter) appear from the cloud bottom
at around *t*=1100 s and that the blue spheres (raindrops about 1 mm in
diameter) that appear just after tend to overtake the red spheres before they
can reach the ground.

The Lagrangian statistics provide information about the surface raindrop growth history. They also show differences for the overall surface raindrops and differences for a fraction of the earliest raindrops. Another set of simulations shows that these differences are sensitive to the evaporation rate. These results suggest that the Lagrangian statistics of the surface raindrops have the potential to provide useful information for estimating atmospheric conditions aloft. They also remind us of the Nakaya diagram (Nakaya, 1954), which can translate the shapes of snow crystals into information on atmospheric conditions experienced by those crystals as they fell to the surface. The present Lagrangian simulation results clearly show that statistics of raindrops, as well as snow crystals, contain information about the sky that formed them.

The Lagrangian cloud simulator (LCS; Onishi et al., 2015) adopts a Eulerian–Lagrangian hybrid framework, in which the flow motion and scalar transport are computed using a Eulerian method, and the particle motion and growth are calculated using a Lagrangian tracking method. In that framework, cloud microphysics are based on direct Lagrangian tracking, which tracks the motion and growth of individual particles. In this study, we developed a CCN activation model for direct Lagrangian tracking simulations and implemented it in the LCS. This implementation enables complete tracking of particle growth from a CCN to surface raindrop, including the CCN activation, condensation–evaporation, collisions, and gravitational settling stages. It should also be noted that HI is considered for colliding droplet pairs in our newly developed LCS, which is currently the only model that can consider HI for tracking the simulation of droplet collisional growth at a reasonable computational cost.

The integrated 3-D LCS has previously been applied to a kinematic shallow
convection cloud development in the case of the *warm-1* KiD model
(Shipway and Hill, 2012), which was originally designed to provide an
intercomparison framework for bulk and bin cloud microphysics models. For
the kinematic simulation used in this study, an
extremely vertically elongated quasi-one-dimensional (quasi-1-D) computational
domain has been adopted. Use of this domain enables the representation of
vertical cloud development while keeping the computational domain volume (and
consequently the total number of particles to track) small. The unique
combination of the high-performance LCS and the
extremely vertically elongated domain has allowed us, for the first time, to
successfully complete a tracking simulation of droplet growth from CCN (dry
aerosols) to raindrops in a reasonably realistic manner. For example, the
process of droplet falling through unsaturated air below cloud layers, which
is ignored in the parcel (zero-dimensional) modeling, can be represented. The
air below the cloud layers is unsaturated (particularly for the onset time of
surface precipitations) and the unsaturated air acts as a barrier that makes
it difficult for raindrops to pass through via evaporation. A raindrop
must grow large enough to overcome this evaporation barrier in order to
become a surface raindrop. The present quasi-1-D approach can provide a
practical platform for realistic discussion on the statistical
deviations of the onset of surface precipitations, for example.

The results obtained from the LCS simulations are first compared to those from the MSSG-Bin model (Onishi and Takahashi, 2012), which is a spectral-bin cloud microphysics model based on the conventional Eulerian framework. The comparison results show good agreement, in terms of bulk statistics such as water mixing ratios, between the LCS and MSSG-Bin, except for slight variations due to CCN treatment differences. Furthermore, it has been confirmed that the LCS is free from particle-phase numerical diffusion.

One of the strengths of direct Lagrangian simulations is that they can provide Lagrangian statistics, on histories of individual Lagrangian particles. This study has focused on surface raindrops, which are the raindrops that reached the ground surface. For example, a surface raindrop with a 1 mm diameter consists of approximately 22 000 constituent particles. This means that each constituent particle grew to 36.1 µm in diameter due to condensation and then grew to 1 mm in diameter because of collisions. Additionally, by applying back-trajectory analysis to the surface raindrops, it was possible to identify the particle that attained the maximum height among all the constituent particles, which we define as the mother particle for a surface raindrop. These analysis results also show that the maximum attained height of mother particles has a positive correlation on their surface precipitation time. This can be easily explained by the fact that it takes more time for a particle to fall from higher altitudes. Interestingly, however, if restricted to the first 10 (the earliest 0.8 %) surface raindrops with the earliest surface precipitation time, the correlation was negative. As a plausible mechanism for this, it was suggested that larger drops among the first 10 surface raindrops, which are lifted up higher, pass the cloud bottom later but overtake the slower falling raindrops before they can reach the surface.

Based on the results obtained in this study, we can conclude that the use of Lagrangian statistics for evaluating surface raindrops can provide useful information for estimating atmospheric conditions aloft. These results also remind us of the Nakaya diagram (Nakaya, 1954), which can interpret the shapes of snow crystals into information on the atmospheric conditions experienced by the crystals. Since our present Lagrangian simulation clearly shows that the statistics of surface raindrops, as well as snow crystals, contain information about the conditions aloft, additional discussion becomes possible, thereby indicating that direct Lagrangian models can be powerful and promising meteorological tools.

Access to the simulation results can be granted upon request, under a collaborative framework with JAMSTEC.

The supplement related to this article is available online at: https://doi.org/10.5194/acp-18-16619-2018-supplement.

YK modified the LCS model for the present computational conditions, conducted the simulations and analyzed the results. RO directed this project and wrote the manuscript.

The authors declare that they have no conflict of interest.

This research was supported by the Japanese Ministry of Education, Culture,
Sports Science and Technology (MEXT) as an “Exploratory Challenge on Post-K
computer” (Challenge of Basic Science–Exploring Extremes through
Multi-Physics and Multi-Scale Simulations). The present numerical simulations
were tested on the Earth Simulator of the Japan Agency for Marine-Earth
Science and Technology (JAMSTEC), while the final runs were performed on the
K computer installed at the RIKEN Advanced Institute for Computational
Science.

Edited by: Eric Jensen

Reviewed by: two anonymous referees

Ayala, O., Grabowski, W. W., and Wang, L.-P.: A hybrid approach for simulating turbulent collisions of hydrodynamically-interacting particles, J. Comput. Phys., 225, 51–73, https://doi.org/10.1016/j.jcp.2006.11.016, 2007. a, b

Chen, S., Yau, M. K., and Bartello, P.: Turbulence Effects of Collision Efficiency and Broadening of Droplet Size Distribution in Cumulus Clouds, J. Atmos. Sci., 75, 203–217, https://doi.org/10.1175/JAS-D-17-0123.1, 2018. a

Dziekan, P. and Pawlowska, H.: Stochastic coalescence in Lagrangian cloud microphysics, Atmos. Chem. Phys., 17, 13509–13520, https://doi.org/10.5194/acp-17-13509-2017, 2017. a

Grabowski, W. W. and Wang, L.-P.: Growth of Cloud Droplets in a Turbulent Environment, Annu. Rev. Fluid Mech., 45, 293–324, https://doi.org/10.1146/annurev-fluid-011212-140750, 2013. a

Houze, R. A.: Cloud Dynamics, vol. 53 of International Geophysics, Academic Press, 1993. a

Kostinski, A. B. and Shaw, R. A.: Fluctuations and Luck in Droplet Growth by Coalescence, Bull. Am. Meteor. Soc., 86, 235–244, https://doi.org/10.1175/BAMS-86-2-235, 2005. a

Kumar, B., Schumacher, J., and Shaw, R. A.: Lagrangian Mixing Dynamics at the Cloudy Clear Air Interface, J. Atmos. Sci., 71, 2564–2580, https://doi.org/10.1175/JAS-D-13-0294.1, 2014. a

Matsuda, K., Onishi, R., and Takahashi, K.: Tree-crown-resolving large-eddy simulation coupled with three-dimensional radiative transfer model, J. Wind Eng. Ind. Aerod., 173, 53–66, https://doi.org/10.1016/j.jweia.2017.11.015, 2018. a

Morinishi, Y., Lund, T. S., Vasilyev, O. V., and Moin, P.: Fully Conservative Higher Order Finite Difference Schemes for Incompressible Flow, J. Comput. Phys., 143, 90–124, https://doi.org/10.1006/jcph.1998.5962, 1998. a

Murray, F. W.: On the Computation of Saturation Vapor Pressure, Journal of Applied Meteorology, 6, 203–204, https://doi.org/10.1175/1520-0450(1967)006<0203:OTCOSV>2.0.CO;2, 1967. a

Nakaya, U.: Snow Crystals: Natural and Artificial, Harvard University Press, 1954. a, b

Onishi, R. and Takahashi, K.: A Warm-Bin–Cold-Bulk Hybrid Cloud Microphysical Model, J. Atmos. Sci., 69, 1474–1497, https://doi.org/10.1175/JAS-D-11-0166.1, 2012. a, b, c, d

Onishi, R., Takahashi, K., and Vassilicos, J. C.: An efficient parallel simulation of interacting inertial particles in homogeneous isotropic turbulence, J. Comput. Phys., 242, 809–827, https://doi.org/10.1016/j.jcp.2013.02.027, 2013. a, b

Onishi, R., Matsuda, K., and Takahashi, K.: Lagrangian Tracking Simulation of Droplet Growth in Turbulence-Turbulence Enhancement of Autoconversion Rate, J. Atmos. Sci., 72, 2591–2607, https://doi.org/10.1175/JAS-D-14-0292.1, 2015. a, b, c, d, e

Riechelmann, T., Noh, Y., and Raasch, S.: A new method for large-eddy simulations of clouds with Lagrangian droplets including the effects of turbulent collision, New J. Phys., 14, 065008, https://doi.org/10.1088/1367-2630/14/6/065008, 2012. a

Rowe, P. N. and Henwood, G. A.: Drag forces in a hydraulic model of a fluidised bed – Part I, T. I. Chem. Eng.-Lond., 39, 43–54, 1961. a

Saito, I. and Gotoh, T.: Turbulence and cloud droplets in cumulus clouds, New J. Phys., 20, 023001, https://doi.org/10.1088/1367-2630/aaa229, 2018. a, b

Shima, S., Kusano, K., Kawano, A., Sugiyama, T., and Kawahara, S.: The super-droplet method for the numerical simulation of clouds and precipitation: a particle-based and probabilistic microphysics model coupled with a non-hydrostatic model, Q. J. Roy. Meteor. Soc., 135, 1307–1320, https://doi.org/10.1002/qj.441, 2009. a

Shipway, B. J. and Hill, A. A.: Diagnosis of systematic differences between multiple parametrizations of warm rain microphysics using a kinematic framework, Q. J. R. Meteorol. Soc., 138, 2196–2211, https://doi.org/10.1002/qj.1913, 2012. a, b

Soong, S.-T.: Numerical Simulation of Warm Rain Development in an Axisymmetric Cloud Model, J. Atmos. Sci., 31, 1262–1285, https://doi.org/10.1175/1520-0469(1974)031<1262:NSOWRD>2.0.CO;2, 1974. a

Takahashi, K., Onishi, R., Baba, Y., Kida, S., Matsuda, K., Goto, K., and Fuchigami, H.: Challenge toward the prediction of typhoon behaviour and down pour, J. Phys.-Conf. Ser., 454, 012072, https://doi.org/10.1088/1742-6596/454/1/012072, 2013. a

Twomey, S.: The nuclei of natural cloud formation part II: The supersaturation in natural clouds and the variation of cloud droplet concentration, Geofisica pura e applicata, 43, 243–249, https://doi.org/10.1007/BF01993560, 1959. a, b

Wilkinson, M.: Large Deviation Analysis of Rapid Onset of Rain Showers, Phs. Rev. Lett., 116, 018501, https://doi.org/10.1103/PhysRevLett.116.018501, 2016. a