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 simu5 lation (DNS) of particle-laden homogeneous isotropic turbulence. The results show that the coherence of the cross spectrum is close to unity for small wavenumbers and decreases almost exponentially with increasing wavenumber. This decreasing trend is dependent on the combination of Stokes numbers. A critical wavenumber is introduced to characterize the exponential decrease of the coherence and parametrized using the Stokes number difference. Comparison with DNS results confirms that the proposed model can reproduce the r p-weighted power spectrum, which is proportional to the clustering influence on the 10 radar reflectivity factor, to a sufficiently high accuracy. 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 for the near-top of turbulent clouds.

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 3 p -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 ob-served at the near-top of the clouds, where the liquid water content and the energy dissipation rate are sufficiently large.

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 6 m −3 ), Z, which is obtained with the following radar equation: where P r and P t are the received and transmitted microwave powers, G is the antenna gain, A e is the effective aperture of the antenna, k m is the microwave wave number, K is the dielectric coefficient of a water droplet, V is the measurement volume, and R is the distance between the antenna and the cloud.The relationship between the radar reflectivity factor and cloud microphysical properties is usually expressed based on the assumption of incoherent scattering (Gossard and Strauch, 1983).Incoherent scattering implies random and uniform dispersion of cloud droplets (Bohren and Huffman, 1983).For this case, the factor is proportional to the sum of the scattering intensities of the individual droplets in the measurement volume.In contrast, droplets form a nonuniform spatial distribution in turbulence; i.e., inertial particles concentrate in small-enstrophy regions during turbulence due to the centrifugal motion (Maxey, 1987;Squires and Eaton, 1991;Wang and Maxey, 1993;Chen et al., 2006).
This preferential concentration is often referred to as turbulent clustering.Nonuniform distribution of cloud droplets results in coherent scattering, which is also referred to as particulate Bragg scattering (Kostinski and Jameson, 2000).For this case, the interference of microwaves scattered by nonuniformly distributed droplets increases the scattered microwave intensity; i.e., the radar reflectivity factor increases due to particulate Bragg scattering.It should be noted that Bragg scattering often indicates coherent scattering due to a nonuniform distribution of the refractive index of clear air.
In the troposphere, turbulent mixing of temperature and water vapor causes spatial variation of the refractive index.To distinguish this effect from the particulate Bragg scattering, it is specifically referred to as clear-air Bragg scattering.The radar reflectivity factor for both types of Bragg scattering is dependent on the microwave frequency f m (= k m c/2π , 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 −5/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)  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.

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 Rayleighscattering approximation: where E inc is the electric field amplitude vector of the incident microwave; r p is the droplet radius; x, x t , and x r are the position vectors of the droplet, microwave transmitter, and microwave receiver; k inc and k sca are the wave number vectors of the incident and scattered microwaves which satisfy |k inc | = |k sca | = k m , and χ is the angle between E inc and k sca .Considering the droplet-size dependence of E sca (t, x, r p ), the electric power of the microwave scattered by a group of droplets, P s , is given by where n(x, r p )dr p dx is the number of droplets with radii from r p to r p +dr 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 where λ is the microwave wavelength and σ is the radar cross section (Gossard and Strauch, 1983), which is defined as where P o is the electric power of the incident microwave, which is given by , and (5) into Eq.(4) yields where the wave number vector κ is defined as κ = k inc − k sca .Note that radar remote sensing typically uses backward scattering; i.e., k sca = −k 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 p ) = n(x, r p )+δn(x, r p ).The temporalaverage part, n(x, r p ), contributes to the separated reflection; therefore, the contribution of this part is negligibly small when n(x, r p ) has no fluctuation on a spatial scale of half the wavelength (Erkelens et al., 2001).Thus, we neglect the contribution of n(x, r p ).Then, we obtain where the angular brackets represent a temporal and spatial average in the measurement volume.
In order to decompose the spatial correlation function δn(x, r p )δn(x + r, r p ) , 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 r (r p ) ≡ 1 for r p = r p is discontinuous at r = 0 because the droplet distribution is composed of spatially discrete points.The singularity is given by n(x, r p ) δ(r)δ(r p − r p ), where δ(r) and δ(r p − r p ) are the Dirac delta functions.Thus, the spatial correlation function is given by δn x, r p δn x + r, r p = n p δ(r)q r r p δ r p − r p + r|r p , r p q r r p q r r p , (8) where n p is the average number density ( n p ≡ N p /V ) and (r|r p , r p ) is defined as the continuous part of δn p (x|r p )δn p (x + r|r p ) .Substitution of Eq. (8) into Eq.( 7) and adoption of the isotropic clustering assumption Gossard and Strauch (1983) yield where κ is κ = |κ|, r 6 p is given by r 6 p = ∞ 0 r 6 p q r (r p )dr p , and E r3np (k) is the r 3 p -weighted power spectrum, defined as The first and second terms on the right-hand side of Eq. ( 9) are the incoherent and coherent scattering parts; particulate Bragg scattering is caused by the second term.Equations ( 9) and ( 10) imply that the particulate Bragg-scattering part of Z for an arbitrary droplet size distribution can be calculated using the cross spectrum for bidisperse droplet size conditions.When droplets are distributed randomly and uniformly, the second term equals zero.Thus, the radar reflectivity factor when assuming a random and uniform droplet distribution is given by the first term; i.e., Z incoh = 2 6 r 6 p n 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 r p1 , the PDF of droplet radius is given by q r (r p ) = δ(r p − r p1 ).Then, the radar reflectivity factor Z is given by where E np (k|r p ) is the power spectrum of number density fluctuations, which satisfies where (k|r p ) is the power spectral density function of n p (x|r p ), defined as 3 Computational method

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.Threedimensional incompressible turbulent airflows were calculated by solving the continuity and Navier-Stokes equations: 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 HS-MAC method (Hirt and Cook, 1972) was used for velocitypressure coupling.The external forcing was applied to maintain the intensity of large-scale eddies for wave numbers k in the range |kL 0 | < 2 (Onishi et al., 2011), where L 0 is the representative length scale.Droplet motions were simulated by Lagrangian pointparticle tracking.Here, we assumed that the droplet density ρ p is much larger than ρ a and the drag term is given based on the Stokes law.The droplet motions were tracked by where v i and g i are the particle velocity and gravitational acceleration in the ith direction.τ p is the droplet relaxation time, which is given by The effects of turbulent modulation and droplet collision were neglected for simplicity because these effects were typically small in the timescale of τ η in clouds.The computational domain was set as a cube with edge lengths of 2π L 0 .Periodic boundary conditions were applied in all three directions.A uniform staggered grid was used for discretization.The number of grid points was set to 512 3 .A Taylor microscale-based turbulent Reynolds number of the obtained flow was Re λ = 204, where Re λ is defined as Re λ ≡ l λ u /ν, 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 1.5 × 10 −5 m 2 s −1 , and the ratio of the droplet density to the air density, ρ p /ρ a , was set to 840, assuming 1 atm and 298 K.The total number of droplets, N p , was set to 1.5 × 10 7 .
For this study, we performed the DNS for monodisperse and polydisperse droplets.Table 1 shows the computational settings for turbulence, droplet size, and gravitational acceleration.For monodisperse droplets, the droplet motions in an identical turbulent flow field were calculated for six values of Stokes number, 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 2 s −3 , which can be observed in cumulus and cumulonimbus clouds (Pinsky et al., 2008).The data for the polydisperse droplet cases were used to discuss the reliability of the proposed cross spectrum model.It should be noted that the droplet size distribution for the polydisperse cases were identical but the Stokes number histograms were different; the Stokes numbers corresponding to the modal radius (10.4 µm) for = 100, 400, and 1000 cm 2 s −3 were 0.035, 0.069, and 0.10, respectively.The gravitational acceleration g ≡ √ g i g i was set to zero for the monodisperse cases.The DNS for polydisperse droplets were performed under the conditions with and without gravitational settling.The Froude numbers, F r (≡ a η /g, in which a η ≡ 3/4 ν −1/4 is the Kolmogorov acceleration), for the cases with gravitational settling were 0.0520, 0.145, and 0.289 for = 100, 400, and 1000 cm 2 s −3 , respectively.The influence of gravitational settling on turbulent particle clustering is often discussed using the settling parameter, S v , which is defined as S v ≡ v T /u η (e.g., Wang and Maxey, 1993;Grabowski and Vaillancourt, 1999;Bec et al., 2014;Ireland et al., 2016), where v T = τ p g is the terminal settling velocity and u η ≡ (ν ) 1/4 is the Kolmogorov velocity.Note that S v satisfies S v = St/F r.The settling parameters corresponding to the modal radius for = 100, 400, and 1000 cm 2 s −3 were 0.67, 0.48, and 0.38, respectively.

Computation of power spectrum and cross spectrum
The power spectral density function (k|r p ) and the cross spectral density function (k|r p1 , r p2 ) are calculated from the Lagrangian droplet distribution data as follows: where n p (k|r p ) is the Fourier component of the droplet number density distribution, n p (x|r p ), and the angle brackets denote an ensemble average.The number density distribution for Lagrangian discrete droplets is given by where x p,j is the position vector of the j th droplet with radius r p , and N p is the total number of droplets with radius r p .The Fourier components of Eq. ( 20) are then given by Substitution of Eq. ( 21) into Eq.( 18) yields Similarly, substitution of Eq. ( 21) into Eq.( 19) yields where n p1 and n p2 are the average number density of droplets with radii of r p1 and r p2 (i.e., n p1 ≡ n p (x|r p1 ) and n p2 ≡ n p (x|r p2 ) ), respectively, N p1 and N p2 are the numbers of droplets with radii of r p1 and r p2 , respectively, x p1,j is the position of the j th droplet with a radius of r p1 , and x p2,j is the position of the j th droplet with a radius of r p2 .It should be noted that the imaginary part of (k|r p1 , r p2 ) is omitted in Eq. ( 25) because, statistically, it should be zero.We confirmed that the imaginary part of C np (k|r p1 , r p2 ) calculated from the DNS data is O(10 −4 ), which is caused by the statistical and truncation errors.
The spectral density functions, (k|r p ) and (k|r p1 , r p2 ), were calculated for discrete wave numbers kL 0 = (h 1 , h 2 , h 3 ), where h 1 , h 2 , and h 3 are arbitrary integers that satisfy k − k/2 ≤ |k| < k + k/2, where k was set to 1/L 0 .E np (k|r p ) and C np (k|r p1 , r p2 ) were then obtained by the following equations: The spectra, E np (k|r p ) and C np (k|r p1 , r p2 ), 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 0 = L 0 /U 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.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 η ≡ ν 3/4 −1/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 (0.05 < kl η < 4.0).Thus, this study used the same normalization for E np (k|r p ) and C np (k|r p1 , r p2 ) as follows: where ξ is the normalized wave number defined as ξ ≡ kl η and St, St 1 , and St 2 are the Stokes numbers of water droplets with radii of r p , r p1 , and r p2 .

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 0 < z < 4l η are indicated.The droplet position data were sampled at the same time step; i.e., the background turbulent flow field is identical.Figure 1a is the overall view of the region with a size of 2π L 0 ×2π L 0 , while Fig. 1b is the magnified view for Figure 3 shows the coherence between St 1 = 0.5 and other Stokes numbers.The coherence coh(ξ |St 1 , St 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 1 and St 2 ; 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., coh where ξ c is the critical wave number normalized by the Kolmogorov scale, given by a function of St 1 and St 2 .In this study, the critical wave number was obtained by finding the best-fitting exponential curve to the coherence for each combination of Stokes numbers.Figure 4 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 c /l η ≈ 5.0|St 1 − St 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.

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 (33) 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 * np (ξ |St) was proposed by Matsuda et al. (2014).The model equation is given by where c 1 , c 2 , α, β, and γ are the model parameters given by the functions of St as follows: (35) 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.
In order to evaluate the reliability of the proposed cross spectrum model, the r 3 p -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 3 p -weighted power spectrum, which is normalized as where r 3 p is given by r 3 p = ∞ 0 r 3 p q r (r p )dr p .The dashed lines show E * r3np (ξ ) predicted using the proposed parameterization including Eq. ( 31), while the dashed-dotted lines show those predicted by assuming perfect coherence; i.e., coh(ξ |St 1 , St 2 ) = 1.The parameterization with coh(ξ |St 1 , St 2 ) = 1 overestimates E * r3np (ξ ) 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 * r3np (ξ ) for large wave numbers, and the assumption of coh(ξ |St 1 , St 2 ) = 1 can be applied only for predicting E * r3np (ξ ) for small wave numbers (ξ < O(10 −1 )).On the other hand, E * r3np (ξ ) values predicted by the parameterization using Eq. ( 31) show good agreement with those obtained by the DNS data for overall wave numbers.The error level of the parameterization using Eq. ( 31) is evaluated by the rms error e rms in units of decibels.e rms is defined as where ξ is defined as ξ = ln ξ 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., 0.05 ≤ ξ ≤ 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 0.05 ≤ ξ ≤ 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.

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. (2014Matsuda et al. ( , 2017)).The large gravitational settling can modulate E * np (ξ |St), and that can be a cause of the significant difference in the radar reflectivity factor.However, the influence on E * np (ξ |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 * np (ξ |St) in Eq. ( 33).Ayala et al. (2008a) and Lu et al. (2010) reported that gravitational settling enlarges the crossover length of the bidisperse RDF.Lu et al. (2010) extended the perturbation expansion analysis of Chun et al. (2005) and presented the formulation for the crossover length in the presence of gravity, which is where C Chun is the coefficient derived by Chun et al. (2005), a 0 is the ratio of the acceleration variance a 2 to the square of the Kolmogorov acceleration a 2 η , τ a is the acceleration correlation timescale, and τ g is the correlation timescale for gravitational settling particles.Chun et al. (2005) obtained the values of C Chun ≈ 5.0 and a 0 ≈ 1.545 based on their DNS results.Lu et al. (2010) further assumed S v 1 so that τ g is simply given by τ g = τ a = 1.5τ η .The crossover length l c of Eq. ( 38) becomes equivalent to that of Chun et al. (2005) when gravitational settling is negligibly small, whereas the gravity effect is dominant for the cases of F r < 0.47, which often appears in cloud turbulence.Thus, in this study, we modify Eq. ( 32) to include the settling influence on the coherence coh(ξ |St 1 , St 2 ).Since ξ c is inversely proportional to l c /l η , we propose the following correction based on Eq. ( 38): Note that this study also adopts the value of a 0 obtained by Chun et al. (2005).
The reliability of the modified parameterization for the case with gravitational settling has been evaluated in the same way as the previous subsection.Figure 7b shows E * r3np (ξ ) for the case with gravitational settling.E * r3np (ξ ) values at large wave numbers are smaller than those for the case without gravitational settling, and the difference from the parameterization with coh(ξ |St 1 , St 2 ) = 1 is larger, indicating that the coherence model is more important than the case without gravitational settling.It is also confirmed that E * r3np (ξ ) 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 0.05 ≤ ξ ≤ 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 0.05 ≤ ξ ≤ 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

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 atmosphereocean 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 G(y)dy = n p m r p q r r p dr p , where y = ln r p , and m(r p ) is the mass of droplets with a radius of r p .The mass coordinate m and logarithmic coordinate y are discretized as where dy = ln 2/(3s), 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 G(y)dy; 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 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 x = y = 25 m and z = 20 m.This study used the temporal slice of cloud simulation data at a higher resolution.

Computational method for radar reflectivity factor
The radar reflectivity factor Z, including the influence of particulate and clear-air Bragg scatterings, is given by where Z PB is the particulate Bragg-scattering part of Z in Eq. ( 9), and Z CB is an additional term reflecting clear-air Bragg scattering.In real clouds, Z PB is caused by droplet number density fluctuations due to turbulent droplet clustering and turbulent entrainment of environmental clear air.E r3np (k) for these factors is given by (Matsuda et al., 2014), where E clust r3np (k) and E entr r3np (k) 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 PB = Z PBc + Z PBe , where Z PBc = 2 7 π 2 κ −2 E clust r3np (κ) and Z PBe = 2 7 π 2 κ −2 E entr r3np (κ).The spectrum E clust r3np (k) was calculated using Eq. ( 10) and the parameterization proposed in the previous section.To determine the Stokes number for each droplet size, the energy dissipation rate of the cloud simulation data was calculated based on the Smagorinsky model; i.e., where C s is the Smagorinsky constant (C s = 0.173 in this study), s the representative grid spacing given by s = ( x y z ) 1/3 , and s ij the strain rate tensor, which is given by , where u i is the air velocity in the resolved scale.
The spectrum E entr r3np (k) was calculated using the wellknown scalar concentration spectrum.Erkelens et al. (2001) estimated this contribution based on the −5/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 entr r3np (k) is approximately given by (Hill, 1978) where χ r3np is the scalar dissipation rate for r 3 p n 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 (γ = 1.4).
The contribution of clear-air Bragg scattering was calculated by where E n (k) is the power spectrum of refractive index fluctuations.The refractive index n ref is given approximately by (Balsley and Gage, 1980) where p and e are the atmospheric pressure and partial pressure of water vapor in hPa, and T is the absolute temperature.The contribution of free electrons was omitted because it is negligibly small in the troposphere.The power spectrum E n (k) is given by the scalar concentration spectrum for P r < 1 (Pao, 1964), where P r ≡ ν/D (D is the scalar diffusive coefficient) is the Prandtl number: where χ n is the scalar dissipation rate for n ref .
In this study, the Prandtl number of refractive index was set to 0.7.The scalar dissipation rates for r 3 p n p and n ref were calculated in the same way.That is, the dissipation rate of an arbitrary scalar θ is given by www.atmos-chem-phys.net/19/1785/2019/where ν t is the eddy kinematic viscosity, Sc t is the turbulent Schmidt number, and θ is the scalar value in the resolved scale.ν t was calculated by using the Smagorinsky model; i.e., ν t = (C s s ) 2 2 s ij s ij 1/2 , and Sc t was set to 0.4 (Moin et al., 1991).

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 τ = Q ext π r 2 p n p z , where Q ext is the extinction efficiency for Mie scattering (Q ext = 2.0 in this study), and r 2 p is given by r 2 p = ∞ 0 r 2 p q r (r p )dr p .Note that the optical transmittance of cloud volume is approximately equal to 1 − τ when τ is sufficiently smaller than unity.Figure 8b shows the isosurfaces of the energy dissipation rate and the vertical velocity u 3 .The locations of upward flows correspond to the locations of clouds and large regions are observed around the upward flows.This indicates that strong turbulence is generated by entrainment motions due to updrafts.This study focused on a vertical cross section that slices the cumulus cloud with the largest upward velocity.Figure 9 shows the liquid water content (LWC) in a logarithmic scale, the energy dissipation rate , the radar reflectivity factor Z dB , the increase of Z dB due to particulate Bragg scattering, and particulate Bragg scattering due to turbulent clustering Z dB PBc in the cross section.The microwave frequency was set to f m = 2.8 GHz, which is the representative frequency of S-band radars.The radar reflectivity factor is shown in units of decibels, which is defined as Z dB (dBZ) = 10log 10 Z (mm 6 m −3 ).Large values of Z dB are observed inside and below the clouds.The strong echo below the clouds reflects drizzling regions, where the LWC is smaller than that inside the clouds but the strong echo returns from the drizzling droplets because Z incoh is proportional to r 6 p .The radar echo on the outside of the isoline of Z dB incoh = −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 largescale 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 dB − (Z incoh + Z CB ) 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 dB PBc is plotted in Fig. 9e.Z dB PBc is larger than −10 dBZ inside the turbulent cloud region, where the LWC is larger than 0.1 g m −3 , and the energy dissipation rate is intermittently larger than 100 cm 2 s −3 .Large values of Z dB PBc 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 num-ber (up to 0.05) distributed intermittently corresponding to the distribution of .The main factor of the height dependence of Z dB PBc is the LWC, which is larger than 1 g m −3 at the near-top of the clouds.Note that Z PBc is proportional to square of the LWC as Eqs.( 9) and (36) implies Z PBc = 2 3 3 2 ρ −2 p κ −2 (LWC) 2 l η E * r3np (κl η ).Thus, the significant influence of turbulent clustering is caused by sufficiently large values of the energy dissipation rate and the LWC.

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 1 and St 2 , has values between the power spectra for St 1 and St 2 at small wave numbers, whereas the spectrum decreases more rapidly than the power spectra as the wave number increases.This decreasing trend is related to the scale dependence of the spatial correlation of cluster locations between two different Stokes numbers.The coherence of the cross spectrum is close to unity for small wave numbers and decreases almost exponentially with increasing wave number.This is qualitatively consistent with the visualization results, in which the clustering locations for different Stokes numbers are almost the same on large scales, whereas a discrepancy in clustering locations is observed at small scales.It is also confirmed that the decreasing trend of the coherence is strongly dependent on the combination of Stokes numbers.
In order to develop a cross spectrum model for estimating the clustering influence on the radar reflectivity factor, we have proposed an exponential model for the wave number dependence of the coherence and introduced the critical wave number (i.e., the decay constant for the model) to consider the dependence of the coherence on the Stokes number combination.The coherence data for all combinations of six Stokes numbers ranging from 0.05 to 2.0 reveal that the critical wave number is inversely proportional to the Stokes number difference, |St 1 − St 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 3 p -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 highresolution 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.
Data availability.Access to the simulation and analysis results can be granted upon request under a collaborative framework with JAM-STEC.

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

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

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

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

Figure 7 .
Figure 7. Comparisons of r 3 p -weighted power spectrum E * r3np (ξ ) 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 .

Figure 8 .
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 2 s −3 and (yellow) the vertical velocity u 3 = 3 m s −1 .
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 ≡ τ p /τ η (τ 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.
)/q r (r p ) so that n(x, r p ) is given by n(x, r p ) = n p (x|r p )q r (r p ).The number density distribution function n p (x|r p ) satisfies x∈V n p (x|r p )dx = N p for arbitrary r p .Note that the spatial correlation function δn(x, r p )δn(x + r, r p ) r p dr p dr p ,(10)where C np (k|r p , r p ) is the cross spectrum of number density fluctuations for n p (x|r p ) and n p (x|r p ): the cross spectrum C np (k|r p , r p ) is defined as C np (k|r p , r p ) ≡ |k|=k (k|r p , r p )dσ k ; i.e., the integration of (k|r p , r p ) over the spherical shell, σ k , at |k| = k, in which (k|r p , r p ) is the cross spectral density function, defined as

Table 1 .
Computational settings of the DNS.
shows the critical wave numbers ξ c for all combinations of St 1 and St 2 , where St 2 > St 1 .ξ c for the same St 1 increases as St 2 decreases, whereas ξ c for the same St 2 increases as St 1 increases.This indicates that ξ c increases as St 1 and St 2 becomes closer to each other.It should be noted that several studies discuss the spatial correlation of bidisperse clustering particles.
cross spectrum of number density fluctuations is a Fourier transform of the bidisperse RDF, the Stokes number ratio, St 2 /St 1 , and the Stokes number difference, St 2 -St 1 , are candidates for the dominant parameter for ξ c . Figure 5 shows ξ c plots against St 2 /St 1 and St 2 -St 1 .These figures clearly inwww.atmos-chem-phys.net/19/1785/2019/ ://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 12.8×12.8×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., x = y = 100 m and z = 40 m.