Effects of turbulence modulation and gravity on particle collision statistics

Dynamics of inertial particles in homogeneous isotropic turbulence is investigated by means of numerical simulations that incorporate the effect of two-way interphase momentum transfer. The continuous phase is solved in the Eulerian approach employing Direct Numerical Simulations (DNS). The dispersed phase is treated using the Lagrangian approach along with the point-particle assumption. The main focus is on computing collision statistics of inertial particles relevant to cloud droplets in typical atmospheric conditions. The vast majority of previous DNS were performed assuming one-way momentum coupling between continuous and dispersed phases. Such simplified approach is adequate only for dilute systems with relatively low mass loading. In this study we investigate the effect of two-way momentum coupling on the kinematic and dynamic collision statistics of the dispersed phase. A number of simulations have been performed at different droplet radii (inertia), mass loading, viscosity and energy dissipation rate. To assess the accuracy of numerical approach the coupling force (exerted by particles on the fluid) was computed using two different techniques, namely particle in cell and projection onto neighboring node. To address the effect of gravity, the simulations have been carried out simultaneously both with and without gravitational acceleration. It has been found that the effect of two-way coupling is significant both for droplet clustering and the radial relative velocity. It turns out that the collision kernel is more sensitive to the particle mass loading when the gravitational acceleration is considered. The collision kernel of settling droplets increases as the droplet mass loading increases. This is direct consequence of larger radial relative velocity. For non-settling droplets the effect of mass loading is opposite, namely, we observe a minor reduction of the collision kernel as the number of droplets increases. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


a b s t r a c t
Dynamics of inertial particles in homogeneous isotropic turbulence is investigated by means of numerical simulations that incorporate the effect of two-way interphase momentum transfer. The continuous phase is solved in the Eulerian approach employing Direct Numerical Simulations (DNS). The dispersed phase is treated using the Lagrangian approach along with the point-particle assumption. The main focus is on computing collision statistics of inertial particles relevant to cloud droplets in typical atmospheric conditions. The vast majority of previous DNS were performed assuming one-way momentum coupling between continuous and dispersed phases. Such simplified approach is adequate only for dilute systems with relatively low mass loading. In this study we investigate the effect of two-way momentum coupling on the kinematic and dynamic collision statistics of the dispersed phase. A number of simulations have been performed at different droplet radii (inertia), mass loading, viscosity and energy dissipation rate. To assess the accuracy of numerical approach the coupling force (exerted by particles on the fluid) was computed using two different techniques, namely particle in cell and projection onto neighboring node. To address the effect of gravity, the simulations have been carried out simultaneously both with and without gravitational acceleration. It has been found that the effect of two-way coupling is significant both for droplet clustering and the radial relative velocity. It turns out that the collision kernel is more sensitive to the particle mass loading when the gravitational acceleration is considered. The collision kernel of settling droplets increases as the droplet mass loading increases. This is direct consequence of larger radial relative velocity. For non-settling droplets the effect of mass loading is opposite, namely, we observe a minor reduction of the collision kernel as the number of droplets increases.

Introduction
Interaction of inertial particles with turbulent flows is a key factor of many environmental processes including the microphysics of clouds ( Devenish et al., 2012;Grabowski and Wang, 2013 ), the dispersion of pollutants, or motion of dispersed particles in water and in the atmosphere ( Eames, 2008;Ruiz et al., 2004;Rosa et al., 2013 ). Turbulent transport is also a driving mechanism in many technological processes in industry ( Brennen, 2009;Crowe et al., 2015 ), including the fields of power generation or process engineering. Examples of these are: efficient combustion of pulverized coal in boilers, removal of particulate matter from flue gases, pneumatic transport in pipelines, combustion of liquid fuel in engines or large scale boilers, spraying of fertilizers and plant protection agents.
It is worth noting that detailed knowledge of the mechanism of turbulence-particle interactions has numerous applications in meteorology. These include: modeling of sandstorms, volcanic ash transport and precipitation formation from cloud droplets and ice crystals. Investigation of the size-broadening and growth of cloud droplets in a turbulent environment is crucial for quantifying rain initiation and rain amount in warm clouds ( Lau and Wu, 2003 ).
In this study the focus is on modeling systems akin to clouds with similar or larger liquid water content. Accurate description of the droplet collision-coalescence and precipitation formation is essential for reliable weather and climate predictions on Earth. These processes cannot be fully resolved in Numerical Weather Prediction  (NWP) simulations because their characteristic length scales are significantly smaller than those defining large-scale atmospheric flows. The detailed examination of the cloud microphysical processes is a necessary step in the development of more accurate parameterizations for NWP models. Due to their practical importance, the cloud microphysical processes have been extensively investigated. Most previous numerical studies were based on point-particle approximation and oneway momentum coupling, for example ( Rosa et al., 2016;Parishani et al., 2015 ). This simplified approach is adequate only for dilute systems with low mass loading, i.e. when the concentration of the disperse phase in a flow is small, with the volume fraction below 10 −6 , see ( Elghobashi, 1994 ). For such suspensions, the particles do not affect the motion of the continuous phase. For larger concentrations (mass loading), the effect of the turbulence modulation by particles becomes important.
The goal of the present study is to investigate the impact of droplet mass loading on turbulent flows and collision statistics under conditions similar to the atmospheric clouds. To model the cloud processes, we use the standard Eulerian-Lagrangian approach along with the point-particle approximation. The method combines fully resolved DNS of turbulent flow with the Lagrangian tracking of individual droplets. The key element of novelty of the present study is the inclusion of two-way momentum coupling effects (mutual interaction between the flow and particles) and gravity. We point out that the radii of considered droplets are at least one order of magnitude smaller than the grid spacing of computational meshes used for modeling ambient turbulent flow (see Table 2 ). This implies that in one computational cell there may be a large number of droplets. Therefore, the numerical method is suitable to represent only the bulk effect of the droplets on the background flow. It is important to note that in simulations under two-way coupling the interaction force between widely separated particles is partially represented. In the case of slowly approaching particles, the interaction effect results from modulation of the background turbulent flow. The fine structures of the flow between nearly touching droplets are not captured because in the point-particle approximation the flow around the individual droplets is not resolved. The size of the droplets is considered implicitly when computing the coupling force between droplets and turbulent air ( Akiki et al., 2017 ). On the other hand, the size of the particles is treated explicitly for detecting collisions between droplets. In this study, we are only interested in geometric collisions without considering local droplet-droplet aerodynamic interactions and surface-surface contact forces. Namely, droplets are allowed to overlap at the beginning of a time step and are not removed from the system after collision. This is the first step and other quantities can be modeled separately as the collision effi- Table 2 Basic characteristics of particles of radii a with respect to turbulent scales. The estimates are based on two different turbulent flows A and C (see Table 1 ) at zero mass loading. τ e is the ratio of the longitudinal integral length scale to the rms fluctuating velocity. The Froude number F r = St * S 2 V is defined as the ratio of particle response time to the residence time of the particles in a Kolmogorov eddy. In all cases the energy dissipation rate was set ε= 400 cm 2 /s 3 . ciency and coalescence efficiency ( Wang et al., 2005 ). Unlike solid particles where the effect of material elasticity has to be considered ( Goswami and Kumaran, 2010a, 2010b, 2011, the coalescence efficiency of small liquid droplets is relatively easy to handle due to the effect of the strong van der Waals force. Most of the simulations in the present study have been performed assuming the ambient energy dissipation rate = 400 cm 2 /s 3 . This value is representative for moderate to strong convection in clouds. In Section 7 we also analyze several simulations performed at lower values of . Our simulations are limited to monodisperse systems, i.e. all particles in domain have the same size. We demonstrate that the turbulence modulation by settling droplets significantly affects their radial relative velocity and consequently alters the collision rate. This important effect should be parameterized in future numerical weather prediction models.

Related studies
There are several studies aimed at quantifying the effects of two-way momentum coupling on turbulence modulation and dynamics of inertial particles. The first rigorous numerical simulations were performed by Squires and Eaton (1990) . Due to computational constraints, their DNS of stationary isotropic turbulence were limited to relatively coarse grid resolutions (meshes with 32 3 and 64 3 nodes) and, consequently, low Taylor microscale Reynolds numbers ( R λ < 38). Considered systems were dilute (low volume fraction) but the particle mass loading was substantial (varied from 0.1 to 1) to have a significant impact on the turbulent flow. The particle inertia, in terms of τ p / τ e , was relatively large and varied from 0.075 to 1.5 (see Table 2 for comparison with the present study). Here, τ p is the Stokes inertial response time of the particle, while τ e is the ratio of the longitudinal integral length scale L s to the rms fluctuating turbulence velocity u . In that pioneering work, the gravitational settling of the particles was not considered.
Several important conclusions result from that study. First, attenuation of the turbulent kinetic energy depends on the mass loading but is weakly sensitive to the particle response time. Second, the energy and dissipation spectra in two-way coupled flows are larger at higher wave numbers and increase with the particle mass loading. It is worth noting that this effect was observed only in simulations performed on the mesh with 64 3 grid nodes. Third, the turbulent flow is modulated more homogeneously by particles of larger inertia. This is a natural consequence of more uniform spatial distribution of large-inertia particles.
The results of Squires and Eaton (1990) were confirmed and extended in a follow-up study by Elghobashi and Truesdell (1993) . The succeeding DNS were performed at larger mesh size 96 3 with particles of inertia St = O (1) . Here, the Stokes number is defined as St = τ p /τ K , where τ K is the Kolmogorov time scale. In contrast to the previous study, the new simulations were limited to decaying turbulence, but the gravitational acceleration was included in the particle equation of motion. Elghobashi and Truesdell (1993) claimed that the increase of energy (at high wave numbers) is the effect of the larger viscous dissipation rate. The increase in dissipation rate results in an increase in the rate of energy transfer from the large-scale motion. They also observed that gravity changes the mechanism of the momentum transfer. The settling particles transfer their momentum to the small-scale motion in a non-isotropic manner. The effect of anisotropy is mitigated by the redistribution of energy (at the same wavenumber) by pressure-strain correlation. Boivin et al., 1998 extended the conclusions formulated by Elghobashi and Truesdell (1993) and Squires and Eaton (1990) . Their study focused on quantifying the effects of Stokes number ( St = 1.26, 4.49 and 11.38) on turbulence modulation. Turbulence was simulated using DNS on mesh with 96 3 nodes ( R max λ = 62 ). The particle equation of motion did not include the gravitational term. Boivin et al., 1998 showed that viscous dissipation in the fluid decreases with increasing mass fraction and is larger for particles with smaller inertia.
In a subsequent study, Sundaram and Collins (1999) analyzed different mechanisms responsible for energy exchange between continuous and dispersed phases. Their DNS were performed for decaying turbulence on a uniform mesh with 128 3 nodes. The number of tracked particles reached 64 3 . Particle gravitational settling was neglected. Interestingly, the numerical model included the effects of particle collision. The collisions were treated as elastic or momentum and energy conserving. They claimed that in the absence of gravity (and other source terms), particles reduce the turbulent kinetic energy of the fluid by increasing the viscous and drag dissipation. Moreover, in simulations with two-way coupling the particle velocities remained correlated at longer distances. Vermorel et al. (2003) used DNS to investigate the turbulence modulation in a particle laden slab flow. The inertial particles were injected at high velocity into a cube with freely decaying isotropic turbulence. Transfer of momentum from the particles resulted in strong fluid acceleration in the slab region. At the borders of the slab turbulence was enhanced due to the production by mean fluid velocity gradients. The opposite effect was observed in the slab core, namely the turbulence was suppressed by particles. Their study was limited to relatively low initial Reynolds number, i.e. R max λ = 35 (mesh size 128 3 ). Gravity forces were neglected due to large relative velocity of the particles and short duration of the simulations. Although valuable per se as a numerical experiment on droplets in turbulence, the geometrical configuration of the problem considered there is not closely related to the problem of cloud microphysics. Bosse et al. (2006) examined the effect of the two-way momentum coupling on the settling rate of the disperse phase in stationary homogeneous isotropic turbulence. They found that particles exert a collective force on the carrier fluid, especially in regions of high concentration, causing local fluid acceleration in the direction aligned with gravity. The enhanced downward fluid motion leads to larger settling velocity in these regions, thus increasing the overall mean particle settling velocity. This study was limited to low Taylor microscale Reynolds number ( ~40).
In a recent study, Monchaux and Dejoan (2017) performed a similar two-way coupled DNS at R λ = 40 on 64 3 grid. Their numerical results were in good agreement with the experiments of Aliseda et al. (2002) . They confirmed earlier observations that the settling velocity increases with increasing volume fraction and local concentration. They also considered the impact of two-way coupling on the particle preferential concentration. It has been shown that clustering of small inertia particles is weaker for larger volume fraction and larger gravity. A reverse tendency is observed for large inertia particles. This behavior is related to an attenuation of the centrifuge effects and to an increase of particle accumulation along the gravity direction.
The simulations described above have been performed under the point-particle assumption. As for the fully resolved simulations of turbulence with finite-size particles, the two-way momentum coupling is automatically accounted for there. Yet, due to high numerical complexity, such studies have become possible only in recent years ( Garcia-Villalba et al., 2012;Hui et al., 2013;Maxey, 2017;Peng et al., 2019 ). With nowadays computational resources it is feasible to simulate systems with up to O (10 5 ) particles.
Most previous studies were focused mainly on investigating the effect of two-way momentum coupling on the continuous phase. In the present work we address the collision statistics of the dispersed phase. In Section 3 we describe the essentials of the numerical method. Detailed characteristics of the turbulent flows are presented in Section 4 . A thorough analysis of the impact of particles on turbulence modulation is provided in Section 5 . The twopoint statistics from the DNS are discussed in Section 6 . The results include sensitivity of the radial distribution function, radial relative velocity and collision kernel to the droplet mass loading. In Section 7 we briefly address the effect of energy dissipation rate on kinematic collision statistics. Section 8 contains a summary and main conclusions.

The numerical method
To simulate homogeneous isotropic turbulence a standard pseudo-spectral method ( Orszag and Patterson, 1972 ) was used. In the method, the Navier-Stokes (N-S) equations are solved on a 3D uniform mesh with N equally spaced grid points in each spatial direction. The flow domain is a cube with size 2 π . Periodic boundary conditions are naturally imposed which is consistent with the 3D discrete Fourier transform applied to the fluid velocity field. The fluid velocity U is found from numerical integration of the N-S equation in rotational form for an incompressible fluid satisfying the continuity equation: Here ω ≡ ∇ × U is the vorticity, P is the pressure, ρ is fluid density and ν is fluid kinematic viscosity. To obtain statistically stationary turbulence, we used a spectral forcing scheme similar to that of Sullivan et al. (1994) . The forcing term f ( x , t ) is nonzero only for a few low wavenumber modes (| k | < 2.5) in the Fourier space. The energy of the first two wavenumber shells (0.5 < | k | < 1.5 and 1.5 < | k | < 2.5) is specified to be a constant. The values of energy in these shells are preset so that their ratio satisfies k −5 / 3 energy spectrum, accounting for the number of modes forced in each shell volume. A total number of 80 modes are forced in this scheme.
The last term f ( p ) in Eq. 1 refers to the cumulative force per unit mass exerted by particles (see Eq. (5) below) on the fluid, as detailed in ( Ireland and Desjardins, 2017 ;Kasbaoui et al., 2019 ) (3) Fig. 1. Schemes of two different approaches for computing mean interphase momentum transfer (as in Garg et al., 2007 ). (a) Particles in cell (PIC) -momentum at grid node 1 is modulated by the motion of all particles located in the square area (or cube in 3D) bounded by a dashed line. Typically, size of the square/cube corresponds to size of grid box. (b) Projection onto neighboring nodes (PNN) -contribution of the particle momentum to the fluid momentum at a grid node depends on the separation distance (particle -node), as for example in Elghobashi and Truesdell (1993) . Alternatively, the particle force is projected onto neighboring nodes using weights which are proportional to cell volumes ( Squires and Eaton, 1990 ). For example, the fluid flow at grid node 1 is affected by a fraction of the particle Stokes drag proportional to the area with marker 1.
where δ is the Dirac delta (suitably mollified in computations), φ m is the local mass loading, V E (x , t ) is the Eulerian particle velocity at the location x and in the absence of gravity. Consistently, V ET (x ) is the particle terminal velocity defined at nodes of the regular mesh. In numerical simulation V E (x , t ) can be computed directly from the Lagrangian particle velocities where σ is an interpolation kernel ( Kasbaoui et al., 2019 ) whose width is equal to the grid spacing. The same method can be used to determine V i ET (x ) . Parameter M is a weighting factor ( Elghobashi, 1994 ) defined as the ratio of N r / N c . Here, N r and N c represent the number of real and computational particles, respectively. Simulations with M > 1 are based on an assumption in which one super-particle (computational particle or parcel) represents a distribution of several smaller (real) particles. This approach allows to reduce the computational cost of simulations with a prohibitively large number of particles, in particular at large mass loadings. The force term f ( p ) in the simplified approach is evaluated for all computational particles and multiplied by the parameter M . As noted in ( Garg et al., 2009 ), the average momentum transfer from particles to the fluid, in traditional Eulerian/Lagrangian simulations with M = 1, strictly depends on the grid resolution. The numerical error associated with the mean interphase momentum transfer increases with the grid refinement. Moreover, the error is non-uniform in space and depends indirectly on the particle inertia. The remedy proposed by ( Garg et al., 2009 ) consists in updating (in time) the statistical weight in a way that the number density of computational particles remains nearly uniform. Although the method yields a numerically convergent solution in terms of mean momentum transfer, the consequences of particles annihilation and cloning on the collision statistics have not been explored in detail. Therefore, most simulations presented in this study were performed for M = 1 . Simulations with M > 1 were necessary to address the problem of turbulence modulation at large mass loadings, but in these simulations, we tried to keep the parameter M possibly close to unity. The Navier-Stokes equations have to be solved together with the particle equation of motion ( Maxey and Riley, 1983 ) where i is the particle number, τ p is the Stokes inertial response time, V i ( t ) is actual particle velocity, U ( Y i ( t ), t ) denotes the fluid velocity at the particle location Y i ( t ) and g is the gravitational acceleration. Re p = 2 aV rel /ν is the particle Reynolds number, where V rel is the particle-fluid relative velocity; f is the drag correction factor (beyond the Stokes regime) which in this study was set to f = 1 .
To incorporate the effect of the particle motion on the background turbulent flow, the coupling force f ( p ) needs to be computed at every time instant. It should be noted, however, that adding this coupling term to the momentum equation ( Eq. 1 ) may have an indirect effect on the continuity equation. As a result, the fluid velocity field will not be divergent free. However, according to Elghobashi and Truesdell (1993) , this undesirable effect can be safely neglected if the volume fraction of the particles is relatively low. Let us notice that Pakseresht and Apte (2019) proposed an update to the map of particle-turbulence interaction regimes ( Elghobashi, 1994 ), extending it with the volumetric two-way coupling. They found that the so-called volumetric displacement effects in the governing equations of the carrier phase become significant for the volume fractions above 5%. In the vast majority of simulations performed in the present study the volume fraction did not exceed 10 −3 . Hence, it can be safely taken that the continuity equation is practically unchanged by the presence of droplets. In the pseudo-spectral code, the divergence of fluid velocity is maintained to be zero by projecting the velocity vector on the plane normal to k , at every time step. The second problem concerns the non-zero mean fluid velocity in simulations with settling particles. This undesirable effect may lead to instability of computations. Therefore, we apply a procedure to fix the mean flow velocity. In the Fourier space the mean flow is represented by the mode at | k | = 0 . We set the amplitude for this mean flow mode to zero. Essentially, this amounts to applying a vertical pressure gradient to counter-balance the weight of the solid particles. In case of simulations with a large number of settling particles, especially those of high inertia, an important question arises about the loss of energy due to this filtering (high-pass). We evaluated the energy losses based on a simulation with 8 million settling droplets of radii 50 μm. The effect of filtering is very small, and the energy loss, at every time step, is of the order of O(10 −7 ) of the total kinetic energy. In reality, the potential energy loss is compensated by the work associated with the applied vertical pressure gradient.
To integrate in time the N-S equation, the coupling force f ( p ) must be evaluated at all grid nodes of the regular mesh. All droplets that are in the vicinity of a given grid node contribute to f ( p ) . The carrier fluid around each droplet is not resolved because the size of the considered droplets is much smaller than the grid spacing. This implies that the coupling force evaluated at the grid node is a function of the local volume-averaged perturbation velocity generated by droplets. According to 3rd Newton's law, the force exerted by particles on the fluid is opposite to the drag exerted on the particles by the fluid. Here comes a conceptual problem, as the forces acting on the particles are known at the particles' locations only. Under the point-particle assumption, the force is zero everywhere except for a delta function at the location of each particle. Thus, for the purpose of practical computations, it is necessary to introduce a mollified / regularized delta function instead, see Eq. 4 . In other words, the corresponding components of f ( p ) need to be projected/interpolated from particle locations to the nodes of the regular mesh. Typically, the contribution to f ( p ) is restricted to the neighboring particles. In applications where the size of the particles is much smaller than both the grid spacing and the Kolmogorov length such simplification is justified. There are several approaches in the literature ( Garg et al., 2007 ) to compute the source term f ( p ) . The most common are: (1) the particle-in-cell, (2) the projection onto neighboring nodes, and (3) the projection onto identical stencil. In this study, we tested two of the methods, namely, particle-in-cell (PIC) and projection onto neighboring node (PNN), see Fig 1 . In the PIC approach, the coupling force is computed as the summation of forces exerted on the fluid by each particle in the control volume surrounding a grid node. In other words, the interpolation kernel σ has a uniform top-hat shape. The PNN method takes into account separation distance between the grid node and the nearby particle. Thus, σ is a standard bi-or tri-linear function. In the present study the weights were computed based on the cell volume partition, see Fig. 1 To integrate the equation of motion, Eq. 5 , the Stokes drag force needs to be evaluated for every particle and at every time step. In the numerical (Eulerian-Lagrangian) approach the drag force is proportional to the difference between the actual particle velocity and the fluid velocity at the particle location. Since the fluid velocity is solved on a regular grid and has discrete representation, another interpolation method is needed to evaluate its value exactly at the particle location. In all simulations performed in the present study the standard 6-point Lagrangian interpolation scheme in each direction was employed ( Ayala et al., 2014 ). This method has been extensively and successfully used in many earlier studies concerning modeling of two-phase flows under oneway momentum coupling. Recently, several alternative methods for computing particle drag have been developed. Ireland and Desjardins (2017) proposed an improved formulations of the drag that provide accurate and grid-independent predictions of particle settling in two-way coupled flows at low particle Reynolds numbers. In turn, Akiki et al. (2017) extended the point-particle model in a way that the drag force includes also effects of hydrodynamic interactions between neighboring particles. In another study, Horwitz and Mani (2016) showed that in two-way coupled point-particle simulations the Stokes drag acting on the particle may be underestimated if evaluated based on the disturbed fluid velocity (disturbed by the particle itself and all neighboring particles). To predict the Stokes drag more accurately a new method was proposed that allows to estimate the undisturbed fluid velocity from the neighboring disturbed fluid velocity information. This improved method was tested and analyzed in the follow up study by Horwitz and Mani (2018) . The important conclusion resulting from these analyses is that the correction to the drag force is required if 2 a/η ≥ 10 −1 (see the regime diagram, Fig. 10 therein). In all simulations performed in the present study the droplet radii were much smaller than the Kolmogorov length Table 3 Correspondence between the number of particles/droplets and the mass fraction in turbulent flows of different statistical characteristics. The estimates are based on data listed in Table 1 . Five different droplet sizes were considered, i.e. radii from 20 to 60 μm. The energy dissipation rate was assumed to be ε = 400 cm 2 /s 3 .

Particle mass loading # particles
Estimates for Flow A ( scale, so according to Horwitz and Mani (2018) the correction is not required. Nevertheless, this aspect is worth checking in future studies.
It is important to note that the PIC/PNN methods and the 6point Lagrangian method differ in the size of the stencil. That one for PIC / PNN methods is 6 times smaller in each direction. Therefore, the question arises about accuracy of these two methods and consequences of spatial discretization. This problem was thoroughly examined by Sundaram and Collins (1996) . They derived the analytical formula for the maximum error associated with inconsistent interpolations. Based on the formula and numerical tests they recommended to use the same high order method for forward and backward interpolations. However, in realistic simulations the numerical error depends mainly on the velocity fluctuations at fine scales. Therefore, it can be assumed that if the velocity gradients of the flow are small (flow cases with higher viscosities, see Table 1 ) the order of interpolation is of secondary importance. In Appendix we show that the statistics computed using the PIC and PNN methods are in good quantitative agreement.
To perform the two-way coupled simulations a massively parallel application was employed. The MPI code was designed to run on supercomputers with distributed memory. A complete description of the code, along with results of former numerical experiments and scalability analysis can be found in ( Ayala et al., 2014;Parishani et al., 2015;Rosa et al., 2015 ). The two-point particle collision statistics such as the radial distribution function and radial relative velocity are handled using a specially designed parallel algorithm with optimized data communication between processes. The algorithm employs the cell-index method and the linked lists concept ( Allen and Tildesley, 1987 ) for efficient detection of closely spaced particles.
In simulations with M = 1 the particle collisions are computed in a deterministic way, that means all collisions between individual particles are treated explicitly. This method cannot be applied when M 1 because collisions of two parcels of particles may bring M or even more (self-collisions of droplets within a parcel) realistic collisions. The problem of the collision detection has been analyzed in several previous studies. In terms of numerical performance, the most promising are methods based on the stochastic approach e.g. ( Sommerfeld, 2001; O'Rourke et al., 2009 ). However, Fig. 2. Energy spectra of turbulent flows simulated using DNS at different parameters of viscosity (marked with different colors of lines). Solid lines -simulations under one-way momentum coupling. Dashed lines -simulations with particles (40 μm droplets) and two-way momentum coupling. (a) Simulations with non-settling particles, (b) effect of gravitational settling included. In each simulation under two-way coupling the particle mass loading was fixed and equal to m = 0.24. the number of assumptions and parametrizations of different physical processes in these models question their accuracy. Recently, ( Johnson, 2019 ) proposed a novel deterministic method to address particle-particle collisions when M = 1. The method involves artificially enhancing the collisional radius depending on the inertia of the tracked particles. For two extreme cases, i.e. low inertia particles ( τ p → 0) and high inertia particles ( τ p → ∞ ) the collisional radius should be increased to 3 √ M (2 a ) and √ M (2 a ) correspondingly. This approach has potential to be used in simulations in large computational domains (grids). It should be added, however, that the method has been developed neglecting two-way momentum coupling effects. Since the currently available methods do not guarantee the required accuracy, the collision statistics of the systems at larger mass loading have been evaluated only coarsely. To reduce the statistical error, we carried out simulations with possibly large number of computational droplets. In Section 6 we show that the statistics computed using the approximate method, at low mass loading, are in quantitative agreement with the statistics evaluated using the exact method. In simulations with larger mass loading we observe some discrepancies, but their absolute value is relatively small. This positive effect may be due to the weaker clustering (more uniform distribution) of particles in two-way coupled systems.

Results on flow statistics
In turbulent clouds the transfer of momentum from microdroplets to the air occurs mainly at fine turbulent scales and ε is the average energy dissipation rate in the whole computational domain, computed in DNS without droplets.
is largely limited to the dissipation range of the energy spectra. This can be explained by the fact that the characteristic length scale of droplet (e.g., the radius) is much smaller than the Kolmogorov length scale, namely a η. In DNS of turbulent flows η is typically smaller or equal to the grid spacing x (that depends on viscosity, see Table 2 ). Based on this dimensional analysis, we confirm that the numerical approach in which the coupling force is projected on the nearest grid nodes only, is accurate enough for modeling the cloud processes. Although the characteristic length scale of the momentum transfer is relatively small, the effect of two-way coupling may be important also on dynamics of the systems at larger scales, corresponding for example to the integral length scale. This stems from the fact that the disperse phase alters the energy cascade. The transfer of kinetic energy among different scales depends on non-trivial triadic interaction of wave numbers. According to Ferrante and Elghobashi (2003) the spectral nonlinear energy-transfer rate to wave number k is given where P ij is the projection tensor and stands for the imaginary part. The superscripts "ˆ "and " * " denote, correspondingly, the Fourier transform and the complex conjugate.
In that regard, it is important to investigate how the mechanism of momentum transfer between droplets and small-scale vortical structures affects the dynamics of the system in the entire range of energy spectra. A typical measure of flow resolution (in terms of small-scale structures) in pseudo-spectral DNS is the parameter k max η. The parameter must be greater than unity for fine scales to be resolved. Here k max is the maximum wave number of computations. In all DNS analyzed in the present study k max = int (N/ 2 − 1 . 5) was fixed and equal to 62. The main factor that determines the flow resolution ( k max η) is fluid viscosity. To quantify this mutual relation six consecutive simulations of homogeneous isotropic turbulence, each with a different value of viscosity, have been performed. The basic parameters and flow statistics at the stationary stage of these flows (without particles) are listed in Table 1 . In addition to above mentioned quantities, Table 1 contains: the energy dissipation rate ε, the r.m.s. fluctuating velocity u , the Taylor microscale Reynolds number R λ = u λ/ν, the integral length scale L s , the transverse Taylor microscale λ, the eddy turnover time T e , the skewness S and flatness F of the fluid velocity gradient.
The data in Table 1 reveal the strict dependence between the resolution parameter ( k max η) and the numerical viscosity. As expected, larger viscous dissipation causes stronger suppression of small-scale motions. This effect can be quantified in terms of characteristic scales of the turbulent flow. Both the Kolmogorov length scale and the time scale increase for the increasing viscosity. Since the size of the computational domain is fixed (2 π in DNS units), larger Kolmogorov scales result in narrower energy spectra and consequently lower Reynolds numbers. It should be noted that for numerical modeling of cloud processes, the effect of Reynolds number may be important ( Rosa et al., 2013 ), especially for droplets with larger inertia. Dynamics of small-size droplets is dominated mainly by fine turbulent structures. In atmospheric clouds R λ is of the order of 10 4 and such value is a few orders of magnitude larger than in DNS ( ~10 2 ). Due to computational cost, achieving such high Reynolds numbers in simulations is not feasible. To maximize R λ in DNS it is necessary to reduce the viscosity parameter to a value that allows to maintain the stability of the numerical method. From this perspective, it is justified to keep k max η close to unity. On the other hand, k max η ≈ 1 results in large velocity gradients at small spatial scales. This in turn may have negative impact on the accuracy of the interpolation in PIC and PNN schemes. Therefore, in most of the simulations we used Flow C (see Table 1 ), for which k max η = 1 . 64 . In order to assess the accuracy of the interpolation methods (i.e. PIC and PNN) an additional set of two-way coupled simulations was performed at maximal value of R λ = 121 . Results from these simulations are presented in Appendix. We conclude that even at the highest R λ , both series of results are in good quantitative agreement.

Impact of particles on turbulence modulation
Motion of inertial particles in turbulent flow alters the energy transfer between different flow scales and thus affects the dynamics of the entire system. Because the dynamical features of the carrier fluid have a major impact on the collision statistics of the disperse phase, it is important to gain a closer insight into the statistical properties of the modeled flows. To quantify the strength of turbulence modulation by the droplets in two-way coupled systems a number of simulations have been carried out. The radii of tracked droplets varied between 20 and 60 μm. Basic properties of the droplets are specified in Table 2 . For converting physical units to spectral units, the kinematic viscosity was assumed equal to 0.17 cm 2 /s. The initial conditions in each DNS were set based on the energy dissipation rate from Table 1 . It should be noted that in simulations under two-way momentum coupling may depend on the particle mass loading. The parameter S V is defined as the ratio of particle still-fluid terminal velocity to the Kolmogorov ve-locity. Further, Table 3 contains data showing the mutual relation between the number of droplets, their radii, and mass loading for flows modeled with two different parameters of viscosity, corresponding to flow cases A and C. As it transpires from Table 3 , in all considered cases the mass loading of the droplets does not exceed 2.5. Equivalently, the maximal volume fraction of the water droplets is of the order of 1%. Therefore, we may assume that for most of the time particles remain far apart, so that including the exact representation of aerodynamic interactions (so-called fourway coupling) among them is not necessary. An important part of this analysis concerns gravitational effects. We address the role of gravitational settling by comparing results of simulations performed with and without gravity.
The first series of simulations was performed using different values of the kinematic viscosity, and identical liquid water content ( m = 0.24). The weighting factor was set to M = 1 . It is worth recalling that different viscosities in DNS yield different values of Kolmogorov scales (see Table 1 ). This, in turn, results in different size of computational domains in physical units. The reason for that is the translation of DNS units to physical units by matching the Kolmogorov scales (length and time). Since the actual domain size is different in each simulation the number of droplets must be suitably adjusted to obtain the same mass loading. Fig. 2 shows the energy spectra of turbulent flows computed in simulations at four different values of viscosity (marked with different line colors). As expected, the energy decreases with viscosity at wavenumbers greater than 2, but with Kolmogorov scaling ( E /( ν 5 ) 1/4 ), all the spectra collapse to one curve (see Fig. 3 ). Results from DNS at zero mass loading, i.e. simulations without dispersed phase, are plotted using solid lines. Dashed lines represent the energy spectra of flows modulated by droplets. The systems are monodisperse, which means that all droplets have the same size ( a = 40 μm). Two panels (a) and (b) correspond to simulations without and with gravity. The data were collected during the sta-tistically stationary stage, i.e. after at least 10 T e and then averaged over time.
Several important conclusions emerge from Figs. 2 and 3 . The kinetic energy at two lowest wavenumbers is identical in all simulations and does not depend on droplets mass loading or gravity. These values were preset in the algorithm for enforcing turbulent flow. As explained in Section 3 , the energy is supplied to the system at every time step to maintain constant level of kinetic energy at two first wavenumber shells. At higher wavenumbers the amount of energy largely depends on gravity. For nonsettling droplets a noticeable reduction of the kinetic energy occurs in the range of intermediate wavenumbers and enhancement is seen at high wavenumbers. The suppression of kinetic energy is a consequence of larger effective dissipation, while the increase is a combined effect of the larger viscous dissipation (at smaller scales) and transfer of momentum from the droplets to the fluid. Such phenomenon was observed in several previous studies, e.g.  ( Squires and Eaton, 1990;Elghobashi and Truesdell, 1993;Bosse et al., 2006 ) and is known as "pivoting".
Interestingly, for very small viscosity, the enhancement of energy at very high wavenumbers is not observed. It should be pointed out, however, that for this particular case the level of kinetic energy of the particle-free turbulent flow is larger than in other simulations. This observation may be a hint to understand why Squires and Eaton (1990) noticed the augmentation of the kinetic energy only in simulations at mesh 64 3 but not at coarser mesh 32 3 . The authors hypothesized that this may be related to the Reynolds number. In the light of new results this inconsistency can be explained as a peculiar effect of different settings of viscosity.
The energy spectra are significantly different if the gravitational settling is considered. In such a case, a large increase of kinetic energy takes place for both medium and high wavenumbers. This increase is due to larger transfer of momentum from particles to the fluid. Furthermore, the settling droplets induce larger velocity gradients in the fluid and thus act as an additional mechanism for enforcing turbulence. It is worth noting that in simulations with gravity, the amount of kinetic energy in the shell corresponding to k = 3 is slightly greater than that at k = 2. In this configuration, the forcing scheme acts as an absorber, which reduces the kinetic energy at larger turbulent scales.
To quantify the effect of turbulence modulation at different mass loadings, an additional set of simulations has been performed. For all these runs we used identical settings for the fluid (i.e. flow C, see Table 1 ), while the number of droplets varied from 1 to 8 millions (consistently M = 1 ). Since the particle inertia has significant impact on the preferential concentration and consequently turbulence modulation, the simulations were performed for different droplet sizes. The quantitative analysis of the obtained results is based on comparison of both the energy and dissipation spectra. Fig. 4 shows the normalized energy spectra of turbulent flows modulated by droplets of radii 30 μm and 40 μm. The simulations were performed for both non-settling droplets (panels (a) and (b)) and including the effects of gravity (panels (c) and (d)). In all considered cases we observe clear dependence of the spectra on the mass loading but in simulations with gravity this effect is more pronounced. It is also noteworthy that the effect of droplet mass loading is different at low (or moderate) and large wavenumbers. In the range of large wavenumbers the kinetic energy increases with the mass loading. This is due to larger momentum transfer from particles to the fluid. At moderate wavenumbers an opposite trend is observed, which can be attributed to a stronger effective dissipation caused by larger concentration of the disperse phase. Furthermore, the spectra are more sensitive to droplets of larger inertia. It should be emphasized that similar simulations but with even smaller droplets ( a < 30 μm) and in absence of gravity yield energy spectra which are not sensitive to droplet concentration (tested up to ~10 7 droplets).
The mass of considered droplets is relatively large, and therefore their motion is largely dominated by gravitational accelera-tion. Rosa et al. (2015) showed that gravity significantly affects the structure of particle clusters. The heavy droplets accumulate in the downward flow regions forming elongated (filament-like) structures. Thus, it is expected that the process of turbulence modulation is no longer isotropic. The 3D spectra plotted in Fig. 4 c and 4 d reveal a strong monotonic relation between the kinetic energy and the particle concentration. To better illustrate the effect of particle settling on the turbulent fluid also dissipation spectra were computed and visualized in Fig. 5 . A tremendous increase of the dissipation is observed in two-way coupled simulations and its magnitude again depends on the droplets mass loading. The results allow us to conclude that droplets transfer a significant amount of momentum to the system, and increase dissipation by enforcing large velocity gradients in the fluid.
To address more broadly the effect of droplet size on the turbulence modulation, an additional analysis has been performed. Fig. 6 shows the normalized spectra of energy and dissipation computed in simulations at different droplet sizes but at the same (rel- atively low) mass loading equal 0.1. The results prove that the size of the droplets (equivalently particle inertia) is of secondary importance to the energy of the system. Alternatively, it means that in terms of large-scale eddy turnover time, all particles have relatively small inertia. This observation is consistent with conclusions formulated by Squires and Eaton (1990) . In simulations with gravity we observe large difference in dissipation spectra between simulations performed under one-way and two-way coupling. The difference is more pronounced at larger wavenumbers, so we can conclude that this is effect of short-range interaction of relatively fast settling droplets with the fluid. The characteristic scales of the particle settling speed with respect to the rms fluctuating velocity are given in Table 2 .
In the subsequent steps of this analysis, the effects of two-way momentum coupling on the local particle distribution will be considered first (see Figs. 7-9 ). Then we address the effects of gravity on turbulence, the rms velocity of the particles, and particle energy budgets. The above analysis confirms that droplets may both enhance and suppress turbulent flows and thus affect their structure. The modulation of turbulence primarily depends on the droplet mass loading and gravity. In simulations without gravity the effect of two-way coupling results in reduction of the kinetic energy at larger scales and enhancement in the dissipation range. If gravity is considered the increase of energy takes place in the entire range of the spectra. In order to gain a deeper insight into the structure of the two-way coupled systems, detailed analysis of instantaneous flow fields is necessary. Therefore, 2D visualizations (cross-sections through domains) of the modeled flows along with locations of droplets have been performed. The used data were taken at the statistically stationary state, mostly at the end of simulations. The flow field is represented by the second invariant II of the velocity gradient tensor i j = ∂u i / ∂x j ( Squires and Eaton, 1990 ) where S ij is the rate of strain, and ω i -vorticity. Fig. 7 shows spatial distributions of II (in DNS units) computed in two-way coupled simulations with droplets of four different radii (20, 30, 40 and 50 μm). Locations of the droplets are marked by tiny black dots. In all simulations we used the same number of droplets (8 millions), which means the mass loading is different in each case. Gravitational settling was not considered. We avoided normalization of II to compare the relative differences between simulations with and without gravity. The blue regions in Fig. 7 corresponding to large and negative values of II indicate areas of high strain rate. The red regions where II is large and positive are regions of high vorticity. This qualitative comparison shows that the location of droplets is strictly correlated with structures of the turbulent flow. The spatial distribution of 20 μm droplets seems to be more uniform than distribution of 50 μm droplets. This is a combined effect of the Stokes number and two-way momentum coupling. As the droplet mass loading increases, more pronounced changes in the flow structure are observed. The extreme values of II remain similar, but the size of the smallest eddies becomes larger. This can be explained as an effect of the momentum transfer from droplets to fluid.
In absence of gravity, the motion of droplets does not have any distinctive direction. Therefore, the effect of turbulence modulation is largely isotropic. On the contrary, in two-way coupled simulations with gravity, the flow is modulated differently in the direction aligned with gravity and in the plane perpendicular to gravity. Therefore, in further analysis we address the differences in vortical structures formed along the vertical and horizontal directions. Fig. 8 shows a qualitative comparison of the flow fields (as in Fig. 7 ) in a vertical cross-section through the computational domain. Gravity is pointed down and aligned vertically. Next, Fig. 9 presents turbulent flows and locations of droplets in the plane perpendicular to gravity. The simulations with gravity were performed using the same settings (number of droplets, fluid viscosity, etc.) as the simulations without gravity. It should be pointed out that the values of II in Figs. 8 and 9 are significantly larger than in Fig. 7 . This is because settling droplets induce larger velocity gradients in the flow. Since the settling velocity depends on droplet inertia the increase of II strictly depends on droplet radii. It is worth noting that the size of vortical structures (in the sta-  tistical sense) is significantly smaller in simulations with gravity than in simulations without gravity. Moreover, the eddies become smaller as the droplet inertia increases. Interestingly, the pattern of turbulent flows in Fig. 8 reveals clear anisotropy. In simulations at large mass loading, eddies are elongated in the vertical direction. Such anisotropy is not present in Fig. 9 because gravity is directed perpendicular to the plane of this cross-section.
A more quantitative measure of turbulence modulation can be obtained from the flow statistics. Our attention is directed to characteristic time scales, energy dissipation rate and fluctuating velocity of the fluid. The statistics were computed in simulations with droplets of different radii and at different mass loadings. The data were averaged over time and normalized by the corresponding values computed in simulations under one-way momentum coupling. Fig. 10 shows (a) the Kolmogorov time scale and (b) the eddy turnover time obtained from simulations with and without gravitational settling. Although τ K and T e characterize turbulence features at different scales, the observed trends in Figs. 10 a and Fig. 10 b are rather similar. In simulations without gravity τ K increases as the mass loading and droplet inertia increase. This is consistent with the qualitative information presented in Fig. 7 , namely, larger τ K results in larger vortical structures. The increase of τ K is a consequence of larger viscous dissipation scales which is an effect induced by droplets. In simulations with gravity the trends are opposite, i.e. both τ K and T e decrease as the droplet mass loading increases.
The energy dissipation rate and the rms fluctuating velocity are presented in Fig. 11 a and Fig. 11 b correspondingly. Here, we observe a tremendous increase of both quantities in simulations with gravity. The enhancement of depends on droplet radii and droplet mass loading. This effect can be explained with larger velocity gradients caused by the fast settling droplets. This is also directly linked to larger values of u .
In order to assess how the fluid anisotropy develops when gravity is included the rms of droplets velocity V' has been computed. For comparison, similar computations were made using data ob- tained in simulations without gravity. The results are presented in Fig. 12 . If gravitational settling is not considered V' decreases with the mass loading and is lower for droplets of larger inertia. This is a direct effect of the suppression of fine turbulent structures by inertial particles. When the gravity is included, similar trend is observed but only for the horizontal component of V' ( ⊥ g ). The vertical component of V' ( g ) significantly increases with mass loading. This is a combined effect of droplet settling and stronger fluid vorticity at fine turbulent scales.
The large increase of the rms fluctuating velocity of the turbulent fluid ( u ∼ √ E ) in simulations with gravity (see Fig. 11 b) is likely due to larger energy transfer from the particles. Since the enhancement occurs for the settling droplets only it is expected that the slope of u in Fig. 11 b is related to the input of the potential energy from the particles. To confirm this interdepen-dence, the time change of the total kinetic energy ( E part ) of 40 μm droplets has been computed. We made use of the formula derived by Sundaram and Collins (1996) The computations were performed for both settling and nonsettling droplets. The time averaged values of temporal changes of energy are presented in Fig. 13 . As expected, the trends of time changes of the particle kinetic energy and u are in good agreement.
For a more complete description of the modeled processes, an additional analysis of the overall energy balance has been performed. This effort also aims at confirming the correctness of the numerical simulations. Again, we considered separately the two cases, i.e. with settling and non-settling 40 μm droplets. The results in form of the time evolution of different components of the temporal changes of energy are presented in Fig. 14 . In both cases we obtained good agreement (red and black lines almost overlap) with theory which says that the sum of additional energy from external forcing and kinetic energy from the particles should be equal to the total dissipation.
The simulation results show that gravity (droplet settling) has significant impact on the kinetic energy of the entire system. Furthermore, the motion of droplets can generate and enhance turbulent flow. Hence the question arises about mutual importance of two different mechanisms maintaining turbulent flows. These are: external forcing scheme (implemented at larger scales) and motion of the droplets (efficient at droplet/grid scales). To analyze this problem more thoroughly additional simulations were developed. In the new simulations the forcing algorithm was turned off and droplets were added to the stagnant air. It should be noted that the condition of fluid incompressibility (∇ · U (x , t ) = 0) and zero mean flow were maintained. The results are presented in Fig. 15 in a way consistent with previous visualizations. The sequence of plots shows that droplets are capable to generate turbulent flow even without any external forcing mechanism. The vortical structures of the flow are similar to those from Fig. 8 , but there is a distinct difference in the droplet distribution. In simulations without forcing scheme the spatial distribution of droplets is more uniform. This may be effect of another mechanism: so called preferential concentration (the situation where the spatial distribution  of particles is correlated to the local properties of the flow). In simulations without the large-scale forcing the initial condition for the flow is U (x , t = 0) = 0 so there is no mechanism for enforcing droplet clustering.
To complete the analysis, the spectra of kinetic energy ( Fig. 16 a), and dissipation ( Fig. 16 b) are also presented. In the range of larger wavenumbers, the values are in quantitative agreement with the results obtained from simulations with external forcing ( Fig. 4 d and Fig. 5 d). However, the difference is significant at low wave numbers, namely, there is no characteristic energy damping. This effect is expected because the level of energy at low wavenumbers is fixed only in simulations with forcing scheme.

Kinematic and dynamic collision statistics for inertial particles
In this section we discuss various effects of two-way momentum coupling on the particle collision statistics. The main focus is on the radial distribution function (RDF) and the radial relative velocity (RRV), see ( Rosa et al., 2013 ). The RDF is a local measure of the effect of preferential concentration of particles on the collision rate. The common method for computing the RDF( r ) involves counting the number of particle pairs at a given separation distance r . In our approach we considered a set of discrete values of r in the range R ≤ r ≤ 10 R , where R = 2 a is the collision radius. Then the RDF can be obtained by dividing the number of pairs at a given separation distance by the number of pairs characteristic for a nominally uniform distribution. According to definition, for monodisperse systems the RDF takes the form where n is the total number of particles in the computational box of volume V box . Then, n pairs is the total number of pairs detected at a separation distance r , falling in a spherical shell of inner radius equal to r − δ and outer radius equal to r + δ. δ is a small fraction ( ~1%) of collision radius and V s is the volume of the spherical shell. The discrete values of RDF( r ; t ), denoted henceforth by g 11 , are averaged over time and the best power-law fit allows to evaluate the RDF at contact g 11 ( R ). The RRV of two nearly touching particles is defined in terms of the relative velocity w in the limit r → 0 as w r (r ) = w · r / | r | . The methodology for computing w r between particles is similar to that for computing the RDF. More details on this method can be found in ( Rosa et al., 2013 ). The RRV and the RDF are directly proportional to the kinematic collision kernel These parameters are key statistical characteristics commonly employed to quantify the collision rate of the inertial particles in turbulent flows. In particular, they are often used to characterize the effects of air turbulence on the growth of cloud droplets during warm rain initiation.
Due to relatively low inertia of the cloud droplets the collision statistics are sensitive to the small-scale vortical structures. In simulations under two-way momentum coupling the dynamics of turbulent flow at small-scales, i.e. corresponding to the dissipation range of the energy spectrum, depends on both fluid viscosity and momentum transfer (from particles to the fluid). Therefore, in the first step, we compare the values of RDF and RRV computed in simulations with different viscosities (equivalently: different resolution parameter k max η). The two series of simulations were limited to droplets of radii 30 and 40 μm. Gravitational settling was not considered. In every simulation we tracked trajectories of 8 millions droplets. The collision statistics were computed "on the fly" at each time step and then averaged over time, as in ( Rosa et al., 2013 ). The postprocessed data are presented in Fig. 17 . The results show that the kinematic collision statistics indeed depend on the resolution parameter. The RDF of nearly touching droplets decreases as the fluid viscosity increases. This is owing to the fact that larger viscosity suppresses small eddies which have a key impact on the droplet clustering. An opposite trend is observed for the radial relative velocity. This, in turn, is the effect of the larger inertia of the droplets (droplet radii are the same). Such an apparent contradiction can be explained by referring to the settings of the numerical simulations. The characteristic scales of particles (e.g. radii) in DNS code are set based on matching the Kolmogorov scales. Since η depends on fluid viscosity the size of Fig. 17. Effect of viscosity (or, equivalently, the resolution parameter) and Reynolds number on the kinematic collision statistics of cloud droplets. Gravitational settling was not considered. All simulations were performed using 8 millions droplets of the same radii. the particles in spectral units is no longer the same. This leads to another conclusion, namely, for the same number of droplets the mass loading increases as the viscosity (and consequently η) increases.
In order to compare the RDF and RRV at the same mass loading and different viscosities additional simulations were performed. In these simulations the number of droplets was adjusted to obtain the same . The results are presented in Fig. 18 . We find that the collision statistics depend on k max η much less than in Fig. 17 . The little reduction of the RDF and increase of the RRV are observed only for droplets of larger radii. This stems from the fact that particles of lower inertia have little effect on turbulence modulation.  Large droplets alter turbulence more strongly, which in turn affects their dynamics and spatial distribution.
In the next step, we address more specifically the effect of droplet inertia, mass loading and gravity on the kinematic collision statistics. The analysis concerns monodisperse systems characterized by the same physical parameters, such as ε= 400 cm 2 /s 3 and R λ = 95 . To assess the effect of droplet inertia on the RDF we performed a number of two-way coupled simulations for different droplet radii and the same mass loading in each series. Fig. 19 shows the RDF of nearly touching droplets computed in simulations (a) without and (b) with gravity. Several important conclusions can be formulated based on these data. First, the RDF is more sensitive to the mass loading in simulations with gravity. Second, if the gravitational settling is not considered the effect of twoway coupling on the RDF is negligible for small (20 μm) and large (60 μm) droplets. This can be explained as follows: the particles with low inertia have little effect on the turbulence modulation, while the motion of particles with large inertia is not sensitive to the small velocity perturbation generated by other particles. It is worthwhile to emphasize that, at given mass loadings there are less larger particles in the computational domain. Third, in the intermediate range of 30 − 50 μm, a reduction of RDF is observed, and its magnitude is proportional to the mass loading. The reduction of RDF is a consequence of vortex suppression, which is the main mechanism causing inhomogeneity of particle distribution. Fourth, if the gravitational settling is considered there is little increase of the RDF with increasing for low inertia droplets (20 μm). This is due to formation of additional vortical structures by settling droplets. Fifth, the RDF of large inertia droplets is significantly reduced. The reason for the homogenization is the high vorticity of the background flow generated by rapidly falling drops.
The relation between the RDF and droplet number is presented more closely in Fig. 20 . The zero value on the X axis corresponds to simulations under one-way momentum coupling. Here, the ef- fect of gravity is clearly discernible. The RDF computed in simulations without gravity ( Fig. 20 a) decreases monotonically (for each droplet size/radius) as the droplet mass loading increases. The largest reduction of the RDF is observed for droplets of radii 40 μm. The difference between simulations under one-way coupling and simulation with 8 millions droplets exceeds 60%. When the gravitational settling is considered ( Fig. 20 b) the RDF increases for droplets of low inertia. For medium size droplets ( 25 − 30 μm) the RDF is not monotonic and reaches a maximum at the intermediate range of droplet number, while for heavy drops a strong reduction of RDF occurs even at very low droplet concentrations. This reduction is significantly larger than in simulations without gravity and for 50 μm reaches 90% (the difference between simulations under one-way coupling and simulation with maximal considered droplet concentration). This effect is due to increase of the relative velocity between droplets and will be further analyzed in detail.
Figs. 20 c and 20 d show a comparison of the RDF computed in simulations with different resolution parameters ( k max η). In other words, we analyze the sensitivity of the collision statistics to the range of turbulent scales ( R λ ). However, it should be emphasized that the droplets mass loading in the corresponding simulations (i.e. with the same droplet number) is not the same. The domain size in simulations at k max η = 1 . 13 is larger, which is a consequence of shorter Kolmogorov length scale. Fig. 20 c shows the RDF of non-settling droplets of four different radii. It turns out that the RDF is less sensitive to the droplet number in simulations at k max η= 1.13. This result is in line with expectations because larger resolution ( k max η) corresponds to smaller domain and consequently greater droplet mass loading. The effect of the flow resolution in simulation with the gravitational settling depends on the droplet radii. The RDF of 20 μm droplets is larger in simulations at k max η = 1 . 64 . This can be explained as a result of the turbulence enhancement. In turn for larger droplets the trend is opposite, namely, the RDF of 50 μm droplets is lower in simulations at A more intuitively appealing comparison of the RDF computed for droplets of different size and the same mass loading is shown in Fig. 21 . Again we consider two cases, i.e. simulations without gravity ( Fig. 21 a) and including gravitational settling ( Fig. 21 b). Due to high numerical complexity some of the simulations were performed using approximate model (i.e. with the weighting factor M = 1). The simplified method was used for systems with large and droplets of small radii. The simulations without gravity were performed for larger values of , up to 1.5. It should be recalled that according to the literature the approximation based on the two-way momentum coupling is accurate for < 1 only. The simulations with gravity were performed for a narrower range of , up to 0.5.
Based on these results, we conclude that the effect of two-way momentum coupling is important for both non-settling and settling droplets. If the gravitational settling is neglected the effect of turbulence modulation on the RDF is significant for medium size droplets and for a quite wide range of the mass loading up to 1. If the gravity is considered the effect of two-way coupling is stronger, but the RDF remains almost constant above > 0.3.
An analogous (to the RDF) analysis has been carried out for the radial relative velocity. Fig. 22 shows the RRV of nearly touching droplets normalized by the Kolmogorov velocity (evaluated for the turbulent flow without particles) for different . The data are consistent with those presented in Fig. 19 . In absence of gravity the differences in the RRV are rather small. The little enhancement is observed for the medium size droplets and its value increases with the mass loading. This is a direct effect of two-way momentum coupling and consequence of the turbulence modulation by moving droplets. The transfer of momentum from particles to the fluid occurs at grid-scale and increases the local fluid velocity. This in turn affects the motion of the neighboring droplets and causes a greater decorrelation of their velocity. If the gravitational settling is considered, we observe a tremendous increase of the normalized RRV. For 50 μm droplets the difference in the RRV between simulations under one-way coupling and these at = 0.1 is one order of magnitude. The reason for that is strong modulation of turbulent flow by fast settling droplets. Formation of small-scale vortical structures (as those in Fig. 8 ) of high angular velocity (vorticity) has a strong impact on the droplets relative motion and consequently alters the RRV. We point out that the increase of RRV is correlated with reduction of the RDF ( Figs. 20 c and 21 b).
Next, we address the relation between the RRV and droplet number. In absence of gravity ( Fig. 23 a) the RRV is weakly sensitive to the droplet concentration. A little increase is observed for droplets of radii 30 − 40 μm. For large droplets (60 μm) the trend is opposite and the RRV decreases with mass loading. For settling droplets, the RRV depends more strongly on the droplet concentration (note the log-scale on the Y axis). There is a continuous and monotonic increase with the mass loading and inertia.
Further, we analyze sensitivity of the RRV to the resolution parameter. Based on data presented in Figs. 23 c and 23 d we conclude that larger k max η results in larger RRV. This applies to the systems with both settling and non-settling droplets. However, in simulations with gravity the differences are greater. The results are somewhat counterintuitive because one may expect that a wider energy spectrum and the presence of small-scale eddies should enhance the RRV. In the present study such regularity is not observed. This is because the corresponding series of simulations were performed at different mass loadings. Fig. 24 shows additional comparison of the RRV for droplets of different radii and the same range of the mass loading. The twopoint collision statistics at large were computed using the approximate model. Interestingly, the results are in good quantitative agreement with these from exact simulations.
Another important aspect worth addressing is the magnitude of the radial relative velocity. In simulations with gravity, at ~0.5, the RRV of 50 μm droplets is about 100 times larger than v K . It should be noted, however, that the v K used for normalization in Figs. 20, 22, 23 and 24 is related to the single-phase turbulent flow. As the mass loading increases the Kolmogorov velocity increases, so that the ratio of the RRV to the actual v K should be smaller (approximately 2 times). But even with this scaling, the normalized values of the RRV are relatively large and their magnitudes exceed 2 u . Here we find an analogy with the variance of droplet velocities presented in Fig. 14 . We hypothesize that this increase of the RRV is a combined effect of large settling velocity and stronger vorticity of fine turbulent scales. It is expected that the large variation of the RRV and sensitivity to in simulations with gravity ( Fig. 24 b) will have a significant impact on the droplet collision rate.
The dynamical collision statistics such as dynamic collision kernel can be obtained from the simulation by detecting (directly) all collision events for a given time period. According to the definition D = ˙ n c /n 1 n 2 , the dynamic collision kernel is the ratio of collision rate to the average number densities of the two size groups of particles (for monodisperse systems n 1 = n 2 ). In absence of aerody-   namic interaction, the dynamic collision kernel matches the kinematic kernel (within statistical uncertainty) ( Rosa et al., 2013 ).
In Fig. 25 we compare the dynamic collision kernels computed in simulations at several values of mass loading and different droplet radii, while the weighting factor was set to 1. There is a notable difference between D computed for settling and nonsettling droplets. In the absence of gravity, the collision kernel decreases with the mass loading and droplet inertia. This reduction is due to lower values of the RDF or equivalently more uniform droplet distribution. If gravity is considered, D increases with the mass loading. Moreover, there is a noteworthy difference between simulations under one-way coupling and simulations at relatively low = 0.015. The increase is mainly due to enhancement of the RRV. Concurrently, the RDF decreases with . The difference be- tween simulations under one-way coupling and two-way coupling is particularly large for heavy droplets. Therefore, we can hypothesize that this is an effect of aerodynamic interaction between droplets. Namely, large drops settling under gravity affect strongly the fluid velocity and indirectly the motion of neighboring droplets (located closer than one grid spacing). This perturbation is more pronounced for larger droplets and consequently alters the twopoint collision statistics. In simulations under one-way coupling the aerodynamic interaction is not considered so the results are much different than results computed at low . In other words, there are different mathematical formulations to model these dynamical systems.
The relations between the collision kernel and the droplet number and mass loading are presented in Fig. 26 . The results confirm that in absence of gravity D , and consequently the collision rate, decrease with mass loading. This is due to more uniform droplet distribution (lower RDF). In simulations with gravity the RDF of large droplets also decreases with . However, in this case the dominant role plays the radial relative velocity. Enhancement of the RRV results in the increase of D up to one order of magnitude.
Finally, we examine the collision kernels of droplets in flows with different range of turbulent scales, equivalently different Reynolds numbers. Two sets of simulations have been performed using, as initial conditions, the flows A and C (see Table 1 ). The R λ values of these particle-free flows are 121 and 95 respectively. It should be added that R λ does not remain constant during simulations, because particles may generate additional vortical structures or suppress them. Moreover, different flows (A and C) are simulated in domains of different sizes so that systems with the same number of droplets have a different mass loading. This time the comparison is restricted to the settling droplets. Fig. 27 a shows D as a function of particle number for these two flows. In all simulations the weighting factor was set to 1. To expose more clearly the differences at the small size droplets the same data are presented in Fig. 27 b but using logarithmic scale. It is observed that D is consistently larger at smaller R λ . In the context of previous results,  it can be concluded that observed effect is mainly due to larger mass loading. The systems with heavier particles are characterized by larger relative velocities and this is reflected in larger D .

Effect of energy dissipation rate
All the simulations discussed in the previous sections were performed at the same value of energy dissipation rate equal to ε= 400 cm 2 /s 3 . Here, we extend the analysis and compare the kinematic collision statistics computed at different ε, in the range 100 − 400 cm 2 /s 3 . This range is typical of the cloud microphysical processes and particularly important for the rate of the precipitation formation. To simplify further analysis the gravitational settling was not considered. The computations were reduced to 30 μm and 40 μm droplets in radii, because the effect of two-way momentum coupling for these droplets was most meaningful.
The results of the above analyzed simulations prove that the main factor that determines the RDF is droplet inertia. The inertia is quantified by a non-dimensional parameter that is the Stokes number. In turn, St depends on the energy dissipation rate as follows St ∼ √ ε . Therefore, we expect reduction of the RDF at smaller values of ε. Results obtained in numerical simulations presented in Fig. 28 confirm these theoretical predictions. There is a systematic reduction of the RDF along with decreasing of ε. For particles of lower size (30 μm) the RDF linearly depends on the mass loading for all considered ε. For larger droplets (40 μm) this dependence is linear only for small St corresponding to ε = 100 cm 2 /s 3 . Similar verification has been done for the radial relative velocity of the nearly touching droplets. It is well known that the velocity decorrelation between fluid and particles depends on their relative inertia. Since the inertia of the fluid is fixed, we may expect larger RRV between particles with greater St . The numerical  results shown in Fig. 29 are in line with these theoretical considerations. It should be noted that the RRV in Fig. 29 is normalized by the Kolmogorov velocity scale which also depends on the energy dissipation rate v K ~( ε) 1/4 . Nevertheless, the increase of the RRV is dominant so the nondimensional quantity increases with ε.

Conclusions
The effects of two-way momentum coupling on the collisioncoalescence of water droplets have been examined using the combined Eulerian-Lagrangian numerical approach. The simulations have been performed for both sedimenting droplets and droplets without sedimentation. The main focus was on modeling kinematic and dynamic collision statistics at different droplet mass loadings.
The simulations have been carried out for droplets of radii in the range 20 − 60 μm. Moreover, we considered different turbulent Reynolds numbers, viscosity and energy dissipation rate. Several important conclusions can be drawn from this study. First, the two-way momentum coupling affects more strongly the dynamics of the systems with the settling droplets. This is mainly reflected in significant increase of the radial relative velocity at larger mass loadings. Second, the effect of two-way coupling on the RDF is rather complex and depends on particle inertia. For smaller droplets we observed little enhancement of the RDF which is a consequence of formation of additional vortical structures by settling droplets. The RDF of large inertia droplets is significantly reduced. Fourth, if gravity is not included, the RRV grows with mass loading and reaches a plateau for droplets of radii 60 μm. This may be due to strong flow perturbations at scales corresponding to the highest wave numbers. Fifth, the RDF of small (20 μm) and large droplets (60 μm) in simulations without gravity is not sensitive to the mass loading because droplets with low inertia have very little effect on the flow also because their mass fraction is lower for a given total number of droplets. In turn the motion of large droplets is weakly sensitive to flow perturbations generated by neighboring droplets. The effect of mass loading on the RDF is important mainly for medium size droplets ( 30 − 50 μm), where monotonic reductions systematically occur. Finally, we computed dynamic collision kernels at different mass loadings. It turned out that there is a fundamental difference between D of settling and non-settling droplets. In the absence of gravity, the collision kernel decreases with the mass loading and droplet inertia, which is mainly the effect of lower RDF. If gravity is considered, D increases with the mass loading.
At the end we outline possible perspectives for further research in this field. The present results are a step forward in quantifying the cloud processes but are limited to the monodisperse systems. Therefore, one potentially significant direction of future research is to consider a polydisperse systems under two-way coupling. Such systems more realistically describe cloud processes. In most former studies, as for example ( Saw et al., 2012a;2012b ) bidisperse systems were simulated assuming one-way coupling. Furthermore, new massively parallel supercomputers open new perspectives for performing such simulations at significantly larger resolutions. In the context of large-eddy simulations, following an earlier work Fig. A1. Comparison of (a) RDF and (b) RRV at contact computed using two different interpolation methods, namely, particle in cell -PIC and projection onto neighboring node -PNN. All simulations were performed at R λ = 121 . limited to one-way coupling ( Rosa and Pozorski, 2017 ), it may be worthwhile to examine the impact of filtering on the collision statistics. An important topic is also investigation of the settling velocity of small heavy particles under two-way coupling, especially for the development of realistic parameterization of thermodynamic processes in NWP models. locity gradients, we used, as the initial condition the flow with minimal value of viscosity (flow A, see Table 1 ).
Based on results presented in Figs. A.30 and A.31 we conclude that both interpolation approaches, i.e. PIC and PNN produce qualitatively similar collision statistics. The largest differences in the RDF are observed for the medium size droplets ( 30 − 40 μm). The RDF computed using the PIC method is about 5% lower comparing to results from simulations with the PNN method. Sensitivity of the RRV to the interpolation formulae is relatively low but slightly increases with the droplet inertia. Kinematic and dynamic collision kernels computed using different interpolations methods are in perfect quantitative agreement.