Journal cover Journal topic
Atmospheric Chemistry and Physics An interactive open-access journal of the European Geosciences Union
Journal topic
Atmos. Chem. Phys., 19, 1785-1799, 2019
https://doi.org/10.5194/acp-19-1785-2019
Atmos. Chem. Phys., 19, 1785-1799, 2019
https://doi.org/10.5194/acp-19-1785-2019

Research article 08 Feb 2019

Research article | 08 Feb 2019

# Turbulent enhancement of radar reflectivity factor for polydisperse cloud droplets

Turbulent enhancement of radar reflectivity factor for polydisperse cloud droplets
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
Abstract

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.

1 Introduction

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 (mm6 m−3), Z, which is obtained with the following radar equation:

$\begin{array}{}\text{(1)}& {P}_{\mathrm{r}}=\frac{{P}_{\mathrm{t}}G{A}_{\mathrm{e}}{k}_{\mathrm{m}}^{\mathrm{4}}|K{|}^{\mathrm{2}}V}{{\mathrm{4}}^{\mathrm{4}}{R}^{\mathrm{4}}}Z,\end{array}$

where Pr and Pt are the received and transmitted microwave powers, G is the antenna gain, Ae is the effective aperture of the antenna, km 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 . Incoherent scattering implies random and uniform dispersion of cloud droplets . 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 . 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 . 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 fm ($={k}_{\mathrm{m}}c/\mathrm{2}\mathit{\pi }$, where c is the speed of light), while the factor for incoherent scattering is independent of fm. 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. 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. 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. first pointed out the possibility that particulate Bragg scattering due to turbulent clustering leads to the frequency dependence reported by . 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 , 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 . 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. 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, Enp(k|rp), where k is the wave number and rp is the droplet radius. The power spectrum Enp(k|rp) is strongly dependent on the droplet size: more specifically, Enp(k|rp) is dependent on the Stokes number, St, 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.

2 Theory

Here, we aim to formulate the radar reflectivity factor Z for a nonuniform distribution of polydisperse cloud droplets based on the discussion of , 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, Esca(t, x,rp), is given by the Rayleigh-scattering approximation:

$\begin{array}{ll}{\mathbit{E}}_{\mathrm{sca}}\left(t,\mathbit{x},{r}_{\mathrm{p}}\right)& ={\mathbit{E}}_{\mathrm{inc}}\frac{{k}_{\mathrm{m}}^{\mathrm{2}}K{r}_{\mathrm{p}}^{\mathrm{3}}}{R}\mathrm{sin}\mathit{\chi }\mathrm{exp}\left[i\left\{\mathit{\omega }t-{\mathbit{k}}_{\mathrm{sca}}\right\\right\\ \text{(2)}& & \cdot \left({\mathbit{x}}_{\mathrm{r}}-\mathbit{x}\right)-{\mathbit{k}}_{\mathrm{inc}}\cdot \left(\mathbit{x}-{\mathbit{x}}_{\mathrm{t}}\right)}],\end{array}$

where Einc is the electric field amplitude vector of the incident microwave; rp is the droplet radius; x, xt, and xr are the position vectors of the droplet, microwave transmitter, and microwave receiver; kinc and ksca are the wave number vectors of the incident and scattered microwaves which satisfy $|{\mathbit{k}}_{\mathrm{inc}}|=|{\mathbit{k}}_{\mathrm{sca}}|={k}_{\mathrm{m}}$, and χ is the angle between Einc and ksca. Considering the droplet-size dependence of Esca(t, x, rp), the electric power of the microwave scattered by a group of droplets, Ps, is given by

$\begin{array}{}\text{(3)}& {P}_{\mathrm{s}}=\stackrel{\mathrm{‾}}{{\left|\underset{\mathbit{x}\in V}{\int }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathbit{E}}_{\mathrm{sca}}\left(t,\mathbit{x},{r}_{\mathrm{p}}\right)n\left(\mathbit{x},{r}_{\mathrm{p}}\right)\mathrm{d}{r}_{\mathrm{p}}\mathrm{d}\mathbit{x}\right|}^{\mathrm{2}}}/\mathit{\zeta },\end{array}$

where n(x, rp)drpdx is the number of droplets with radii from rp to ${r}_{{}_{\mathrm{p}}}+\mathrm{d}{r}_{{}_{\mathrm{p}}}$ in an infinitesimal volume dx at position x, and ζ 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

$\begin{array}{}\text{(4)}& Z=\frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K{|}^{\mathrm{2}}V{\mathrm{sin}}^{\mathrm{2}}\mathit{\chi }}\mathit{\sigma },\end{array}$

where λ is the microwave wavelength and σ is the radar cross section , which is defined as

$\begin{array}{}\text{(5)}& \mathit{\sigma }=\mathrm{4}\mathit{\pi }{R}^{\mathrm{2}}\frac{{P}_{\mathrm{s}}}{{P}_{\mathrm{o}}},\end{array}$

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

$\begin{array}{}\text{(6)}& Z=\frac{{\mathrm{2}}^{\mathrm{6}}}{V}\stackrel{\mathrm{‾}}{{\left|\underset{\mathbit{x}\in V}{\int }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{r}_{\mathrm{p}}^{\mathrm{3}}n\left(\mathbit{x},{r}_{\mathrm{p}}\right)\mathrm{exp}\left(-i\mathbit{\kappa }\cdot \mathbit{x}\right)\mathrm{d}{r}_{\mathrm{p}}\mathrm{d}\mathbit{x}\right|}^{\mathrm{2}}},\end{array}$

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

Similarly to , we assume n(x, rp) to be composed of the temporal-average and fluctuation parts; i.e., n(x, ${r}_{\mathrm{p}}\right)=\stackrel{\mathrm{‾}}{n\left(\mathbit{x},{r}_{\mathrm{p}}\right)}+\mathit{\delta }n\left(\mathbit{x},{r}_{\mathrm{p}}\right)$. The temporal-average part, $\stackrel{\mathrm{‾}}{n\left(\mathbit{x},{r}_{\mathrm{p}}\right)}$, contributes to the separated reflection; therefore, the contribution of this part is negligibly small when $\stackrel{\mathrm{‾}}{n\left(\mathbit{x},{r}_{\mathrm{p}}\right)}$ has no fluctuation on a spatial scale of half the wavelength . Thus, we neglect the contribution of $\stackrel{\mathrm{‾}}{n\left(\mathbit{x},{r}_{\mathrm{p}}\right)}$. Then, we obtain

$\begin{array}{ll}Z=& {\mathrm{2}}^{\mathrm{6}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left\{{r}_{\mathrm{p}}^{\mathrm{3}}{{r}_{\mathrm{p}}^{\prime }}^{\mathrm{3}}\underset{\mathbit{r}}{\int }〈\mathit{\delta }n\left(\mathbit{x},{r}_{\mathrm{p}}\right)\mathit{\delta }n\left(\mathbit{x}+\mathbit{r},{r}_{\mathrm{p}}^{\prime }\right)〉\right\\\ \text{(7)}& & \mathrm{exp}\left(-i\mathbit{\kappa }\cdot \mathbit{r}\right)\mathrm{d}\mathbit{r}}\mathrm{d}{r}_{\mathrm{p}}\mathrm{d}{r}_{\mathrm{p}}^{\prime },\end{array}$

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

In order to decompose the spatial correlation function $〈\mathit{\delta }n\left(\mathbit{x},{r}_{\mathrm{p}}\right)\mathit{\delta }n\left(\mathbit{x}+\mathbit{r},{r}_{\mathrm{p}}^{\prime }\right)〉$, we introduce the probability density function (PDF) of droplet radius rp to the measurement volume, qr(rp), and the number density distribution function for monodisperse droplets with a radius of rp, np(x|rp). The PDF is defined as ${q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)\equiv \frac{\mathrm{1}}{{N}_{\mathrm{p}}}\underset{\mathbit{x}\in V}{\int }\stackrel{\mathrm{‾}}{n\left(\mathbit{x},{r}_{\mathrm{p}}\right)}\mathrm{d}\mathbit{x}$, where Np is the total number of droplets in the measurement volume; i.e., ${N}_{\mathrm{p}}\equiv \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathbit{x}\in V}{\int }\stackrel{\mathrm{‾}}{n\left(\mathbit{x},{r}_{\mathrm{p}}\right)}\mathrm{d}\mathbit{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(\mathbit{x}|{r}_{\mathrm{p}}\right)\equiv n\left(\mathbit{x},{r}_{\mathrm{p}}\right)/{q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)$ so that n(x,rp) is given by $n\left(\mathbit{x},{r}_{\mathrm{p}}\right)={n}_{\mathrm{p}}\left(\mathbit{x}|{r}_{\mathrm{p}}\right){q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)$. The number density distribution function np(x|rp) satisfies $\underset{\mathbit{x}\in V}{\int }\stackrel{\mathrm{‾}}{{n}_{\mathrm{p}}\left(\mathbit{x}|{r}_{\mathrm{p}}\right)}\mathrm{d}\mathbit{x}={N}_{\mathrm{p}}$ for arbitrary rp. Note that the spatial correlation function $〈\mathit{\delta }n\left(\mathbit{x},{r}_{\mathrm{p}}\right)\mathit{\delta }n\left(\mathbit{x}+\mathbit{r},{r}_{\mathrm{p}}^{\prime }\right)〉$ 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 $〈n\left(\mathbit{x},{r}_{\mathrm{p}}\right)〉\mathit{\delta }\left(\mathbit{r}\right)\mathit{\delta }\left({r}_{\mathrm{p}}^{\prime }-{r}_{\mathrm{p}}\right)$, where δ(r) and $\mathit{\delta }\left({r}_{\mathrm{p}}^{\prime }-{r}_{\mathrm{p}}\right)$ are the Dirac delta functions. Thus, the spatial correlation function is given by

$\begin{array}{ll}〈\mathit{\delta }n\left(\mathbit{x},{r}_{\mathrm{p}}\right)\mathit{\delta }n\left(\mathbit{x}+\mathbit{r},{r}_{\mathrm{p}}^{\prime }\right)〉& =〈{n}_{\mathrm{p}}〉\mathit{\delta }\left(\mathbit{r}\right){q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)\mathit{\delta }\left({r}_{\mathrm{p}}^{\prime }-{r}_{\mathrm{p}}\right)\\ \text{(8)}& & +\mathrm{\Psi }\left(\mathbit{r}|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right){q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right){q}_{\mathrm{r}}\left({r}_{\mathrm{p}}^{\prime }\right),\end{array}$

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

$\begin{array}{}\text{(9)}& Z={\mathrm{2}}^{\mathrm{6}}〈{r}_{\mathrm{p}}^{\mathrm{6}}〉〈{n}_{\mathrm{p}}〉+{\mathrm{2}}^{\mathrm{7}}{\mathit{\pi }}^{\mathrm{2}}{\mathit{\kappa }}^{-\mathrm{2}}{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}\left(\mathit{\kappa }\right),\end{array}$

where κ is $\mathit{\kappa }=|\mathbit{\kappa }|$, $〈{r}_{\mathrm{p}}^{\mathrm{6}}〉$ is given by $〈{r}_{\mathrm{p}}^{\mathrm{6}}〉=\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 Er3np(k) is the ${r}_{\mathrm{p}}^{\mathrm{3}}$-weighted power spectrum, defined as

$\begin{array}{}\text{(10)}& {E}_{\mathrm{r}\mathrm{3}\mathrm{np}}\left(k\right)\equiv \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{r}_{\mathrm{p}}^{\mathrm{3}}{r}_{\mathrm{p}}^{{\mathrm{3}}^{\prime }}{q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right){q}_{\mathrm{r}}\left({r}_{\mathrm{p}}^{\prime }\right){C}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)\mathrm{d}{r}_{\mathrm{p}}\mathrm{d}{r}_{\mathrm{p}}^{\prime },\end{array}$

where ${C}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)$ is the cross spectrum of number density fluctuations for np(x|rp) and ${n}_{\mathrm{p}}\left(\mathbit{x}|{r}_{\mathrm{p}}^{\prime }\right)$: the cross spectrum ${C}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)$ is defined as ${C}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)\equiv \underset{|\mathbit{k}|=k}{\int }\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathbit{k}|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)\mathrm{d}{\mathit{\sigma }}_{k}$; i.e., the integration of $\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathbit{k}|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)$ over the spherical shell, σk, at $|\mathbit{k}|=k$, in which $\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathbit{k}|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)$ is the cross spectral density function, defined as

$\begin{array}{}\text{(11)}& \stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathbit{k}|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)=\frac{\mathrm{1}}{\left(\mathrm{2}\mathit{\pi }{\right)}^{\mathrm{3}}}\underset{\mathbit{r}}{\int }\mathrm{\Psi }\left(\mathbit{r}|{r}_{\mathrm{p}},{r}_{\mathrm{p}}^{\prime }\right)\mathrm{exp}\left(-i\mathbit{k}\cdot \mathbit{r}\right)\mathrm{d}\mathbit{r}.\end{array}$

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}}〈{r}_{\mathrm{p}}^{\mathrm{6}}〉〈{n}_{\mathrm{p}}〉$.

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 rp1, the PDF of droplet radius is given by ${q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)=\mathit{\delta }\left({r}_{\mathrm{p}}-{r}_{\mathrm{p}\mathrm{1}}\right)$. Then, the radar reflectivity factor Z is given by

$\begin{array}{}\text{(12)}& Z={\mathrm{2}}^{\mathrm{6}}{r}_{\mathrm{p}\mathrm{1}}^{\mathrm{6}}〈{n}_{\mathrm{p}}〉+{\mathrm{2}}^{\mathrm{7}}{\mathit{\pi }}^{\mathrm{2}}{\mathit{\kappa }}^{-\mathrm{2}}{r}_{\mathrm{p}\mathrm{1}}^{\mathrm{6}}{E}_{\mathrm{np}}\left(\mathit{\kappa }|{r}_{\mathrm{p}\mathrm{1}}\right),\end{array}$

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

$\begin{array}{}\text{(13)}& \stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathbit{k}|{r}_{\mathrm{p}}\right)=\frac{\mathrm{1}}{\left(\mathrm{2}\mathit{\pi }{\right)}^{\mathrm{3}}}\underset{\mathbit{r}}{\int }\mathrm{\Psi }\left(\mathbit{r}|{r}_{\mathrm{p}},{r}_{\mathrm{p}}\right)\mathrm{exp}\left(-i\mathbit{k}\cdot \mathbit{r}\right)\mathrm{d}\mathbit{r}.\end{array}$
3 Computational method

## 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:

$\begin{array}{}\text{(14)}& & \frac{\partial {u}_{i}}{\partial {x}_{i}}=\mathrm{0},\text{(15)}& & \frac{\partial {u}_{i}}{\partial t}+\frac{\partial {u}_{i}{u}_{j}}{\partial {x}_{j}}=-\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}\frac{\partial p}{\partial {x}_{i}}+\mathit{\nu }\frac{{\partial }^{\mathrm{2}}{u}_{i}}{\partial {x}_{j}\partial {x}_{j}}+{F}_{i},\end{array}$

where ui is the flow velocity in the ith direction, p is the pressure, ρa is the air density, ν is the kinematic viscosity, and Fi is the external forcing term. The nonlinear term was discretized by the fourth-order central difference scheme . The time integration was calculated by the second-order Runge–Kutta scheme. The HSMAC method 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 $|\mathbit{k}{L}_{\mathrm{0}}|<\mathrm{2}$ , where L0 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

$\begin{array}{}\text{(16)}& \frac{\mathrm{d}{v}_{i}}{\mathrm{d}t}=-\frac{{v}_{i}-{u}_{i}}{{\mathit{\tau }}_{\mathrm{p}}}+{g}_{i},\end{array}$

where vi and gi are the particle velocity and gravitational acceleration in the ith direction. τp is the droplet relaxation time, which is given by

$\begin{array}{}\text{(17)}& {\mathit{\tau }}_{\mathrm{p}}=\frac{{\mathit{\rho }}_{\mathrm{p}}}{{\mathit{\rho }}_{\mathrm{a}}}\frac{\mathrm{2}{r}_{\mathrm{p}}^{\mathrm{2}}}{\mathrm{9}\mathit{\nu }}.\end{array}$

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πL0. 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 5123. A Taylor microscale-based turbulent Reynolds number of the obtained flow was Reλ=204, where Reλ 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 Reλ is sufficiently large to obtain turbulent clustering data for high Reynolds number turbulence in the wave number range relevant to radar observations (see Sect. 3.2). The kinematic viscosity, ν, was set to $\mathrm{1.5}×{\mathrm{10}}^{-\mathrm{5}}$ m2 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, Np, was set to 1.5×107.

Table 1Computational settings of the DNS.

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, St. The clustering data for the monodisperse cases were used for calculating the cross spectrum of number density fluctuations for any combinations of these St values. For polydisperse droplets, a typical droplet size distribution for maritime cumulus clouds (the size distribution data named “CUMA” in ) 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 cm2 s−3, which can be observed in cumulus and cumulonimbus clouds . 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 cm2 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, Fr ($\equiv {a}_{\mathit{\eta }}/g$, in which ${a}_{\mathit{\eta }}\equiv {\mathit{ϵ}}^{\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 cm2 s−3, respectively. The influence of gravitational settling on turbulent particle clustering is often discussed using the settling parameter, Sv, which is defined as ${S}_{\mathrm{v}}\equiv {v}_{\mathrm{T}}/{u}_{\mathit{\eta }}$ , where vT=τpg is the terminal settling velocity and ${u}_{\mathit{\eta }}\equiv \left(\mathit{\nu }\mathit{ϵ}{\right)}^{\mathrm{1}/\mathrm{4}}$ is the Kolmogorov velocity. Note that Sv satisfies ${S}_{\mathrm{v}}=St/Fr$. The settling parameters corresponding to the modal radius for ϵ=100, 400, and 1000 cm2 s−3 were 0.67, 0.48, and 0.38, respectively.

Figure 1Spatial distribution of droplets for (orange) St=0.2, (green) 0.5, and (blue) 1.0 in the regions of (a) 2πL0×2πL0 and (b) 0.5πL0×0.5πL0. Droplets located in the range of $\mathrm{0} are shown. The square frame in (a) corresponds to the region of (b).

## 3.2 Computation of power spectrum and cross spectrum

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

$\begin{array}{}\text{(18)}& & \stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathbit{k}|{r}_{\mathrm{p}}\right)={L}_{\mathrm{0}}^{-\mathrm{3}}〈\stackrel{\mathrm{̃}}{{n}_{\mathrm{p}}}\left(\mathbit{k}|{r}_{\mathrm{p}}\right)\stackrel{\mathrm{̃}}{{n}_{\mathrm{p}}}\left(-\mathbit{k}|{r}_{\mathrm{p}}\right)〉,\text{(19)}& & \stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathbit{k}|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}}\right)={L}_{\mathrm{0}}^{-\mathrm{3}}〈\stackrel{\mathrm{̃}}{{n}_{\mathrm{p}}}\left(\mathbit{k}|{r}_{\mathrm{p}\mathrm{1}}\right)\stackrel{\mathrm{̃}}{{n}_{\mathrm{p}}}\left(-\mathbit{k}|{r}_{\mathrm{p}\mathrm{2}}\right)〉,\end{array}$

where $\stackrel{\mathrm{̃}}{{n}_{\mathrm{p}}}\left(\mathbit{k}|{r}_{\mathrm{p}}\right)$ is the Fourier component of the droplet number density distribution, np(x|rp), and the angle brackets denote an ensemble average. The number density distribution for Lagrangian discrete droplets is given by

$\begin{array}{}\text{(20)}& {n}_{\mathrm{p}}\left(\mathbit{x}|{r}_{\mathrm{p}}\right)=\sum _{j=\mathrm{1}}^{{N}_{\mathrm{p}}}\mathit{\delta }\left(\mathbit{x}-{\mathbit{x}}_{\mathrm{p},j}\right),\end{array}$

where xp,j is the position vector of the jth droplet with radius rp, and Np is the total number of droplets with radius rp. The Fourier components of Eq. (20) are then given by

$\begin{array}{}\text{(21)}& \stackrel{\mathrm{̃}}{{n}_{\mathrm{p}}}\left(\mathbit{k}|{r}_{\mathrm{p}}\right)=\frac{\mathrm{1}}{\left(\mathrm{2}\mathit{\pi }{\right)}^{\mathrm{3}}}\sum _{j=\mathrm{1}}^{{N}_{\mathrm{p}}}\mathrm{exp}\left(-i\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p},j}\right).\end{array}$

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

$\begin{array}{}\text{(22)}& & \frac{\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathbit{k}|{r}_{\mathrm{p}}\right)}{〈{n}_{\mathrm{p}}{〉}^{\mathrm{2}}{L}_{\mathrm{0}}^{\mathrm{3}}}〉=\frac{\mathrm{1}}{{N}_{\mathrm{p}}^{\mathrm{2}}}〈\sum _{j=\mathrm{1}}^{{N}_{\mathrm{p}}}\mathrm{exp}\left(-i\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p},j}\right)\sum _{{j}^{\prime }=\mathrm{1},{j}^{\prime }\ne j}^{{N}_{\mathrm{p}}}\mathrm{exp}\left(i\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p},{j}^{\prime }}\right)〉\text{(23)}& & =\frac{\mathrm{1}}{{N}_{\mathrm{p}}^{\mathrm{2}}}\left[〈{\left\{\sum _{j=\mathrm{1}}^{{N}_{\mathrm{p}}}\mathrm{cos}\left(\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p},j}\right)\right\}}^{\mathrm{2}}〉+〈{\left\{\sum _{j=\mathrm{1}}^{{N}_{\mathrm{p}}}\mathrm{sin}\left(\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p},j}\right)\right\}}^{\mathrm{2}}〉\right]-\frac{\mathrm{1}}{{N}_{\mathrm{p}}}.\end{array}$

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

$\begin{array}{ll}\frac{\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathbit{k}|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}}\right)}{〈{n}_{\mathrm{p}\mathrm{1}}〉〈{n}_{\mathrm{p}\mathrm{2}}〉{L}_{\mathrm{0}}^{\mathrm{3}}}=& \frac{\mathrm{1}}{{N}_{\mathrm{p}\mathrm{1}}{N}_{\mathrm{p}\mathrm{2}}}〈\sum _{j=\mathrm{1}}^{{N}_{\mathrm{p}\mathrm{1}}}\mathrm{exp}\left(-i\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p}\mathrm{1},j}\right)\\ \text{(24)}& & \sum _{{j}^{\prime }=\mathrm{1}}^{{N}_{\mathrm{p}\mathrm{2}}}\mathrm{exp}\left(i\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p}\mathrm{2},{j}^{\prime }}\right)〉\end{array}$

$\begin{array}{ll}=\frac{\mathrm{1}}{{N}_{\mathrm{p}\mathrm{1}}{N}_{\mathrm{p}\mathrm{2}}}& \left\{〈\sum _{j=\mathrm{1}}^{{N}_{\mathrm{p}\mathrm{1}}}\mathrm{cos}\left(\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p}\mathrm{1},j}\right)\sum _{{j}^{\prime }=\mathrm{1}}^{{N}_{\mathrm{p}\mathrm{2}}}\mathrm{cos}\left(\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p}\mathrm{2},{j}^{\prime }}\right)〉\right\\\ \text{(25)}& & +〈\sum _{j=\mathrm{1}}^{{N}_{\mathrm{p}\mathrm{1}}}\mathrm{sin}\left(\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p}\mathrm{1},j}\right)\sum _{{j}^{\prime }=\mathrm{1}}^{{N}_{\mathrm{p}\mathrm{2}}}\mathrm{sin}\left(\mathbit{k}\cdot {\mathbit{x}}_{\mathrm{p}\mathrm{2},{j}^{\prime }}\right)〉},\end{array}$

where np1 and np2 are the average number density of droplets with radii of rp1 and rp2 (i.e., $〈{n}_{\mathrm{p}\mathrm{1}}〉\equiv 〈{n}_{\mathrm{p}}\left(\mathbit{x}|{r}_{\mathrm{p}\mathrm{1}}\right)〉$ and $〈{n}_{\mathrm{p}\mathrm{2}}〉\equiv 〈{n}_{\mathrm{p}}\left(\mathbit{x}|{r}_{\mathrm{p}\mathrm{2}}\right)〉$), respectively, Np1 and Np2 are the numbers of droplets with radii of rp1 and rp2, respectively, xp1,j is the position of the jth droplet with a radius of rp1, and ${\mathbit{x}}_{\mathrm{p}\mathrm{2},{j}^{\prime }}$ is the position of the jth droplet with a radius of rp2. It should be noted that the imaginary part of $\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathbit{k}|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}}\right)$ is omitted in Eq. (25) because, statistically, it should be zero. We confirmed that the imaginary part of ${C}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}}\right)$ calculated from the DNS data is O(10−4), which is caused by the statistical and truncation errors.

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

$\begin{array}{}\text{(26)}& & {E}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}}\right)=\frac{\mathrm{1}}{\mathrm{\Delta }k}\sum _{k-\mathrm{\Delta }k/\mathrm{2}\le |\mathbit{k}|

The spectra, Enp(k|rp) and ${C}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}}\right)$, were calculated for 19 representative wave numbers in the same way as . 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 St cases.

Figure 2Normalized cross spectra ${C}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}}$, St2) of droplet number density fluctuations compared to normalized power spectra ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|St\right)$.

normalized the power spectrum Enp(k|rp) by using np and the Kolmogorov scale, lη, defined as ${l}_{\mathit{\eta }}\equiv {\mathit{\nu }}^{\mathrm{3}/\mathrm{4}}{\mathit{ϵ}}^{-\mathrm{1}/\mathrm{4}}$, and confirmed that, for Reλ≥204, the Reλ dependence of the normalized power spectrum is negligible at the wave number range relevant to radar observations ($\mathrm{0.05}). Thus, this study used the same normalization for Enp(k|rp) and ${C}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}}\right)$ as follows:

$\begin{array}{}\text{(28)}& {E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|St\right)={E}_{\mathrm{np}}^{*}\left(k{l}_{\mathit{\eta }}|St\right)\equiv \frac{{E}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}}\right)}{〈{n}_{\mathrm{p}}{〉}^{\mathrm{2}}{l}_{\mathit{\eta }}},\end{array}$

$\begin{array}{}\text{(29)}& {C}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)={C}_{\mathrm{np}}^{*}\left(k{l}_{\mathit{\eta }}|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)\equiv \frac{{C}_{\mathrm{np}}\left(k|{r}_{\mathrm{p}\mathrm{1}},{r}_{\mathrm{p}\mathrm{2}}\right)}{〈{n}_{\mathrm{p}\mathrm{1}}〉〈{n}_{\mathrm{p}\mathrm{2}}〉{l}_{\mathit{\eta }}},\end{array}$

where ξ is the normalized wave number defined as ξklη and St, St1, and St2 are the Stokes numbers of water droplets with radii of rp, rp1, and rp2.

Figure 3Coherence of cross spectra for combinations of St1=0.5 and other Stokes numbers. The same coherence is plotted in (a) a double logarithmic chart and (b) a vertically logarithmic chart.

4 DNS results and parameterization

## 4.1 Spatial droplet distribution and cross spectrum

Figure 1 shows the spatial distribution of droplets for St=0.2, 0.5, and 1.0. Droplets located in the range of $\mathrm{0} 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πL0×2πL0, while Fig. 1b is the magnified view for the region with a size of 0.5πL0×0.5πL0. The clusters and void areas are clearly observed for all St cases. Their location for St=0.5 is almost the same as that for the other St. This is because the locations of clusters and void areas for small St 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 St increases. Figure 2 shows the power spectra ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|St\right)$ and the cross spectra ${C}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)$ 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 }|St\right)$ 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 . 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 St1=0.5 and St2=0.2 is located between the power spectra ${E}_{\mathrm{np}}^{*}$ for St=0.5 and St=0.2. Similarly, the curve of ${C}_{\mathrm{np}}^{*}$ for St1=0.5 and St2=1.0 is located between the power spectra ${E}_{\mathrm{np}}^{*}$ for St=0.5 and St=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 St. In order to focus on the influence of the spatial correlation of clusters, we have evaluated the coherence $\mathrm{coh}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)$, which is defined as

$\begin{array}{}\text{(30)}& \mathrm{coh}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)=\frac{|{C}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)|}{\sqrt{{E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}}\right){E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{2}}\right)}}.\end{array}$

Figure 3 shows the coherence between St1=0.5 and other Stokes numbers. The coherence $\mathrm{coh}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)$ 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 St 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 St1 and St2; the slope becomes steeper as the difference in St increases. These results indicate that the decreasing trend of coherence can be approximated by an exponential function; i.e.,

$\begin{array}{}\text{(31)}& \mathrm{coh}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)\approx \mathrm{exp}\left(-\mathit{\xi }/{\mathit{\xi }}_{\mathrm{c}}\right),\end{array}$

where ξc is the critical wave number normalized by the Kolmogorov scale, given by a function of St1 and St2. 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 St1 and St2, where St2>St1. ξc for the same St1 increases as St2 decreases, whereas ξc for the same St2 increases as St1 increases. This indicates that ξc increases as St1 and St2 becomes closer to each other. It should be noted that several studies discuss the spatial correlation of bidisperse clustering particles. For example, 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. 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, lc, 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, St2St1, and the Stokes number difference, St2St1, are candidates for the dominant parameter for ξc. Figure 5 shows ξc plots against St2St1 and St2St1. These figures clearly indicate that the Stokes number dependence of ξc is explained better by St2St1 than St2St1. This would be because both the critical wave number ξc and crossover length lc 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 lc. 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

$\begin{array}{}\text{(32)}& {\mathit{\xi }}_{\mathrm{c}}\left(S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)\approx \frac{\mathrm{0.191}}{|S{t}_{\mathrm{1}}-S{t}_{\mathrm{2}}|}.\end{array}$

This implies that ξc is closely related to the inverse of lclη because the crossover length of 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 are valid for St≪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 St≤2.0.

Figure 4Critical wave number ξc for a combination of Stokes numbers, St1 and St2. The symbol color and type indicate the combination of Stokes numbers. The black, red, green, blue, and light blue symbols are St1=0.05, 0.1, 0.2, 0.5, and 1.0, respectively. The square, circle, triangle, inverse triangle, and diamond symbols are St2=0.1, 0.2, 0.5, 1.0, and 2.0.

Figure 5(a) Critical wave number ξc against the Stokes number ratio St2St1, and (b) against the Stokes number difference, St2St1. Notations are as in Fig. 4.

## 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 qr(rp) and np are given. That is, the normalized cross spectrum ${C}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)$ is given by

$\begin{array}{}\text{(33)}& {C}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)=\mathrm{coh}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)\sqrt{{E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}}\right){E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{2}}\right)}.\end{array}$

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 }|St\right)$ was proposed by . The model equation is given by

$\begin{array}{}\text{(34)}& {E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|St\right)=\frac{{c}_{\mathrm{1}}{\mathit{\xi }}^{\mathit{\alpha }}}{{\left\{\mathrm{1}+{\left({c}_{\mathrm{1}}/{c}_{\mathrm{2}}\right)}^{\mathrm{2}\mathit{\gamma }/\left(\mathit{\alpha }-\mathit{\beta }\right)}{\mathit{\xi }}^{\mathrm{2}\mathit{\gamma }}\right\}}^{\left(\mathit{\alpha }-\mathit{\beta }\right)/\mathrm{2}\mathit{\gamma }}},\end{array}$

where c1, c2, α, β, and γ are the model parameters given by the functions of St as follows:

$\begin{array}{}\text{(35)}& \left\{\begin{array}{rcl}{c}_{\mathrm{1}}& =& \mathrm{13.4}/\left[\mathrm{1}+\left(St/\mathrm{0.29}{\right)}^{-\mathrm{1.25}}\right],\\ {c}_{\mathrm{2}}& =& \mathrm{6.7}S{t}^{\mathrm{1.6}}/\left[\mathrm{1}+\mathrm{0.68}S{t}^{\mathrm{3.7}}\right],\\ \mathit{\alpha }& =& \mathrm{0.44}-\mathrm{0.20}\mathrm{ln}St,\\ \mathit{\beta }& =& -\mathrm{1}+\mathrm{0.77}S{t}^{-\mathrm{1}}\mathrm{exp}\left[-\left(\mathrm{ln}St-\mathrm{0.55}{\right)}^{\mathrm{2}}/\mathrm{2.0}\right],\\ \mathit{\gamma }& =& \mathrm{1.6}\end{array}\right\.\end{array}$

confirmed that this parameterization is reliable for St≤2.0: in this range, the error of the parameterization is smaller than 1 dB.

Figure 6Critical wave number ξc against the inverse of Stokes number difference. The dashed line is the best-fitting curve to the ξc data. Other notations are as in Fig. 4.

In order to evaluate the reliability of the proposed cross spectrum model, the ${r}_{\mathrm{p}}^{\mathrm{3}}$-weighted power spectrum Er3np(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

$\begin{array}{}\text{(36)}& {E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi }\right)=\frac{{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}\left(k\right)}{〈{r}_{\mathrm{p}}^{\mathrm{3}}{〉}^{\mathrm{2}}〈{n}_{\mathrm{p}}{〉}^{\mathrm{2}}{l}_{\mathit{\eta }}},\end{array}$

where $〈{r}_{\mathrm{p}}^{\mathrm{3}}〉$ is given by $〈{r}_{\mathrm{p}}^{\mathrm{3}}〉=\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 }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)=\mathrm{1}$. The parameterization with $\mathrm{coh}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)=\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 }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)=\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 }). 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 erms in units of decibels. erms is defined as

$\begin{array}{}\text{(37)}& {e}_{\mathrm{rms}}=\frac{\mathrm{1}}{{\mathit{\xi }}_{\mathrm{max}}^{\prime }-{\mathit{\xi }}_{\mathrm{min}}^{\prime }}\underset{{\mathit{\xi }}_{\mathrm{min}}^{\prime }}{\overset{{\mathit{\xi }}_{\mathrm{max}}^{\prime }}{\int }}{\left\{{E}_{\mathrm{r}\mathrm{3}\mathrm{np},\mathrm{model}}^{*\mathrm{dB}}\left(\mathit{\xi }\right)-{E}_{\mathrm{r}\mathrm{3}\mathrm{np},\mathrm{DNS}}^{*\mathrm{dB}}\left(\mathit{\xi }\right)\right\}}^{\mathrm{2}}\mathrm{d}{\mathit{\xi }}^{\prime },\end{array}$

where ξ is defined as ${\mathit{\xi }}^{\prime }=\mathrm{ln}\mathit{\xi }$ and superscript dB denotes a value in units of decibels. erms was calculated for the wave number range relevant to radar observations; i.e., $\mathrm{0.05}\le \mathit{\xi }\le \mathrm{4.0}$. erms 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 , erms values for CUMA_eps400 and CUMA_eps1000 are negligibly small. erms 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, erms 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.

Figure 7Comparisons of ${r}_{\mathrm{p}}^{\mathrm{3}}$-weighted power spectrum ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi }\right)$ obtained from DNS data and estimated by the proposed cross spectrum model for the cases of (a) g=0.0 and (b) 9.8 m s−2.

## 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 . The large gravitational settling can modulate ${E}_{\mathrm{np}}^{*}\left(\mathit{\xi }|St\right)$, 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 }|St\right)$ is insignificant for Sv<3 .

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 }|St\right)$ in Eq. (33). and reported that gravitational settling enlarges the crossover length of the bidisperse RDF. extended the perturbation expansion analysis of and presented the formulation for the crossover length in the presence of gravity, which is

$\begin{array}{}\text{(38)}& \frac{{l}_{\mathrm{c}}}{{l}_{\mathit{\eta }}}={C}_{\mathrm{Chun}}\left|S{t}_{\mathrm{1}}-S{t}_{\mathrm{2}}\right|{\left[\mathrm{1}+\frac{\mathrm{1}}{\mathrm{3}{a}_{\mathrm{0}}}\left(\frac{{\mathit{\tau }}_{\mathrm{g}}}{{\mathit{\tau }}_{\mathrm{a}}}\right)F{r}^{-\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}},\end{array}$

where CChun is the coefficient derived by , a0 is the ratio of the acceleration variance a2 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. obtained the values of CChun≈5.0 and a0≈1.545 based on their DNS results. further assumed Sv1 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 lc of Eq. (38) becomes equivalent to that of when gravitational settling is negligibly small, whereas the gravity effect is dominant for the cases of Fr<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 }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)$. Since ξc is inversely proportional to lclη, we propose the following correction based on Eq. (38):

$\begin{array}{}\text{(39)}& {\mathit{\xi }}_{\mathrm{c}}\left(S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)=\frac{\mathrm{0.191}}{|S{t}_{\mathrm{1}}-S{t}_{\mathrm{2}}|}{\left[\mathrm{1}+\frac{\mathrm{1}}{\mathrm{3}{a}_{\mathrm{0}}}F{r}^{-\mathrm{2}}\right]}^{-\mathrm{1}/\mathrm{2}}.\end{array}$

Note that this study also adopts the value of a0 obtained by .

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 }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)=\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. erms 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, erms 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, Sv for the modal radius is smaller than unity. For the cases of Sv>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 . When the clustering structure is anisotropic, the influence on the radar reflectivity factor theoretically depends on the direction of microwave propagation.

5 Application to cloud simulation data

## 5.1 Cloud simulation data

We have applied the proposed model to the high-resolution cloud-simulation data of 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. 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

$\begin{array}{}\text{(40)}& G\left(y\right)\mathrm{d}y={n}_{\mathrm{p}}m\left({r}_{\mathrm{p}}\right){q}_{\mathrm{r}}\left({r}_{\mathrm{p}}\right)\mathrm{d}{r}_{\mathrm{p}},\end{array}$

where y=ln rp, and m(rp) is the mass of droplets with a radius of rp. The mass coordinate m and logarithmic coordinate y are discretized as

$\begin{array}{}\text{(41)}& & {m}_{k}={\mathrm{2}}^{\mathrm{1}/s}{m}_{k-\mathrm{1}}\text{(42)}& & {y}_{k}={y}_{k-\mathrm{1}}+\mathrm{d}y,\end{array}$

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, rp1, was 3 µm; thus, the representative radius of the 33rd bin (the largest droplet class) was rp33=4.9 mm. The prognostic variable for liquid water is the water mass content, Mk, 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 Mk 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 (Twomey1959). The activated droplets were added to the bins using the prescribed spectrum method (Soong1974). Details of the model configuration are described in . The model settings and computational conditions were based on the protocol of the RICO model intercomparison project (van Zanten et al.2006http://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}×\mathrm{12.8}×\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. 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

$\begin{array}{}\text{(43)}& Z={Z}_{\mathrm{incoh}}+{Z}_{\mathrm{PB}}+{Z}_{\mathrm{CB}},\end{array}$

where ZPB is the particulate Bragg-scattering part of Z in Eq. (9), and ZCB is an additional term reflecting clear-air Bragg scattering. In real clouds, ZPB is caused by droplet number density fluctuations due to turbulent droplet clustering and turbulent entrainment of environmental clear air. Er3np(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)$ , 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, ZPB 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.,

$\begin{array}{}\text{(44)}& \mathit{ϵ}={\left({C}_{\mathrm{s}}{\mathrm{\Delta }}_{\mathrm{s}}\right)}^{\mathrm{2}}{\left(\mathrm{2}\stackrel{\mathrm{^}}{{s}_{ij}}\stackrel{\mathrm{^}}{{s}_{ij}}\right)}^{\mathrm{3}/\mathrm{2}},\end{array}$

where Cs is the Smagorinsky constant (Cs=0.173 in this study), Δs the representative grid spacing given by ${\mathrm{\Delta }}_{\mathrm{s}}=\left({\mathrm{\Delta }}_{x}{\mathrm{\Delta }}_{y}{\mathrm{\Delta }}_{z}{\right)}^{\mathrm{1}/\mathrm{3}}$, and $\stackrel{\mathrm{^}}{{s}_{ij}}$ the strain rate tensor, which is given by $\stackrel{\mathrm{^}}{{s}_{ij}}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial \stackrel{\mathrm{^}}{{u}_{i}}}{\partial {x}_{j}}+\frac{\partial \stackrel{\mathrm{^}}{{u}_{j}}}{\partial {x}_{i}}\right)$, where $\stackrel{\mathrm{^}}{{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. 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 (klη<0.1) is also considered, since the diffusive coefficient Dnp for droplet number density is much smaller than ν. ${E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{entr}}\left(k\right)$ is approximately given by (Hill1978)

$\begin{array}{}\text{(45)}& \frac{{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{\mathrm{entr}}\left(k\right)}{{\mathit{\chi }}_{\mathrm{r}\mathrm{3}\mathrm{np}}{\mathit{ϵ}}^{-\mathrm{3}/\mathrm{4}}{\mathit{\nu }}^{\mathrm{5}/\mathrm{4}}}={C}_{\mathrm{c}}{\left(k{l}_{\mathit{\eta }}\right)}^{-\mathrm{5}/\mathrm{3}}{\left[{\left\{k{l}_{\mathit{\eta }}{\left({C}_{\mathrm{b}}/{C}_{\mathrm{c}}\right)}^{\mathrm{3}/\mathrm{2}}\right\}}^{\mathrm{2}{\mathit{\gamma }}^{\prime }}+\mathrm{1}\right]}^{\mathrm{1}/\mathrm{3}{\mathit{\gamma }}^{\prime }},\end{array}$

where χr3np is the scalar dissipation rate for $〈{r}_{\mathrm{p}}^{\mathrm{3}}〉{n}_{\mathrm{p}}$, Cc is the Obukhov–Corrsin constant (Cc=0.67) , Cb is the Batchelor constant (Cb=3.7) , and γ is the model parameter (${\mathit{\gamma }}^{\prime }=\mathrm{1.4}$).

Figure 8Three-dimensional visualization of cloud simulation data: (a) volume rendering of optical depth and (b) isosurfaces of (blue) the energy dissipation rate ϵ=100 cm2 s−3 and (yellow) the vertical velocity u3=3 m s−1.

The contribution of clear-air Bragg scattering was calculated by

$\begin{array}{}\text{(46)}& {Z}_{\mathrm{CB}}={\mathrm{2}}^{\mathrm{5}}|K{|}^{-\mathrm{2}}{\mathit{\kappa }}^{-\mathrm{2}}{E}_{n}\left(\mathit{\kappa }\right),\end{array}$

where En(k) is the power spectrum of refractive index fluctuations. The refractive index nref is given approximately by

$\begin{array}{}\text{(47)}& {n}_{\mathrm{ref}}=\mathrm{1}+\mathrm{3.73}×{\mathrm{10}}^{-\mathrm{1}}\frac{e}{{T}^{\mathrm{2}}}+\mathrm{7.76}×{\mathrm{10}}^{-\mathrm{5}}\frac{p}{T},\end{array}$

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 En(k) is given by the scalar concentration spectrum for Pr<1 (Pao1964), where $Pr\equiv \mathit{\nu }/D$ (D is the scalar diffusive coefficient) is the Prandtl number:

$\begin{array}{}\text{(48)}& \frac{{E}_{n}\left(k\right)}{{\mathit{\chi }}_{n}{\mathit{ϵ}}^{-\mathrm{3}/\mathrm{4}}{\mathit{\nu }}^{\mathrm{5}/\mathrm{4}}}={C}_{\mathrm{c}}{\left(k{l}_{\mathit{\eta }}\right)}^{-\mathrm{5}/\mathrm{3}}\mathrm{exp}\left[-\mathrm{1.5}\frac{{C}_{\mathrm{c}}}{Pr}{\left(k{l}_{\mathit{\eta }}\right)}^{\mathrm{4}/\mathrm{3}}\right],\end{array}$

where χn is the scalar dissipation rate for nref. In this study, the Prandtl number of refractive index was set to 0.7.

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

$\begin{array}{}\text{(49)}& {\mathit{\chi }}_{\mathit{\theta }}=\mathrm{2}\frac{{\mathit{\nu }}_{\mathrm{t}}}{S{c}_{\mathrm{t}}}\frac{\partial \stackrel{\mathrm{^}}{\mathit{\theta }}}{\partial {x}_{i}}\frac{\partial \stackrel{\mathrm{^}}{\mathit{\theta }}}{\partial {x}_{i}},\end{array}$

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

## 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 }〈{r}_{\mathrm{p}}^{\mathrm{2}}〉{n}_{\mathrm{p}}{\mathrm{\Delta }}_{z}$, where Qext is the extinction efficiency for Mie scattering (Qext=2.0 in this study), and $〈{r}_{\mathrm{p}}^{\mathrm{2}}〉$ is given by $〈{r}_{\mathrm{p}}^{\mathrm{2}}〉=\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 u3. 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.

Figure 9Liquid water content (LWC), energy dissipation rate ϵ, and radar reflectivity factors for S-band microwaves in the vertical cross section: (a) LWC in a logarithmic scale (i.e., log 10LWC – g m−3), (b) ϵ (m2 s−3), (c) ZdB (dBZ), (d) ${Z}^{\mathrm{dB}}-\left({Z}_{\mathrm{incoh}}+{Z}_{\mathrm{CB}}{\right)}^{\mathrm{dB}}$ (dB), and (e) ${Z}_{\mathrm{PBc}}^{\mathrm{dB}}$ (dBZ). The solid lines in (c)(e) indicate the isoline of ${Z}_{\mathrm{incoh}}^{\mathrm{dB}}=-\mathrm{18}$ dBZ.

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 ZdB, the increase of ZdB 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 fm=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 ZdB (dBZ)=10log 10Z (mm6 m−3). Large values of ZdB 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 Zincoh is proportional to $〈{r}_{\mathrm{p}}^{\mathrm{6}}〉$. 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., ZCB. The radar echo layer at the height from 2.2 to 2.5 km is caused by ZCB due to a large humidity gap in the inversion layer.

Figure 9c does not show a clear sign of the mantle echo reported by . 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., ZPBe) because the large-scale cloud water inhomogeneity at cloud edges produces small-scale fluctuations due to the turbulent cascade . The influence of turbulent entrainment on ZdB 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}}-\left({Z}_{\mathrm{incoh}}+{Z}_{\mathrm{CB}}{\right)}^{\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 cm2 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 ZPBc 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}{\right)}^{\mathrm{2}}{l}_{\mathit{\eta }}{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\kappa }{l}_{\mathit{\eta }}\right)$. Thus, the significant influence of turbulent clustering is caused by sufficiently large values of the energy dissipation rate and the LWC.

6 Conclusions

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, St1 and St2, has values between the power spectra for St1 and St2 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 . 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 Sv1, where Sv 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 . 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 . 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.

Data availability
Data availability.

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

Author contributions
Author contributions.

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.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

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

References

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