**Research article**
08 Feb 2019

**Research article** | 08 Feb 2019

# Turbulent enhancement of radar reflectivity factor for polydisperse cloud droplets

Keigo Matsuda and Ryo Onishi

**Keigo Matsuda and Ryo Onishi**Keigo Matsuda and Ryo Onishi

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

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

**Correspondence**: Keigo Matsuda (k.matsuda@jamstec.go.jp)

**Correspondence**: Keigo Matsuda (k.matsuda@jamstec.go.jp)

Received: 14 May 2018 – Discussion started: 13 Jun 2018 – Revised: 05 Nov 2018 – Accepted: 07 Jan 2019 – Published: 08 Feb 2019

The radar reflectivity factor is important for estimating cloud microphysical properties; thus, in this study, we determine the quantitative influence of microscale turbulent clustering of polydisperse droplets on the radar reflectivity factor. The theoretical solution for particulate Bragg scattering is obtained without assuming monodisperse droplet sizes. The scattering intensity is given by an integral function including the cross spectrum of number density fluctuations for two different droplet sizes. We calculate the cross spectrum based on turbulent clustering data, which are obtained by the direct numerical simulation (DNS) of particle-laden homogeneous isotropic turbulence. The results show that the coherence of the cross spectrum is close to unity for small wave numbers and decreases almost exponentially with increasing wave number. This decreasing trend is dependent on the combination of Stokes numbers. A critical wave number is introduced to characterize the exponential decrease of the coherence and parameterized using the Stokes number difference. Comparison with DNS results confirms that the proposed model can reproduce the ${r}_{\mathrm{p}}^{\mathrm{3}}$-weighted power spectrum, which is proportional to the clustering influence on the radar reflectivity factor to a sufficiently high accuracy. Furthermore, the proposed model is extended to incorporate the gravitational settling influence by modifying the critical wave number based on the analytical equation derived for the bidisperse radial distribution function. The estimate of the modified model also shows good agreement with the DNS results for the case with gravitational droplet settling. The model is then applied to high-resolution cloud-simulation data obtained from a spectral-bin cloud simulation. The result shows that the influence of turbulent clustering can be significant inside turbulent clouds. The large influence is observed at the near-top of the clouds, where the liquid water content and the energy dissipation rate are sufficiently large.

Radar remote sensing is widely used for observing a spatial distribution of
cloud and precipitation particles because it can also provide estimates of
cloud microphysical properties. The remote-sensing data are analyzed and
displayed using the radar reflectivity factor (mm^{6} m^{−3}), *Z*, which
is obtained with the following radar equation:

where *P*_{r} and *P*_{t} are the received and transmitted microwave
powers, *G* is the antenna gain, *A*_{e} is the effective
aperture of the antenna, *k*_{m} is the microwave wave number, *K* is the
dielectric coefficient of a water droplet, *V* is the measurement volume, and
*R* is the distance between the antenna and the cloud. The relationship
between the radar reflectivity factor and cloud microphysical properties is
usually expressed based on the assumption of incoherent scattering
(Gossard and Strauch, 1983). Incoherent scattering implies random and
uniform dispersion of cloud droplets (Bohren and Huffman, 1983). For this
case, the factor is proportional to the sum of the scattering intensities of
the individual droplets in the measurement volume. In contrast, droplets form
a nonuniform spatial distribution in turbulence; i.e., inertial particles
concentrate in small-enstrophy regions during turbulence due to the
centrifugal motion (Maxey, 1987; Squires and Eaton, 1991; Wang and Maxey, 1993; Chen et al., 2006). This
preferential concentration is often referred to as turbulent clustering.
Nonuniform distribution of cloud droplets results in coherent scattering,
which is also referred to as particulate Bragg scattering
(Kostinski and Jameson, 2000). For this case, the interference of
microwaves scattered by nonuniformly distributed droplets increases the
scattered microwave intensity; i.e., the radar reflectivity factor increases
due to particulate Bragg scattering. It should be noted that Bragg scattering
often indicates coherent scattering due to a nonuniform distribution of the
refractive index of clear air. In the troposphere, turbulent mixing of
temperature and water vapor causes spatial variation of the refractive index.
To distinguish this effect from the particulate Bragg scattering, it is
specifically referred to as clear-air Bragg scattering. The radar
reflectivity factor for both types of Bragg scattering is dependent on the microwave
frequency *f*_{m} ($={k}_{\mathrm{m}}c/\mathrm{2}\mathit{\pi}$, where *c* is the speed of light),
while the factor for incoherent scattering is independent of *f*_{m}.
Knight and Miller (1998) reported radar frequency dependence of their
observation results for small warm cumulus clouds using S- and X-band
microwaves, which have wavelengths of 10 and 3 cm, respectively. Their
observation data of S-band radar show a characteristic echo pattern of the
mantle echo, in which a strong radar echo was observed in cloud edges, while it
is relatively weak in cloud core regions. The mantle echo can be explained by
clear-air Bragg scattering, since the radar reflectivity factor difference is
about 19 dB as expected for ideal Bragg scattering. They also reported
that there are many cases in which frequency dependence is not explained by
clear-air Bragg scattering. In such cases, the S-band radar echo is about
10 dB stronger than the X-band radar echo, and the difference is also observed
in the cloud core regions. Rogers and Brown (1997) also reported a similar
frequency dependence of radar returns during their observation of a smoke
plume from intense fire, using a UHF wind profiler and an X-band radar.
Erkelens et al. (2001) analyzed the observation data and estimated the
influence of coherent scattering by cloud droplets using the $-\mathrm{5}/\mathrm{3}$ power law
of turbulent spectrum, which represents turbulent mixing of cloud water with
environmental clear air, i.e., turbulent entrainment. They concluded that
coherent scattering by cloud droplets can be more significant than clear-air
Bragg scattering, whereas turbulent entrainment is not the only factor
relevant to the frequency dependence in the observation data.
Kostinski and Jameson (2000) first pointed out the possibility that
particulate Bragg scattering due to turbulent clustering leads to the
frequency dependence reported by Knight and Miller (1998). To evaluate the
quantitative influence of turbulent clustering on the radar reflectivity
factor, it is crucial to understand the spatial structure of turbulent
clustering. Turbulent clustering has been discussed in much of the literature
because it can enhance the collision growth of cloud droplets
(e.g., Sundaram and Collins, 1997; Reade and Collins, 2000; Ayala et al., 2008a, b; Onishi et al., 2009; Wang et al., 2009; Onishi and Vassilicos, 2014), and
statistical data on turbulent clustering have been obtained for scales
relevant to droplet collisions. However, these data cannot be adopted for
particulate Bragg scattering because the clustering scales relevant to
particulate Bragg scattering are on the microwave wavelength, which is larger
than droplet collision scales. A quantitative estimate of particulate Bragg
scattering due to turbulent clustering is first provided by
Dombrovsky and Zaichik (2010). Their analytical estimate was based on a
clustering model for droplet collision scales but indicated that turbulent
clustering can lead to a considerable increase in the radar reflectivity
factor. Matsuda et al. (2014) clarified the quantitative influence based on
turbulent clustering data obtained by a three-dimensional direct numerical
simulation (DNS), which covered the clustering scales on the microwave
wavelength. They estimated the increase in the radar reflectivity factor due
to turbulent clustering by calculating the power spectrum of number density
fluctuations, *E*_{np}(*k*|*r*_{p}), where *k* is the wave number and
*r*_{p} is the droplet radius. The power spectrum *E*_{np}(*k*|*r*_{p})
is strongly dependent on the droplet size: more specifically, *E*_{np}(*k*|*r*_{p})
is dependent on the Stokes number, *S**t*, which is defined as $St\equiv {\mathit{\tau}}_{\mathrm{p}}/{\mathit{\tau}}_{\mathit{\eta}}$ (*τ*_{p} is the
relaxation time of droplet motion and *τ*_{η} is the Kolmogorov time).
However, the discussion of radar reflectivity factor increases is limited to
cases of monodisperse particles. Thus, the results are not directly
applicable to particulate Bragg scattering for real cloud systems, in which
cloud droplets have broad droplet size distributions.

Therefore, this study aims to investigate the influence of turbulent clustering of polydisperse droplets on particulate Bragg scattering. Firstly, the theoretical formulation of particulate Bragg scattering is extended for polydisperse particles and expressed using the cross spectrum of number density fluctuations for two different droplet sizes. Secondly, the three-dimensional DNS of particle-laden homogeneous isotropic turbulence is performed to obtain turbulent droplet clustering data, which are used to calculate the power spectrum and the cross spectrum of number density fluctuations. A parameterization for the cross spectrum is then proposed considering the dependence of the Stokes number combination, and the influence of gravitational settling is discussed and incorporated. Finally, in order to investigate the impact of turbulent clustering on radar observations of realistic clouds, the proposed model is applied to high-resolution cloud-simulation data obtained by a spectral-bin cloud microphysics simulation.

Here, we aim to formulate the radar reflectivity factor *Z* for a nonuniform
distribution of polydisperse cloud droplets based on the discussion of
Gossard and Strauch (1983), but without assuming monodisperse droplet
sizes. Because the radii of cloud droplets are much smaller than the
microwave wavelength, the electric field vector of the microwaves scattered
by a single droplet, *E*_{sca}(*t*, ** x**,

*r*

_{p}), is given by the Rayleigh-scattering approximation:

where *E*_{inc} is the electric field amplitude vector of the
incident microwave; *r*_{p} is the droplet radius; ** x**,

*x*_{t}, and

*x*_{r}are the position vectors of the droplet, microwave transmitter, and microwave receiver;

*k*_{inc}and

*k*_{sca}are the wave number vectors of the incident and scattered microwaves which satisfy $\left|{\mathit{k}}_{\mathrm{inc}}\right|=\left|{\mathit{k}}_{\mathrm{sca}}\right|={k}_{\mathrm{m}}$, and

*χ*is the angle between

*E*_{inc}and

*k*_{sca}. Considering the droplet-size dependence of

*E*_{sca}(

*t*,

**,**

*x**r*

_{p}), the electric power of the microwave scattered by a group of droplets,

*P*

_{s}, is given by

where *n*(** x**,

*r*

_{p})d

*r*

_{p}d

**is the number of droplets with radii from**

*x**r*

_{p}to ${r}_{{}_{\mathrm{p}}}+\mathrm{d}{r}_{{}_{\mathrm{p}}}$ in an infinitesimal volume d

**at position**

*x***, and**

*x**ζ*is the intrinsic impedance. The overbar denotes the temporal average. The relationship between the radar reflectivity factor

*Z*and the scattering properties of target clouds is given by

where *λ* is the microwave wavelength and *σ* is the radar cross
section (Gossard and Strauch, 1983), which is defined as

where *P*_{o} is the electric power of the incident microwave, which is
given by ${P}_{\mathrm{o}}=|{\mathit{E}}_{\mathrm{inc}}{|}^{\mathrm{2}}/\mathit{\zeta}$. Substitution of
Eqs. (2), (3), and (5)
into Eq. (4) yields

where the wave number vector ** κ** is defined as
$\mathit{\kappa}={\mathit{k}}_{\mathrm{inc}}-{\mathit{k}}_{\mathrm{sca}}$. Note that
radar remote sensing typically uses backward scattering; i.e.,
${\mathit{k}}_{\mathrm{sca}}=-{\mathit{k}}_{\mathrm{inc}}$. Thus,

**=2**

*κ*

*k*_{inc}.

Similarly to Gossard and Strauch (1983), we assume *n*(** x**,

*r*

_{p}) to be composed of the temporal-average and fluctuation parts; i.e.,

*n*(

**, ${r}_{\mathrm{p}})=\stackrel{\mathrm{\u203e}}{n(\mathit{x},{r}_{\mathrm{p}})}+\mathit{\delta}n(\mathit{x},{r}_{\mathrm{p}})$. The temporal-average part, $\stackrel{\mathrm{\u203e}}{n(\mathit{x},{r}_{\mathrm{p}})}$, contributes to the separated reflection; therefore, the contribution of this part is negligibly small when $\stackrel{\mathrm{\u203e}}{n(\mathit{x},{r}_{\mathrm{p}})}$ has no fluctuation on a spatial scale of half the wavelength (Erkelens et al., 2001). Thus, we neglect the contribution of $\stackrel{\mathrm{\u203e}}{n(\mathit{x},{r}_{\mathrm{p}})}$. Then, we obtain**

*x*where the angular brackets represent a temporal and spatial average in the measurement volume.

In order to decompose the spatial correlation function $\langle \mathit{\delta}n(\mathit{x},{r}_{\mathrm{p}})\mathit{\delta}n(\mathit{x}+\mathit{r},{r}_{\mathrm{p}}^{\prime})\rangle $, we introduce the
probability density function (PDF) of droplet radius *r*_{p} to the
measurement volume, *q*_{r}(*r*_{p}), and the number density
distribution function for monodisperse droplets with a radius of *r*_{p},
*n*_{p}(** x**|

*r*

_{p}). The PDF is defined as ${q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)\equiv \frac{\mathrm{1}}{{N}_{\mathrm{p}}}\underset{\mathit{x}\in V}{\int}\stackrel{\mathrm{\u203e}}{n(\mathit{x},{r}_{\mathrm{p}})}\mathrm{d}\mathit{x}$, where

*N*

_{p}is the total number of droplets in the measurement volume; i.e., ${N}_{\mathrm{p}}\equiv \underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}\underset{\mathit{x}\in V}{\int}\stackrel{\mathrm{\u203e}}{n(\mathit{x},{r}_{\mathrm{p}})}\mathrm{d}\mathit{x}\mathrm{d}{r}_{\mathrm{p}}$. The PDF satisfies $\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)\mathrm{d}{r}_{\mathrm{p}}=\mathrm{1}$. The number density distribution function for monodisperse droplets is then defined as ${n}_{\mathrm{p}}\left(\mathit{x}\right|{r}_{\mathrm{p}})\equiv n(\mathit{x},{r}_{\mathrm{p}})/{q}_{\mathrm{r}}({r}_{\mathrm{p}})$ so that

*n*(

**,**

*x**r*

_{p}) is given by $n(\mathit{x},{r}_{\mathrm{p}})={n}_{\mathrm{p}}\left(\mathit{x}\right|{r}_{\mathrm{p}}\left){q}_{\mathrm{r}}\right({r}_{\mathrm{p}})$. The number density distribution function

*n*

_{p}(

**|**

*x**r*

_{p}) satisfies $\underset{\mathit{x}\in V}{\int}\stackrel{\mathrm{\u203e}}{{n}_{\mathrm{p}}\left(\mathit{x}\right|{r}_{\mathrm{p}})}\mathrm{d}\mathit{x}={N}_{\mathrm{p}}$ for arbitrary

*r*

_{p}. Note that the spatial correlation function $\langle \mathit{\delta}n(\mathit{x},{r}_{\mathrm{p}})\mathit{\delta}n(\mathit{x}+\mathit{r},{r}_{\mathrm{p}}^{\prime})\rangle $ for ${r}_{\mathrm{p}}^{\prime}={r}_{\mathrm{p}}$ is discontinuous at

**=**

*r***0**because the droplet distribution is composed of spatially discrete points. The singularity is given by $\langle n(\mathit{x},{r}_{\mathrm{p}})\rangle \mathit{\delta}\left(\mathit{r}\right)\mathit{\delta}({r}_{\mathrm{p}}^{\prime}-{r}_{\mathrm{p}})$, where

*δ*(

**) and $\mathit{\delta}({r}_{\mathrm{p}}^{\prime}-{r}_{\mathrm{p}})$ are the Dirac delta functions. Thus, the spatial correlation function is given by**

*r*
where 〈*n*_{p}〉 is the average number density ($\langle {n}_{\mathrm{p}}\rangle \equiv {N}_{\mathrm{p}}/V$)
and $\mathrm{\Psi}\left(\mathit{r}\right|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime})$ is defined as the continuous part of
$\langle \mathit{\delta}{n}_{\mathrm{p}}\left(\mathit{x}\right|{r}_{\mathrm{p}}\left)\mathit{\delta}{n}_{\mathrm{p}}\right(\mathit{x}+\mathit{r}\left|{r}_{\mathrm{p}}^{\prime}\right)\rangle $.
Substitution of Eq. (8) into Eq. (7) and
adoption of the isotropic clustering assumption Gossard and Strauch (1983) yield

where *κ* is $\mathit{\kappa}=\left|\mathit{\kappa}\right|$, $\langle {r}_{\mathrm{p}}^{\mathrm{6}}\rangle $
is given by $\langle {r}_{\mathrm{p}}^{\mathrm{6}}\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{r}_{\mathrm{p}}^{\mathrm{6}}{q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)\mathrm{d}{r}_{\mathrm{p}}$,
and *E*_{r3np}(*k*) is the ${r}_{\mathrm{p}}^{\mathrm{3}}$-weighted power spectrum,
defined as

where ${C}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime})$ is the cross spectrum of number
density fluctuations for *n*_{p}(** x**|

*r*

_{p}) and ${n}_{\mathrm{p}}\left(\mathit{x}\right|{r}_{\mathrm{p}}^{\prime})$: the cross spectrum ${C}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime})$ is defined as ${C}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime})\equiv \underset{\left|\mathit{k}\right|=k}{\int}\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}(\mathit{k}|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime})\mathrm{d}{\mathit{\sigma}}_{k}$; i.e., the integration of $\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}\left(\mathit{k}\right|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime})$ over the spherical shell,

*σ*

_{k}, at $\left|\mathit{k}\right|=k$, in which $\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}\left(\mathit{k}\right|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime})$ is the cross spectral density function, defined as

The first and second terms on the right-hand side of Eq. (9)
are the incoherent and coherent scattering parts; particulate
Bragg scattering is caused by the second term. Equations (9)
and (10) imply that the particulate Bragg-scattering part of *Z* for
an arbitrary droplet size distribution can be calculated using the cross
spectrum for bidisperse droplet size conditions. When droplets are
distributed randomly and uniformly, the second term equals zero. Thus, the
radar reflectivity factor when assuming a random and uniform droplet
distribution is given by the first term; i.e., ${Z}_{\mathrm{incoh}}={\mathrm{2}}^{\mathrm{6}}\langle {r}_{\mathrm{p}}^{\mathrm{6}}\rangle \langle {n}_{\mathrm{p}}\rangle $.

It should be noted that Eq. (9) satisfies the theoretical
solution for particulate Bragg scattering of monodisperse droplets: for the
case of monodisperse droplets with radii of *r*_{p1}, the PDF of droplet
radius is given by ${q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)=\mathit{\delta}({r}_{\mathrm{p}}-{r}_{\mathrm{p}\mathrm{1}})$. Then,
the radar reflectivity factor *Z* is given by

where *E*_{np}(*k*|*r*_{p}) is the power spectrum of number density
fluctuations, which satisfies ${E}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}})={C}_{\mathrm{np}}(k|{r}_{\mathrm{p}},{r}_{\mathrm{p}})$.
Note that *E*_{np}(*k*|*r*_{p}) is defined as ${E}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}})\equiv \underset{\left|\mathit{k}\right|=k}{\int}\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}(\mathit{k}\left|{r}_{\mathrm{p}}\right)\mathrm{d}{\mathit{\sigma}}_{k}$,
where $\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}\left(\mathit{k}\right|{r}_{\mathrm{p}})$ is the power
spectral density function of *n*_{p}(** x**|

*r*

_{p}), defined as

## 3.1 Direct numerical simulation

In order to obtain turbulent clustering data for calculating the cross spectrum, we have performed a three-dimensional DNS for particle-laden homogeneous isotropic turbulence. Three-dimensional incompressible turbulent airflows were calculated by solving the continuity and Navier–Stokes equations:

where *u*_{i} is the flow velocity in the *i*th direction, *p* is the pressure,
*ρ*_{a} is the air density, *ν* is the kinematic viscosity, and
*F*_{i} is the external forcing term. The nonlinear term was discretized by the
fourth-order central difference scheme (Morinishi et al., 1998). The time
integration was calculated by the second-order Runge–Kutta scheme. The
HSMAC method (Hirt and Cook, 1972) was used for velocity-pressure coupling. The
external forcing was applied to maintain the intensity of large-scale eddies
for wave numbers ** k** in the range $\left|\mathit{k}{L}_{\mathrm{0}}\right|<\mathrm{2}$
(Onishi et al., 2011), where

*L*

_{0}is the representative length scale.

Droplet motions were simulated by Lagrangian point-particle tracking. Here,
we assumed that the droplet density *ρ*_{p} is much larger than *ρ*_{a}
and the drag term is given based on the Stokes law. The droplet motions were tracked by

where *v*_{i} and *g*_{i} are the particle velocity and gravitational
acceleration in the *i*th direction. *τ*_{p} is the
droplet relaxation time, which is given by

The effects of turbulent modulation and droplet collision were neglected for
simplicity because these effects were typically small in the timescale of *τ*_{η} in clouds.

The computational domain was set as a cube with edge lengths of 2*π**L*_{0}.
Periodic boundary conditions were applied in all three directions. A uniform
staggered grid was used for discretization. The number of grid points was set
to 512^{3}. A Taylor microscale-based turbulent Reynolds number of the
obtained flow was *R**e*_{λ}=204, where *R**e*_{λ} is
defined as $R{e}_{\mathit{\lambda}}\equiv {l}_{\mathit{\lambda}}{u}^{\prime}/\mathit{\nu}$, where
*l*_{λ} is the Taylor microscale and *u*^{′} is the rms value of the velocity
fluctuation. Note that this value of *R**e*_{λ} is sufficiently
large to obtain turbulent clustering data for high Reynolds number turbulence
in the wave number range relevant to radar observations (Matsuda et al., 2014)
(see Sect. 3.2). The kinematic viscosity, *ν*, was set
to $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m^{2} s^{−1}, and the ratio of the droplet density to the
air density, *ρ*_{p}∕*ρ*_{a}, was set to 840, assuming 1 atm and
298 K. The total number of droplets, *N*_{p}, was set to 1.5×10^{7}.

For this study, we performed the DNS for monodisperse and polydisperse
droplets. Table 1 shows the computational settings for
turbulence, droplet size, and gravitational acceleration. For monodisperse
droplets, the droplet motions in an identical turbulent flow field were
calculated for six values of Stokes number, *S**t*. The clustering data
for the monodisperse cases were used for calculating the cross spectrum of
number density fluctuations for any combinations of these *S**t* values. For
polydisperse droplets, a typical droplet size distribution for maritime
cumulus clouds (the size distribution data named “CUMA” in
Hess et al., 1998) was applied, and the droplets were tracked in turbulent
flows using three different energy dissipation rates *ϵ*.
*ϵ* values of the obtained turbulent flows were approximately 100, 400,
and 1000 cm^{2} s^{−3}, which can be observed in cumulus and cumulonimbus clouds
(Pinsky et al., 2008). The data for the polydisperse droplet cases were used
to discuss the reliability of the proposed cross spectrum model. It should be
noted that the droplet size distribution for the polydisperse cases were
identical but the Stokes number histograms were different; the Stokes numbers
corresponding to the modal radius (10.4 *μ*m) for *ϵ*=100, 400,
and 1000 cm^{2} s^{−3} were 0.035, 0.069, and 0.10, respectively. The
gravitational acceleration $g\equiv \sqrt{{g}_{i}{g}_{i}}$ was set to zero for the
monodisperse cases. The DNS for polydisperse droplets were performed under
the conditions with and without gravitational settling. The Froude numbers,
*F**r* ($\equiv {a}_{\mathit{\eta}}/g$, in which ${a}_{\mathit{\eta}}\equiv {\mathit{\u03f5}}^{\mathrm{3}/\mathrm{4}}{\mathit{\nu}}^{-\mathrm{1}/\mathrm{4}}$
is the Kolmogorov acceleration), for the cases with gravitational
settling were 0.0520, 0.145, and 0.289 for *ϵ*=100, 400, and
1000 cm^{2} s^{−3}, respectively. The influence of gravitational settling on
turbulent particle clustering is often discussed using the settling
parameter, *S*_{v}, which is defined as ${S}_{\mathrm{v}}\equiv {v}_{\mathrm{T}}/{u}_{\mathit{\eta}}$
(e.g., Wang and Maxey, 1993; Grabowski and Vaillancourt, 1999; Bec et al., 2014; Ireland et al., 2016),
where *v*_{T}=*τ*_{p}*g* is the terminal settling velocity and
${u}_{\mathit{\eta}}\equiv (\mathit{\nu}\mathit{\u03f5}{)}^{\mathrm{1}/\mathrm{4}}$ is the Kolmogorov velocity. Note that
*S*_{v} satisfies ${S}_{\mathrm{v}}=St/Fr$. The settling parameters
corresponding to the modal radius for *ϵ*=100, 400, and 1000 cm^{2} s^{−3}
were 0.67, 0.48, and 0.38, respectively.

## 3.2 Computation of power spectrum and cross spectrum

The power spectral density function $\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}\left(\mathit{k}\right|{r}_{\mathrm{p}})$ and the cross spectral density function $\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}\left(\mathit{k}\right|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}})$ are calculated from the Lagrangian droplet distribution data as follows:

where $\stackrel{\mathrm{\u0303}}{{n}_{\mathrm{p}}}\left(\mathit{k}\right|{r}_{\mathrm{p}})$ is the Fourier component of
the droplet number density distribution, *n*_{p}(** x**|

*r*

_{p}), and the angle brackets denote an ensemble average. The number density distribution for Lagrangian discrete droplets is given by

where *x*_{p,j} is the position vector of the *j*th droplet with
radius *r*_{p}, and *N*_{p} is the total number of droplets with
radius *r*_{p}. The Fourier components of Eq. (20) are
then given by

Substitution of Eq. (21) into Eq. (18) yields

Similarly, substitution of Eq. (21) into Eq. (19) yields

where 〈*n*_{p1}〉 and 〈*n*_{p2}〉 are the
average number density of droplets with radii of *r*_{p1} and
*r*_{p2} (i.e., $\langle {n}_{\mathrm{p}\mathrm{1}}\rangle \equiv \langle {n}_{\mathrm{p}}\left(\mathit{x}\right|{r}_{\mathrm{p}\mathrm{1}})\rangle $
and $\langle {n}_{\mathrm{p}\mathrm{2}}\rangle \equiv \langle {n}_{\mathrm{p}}\left(\mathit{x}\right|{r}_{\mathrm{p}\mathrm{2}})\rangle $),
respectively, *N*_{p1} and *N*_{p2} are the numbers
of droplets with radii of *r*_{p1} and *r*_{p2},
respectively, *x*_{p1,j} is the position of the *j*th droplet
with a radius of *r*_{p1}, and ${\mathit{x}}_{\mathrm{p}\mathrm{2},{j}^{\prime}}$ is the position of
the *j*^{′}th droplet with a radius of *r*_{p2}. It should be noted that the
imaginary part of $\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}\left(\mathit{k}\right|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}})$ is
omitted in Eq. (25) because, statistically, it should be
zero. We confirmed that the imaginary part of
${C}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}})$ calculated from the DNS data is
*O*(10^{−4}), which is caused by the statistical and truncation errors.

The spectral density functions, $\stackrel{\mathrm{\u0303}}{\mathrm{\Phi}}\left(\mathit{k}\right|{r}_{\mathrm{p}})$ and
$\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}\left(\mathit{k}\right|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}})$, were calculated for
discrete wave numbers $\mathit{k}{L}_{\mathrm{0}}=({h}_{\mathrm{1}},{h}_{\mathrm{2}},{h}_{\mathrm{3}})$, where *h*_{1}, *h*_{2}, and
*h*_{3} are arbitrary integers that satisfy $k-\mathrm{\Delta}k/\mathrm{2}\le \left|\mathit{k}\right|<k+\mathrm{\Delta}k/\mathrm{2}$,
where Δ*k* was set to 1∕*L*_{0}. *E*_{np}(*k*|*r*_{p})
and ${C}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}})$ were then obtained by the following equations:

The spectra, *E*_{np}(*k*|*r*_{p}) and
${C}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}})$, were calculated for 19 representative wave numbers in the same way as
Matsuda et al. (2014). The ensemble average in Eq. (23)
was obtained by averaging 10 temporal slices of the droplet distributions,
which were sampled for intervals of ${T}_{\mathrm{0}}={L}_{\mathrm{0}}/{U}_{\mathrm{0}}$. The ensemble average in
Eq. (25) was also obtained for 10 pairs of temporal
slices. Each pair was composed of temporal slices at the same time step for
different *S**t* cases.

Matsuda et al. (2014) normalized the power spectrum *E*_{np}(*k*|*r*_{p})
by using 〈*n*_{p}〉 and the Kolmogorov scale, *l*_{η},
defined as ${l}_{\mathit{\eta}}\equiv {\mathit{\nu}}^{\mathrm{3}/\mathrm{4}}{\mathit{\u03f5}}^{-\mathrm{1}/\mathrm{4}}$, and confirmed that,
for *R**e*_{λ}≥204, the *R**e*_{λ} dependence of the
normalized power spectrum is negligible at the wave number range relevant to
radar observations ($\mathrm{0.05}<k{l}_{\mathit{\eta}}<\mathrm{4.0}$). Thus, this study used the same
normalization for *E*_{np}(*k*|*r*_{p}) and
${C}_{\mathrm{np}}\left(k\right|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}})$ as follows:

where *ξ* is the normalized wave number defined as *ξ*≡*k**l*_{η}
and *S**t*, *S**t*_{1}, and *S**t*_{2} are the Stokes numbers of
water droplets with radii of *r*_{p}, *r*_{p1}, and *r*_{p2}.

## 4.1 Spatial droplet distribution and cross spectrum

Figure 1 shows the spatial distribution of droplets for
*S**t*=0.2, 0.5, and 1.0. Droplets located in the range of
$\mathrm{0}<z<\mathrm{4}{l}_{\mathit{\eta}}$ are indicated. The droplet position data were sampled at the
same time step; i.e., the background turbulent flow field is identical.
Figure 1a is the overall view of the region with a size
of 2*π**L*_{0}×2*π**L*_{0}, while Fig. 1b is the
magnified view for the region with a size of 0.5*π**L*_{0}×0.5*π**L*_{0}.
The clusters and void areas are clearly observed for all *S**t* cases.
Their location for *S**t*=0.5 is almost the same as that for the other *S**t*.
This is because the locations of clusters and void areas for
small *S**t* are strongly dependent on the instantaneous turbulent flow
field. However, the small-scale structure of clusters is not exactly the
same; i.e., the droplets become more concentrated in clusters as *S**t* increases.
Figure 2 shows the power spectra ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi}\right|St)$
and the cross spectra ${C}_{\mathrm{np}}^{*}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})$
obtained from the DNS data. Note that the high wave number portion is
omitted when the value of the cross spectrum is smaller than the
computational error level (10^{−4}). The power spectra ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi}\right|St)$
show power-law-type slopes at wave numbers smaller and
larger than the peak location. The peak height and slope of the spectra are
strongly dependent on the Stokes number; this Stokes number dependence was
discussed by Matsuda et al. (2014). In the small wave number region, the cross
spectra ${C}_{\mathrm{np}}^{*}$ also show power-law-type slopes. In this region, the
curve of ${C}_{\mathrm{np}}^{*}$ for *S**t*_{1}=0.5 and *S**t*_{2}=0.2 is located
between the power spectra ${E}_{\mathrm{np}}^{*}$ for *S**t*=0.5 and
*S**t*=0.2. Similarly, the curve of ${C}_{\mathrm{np}}^{*}$ for *S**t*_{1}=0.5 and
*S**t*_{2}=1.0 is located between the power spectra ${E}_{\mathrm{np}}^{*}$ for
*S**t*=0.5 and *S**t*=1.0. On the other hand, in the large wave number
region, both cross spectra ${C}_{\mathrm{np}}^{*}$ become smaller than the power
spectra ${E}_{\mathrm{np}}^{*}$ without showing power-law-type slopes. These trends
imply that the cross spectrum is influenced by not only the Stokes number
dependence of the clustering intensity, but also the spatial correlation of
clusters between different values of *S**t*. In order to focus on the
influence of the spatial correlation of clusters, we have evaluated the
coherence $\mathrm{coh}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})$, which is defined as

Figure 3 shows the coherence between *S**t*_{1}=0.5 and
other Stokes numbers. The coherence $\mathrm{coh}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})$ is close
to unity in the small wave number region and decreases to zero as the
wave number increases. These trends correspond to the spatial correlation
between cluster locations for different *S**t* cases in
Fig. 1. Figure 3b shows that the
coherence decreases with an almost constant slope in the vertically
logarithmic and horizontally linear chart. The slope of the coherence is
dependent on the combination of *S**t*_{1} and *S**t*_{2}; the slope
becomes steeper as the difference in *S**t* increases. These results
indicate that the decreasing trend of coherence can be approximated by an
exponential function; i.e.,

where *ξ*_{c} is the critical wave number normalized by the Kolmogorov
scale, given by a function of *S**t*_{1} and *S**t*_{2}. In this study,
the critical wave number was obtained by finding the best-fitting exponential
curve to the coherence for each combination of Stokes numbers.
Figure 4 shows the critical wave numbers *ξ*_{c} for
all combinations of *S**t*_{1} and *S**t*_{2}, where *S**t*_{2}>*S**t*_{1}. *ξ*_{c} for
the same *S**t*_{1} increases as *S**t*_{2} decreases, whereas
*ξ*_{c} for the same *S**t*_{2} increases as *S**t*_{1} increases.
This indicates that *ξ*_{c} increases as
*S**t*_{1} and *S**t*_{2} becomes closer to each other. It should be noted that several
studies discuss the spatial correlation of bidisperse clustering particles.
For example, Zhou et al. (2001) developed the radial distribution function (RDF)
model at the separation length of the Kolmogorov scale. They reported
that the correlation coefficient of the bidisperse RDF obtained by their DNS
is explained well by the ratio of two Stokes numbers. Chun et al. (2005) also
discussed the bidisperse RDF of clustering particles. The result of their
perturbation expansion analysis indicated that the bidisperse RDF becomes
constant at separation lengths sufficiently smaller than the crossover
length, *l*_{c}, which is proportional to the Stokes number difference.
As the cross spectrum of number density fluctuations is a Fourier transform
of the bidisperse RDF, the Stokes number ratio, *S**t*_{2}∕*S**t*_{1}, and
the Stokes number difference, *S**t*_{2}–*S**t*_{1}, are candidates for the
dominant parameter for *ξ*_{c}. Figure 5 shows
*ξ*_{c} plots against *S**t*_{2}∕*S**t*_{1} and *S**t*_{2}–*S**t*_{1}. These
figures clearly indicate that the Stokes number dependence of *ξ*_{c}
is explained better by *S**t*_{2}–*S**t*_{1} than *S**t*_{2}∕*S**t*_{1}.
This would be because both the critical wave number *ξ*_{c} and
crossover length *l*_{c} represent the critical scale of the spatial
correlation between the clusters for different Stokes numbers; i.e., the
cluster locations are less correlated on a scale smaller than the critical
scale. Based on this insight, we estimate that the critical wave number *ξ*_{c}
is inversely proportional to the crossover length *l*_{c}.
Figure 6 shows *ξ*_{c} against the inverse of the
Stokes number difference, which is generally expressed as $\mathrm{1}/|S{t}_{\mathrm{1}}-S{t}_{\mathrm{2}}|$.
This figure confirms that *ξ*_{c} is approximately
proportional to $\mathrm{1}/|S{t}_{\mathrm{1}}-S{t}_{\mathrm{2}}|$; the least square fitting gives

This implies that *ξ*_{c} is closely related to the inverse
of *l*_{c}∕*l*_{η} because the crossover length of Chun et al. (2005) is
approximately ${l}_{\mathrm{c}}/{l}_{\mathit{\eta}}\approx \mathrm{5.0}|S{t}_{\mathrm{1}}-S{t}_{\mathrm{2}}|$ based
on their DNS data. It should be noted that the analytical results of
Chun et al. (2005) are valid for *S**t*≪1. Thus, the deviation of *ξ*_{c}
from the linear curve is due to the higher-order response of
particle motions to the turbulent flow. However, Fig. 6
confirms that the Stokes number difference is the dominant parameter
for *ξ*_{c}, at least for *S**t*≤2.0.

## 4.2 Modeling the influence of polydisperse clustering droplets on radar reflectivity factor

According to the above discussion, we can estimate the increase in the radar
reflectivity factor due to turbulent clustering of polydisperse droplets in
Eq. (9), provided that *q*_{r}(*r*_{p}) and
〈*n*_{p}〉 are given. That is, the normalized cross spectrum
${C}_{\mathrm{np}}^{*}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})$ is given by

Here, we assume that the cross spectrum is a positive real number. The coherence is estimated using Eqs. (31) and (32). The parameterization for ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi}\right|St)$ was proposed by Matsuda et al. (2014). The model equation is given by

where *c*_{1}, *c*_{2}, *α*, *β*, and *γ* are the model parameters
given by the functions of *S**t* as follows:

Matsuda et al. (2014) confirmed that this parameterization is reliable for
*S**t*≤2.0: in this range, the error of the parameterization is smaller than 1 dB.

In order to evaluate the reliability of the proposed cross spectrum model,
the ${r}_{\mathrm{p}}^{\mathrm{3}}$-weighted power spectrum *E*_{r3np}(*k*) for the droplet
size distribution of CUMA has been estimated using the proposed cross
spectrum model and compared with the spectrum obtained from the DNS data for
the cases of CUMA_eps100, CUMA_eps400, and CUMA_eps1000. Figure 7a
shows the ${r}_{\mathrm{p}}^{\mathrm{3}}$-weighted power spectrum, which is normalized as

where $\langle {r}_{\mathrm{p}}^{\mathrm{3}}\rangle $ is given by $\langle {r}_{\mathrm{p}}^{\mathrm{3}}\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{r}_{\mathrm{p}}^{\mathrm{3}}{q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)\mathrm{d}{r}_{\mathrm{p}}$. The
dashed lines show ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi}\right)$ predicted using the proposed
parameterization including Eq. (31), while the dashed-dotted
lines show those predicted by assuming perfect coherence; i.e., $\mathrm{coh}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})=\mathrm{1}$.
The parameterization with $\mathrm{coh}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})=\mathrm{1}$ overestimates
${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi}\right)$ at large wave numbers and the
difference becomes larger as *ϵ* becomes larger. This indicates that
the influence of the weak spatial correlation of cluster locations between
different Stokes numbers is not negligible for predicting the spectrum ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi}\right)$
for large wave numbers, and the assumption of
$\mathrm{coh}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})=\mathrm{1}$ can be applied only for
predicting ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi}\right)$ for small wave numbers ($\mathit{\xi}<O\left({\mathrm{10}}^{-\mathrm{1}}\right)$). On the other hand,
${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi}\right)$ values predicted by the parameterization using
Eq. (31) show good agreement with those obtained by the DNS
data for overall wave numbers. The error level of the parameterization using
Eq. (31) is evaluated by the rms error *e*_{rms} in units
of decibels. *e*_{rms} is defined as

where *ξ*^{′} is defined as ${\mathit{\xi}}^{\prime}=\mathrm{ln}\mathit{\xi}$ and superscript dB denotes a value
in units of decibels. *e*_{rms} was calculated for the wave number range
relevant to radar observations; i.e., $\mathrm{0.05}\le \mathit{\xi}\le \mathrm{4.0}$. *e*_{rms}
for the cases of CUMA_eps100, CUMA_eps400, and CUMA_eps1000 are 1.41,
0.152, and 0.251 dB, respectively. Because the error level of 1 dB is
unavoidable for radar observations (Bringi et al., 1990; Carey et al., 2000),
*e*_{rms} values for CUMA_eps400 and CUMA_eps1000 are negligibly small.
*e*_{rms} for CUMA_eps100 is slightly larger than the threshold, but this is
caused by the error of calculating the reference spectrum based on the DNS
data at *ξ*>2. We confirm that, for CUMA_eps100, *e*_{rms} evaluated
at the range of $\mathrm{0.05}\le \mathit{\xi}\le \mathrm{2.0}$ is smaller than 1 dB. Thus, the Stokes
number dependence of the cross spectrum in the absence of gravity is
appropriately modeled to predict the influence of turbulent clustering to a
sufficient accuracy.

## 4.3 Influence of gravitational settling on cross spectrum

The parameterization summarized in the previous subsection was obtained under
the condition without gravitational droplet settling. The settling influence for the monodispersed cases was discussed by
Matsuda et al. (2014, 2017). The large gravitational
settling can modulate ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi}\right|St)$, and that can be a cause
of the significant difference in the radar reflectivity factor. However, the
influence on ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi}\right|St)$ is insignificant for *S*_{v}<3
(Matsuda et al., 2014).

For the cases of polydisperse particles, the settling influence on the coherence term must be considered as well as the influence on ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi}\right|St)$ in Eq. (33). Ayala et al. (2008a) and Lu et al. (2010) reported that gravitational settling enlarges the crossover length of the bidisperse RDF. Lu et al. (2010) extended the perturbation expansion analysis of Chun et al. (2005) and presented the formulation for the crossover length in the presence of gravity, which is

where *C*_{Chun} is the coefficient derived by Chun et al. (2005),
*a*_{0} is the ratio of the acceleration variance 〈*a*^{2}〉 to the square of
the Kolmogorov acceleration ${a}_{\mathit{\eta}}^{\mathrm{2}}$, *τ*_{a} is the acceleration
correlation timescale, and *τ*_{g} is the correlation timescale for
gravitational settling particles. Chun et al. (2005) obtained the values of
*C*_{Chun}≈5.0 and *a*_{0}≈1.545 based on their DNS
results. Lu et al. (2010) further assumed *S*_{v}*≲*1 so that
*τ*_{g} is simply given by ${\mathit{\tau}}_{\mathrm{g}}={\mathit{\tau}}_{\mathrm{a}}=\mathrm{1.5}{\mathit{\tau}}_{\mathit{\eta}}$.
The crossover length *l*_{c} of Eq. (38) becomes
equivalent to that of Chun et al. (2005) when gravitational settling is
negligibly small, whereas the gravity effect is dominant for the cases of
*F**r*<0.47, which often appears in cloud turbulence. Thus, in this
study, we modify Eq. (32) to include the settling influence on
the coherence $\mathrm{coh}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})$. Since *ξ*_{c} is
inversely proportional to *l*_{c}∕*l*_{η}, we propose the following
correction based on Eq. (38):

Note that this study also adopts the value of *a*_{0} obtained by Chun et al. (2005).

The reliability of the modified parameterization for the case with
gravitational settling has been evaluated in the same way as the previous
subsection. Figure 7b shows ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi}\right)$ for the
case with gravitational settling. ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi}\right)$ values at large
wave numbers are smaller than those for the case without gravitational
settling, and the difference from the parameterization with $\mathrm{coh}\left(\mathit{\xi}\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})=\mathrm{1}$
is larger, indicating that the coherence model is more
important than the case without gravitational settling. It is also confirmed
that ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi}\right)$ values predicted by the proposed parameterization
show good agreement with those of the DNS results for the case with
gravitational settling. *e*_{rms} evaluated at the range of $\mathrm{0.05}\le \mathit{\xi}\le \mathrm{4.0}$
are 0.93 and 0.31 dB for the cases of CUMA_eps400 and
CUMA_eps1000, respectively, and that for CUMA_eps100 at the range of
$\mathrm{0.05}\le \mathit{\xi}\le \mathrm{2.0}$ is 0.26 dB. Thus, *e*_{rms} remains smaller than
1 dB, even for the case with gravitational settling. These results indicate that
the proposed parameterization can predict the influence of turbulent
clustering for polydisperse droplets considering the gravity effect to a
sufficient accuracy. For the CUMA cases in Fig. 7b,
*S*_{v} for the modal radius is smaller than unity. For the cases of
*S*_{v}>*O*(1), the proposed parameterization would become less reliable.
To improve the parameterization, it is necessary to consider an anisotropic
clustering structure of settling inertial particles. Inertial particles with
a large settling velocity form anisotropic clusters, which are vertically
elongated and horizontally confined (Bec et al., 2014; Ireland et al., 2016; Matsuda et al., 2017). When the clustering structure
is anisotropic, the influence on the radar reflectivity factor theoretically
depends on the direction of microwave propagation.

## 5.1 Cloud simulation data

We have applied the proposed model to the high-resolution cloud-simulation
data of Onishi and Takahashi (2012) to investigate the influence of
turbulent clustering on radar observations. They used the Multi-Scale
Simulator for the Geoenvironment (MSSG), which is a multi-scale
atmosphere–ocean coupled model developed by the Japan Agency for Marine-Earth
Science and Technology. The atmospheric component of MSSG (MSSG-A) solves
nonhydrostatic equations and predicts three wind components, air density,
and pressure, as well as water substance. Finite-difference schemes are used
for calculating spatial derivatives. Turbulent diffusion is calculated using
the static Smagorinsky model. Onishi and Takahashi (2012) used a
spectral-bin scheme for liquid water to explicitly account for the droplet
size distributions. The spectral bin scheme predicts the mass distribution
function *G*(*y*), which is given by

where *y*=ln *r*_{p}, and *m*(*r*_{p}) is the mass of droplets with a radius of *r*_{p}.
The mass coordinate *m* and logarithmic coordinate *y* are discretized as

where $\mathrm{d}y=\mathrm{ln}\mathrm{2}/\left(\mathrm{3}s\right)$, and *s* is a constant; *s*=1 was used. The number
of bins was 33. The representative radius of the first bin, *r*_{p1}, was
3 µm; thus, the representative radius of the 33rd bin (the largest
droplet class) was *r*_{p33}=4.9 mm. The prognostic variable for liquid
water is the water mass content, *M*_{k}, which is defined as
${M}_{k}=\underset{{y}_{k-\mathrm{1}/\mathrm{2}}}{\overset{{y}_{k+\mathrm{1}/\mathrm{2}}}{\int}}G\left(y\right)\mathrm{d}y$; i.e., 33 transport equations for
*M*_{k} were solved in this simulation. The activation process of cloud
condensation nuclei (CCN) was considered based on Twomey's relationship
between the number of activated CCN and the saturation ratio
(Twomey, 1959). The activated droplets were added to the bins using
the prescribed spectrum method (Soong, 1974). Details of the model
configuration are described in Onishi and Takahashi (2012). The model
settings and computational conditions were based on the protocol of the RICO
model intercomparison project (van Zanten et al., 2006, http://www.knmi.nl/samenw/rico/, last access: 6 February 2019). The protocol is based
on the rain in cumulus over the ocean (RICO) field campaign. The domain size
is $\mathrm{12.8}\times \mathrm{12.8}\times \mathrm{4.0}$ km. The resolution setting of the original
RICO protocol is 128×128 points in the horizontal directions and
100 points in the vertical direction; i.e., ${\mathrm{\Delta}}_{x}={\mathrm{\Delta}}_{y}=\mathrm{100}$ m
and Δ_{z}=40 m. Onishi and Takahashi (2012) performed the cloud
simulation for 24 h using the original resolution setting, then continued it
for an additional hour using a higher-resolution setting, generating
512×512 points in the horizontal directions and 200 points in the vertical
direction, giving grid spacings of ${\mathrm{\Delta}}_{x}={\mathrm{\Delta}}_{y}=\mathrm{25}$ m and
Δ_{z}=20 m. This study used the temporal slice of cloud simulation
data at a higher resolution.

## 5.2 Computational method for radar reflectivity factor

The radar reflectivity factor *Z*, including the influence of particulate and
clear-air Bragg scatterings, is given by

where *Z*_{PB} is the particulate Bragg-scattering part of *Z* in
Eq. (9), and *Z*_{CB} is an additional term reflecting
clear-air Bragg scattering. In real clouds, *Z*_{PB} is caused by droplet
number density fluctuations due to turbulent droplet clustering and turbulent
entrainment of environmental clear air. *E*_{r3np}(*k*) for these factors
is given by ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}\left(k\right)={E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{clust}}\left(k\right)+{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{entr}}\left(k\right)$
(Matsuda et al., 2014), where ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{clust}}\left(k\right)$ and
${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{entr}}\left(k\right)$ are the power spectra for turbulent clustering
and entrainment. Note that the correlation term between the cloud water
fluctuations due to these factors is negligibly small because the scales of
the clustering and entrainment sources are typically separated. Thus,
*Z*_{PB} is also given by the linear combination; i.e.,
${Z}_{\mathrm{PB}}={Z}_{\mathrm{PBc}}+{Z}_{\mathrm{PBe}}$, where
${Z}_{\mathrm{PBc}}={\mathrm{2}}^{\mathrm{7}}{\mathit{\pi}}^{\mathrm{2}}{\mathit{\kappa}}^{-\mathrm{2}}{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{clust}}\left(\mathit{\kappa}\right)$
and ${Z}_{\mathrm{PBe}}={\mathrm{2}}^{\mathrm{7}}{\mathit{\pi}}^{\mathrm{2}}{\mathit{\kappa}}^{-\mathrm{2}}{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{entr}}\left(\mathit{\kappa}\right)$.

The spectrum ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{clust}}\left(k\right)$ was calculated using
Eq. (10) and the parameterization proposed in the previous
section. To determine the Stokes number for each droplet size, the energy
dissipation rate *ϵ* of the cloud simulation data was calculated based
on the Smagorinsky model; i.e.,

where *C*_{s} is the Smagorinsky constant (*C*_{s}=0.173 in this
study), Δ_{s} the representative grid spacing given by
${\mathrm{\Delta}}_{\mathrm{s}}=({\mathrm{\Delta}}_{x}{\mathrm{\Delta}}_{y}{\mathrm{\Delta}}_{z}{)}^{\mathrm{1}/\mathrm{3}}$, and
$\widehat{{s}_{ij}}$ the strain rate tensor, which is given by
$\widehat{{s}_{ij}}=\frac{\mathrm{1}}{\mathrm{2}}(\frac{\partial \widehat{{u}_{i}}}{\partial {x}_{j}}+\frac{\partial \widehat{{u}_{j}}}{\partial {x}_{i}})$,
where $\widehat{{u}_{i}}$ is the air velocity in the resolved scale.

The spectrum ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{entr}}\left(k\right)$ was calculated using the well-known
scalar concentration spectrum. Erkelens et al. (2001) estimated this
contribution based on the $-\mathrm{5}/\mathrm{3}$ power law in the inertial-convective range
of the spectrum. In this study, the −1 power law in the viscous-convective
range (*k**l*_{η}<0.1) is also considered, since the diffusive
coefficient *D*_{np} for droplet number density is much smaller than *ν*.
${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{entr}}\left(k\right)$ is approximately given by (Hill, 1978)

where *χ*_{r3np} is the scalar dissipation rate for $\langle {r}_{\mathrm{p}}^{\mathrm{3}}\rangle {n}_{\mathrm{p}}$,
*C*_{c} is the Obukhov–Corrsin constant (*C*_{c}=0.67)
(Sreenivasan, 1996; Goto and Kida, 1999), *C*_{b} is the
Batchelor constant (*C*_{b}=3.7) (Grant et al., 1968; Goto and Kida, 1999),
and *γ*^{′} is the model parameter (${\mathit{\gamma}}^{\prime}=\mathrm{1.4}$).

The contribution of clear-air Bragg scattering was calculated by

where *E*_{n}(*k*) is the power spectrum of refractive index fluctuations.
The refractive index *n*_{ref} is given approximately by
(Balsley and Gage, 1980)

where *p* and *e* are the atmospheric pressure and partial pressure of water
vapor in hPa, and *T* is the absolute temperature. The contribution of free
electrons was omitted because it is negligibly small in the troposphere. The
power spectrum *E*_{n}(*k*) is given by the scalar concentration spectrum
for *P**r*<1 (Pao, 1964), where $Pr\equiv \mathit{\nu}/D$ (*D* is the
scalar diffusive coefficient) is the Prandtl number:

where *χ*_{n} is the scalar dissipation rate for *n*_{ref}.
In this study, the Prandtl number of refractive index was set to 0.7.

The scalar dissipation rates for $\langle {r}_{\mathrm{p}}^{\mathrm{3}}\rangle {n}_{\mathrm{p}}$ and
*n*_{ref} were calculated in the same way. That is, the dissipation rate
of an arbitrary scalar *θ* is given by

where *ν*_{t} is the eddy kinematic viscosity, *S**c*_{t} is
the turbulent Schmidt number, and $\widehat{\mathit{\theta}}$ is the scalar value in
the resolved scale. *ν*_{t} was calculated by using the Smagorinsky
model; i.e., ${\mathit{\nu}}_{\mathrm{t}}=({C}_{\mathrm{s}}{\mathrm{\Delta}}_{\mathrm{s}}{)}^{\mathrm{2}}{\left(\mathrm{2}\widehat{{s}_{ij}}\widehat{{s}_{ij}}\right)}^{\mathrm{1}/\mathrm{2}}$,
and *S**c*_{t} was set to 0.4 (Moin et al., 1991).

## 5.3 Results and discussion

Figure 8a shows the three-dimensional visualization of
liquid water. The optical thickness of each grid cell, *τ*_{Δ}, is
visualized by volume rendering to mimic human-eye observations of clouds.
Here, the optical thickness is defined by ${\mathit{\tau}}_{\mathrm{\Delta}}={Q}_{\mathrm{ext}}\mathit{\pi}\langle {r}_{\mathrm{p}}^{\mathrm{2}}\rangle {n}_{\mathrm{p}}{\mathrm{\Delta}}_{z}$,
where *Q*_{ext} is the extinction efficiency for Mie scattering (*Q*_{ext}=2.0 in this
study), and $\langle {r}_{\mathrm{p}}^{\mathrm{2}}\rangle $ is given by $\langle {r}_{\mathrm{p}}^{\mathrm{2}}\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{r}_{\mathrm{p}}^{\mathrm{2}}{q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)\mathrm{d}{r}_{\mathrm{p}}$. Note
that the optical transmittance of cloud volume is approximately equal to
1−*τ*_{Δ} when *τ*_{Δ} is sufficiently smaller than unity.
Figure 8b shows the isosurfaces of the energy dissipation
rate *ϵ* and the vertical velocity *u*_{3}. The locations of upward
flows correspond to the locations of clouds and large *ϵ* regions are
observed around the upward flows. This indicates that strong turbulence is
generated by entrainment motions due to updrafts.

This study focused on a vertical cross section that slices the cumulus cloud
with the largest upward velocity. Figure 9 shows the liquid
water content (LWC) in a logarithmic scale, the energy dissipation rate *ϵ*,
the radar reflectivity factor *Z*^{dB}, the increase of *Z*^{dB}
due to particulate Bragg scattering, and particulate Bragg
scattering due to turbulent clustering ${Z}_{\mathrm{PBc}}^{\mathrm{dB}}$ in the cross
section. The microwave frequency was set to *f*_{m}=2.8 GHz, which is the
representative frequency of S-band radars. The radar reflectivity factor is
shown in units of decibels, which is defined as *Z*^{dB} (dBZ)=10log _{10}*Z* (mm^{6} m^{−3}).
Large values of *Z*^{dB} are observed
inside and below the clouds. The strong echo below the clouds reflects
drizzling regions, where the LWC is smaller than that inside the clouds but
the strong echo returns from the drizzling droplets because *Z*_{incoh} is
proportional to $\langle {r}_{\mathrm{p}}^{\mathrm{6}}\rangle $. The radar echo on the outside
of the isoline of ${Z}_{\mathrm{incoh}}^{\mathrm{dB}}=-\mathrm{18}$ dBZ is caused by clear-air
Bragg scattering, i.e., *Z*_{CB}. The radar echo layer at the height from
2.2 to 2.5 km is caused by *Z*_{CB} due to a large humidity gap in the
inversion layer.

Figure 9c does not show a clear sign of the mantle echo
reported by Knight and Miller (1998). Knight and Miller (1998) reported
that the predominant mantle echo is observed on dry days, while it is poorly
observed on the most humid day. Our result is in accord with this report, as
the relative humidity of the environmental air is above 80 % at heights below
about 2.2 km in our cloud simulation data. A possible cause of the mantle
echo is particulate Bragg scattering due to turbulent entrainment (i.e.,
*Z*_{PBe}) because the large-scale cloud water inhomogeneity at cloud
edges produces small-scale fluctuations due to the turbulent cascade
(Erkelens et al., 2001). The influence of turbulent entrainment on *Z*^{dB}
turned out to be, however, negligibly small. That is, the fluctuations
caused by the large-scale inhomogeneity were not significantly large at the
scale of the half wavelength in the present simulation.

In Fig. 9d, the increase due to particulate Bragg
scattering, ${Z}^{\mathrm{dB}}-({Z}_{\mathrm{incoh}}+{Z}_{\mathrm{CB}}{)}^{\mathrm{dB}}$, is significant
at the near-top of the clouds. The maximum difference is larger than 5 dB. In
order to discuss the reason for the strong clustering influence at the
near-top of the clouds, the raw value of ${Z}_{\mathrm{PBc}}^{\mathrm{dB}}$ is plotted in
Fig. 9e. ${Z}_{\mathrm{PBc}}^{\mathrm{dB}}$ is larger than −10 dBZ
inside the turbulent cloud region, where the LWC is larger than
0.1 g m^{−3}, and the energy dissipation rate *ϵ* is intermittently larger than
100 cm^{2} s^{−3}. Large values of ${Z}_{\mathrm{PBc}}^{\mathrm{dB}}$ are shown at the
near-top inside this cloud region. We have confirmed that the droplet size in
this cloud region was almost homogeneous: the volume-averaged droplet radius
ranged within 7 to 11 µm. As a result, large values of the Stokes number
(up to 0.05) distributed intermittently corresponding to the distribution
of *ϵ*. The main factor of the height dependence
of ${Z}_{\mathrm{PBc}}^{\mathrm{dB}}$ is the LWC, which is larger than 1 g m^{−3} at the near-top of the
clouds. Note that *Z*_{PBc} is proportional to square of the LWC as
Eqs. (9) and (36) implies
${Z}_{\mathrm{PBc}}={\mathrm{2}}^{\mathrm{3}}{\mathrm{3}}^{\mathrm{2}}{\mathit{\rho}}_{\mathrm{p}}^{-\mathrm{2}}{\mathit{\kappa}}^{-\mathrm{2}}\left(\mathrm{LWC}{)}^{\mathrm{2}}{l}_{\mathit{\eta}}{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\right(\mathit{\kappa}{l}_{\mathit{\eta}})$.
Thus, the significant influence of turbulent clustering is caused
by sufficiently large values of the energy dissipation rate and the LWC.

This study has investigated the influence of microscale
turbulent clustering of polydisperse cloud droplets on the radar reflectivity
factor. Firstly, the theoretical solution for particulate Bragg scattering
for polydisperse droplets has been obtained considering the droplet size
distribution in the measurement volume and the droplet size dependence of
turbulent clustering. The obtained formula shows that the particulate Bragg-scattering part of the radar reflectivity factor is given by a double
integral function including the cross spectrum of number density fluctuations
for bidisperse droplets. Secondly, the wave number and Stokes number
dependence of the cross spectrum has been investigated using the turbulent
droplet clustering data obtained from a direct numerical simulation (DNS) of
particle-laden homogeneous isotropic turbulence without gravitational
settling. The result shows that the cross spectrum for a combination of
Stokes numbers, *S**t*_{1} and *S**t*_{2}, has values between the power
spectra for *S**t*_{1} and *S**t*_{2} at small wave numbers, whereas the
spectrum decreases more rapidly than the power spectra as the wave number
increases. This decreasing trend is related to the scale dependence of the
spatial correlation of cluster locations between two different Stokes
numbers. The coherence of the cross spectrum is close to unity for small
wave numbers and decreases almost exponentially with increasing wave number.
This is qualitatively consistent with the visualization results, in which the
clustering locations for different Stokes numbers are almost the same on
large scales, whereas a discrepancy in clustering locations is observed at
small scales. It is also confirmed that the decreasing trend of the coherence
is strongly dependent on the combination of Stokes numbers.

In order to develop a cross spectrum model for estimating the clustering influence on the radar reflectivity factor, we have proposed an exponential model for the wave number dependence of the coherence and introduced the critical wave number (i.e., the decay constant for the model) to consider the dependence of the coherence on the Stokes number combination. The coherence data for all combinations of six Stokes numbers ranging from 0.05 to 2.0 reveal that the critical wave number is inversely proportional to the Stokes number difference, $|S{t}_{\mathrm{1}}-S{t}_{\mathrm{2}}|$. This implies that the critical wave number is inversely proportional to the crossover length for the bidisperse radial distribution function (RDF). The proposed coherence model enables us to estimate the cross spectrum for arbitrary combinations of Stokes numbers using the power spectrum model proposed by Matsuda et al. (2014). Comparison of the model estimate with the DNS results for a typical droplet size distribution in cumulus clouds confirms the reliability of the Stokes number dependence of the proposed model.

The proposed model has been further extended for the case with gravitational
settling. We have assumed *S*_{v}*≲*1, where *S*_{v} is the
settling parameter, and modified the parameterization for the critical
wave number based on the analytical equation for the crossover length
considering the settling influence (Lu et al., 2010). The ${r}_{\mathrm{p}}^{\mathrm{3}}$-weighted
power spectrum estimated by the modified model shows a good
agreement with that obtained by the DNS data considering the droplet size
distribution in cumulus clouds and gravitational settling, indicating that
the proposed model can estimate the clustering influence on the radar
reflectivity factor to a sufficient accuracy.

Finally, the proposed model has been applied to high-resolution cloud-simulation data of Onishi and Takahashi (2012). The data were obtained using the Multi-Scale Simulator for the Geoenvironment (MSSG), which is a multi-scale nonhydrostatic atmosphere–ocean coupled model. The cloud and rain droplet size distribution was explicitly calculated at each grid using a spectral-bin cloud microphysics scheme. The radar reflectivity factor has been calculated considering particulate Bragg scattering due to turbulent clustering and turbulent entrainment as well as clear-air Bragg scattering caused by temperature and humidity fluctuations. The result shows that the influence of turbulent entrainment is negligibly small in our case, whereas the influence of turbulent clustering can be significant inside turbulent clouds. The large influence is observed at the near-top of the clouds, where the liquid water content (LWC) and the energy dissipation rate are sufficiently large.

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

KM designed the study, conducted the DNSs, developed the parameterization based n the DNS data, and analyzed the cloud simulation data. RO developed the DNS program, and provided the cloud simulation data. Both authors contributed to the discussion, interpretation of the results, and to writing the paper.

The authors declare that they have no conflict of interest.

This research was supported by JSPS KAKENHI grant number JP17K14598. The
numerical simulations presented were carried out on the Earth Simulator
supercomputer system of the Japan Agency for Marine-Earth Science and
Technology (JAMSTEC).

Edited by: Graham Feingold

Reviewed by: two anonymous referees

Ayala, O., Rosa, B., and Wang, L.-P.: Effects of turbulence on the geometric collision rate of sedimenting droplets. Part 2. Theory and parameterization, New J. Phys., 10, 075016, https://doi.org/10.1088/1367-2630/10/7/075016, 2008a. a, b

Ayala, O., Rosa, B., Wang, L.-P., and Grabowski, W.: Effects of turbulence on the geometric collision rate of sedimenting droplets. Part 1. Results from direct numerical simulation, New J. Phys., 10, 075015, https://doi.org/10.1088/1367-2630/10/7/075015, 2008b. a

Balsley, B. and Gage, K.: The MST Radar Technique: Potential for Middle Atmospheric Studies, Pure Appl. Geophys., 118, 452–493, 1980. a

Bec, J., Homann, H., and Ray, S.: Gravity-Driven Enhancement of Heavy Particle Clustering in Turbulent Flow, Phys. Rev. Lett., 112, 184501, https://doi.org/10.1103/PhysRevLett.112.184501, 2014. a, b

Bohren, C. and Huffman, D.: Absorption and Scattering of Light by Small Particles, Wiley, 1983. a

Bringi, V., Chandrasekar, V., Balakrishnan, N., and Zrnić, D.: An examination of propagation effects in rainfall on radar measurements at microwave frequencies, J. Atmos. Ocean. Tech., 7, 829–840, 1990. a

Carey, L., Rutledge, S., Ahijevych, D., and Keenan, T.: Correcting Propagation Effects in C-Band Polarimetric Radar Observations of Tropical Convection Using Differential Propagation Phase, J. Atmos. Sci., 39, 1405–1433, 2000. a

Chen, L., Goto, S., and Vassilicos, J.: Turbulent clustering of stagnation points and inertial particles, J. Fluid Mech., 553, 143–154, 2006. a

Chun, J., Koch, D., Rani, S., Ahluwalia, A., and Collins, L.: Clustering of aerosol particles in isotropic turbulence, J. Fluid Mech., 536, 219–251, 2005. a, b, c, d, e, f, g, h

Dombrovsky, L. and Zaichik, L.: An effect of turbulent clustering on scattering of microwave radiation by small particles in the atmosphere, J. Quant. Spectrosc. Ra., 111, 234–242, 2010. a

Erkelens, J., Venema, V., Russchenberg, H., and Ligthart, L.: Coherent Scattering of Microwave by Particles: Evidence from Clouds and Smoke, J. Atmos. Sci., 58, 1091–1102, 2001. a, b, c, d

Gossard, E. and Strauch, R.: Radar Observation of Clear Air and Clouds, in: vol. 14 of Developments in Atmospheric Science, Elsevir, New York, 1983. a, b, c, d, e

Goto, S. and Kida, S.: Passive scalar spectrum in isotropic turbulence: Prediction by the Lagrangian direct-interaction approximation, Phys. Fluids, 11, 1936–1952, 1999. a, b

Grabowski, W. and Vaillancourt, P.: Comments on “Preferential Concentration of Cloud Droplets by Turbulence: Effects on the Early Evolution of Cumulus Cloud Droplet Spectra”, J. Atmos. Sci., 56, 1433–1436, 1999. a

Grant, H., Hughes, B., Vogel, W., and Moilliet, A.: The spectrum of temperature fluctuations in turbulent flow, J. Fluid Mech., 34, 423–442, 1968. a

Hess, M., Koepke, P., and Schult, I.: Optical Properties of Aerosols and Clouds: The Software Package OPAC, B. Am. Meteorol. Soc., 79, 831–844, 1998. a

Hill, R.: Models of the scalar spectrum for turbulent advection, J. Fluid Mech., 88, 541–562, 1978. a

Hirt, C. and Cook, J.: Calculating three-dimensional flow around structures, J. Comput. Phys., 10, 324–340, 1972. a

Ireland, P., Bragg, A., and Collins, L.: The effect of Reynolds number on inertial particle dynamics in isotropic turbulence. Part 2. Simulations with gravitational effects, J. Fluid Mech., 796, 659–711, 2016. a, b

Knight, C. and Miller, L.: Early radar echoes from small, warm cumulus: Bragg and hydrometeor scattering, J. Atmos. Sci., 55, 2974–2992, 1998. a, b, c, d

Kostinski, A. and Jameson, A.: On the spatial distribution of cloud particles, J. Atmos. Sci., 57, 901–915, 2000. a, b

Lu, J., Nordslek, H., and Shaw, R.: Clustering of settling charged particles in turbulence: theory and experiments, New J. Phys., 12, 123030, https://doi.org/10.1088/1367-2630/12/12/123030, 2010. a, b, c, d

Matsuda, K., Onishi, R., Hirahara, M., Kurose, R., Takahashi, K., and Komori, S.: Influence of microscale turbulent droplet clustering on radar cloud observations, J. Atmos. Sci., 71, 3569–3582, 2014. a, b, c, d, e, f, g, h, i, j, k

Matsuda, K., Onishi, R., and Takahashi, K.: Influence of gravitational settling on turbulent droplet clustering and radar reflectivity factor, Flow Turb. Combust., 98, 327–340, 2017. a, b

Maxey, M.: The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields, J. Fluid Mech., 174, 441–465, 1987. a

Moin, P., Squires, K., Cabot, W., and Lee, S.: A dynamic subgridscale model for compressible turbulence and scalar transport, Phys. Fluids A, 3, 2746–2757, 1991. a

Morinishi, Y., Lundm, T., Vasilyev, O., and Moin, P.: Fully conservative higher order finite difference schemes for incompressible flow, J. Comput. Phys., 143, 90–124, 1998. a

Onishi, R. and Takahashi, K.: A Warm-Bin–Cold-Bulk Cloud Microphysical Model, J. Atmos. Sci., 69, 1474–1497, 2012. a, b, c, d, e

Onishi, R. and Vassilicos, J.: Collision Statistics of Inertial Particles in Two-Dimensional Homogeneous Isotropic Turbulence with an Inverse Cascade, J. Fluid Mech., 745, 279–299, 2014. a

Onishi, R., Takahashi, K., and Komori, S.: Influence of gravity on collisions of monodispersed droplets in homogeneous isotropic turbulence, Phys. Fluids, 21, 125108, https://doi.org/10.1063/1.3276906, 2009. a

Onishi, R., Baba, Y., and Takahashi, K.: Large-Scale Forcing with Less Communication in Finite-Difference Simulations of Steady Isotropic Turbulence, J. Comput. Phys., 230, 4088–4099, 2011. a

Pao, Y.: Statistical Behavior of a Turbulent Multicomponent Mixture with First-Order Reactions, AIAA J., 2, 1550–1559, 1964. a

Pinsky, M., Khain, A., and Krugliak, H.: Collision of Cloud Droplets in a Turbulent Flow. Part V: Application of Detailed Tables of Turbulent Collision Rate Enhancement to Simulation of Droplet Spectra Evolution, J. Atmos. Sci., 65, 357–374, 2008. a

Reade, W. and Collins, L.: Effect of preferential concentration on turbulent collision rates, Phys. Fluids, 12, 2530, https://doi.org/10.1063/1.1288515, 2000. a

Rogers, R. and Brown, W.: Radar Observations of a Major Industrial Fire, B. Am. Meteorol. Soc., 78, 803–814, 1997. a

Soong, S.-T.: Numerical simulation of warm rain development in an axisymmetric cloud model, J. Atmos. Sci., 31, 1262–1285, 1974. a

Squires, K. and Eaton, J.: Preferential concentration of particles by turbulence, Phys. Fluids A, 3, 1169–1178, 1991. a

Sreenivasan, K.: The passive scalar spectrum and the Obukhov–Corrsin constant, Phys. Fluids, 8, 189–196, 1996. a

Sundaram, S. and Collins, L.: Collision statistics in an isotropic particle-laden turbulent suspension. Part 1. Direct numerical simulations, J. Fluid Mech., 335, 75–109, 1997. a

Twomey, S.: The nuclei of natural cloud formation. Part II: The supersaturation in natural clouds and the variation of cloud droplets concentrations, Geofis. Pura Appl., 43, 243–249, 1959. a

van Zanten, M., Nuijens, L., and Siebesma, A.: Precipitating Shallow Cumulus Case1, available at: http://www.knmi.nl/samenw/rico/ (last access: 6 February 2019), 2006. a

Wang, L. and Maxey, M.: Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence, J. Fluid Mach., 256, 27–68, 1993. a, b

Wang, L., Rosa, B., Gao, H., He, G., and Jin, G.: Turbulent collision of inertial particles: Point-particle based, hybrid simulations and beyond, Int. J. Multiphase Flow, 35, 854–867, 2009. a

Zhou, Y., Wexler, A., and Wang, L.: Modelling turbulent collision of bidisperse inertial particles, J. Fluid Mech., 433, 77–104, 2001. a