Turbulent enhancement of radar reflectivity factor for polydisperse cloud droplets

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 (mm⁶ m⁻³), Z, which is obtained with the following radar equation:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{P}_{\mathrm{r}}={\displaystyle \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 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, 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 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:

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

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

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

where n(x, r_p)dr_pdx is the number of droplets with radii from r_p 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)}& {\displaystyle }{\displaystyle }Z={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K{|}^{\mathrm{2}}V{\sin }^{\mathrm{2}}\mathit{\chi }}}\mathit{\sigma },\end{array}$$

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

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

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

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

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, κ=2k_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(x, $r_{\mathrm{p}})=\stackrel{\mathrm{‾}}{n(\mathit{x},r_{\mathrm{p}})}+\mathit{\delta }n(\mathit{x},r_{\mathrm{p}})$. The temporal-average part, $\stackrel{\mathrm{‾}}{n(\mathit{x},r_{\mathrm{p}})}$, contributes to the separated reflection; therefore, the contribution of this part is negligibly small when $\stackrel{\mathrm{‾}}{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{‾}}{n(\mathit{x},r_{\mathrm{p}})}$. Then, we obtain

$$\begin{array}{ll}{\displaystyle }Z=& {\displaystyle }{\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{\mathit{r}}{\int }〈\mathit{\delta }n\left(\mathit{x},r_{\mathrm{p}}\right)\mathit{\delta }n\left(\mathit{x}+\mathit{r},r_{\mathrm{p}}^{\prime }\right)〉\right.\\ \text{(7)}& {\displaystyle }& {\displaystyle }\left.\exp \left(-i\mathit{\kappa }\cdot \mathit{r}\right)\mathrm{d}\mathit{r}\right\}\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 $\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{‾}}{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{‾}}{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{‾}}{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 δ(r) and $\mathit{\delta }(r_{\mathrm{p}}^{\prime }-r_{\mathrm{p}})$ are the Dirac delta functions. Thus, the spatial correlation function is given by

$$\begin{array}{ll}{\displaystyle }〈\mathit{\delta }n\left(\mathit{x},r_{\mathrm{p}}\right)\mathit{\delta }n\left(\mathit{x}+\mathit{r},r_{\mathrm{p}}^{\prime }\right)〉& {\displaystyle }=\langle n_{\mathrm{p}}\rangle \mathit{\delta }\left(\mathit{r}\right)q_{\mathrm{r}}\left(r_{\mathrm{p}}\right)\mathit{\delta }\left(r_{\mathrm{p}}^{\prime }-r_{\mathrm{p}}\right)\\ \text{(8)}& {\displaystyle }& {\displaystyle }+\mathrm{\Psi }\left(\mathit{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 〈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

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }Z={\mathrm{2}}^{\mathrm{6}}\langle r_{\mathrm{p}}^{\mathrm{6}}\rangle \langle n_{\mathrm{p}}\rangle +{\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 }=\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

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{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\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{̃}}{\mathrm{\Psi }}(\mathit{k}|r_{\mathrm{p}},r_{\mathrm{p}}^{\prime })\mathrm{d}{\mathit{\sigma }}_k$; i.e., the integration of $\stackrel{\mathrm{̃}}{\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{̃}}{\mathrm{\Psi }}\left(\mathit{k}\right|r_{\mathrm{p}},r_{\mathrm{p}}^{\prime })$ is the cross spectral density function, defined as

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathit{k}|r_{\mathrm{p}},r_{\mathrm{p}}^{\prime }\right)={\displaystyle \frac{\mathrm{1}}{(\mathrm{2}\mathit{\pi }{)}^{\mathrm{3}}}}\underset{\mathit{r}}{\int }\mathrm{\Psi }\left(\mathit{r}|r_{\mathrm{p}},r_{\mathrm{p}}^{\prime }\right)\exp \left(-i\mathit{k}\cdot \mathit{r}\right)\mathrm{d}\mathit{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}}\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

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }Z={\mathrm{2}}^{\mathrm{6}}{r}_{\mathrm{p}\mathrm{1}}^{\mathrm{6}}\langle n_{\mathrm{p}}\rangle +{\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 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{̃}}{\mathrm{\Phi }}(\mathit{k}\left|r_{\mathrm{p}}\right)\mathrm{d}{\mathit{\sigma }}_k$, where $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathit{k}\right|r_{\mathrm{p}})$ is the power spectral density function of n_p(x|r_p), defined as

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathit{k}|r_{\mathrm{p}}\right)={\displaystyle \frac{\mathrm{1}}{(\mathrm{2}\mathit{\pi }{)}^{\mathrm{3}}}}\underset{\mathit{r}}{\int }\mathrm{\Psi }\left(\mathit{r}|r_{\mathrm{p}},r_{\mathrm{p}}\right)\exp \left(-i\mathit{k}\cdot \mathit{r}\right)\mathrm{d}\mathit{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)}& {\displaystyle }& {\displaystyle \frac{\partial u_i}{\partial x_i}}=\mathrm{0},\text{(15)}& {\displaystyle }& {\displaystyle \frac{\partial u_i}{\partial t}}+{\displaystyle \frac{\partial u_i{u}_j}{\partial x_j}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}}{\displaystyle \frac{\partial p}{\partial x_i}}+\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_i}{\partial x_j\partial x_j}}+F_i,\end{array}$$

where u_i is the flow velocity in the ith 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₀ 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)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{v}_i}{\mathrm{d}t}}=-{\displaystyle \frac{v_i-u_i}{{\mathit{\tau }}_{\mathrm{p}}}}+g_i,\end{array}$$

where v_i and g_i 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)}& {\displaystyle }{\displaystyle }{\mathit{\tau }}_{\mathrm{p}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{p}}}{{\mathit{\rho }}_{\mathrm{a}}}}{\displaystyle \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πL₀. 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³. 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 (Matsuda et al., 2014) (see Sect. 3.2). The kinematic viscosity, ν, was set to $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹, 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⁷.

Table 1 Computational 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 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² s⁻³, 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² s⁻³ 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 cm² s⁻³, 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=τ_pg is the terminal settling velocity and $u_{\mathit{\eta }}\equiv (\mathit{\nu }\mathit{ϵ}{)}^{\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² s⁻³ were 0.67, 0.48, and 0.38, respectively.

Figure 1 Spatial distribution of droplets for (orange) St=0.2, (green) 0.5, and (blue) 1.0 in the regions of (a) 2πL₀×2πL₀ and (b) 0.5πL₀×0.5πL₀. Droplets located in the range of $\mathrm{0}<z<\mathrm{4}{l}_{\mathit{\eta }}$ 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(\mathit{k}\right|r_{\mathrm{p}})$ and the cross spectral density function $\stackrel{\mathrm{̃}}{\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:

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathit{k}|r_{\mathrm{p}}\right)=L_{\mathrm{0}}^{-\mathrm{3}}〈\stackrel{\mathrm{̃}}{n_{\mathrm{p}}}\left(\mathit{k}|r_{\mathrm{p}}\right)\stackrel{\mathrm{̃}}{n_{\mathrm{p}}}\left(-\mathit{k}|r_{\mathrm{p}}\right)〉,\text{(19)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathit{k}|r_{\mathrm{p}\mathrm{1}},r_{\mathrm{p}\mathrm{2}}\right)=L_{\mathrm{0}}^{-\mathrm{3}}〈\stackrel{\mathrm{̃}}{n_{\mathrm{p}}}\left(\mathit{k}|r_{\mathrm{p}\mathrm{1}}\right)\stackrel{\mathrm{̃}}{n_{\mathrm{p}}}\left(-\mathit{k}|r_{\mathrm{p}\mathrm{2}}\right)〉,\end{array}$$

where $\stackrel{\mathrm{̃}}{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

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

where x_p,j is the position vector of the jth 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

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

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

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathit{k}|r_{\mathrm{p}}\right)}{\langle n_{\mathrm{p}}{\rangle }^{\mathrm{2}}{L}_{\mathrm{0}}^{\mathrm{3}}}}\rangle ={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{p}}^{\mathrm{2}}}}〈\sum _{j=\mathrm{1}}^{N_{\mathrm{p}}}\exp \left(-i\mathit{k}\cdot {\mathit{x}}_{\mathrm{p},j}\right)\sum _{j^{\prime }=\mathrm{1},j^{\prime }\ne j}^{N_{\mathrm{p}}}\exp \left(i\mathit{k}\cdot {\mathit{x}}_{\mathrm{p},j^{\prime }}\right)〉\text{(23)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{p}}^{\mathrm{2}}}}\left[〈{\left\{\sum _{j=\mathrm{1}}^{N_{\mathrm{p}}}\cos \left(\mathit{k}\cdot {\mathit{x}}_{\mathrm{p},j}\right)\right\}}^{\mathrm{2}}〉+〈{\left\{\sum _{j=\mathrm{1}}^{N_{\mathrm{p}}}\sin \left(\mathit{k}\cdot {\mathit{x}}_{\mathrm{p},j}\right)\right\}}^{\mathrm{2}}〉\right]-{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{p}}}}.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathit{k}|r_{\mathrm{p}\mathrm{1}},r_{\mathrm{p}\mathrm{2}}\right)}{\langle n_{\mathrm{p}\mathrm{1}}\rangle \langle n_{\mathrm{p}\mathrm{2}}\rangle L_{\mathrm{0}}^{\mathrm{3}}}}=& {\displaystyle \frac{\mathrm{1}}{N_{\mathrm{p}\mathrm{1}}{N}_{\mathrm{p}\mathrm{2}}}}〈\sum _{j=\mathrm{1}}^{N_{\mathrm{p}\mathrm{1}}}\exp \left(-i\mathit{k}\cdot {\mathit{x}}_{\mathrm{p}\mathrm{1},j}\right)\\ \text{(24)}& {\displaystyle }& {\displaystyle }\sum _{j^{\prime }=\mathrm{1}}^{N_{\mathrm{p}\mathrm{2}}}\exp \left(i\mathit{k}\cdot {\mathit{x}}_{\mathrm{p}\mathrm{2},j^{\prime }}\right)〉\end{array}$$

$$\begin{array}{ll}{\displaystyle }={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{p}\mathrm{1}}{N}_{\mathrm{p}\mathrm{2}}}}& {\displaystyle }\left\{〈\sum _{j=\mathrm{1}}^{N_{\mathrm{p}\mathrm{1}}}\cos \left(\mathit{k}\cdot {\mathit{x}}_{\mathrm{p}\mathrm{1},j}\right)\sum _{j^{\prime }=\mathrm{1}}^{N_{\mathrm{p}\mathrm{2}}}\cos \left(\mathit{k}\cdot {\mathit{x}}_{\mathrm{p}\mathrm{2},j^{\prime }}\right)〉\right.\\ \text{(25)}& {\displaystyle }& {\displaystyle }\left.+〈\sum _{j=\mathrm{1}}^{N_{\mathrm{p}\mathrm{1}}}\sin \left(\mathit{k}\cdot {\mathit{x}}_{\mathrm{p}\mathrm{1},j}\right)\sum _{j^{\prime }=\mathrm{1}}^{N_{\mathrm{p}\mathrm{2}}}\sin \left(\mathit{k}\cdot {\mathit{x}}_{\mathrm{p}\mathrm{2},j^{\prime }}\right)〉\right\},\end{array}$$

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 jth 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{̃}}{\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⁻⁴), which is caused by the statistical and truncation errors.

The spectral density functions, $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathit{k}\right|r_{\mathrm{p}})$ and $\stackrel{\mathrm{̃}}{\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₁, h₂, and h₃ 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₀. 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:

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{np}}\left(k|r_{\mathrm{p}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }k}}\sum _{k-\mathrm{\Delta }k/\mathrm{2}\le \left|\mathit{k}\right|<k+\mathrm{\Delta }k/\mathrm{2}}\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathit{k}|r_{\mathrm{p}}\right),\text{(27)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{np}}\left(k|r_{\mathrm{p}\mathrm{1}},r_{\mathrm{p}\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }k}}\sum _{k-\mathrm{\Delta }k/\mathrm{2}\le \left|\mathit{k}\right|<k+\mathrm{\Delta }k/\mathrm{2}}\stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathit{k}|r_{\mathrm{p}\mathrm{1}},r_{\mathrm{p}\mathrm{2}}\right).\end{array}$$

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 St cases.

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

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{ϵ}}^{-\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}<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:

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }{E}_{\mathrm{np}}^{*}\left(\mathit{\xi }\right|St)=E_{\mathrm{np}}^{*}\left(k{l}_{\mathit{\eta }}|St\right)\equiv {\displaystyle \frac{E_{\mathrm{np}}\left(k|r_{\mathrm{p}}\right)}{\langle n_{\mathrm{p}}{\rangle }^{\mathrm{2}}{l}_{\mathit{\eta }}}},\end{array}$$

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }{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 {\displaystyle \frac{C_{\mathrm{np}}\left(k|r_{\mathrm{p}\mathrm{1}},r_{\mathrm{p}\mathrm{2}}\right)}{\langle n_{\mathrm{p}\mathrm{1}}\rangle \langle n_{\mathrm{p}\mathrm{2}}\rangle l_{\mathit{\eta }}}},\end{array}$$

where ξ is the normalized wave number defined as ξ≡kl_η and St, St₁, and St₂ are the Stokes numbers of water droplets with radii of r_p, r_p1, and r_p2.

Figure 3 Coherence of cross spectra for combinations of St₁=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}<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₀×2πL₀, while Fig. 1b is the magnified view for the region with a size of 0.5πL₀×0.5πL₀. 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 }\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⁻⁴). 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 St₁=0.5 and St₂=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 St₁=0.5 and St₂=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 }\right|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}})$, which is defined as

$$\begin{array}{}\text{(30)}& {\displaystyle }{\displaystyle }\mathrm{coh}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)={\displaystyle \frac{\left|C_{\mathrm{np}}^{*}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)\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 St₁=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 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 St₁ and St₂; 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)}& {\displaystyle }{\displaystyle }\mathrm{coh}\left(\mathit{\xi }|S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)\approx \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 St₁ and St₂. 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 St₁ and St₂, where St₂>St₁. ξ_c for the same St₁ increases as St₂ decreases, whereas ξ_c for the same St₂ increases as St₁ increases. This indicates that ξ_c increases as St₁ and St₂ 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, St₂∕St₁, and the Stokes number difference, St₂–St₁, are candidates for the dominant parameter for ξ_c. Figure 5 shows ξ_c plots against St₂∕St₁ and St₂–St₁. These figures clearly indicate that the Stokes number dependence of ξ_c is explained better by St₂–St₁ than St₂∕St₁. 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

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }{\mathit{\xi }}_{\mathrm{c}}\left(S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)\approx {\displaystyle \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 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 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 4 Critical wave number ξ_c for a combination of Stokes numbers, St₁ and St₂. The symbol color and type indicate the combination of Stokes numbers. The black, red, green, blue, and light blue symbols are St₁=0.05, 0.1, 0.2, 0.5, and 1.0, respectively. The square, circle, triangle, inverse triangle, and diamond symbols are St₂=0.1, 0.2, 0.5, 1.0, and 2.0.

Figure 5 (a) Critical wave number ξ_c against the Stokes number ratio St₂∕St₁, and (b) against the Stokes number difference, St₂–St₁. 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 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

$$\begin{array}{}\text{(33)}& {\displaystyle }{\displaystyle }{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 }\right|St)$ was proposed by Matsuda et al. (2014). The model equation is given by

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

where c₁, c₂, α, β, and γ are the model parameters given by the functions of St as follows:

$$\begin{array}{}\text{(35)}& {\displaystyle }{\displaystyle }\left\{\begin{array}{rcl}{c}_{\mathrm{1}}& =& \mathrm{13.4}/\left[\mathrm{1}+(St/\mathrm{0.29}{)}^{-\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}\ln St,\\ \mathit{\beta }& =& -\mathrm{1}+\mathrm{0.77}S{t}^{-\mathrm{1}}\exp \left[-(\ln St-\mathrm{0.55}{)}^{\mathrm{2}}/\mathrm{2.0}\right],\\ \mathit{\gamma }& =& \mathrm{1.6}\end{array}\right..\end{array}$$

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

Figure 6 Critical 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 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

$$\begin{array}{}\text{(36)}& {\displaystyle }{\displaystyle }{E}_{\mathrm{r}\mathrm{3}\mathrm{np}}^{*}\left(\mathit{\xi }\right)={\displaystyle \frac{E_{\mathrm{r}\mathrm{3}\mathrm{np}}\left(k\right)}{\langle r_{\mathrm{p}}^{\mathrm{3}}{\rangle }^{\mathrm{2}}\langle n_{\mathrm{p}}{\rangle }^{\mathrm{2}}{l}_{\mathit{\eta }}}},\end{array}$$

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

$$\begin{array}{}\text{(37)}& {\displaystyle }{\displaystyle }{e}_{\mathrm{rms}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\xi }}_{\max }^{\prime }-{\mathit{\xi }}_{\min }^{\prime }}}\underset{{\mathit{\xi }}_{\min }^{\prime }}{\overset{{\mathit{\xi }}_{\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 }=\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.

Figure 7 Comparisons 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⁻².

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

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

where C_Chun is the coefficient derived by Chun et al. (2005), a₀ is the ratio of the acceleration variance 〈a²〉 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₀≈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 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 }\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):

$$\begin{array}{}\text{(39)}& {\displaystyle }{\displaystyle }{\mathit{\xi }}_{\mathrm{c}}\left(S{t}_{\mathrm{1}},S{t}_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{0.191}}{|S{t}_{\mathrm{1}}-S{t}_{\mathrm{2}}|}}{\left[\mathrm{1}+{\displaystyle \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 a₀ 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 Application to cloud simulation data

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

$$\begin{array}{}\text{(40)}& {\displaystyle }{\displaystyle }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 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

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

where $\mathrm{d}y=\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

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

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.,

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

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 (kl_η<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)

$$\begin{array}{}\text{(45)}& {\displaystyle }{\displaystyle \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 $\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}$).

Figure 8 Three-dimensional visualization of cloud simulation data: (a) volume rendering of optical depth and (b) isosurfaces of (blue) the energy dissipation rate ϵ=100 cm² s⁻³ and (yellow) the vertical velocity u₃=3 m s⁻¹.

The contribution of clear-air Bragg scattering was calculated by

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

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)

$$\begin{array}{}\text{(47)}& {\displaystyle }{\displaystyle }{n}_{\mathrm{ref}}=\mathrm{1}+\mathrm{3.73}\times {\mathrm{10}}^{-\mathrm{1}}{\displaystyle \frac{e}{T^{\mathrm{2}}}}+\mathrm{7.76}\times {\mathrm{10}}^{-\mathrm{5}}{\displaystyle \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 E_n(k) is given by the scalar concentration spectrum for Pr<1 (Pao, 1964), where $Pr\equiv \mathit{\nu }/D$ (D is the scalar diffusive coefficient) is the Prandtl number:

$$\begin{array}{}\text{(48)}& {\displaystyle }{\displaystyle \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}}\exp \left[-\mathrm{1.5}{\displaystyle \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 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

$$\begin{array}{}\text{(49)}& {\displaystyle }{\displaystyle }{\mathit{\chi }}_{\mathit{\theta }}=\mathrm{2}{\displaystyle \frac{{\mathit{\nu }}_{\mathrm{t}}}{S{c}_{\mathrm{t}}}}{\displaystyle \frac{\partial \widehat{\mathit{\theta }}}{\partial x_i}}{\displaystyle \frac{\partial \widehat{\mathit{\theta }}}{\partial x_i}},\end{array}$$

where ν_t is the eddy kinematic viscosity, Sc_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 Sc_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₃. 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 9 Liquid 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 ₁₀LWC – g m⁻³), (b) ϵ (m² s⁻³), (c) Z^dB (dBZ), (d) $Z^{\mathrm{dB}}-(Z_{\mathrm{incoh}}+Z_{\mathrm{CB}}{)}^{\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 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 ₁₀Z (mm⁶ m⁻³). 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⁻³, and the energy dissipation rate ϵ is intermittently larger than 100 cm² s⁻³. 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⁻³ 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.

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, St₁ and St₂, has values between the power spectra for St₁ and St₂ 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.