Effect of Gravity on Particle Clustering and Collisions in Decaying Turbulence

The preferential concentration of sedimenting particles in decaying homogeneous isotropic turbulence is investigated using radial distribution functions (RDF). Direct numerical simulations of polydisperse distributions of non-sedimenting and sedimenting particles of radii 10–55 μm are performed. We see a power law behaviour for the RDF in decaying turbulence and the power-law relation derived by Chun et al. (J Fluid Mech 536:219–251, 2005) for the RDF of non-sedimenting particles holds for sedimenting particles as well. Empirical formulas are generated for the power-law coefficients which are shown to be functions of the Stokes number and the Taylor Reynolds number for sedimenting particles. An in-depth analysis of the turbulent kinematic collision kernel for both non-sedimenting and sedimenting collision kernels confirms that gravity enhances the collision kernel for unequal sized particles and decreases for same-sized particles. Models are created for both monodisperse and bidisperse RDFs which are combined with existing models for the conditional radial relative velocities of colliding particles to predict kinematic collision kernels for both non-sedimenting and sedimenting particles. The effect on the collision kernel due to turbulence is also explored and enhancement of factors of up to three is observed with respect to the gravitational collision kernel.


Introduction
The clustering of intertial particles in incompressible turbulent flow has been established using experimental, numerical and theoretical investigations. Much of the early observation of the preferential concentration of inertial particles in turbulent flow was through observing the particle concentration field from direct numerical simulations (DNS) (Maxey 1987;Squires and Eaton 1991;Wang and Maxey 1993;Eaton and Fessler 1994). Particle clustering was also observed in a variety of experimental studies Wood et al. 2005;Aliseda et al. 2002;Salazar et al. 2008). The DNS and experimental studies have been complemented by numerous theoretical studies that have focused on describing and predicting the mechanisms behind the clustering of particles at sub-Kolmogorov scales. (Chun et al. 2005;Balkovsky et al. 2001;Bec 2003;Zaichik and Alipchenkov 2003;Gustavsson and Mehlig 2011). Understanding the dynamics of particles suspended in turbulent flows is valuable for a wide range of problems, with a majority of studies dedicated to understanding atmospheric cloud formation. In clouds, the implications of particle clustering may be significant; standard microphysics models are unable to explain the rapid cloud droplet growth in the intermediate size range of 15-40 μm where neither condensational growth or growth due to gravitational collision-coalescence are effective . It has been shown that turbulence intermittencies and an interplay between turbulence and gravity increases the inertial effects leading to an accelerated formation of larger droplets (Falkovich et al. 2002). However many open questions still remain, making this an active field of research. The physical mechanism behind particle clustering was first discussed by Maxey (1987) who showed that the velocity divergence of the particle flow field was positive in regions of high vorticity and negative in regions of high strain rate. Maxey (1987) predicted that (weakly) inertial particles would hence be 'centrifuged' out of vortex cores and accumulate preferentially in regions of low vorticity or high strain rate. Studies have since suggested alternative mechanisms by which particle clustering occurs, which includes (among others) the ergodic-non-ergodic mechanism (Duncan et al. 2005;Gustavsson and Mehlig 2011) and the 'sweep-stick' mechanism (Chen et al. 2006;Goto and Vassilicos 2006). More recently, in a series of papers (Bragg and Collins 2014a, b;Bragg et al. 2015a, b), the physical mechanism behind the clustering of inertial particles have been studied in detail across the dissipative and inertial scales. It was shown that the centrifuging mechanism is dominant for St ≪ 1 where St = p ∕ is the Stokes number based on the particle response time p and the Kolmogorov time scale . A non-local, symmetry breaking path-history mechanism takes over for higher Stokes numbers. This was also confirmed by DNS studies of inertial particles in isotropic turbulence by Ireland et al. (2016b) who observed preferential sampling for St < 0.1 and non-local effects for St > 0.2.
Different methods have been used to quantify and characterize inertial particle clustering, with a few popular examples being bin counting (Wang and Maxey 1993;Eaton and Fessler 1994;Wood et al. 2005), radial distribution functions (RDF) (Sundaram and Collins 1997;Reade and Collins 2000), Lyapunov exponents (Bec et al. 2006) and Voronoi analyses (Monchaux et al. 2010;Baker et al. 2017;Huck et al. 2018;Petersen et al. 2019;Momenifar and Bragg 2020). In this work, we will be working extensively with the RDF which is defined as the ratio of the probability of finding particle pairs at a certain separation distance to the probability of finding particle pairs at the same separation in a uniformly distributed particle field, and is one of the most widely used tools to study particle clustering. The RDF is an important tool as it provides a scale dependent quantification of 1 3 clustering and also appears as a correction to the kinematic collision kernel. Sundaram and Collins (1997) quantified the role of particle clustering on droplet collisions by multiplying the collision kernel calculated by Saffman and Turner (1956) with the RDF calculated at the point of collision for a monodisperse distribution of particles in homogeneous isotropic turbulence. DNS studies have shown that turbulence can enhance the collision kernel relative to stagnant flow through an increase in both the relative velocities of the colliding droplets as well as the increase in the RDF (Franklin et al. 2005(Franklin et al. , 2007Ayala et al. 2008b). Reade and Collins (2000) showed that the RDF can have values up to a 100 at separation distances equal to the collision radius further highlighting the importance of clustering towards the collision process. They modelled the monodisperse RDF g 11 using DNS data (at Taylor Reynolds number Re = 54.5 ) of non-sedimenting particles in homogeneous isotropic turbulence and suggested the following function: Here r = r∕ where r is the inter-particle separation distance and is the Kolmogorov length scale. The coefficients c 0 , c 1 and c 2 are all functions of St in principle. The study reveals an inverse power-law of the form r −c 1 at small r , and transitions into an exponentially decaying tail at higher r . Chun et al. (2005) derived an analytical expression for the RDF of a monodisperse collection of particles with St ≪ 1 as which is equivalent to Eq. (1). They showed that in a monodisperse distribution, the powerlaw dependence arises from a balance between an inward drift velocity caused by the particle inertia (inertial particles sample more strain than rotation giving rise to a net inward drift velocity) and a diffusion due to the random nature of the flow. For particles with different Stokes numbers, an additional diffusion arises due to the different particle response times to the local turbulent accelerations. This tends to homogenise the particle pair concentration leading to a reduction in the RDF. The additional acceleration diffusion term is independent of r and the resulting behaviour of the RDF is a power law dependence ( r −c 1 ) at larger scales and a transition to a plateau (independent of r ) at smaller scales. An expression is also derived for the cross-over length scale r c at which this transition occurs.
There are also numerous studies on particle clustering using tools from dynamical systems. These studies reveal that preferential concentration follows from the dissipative nature of particle dynamics (due to the Stokes drag) and in the position-velocity phase space, the particle trajectories converge to a dynamically evolving fractal attractor (Bec 2003;Bec et al. 2006Bec et al. , 2007. In such studies, the correlation dimension ( D 2 ) of the spatial distribution is often used to quantify particle clustering and is also related to the exponent c 1 . This is an important quantity and has been investigated in detail in Gustavsson et al. (2015).
The presence of gravity adds complexity and breaks the spherical symmetry in particle pair concentration as well as relative velocities between two neighbouring particles. The effect of gravity on two-point particle statistics is extremely important for cloud microphysics and has been studied using DNS for monodisperse distributions recently by Bec et al. (2014); ; Ireland et al. (2016a); Baker et al. (2017), and for bidisperse distributions by Ayala et al. (2008a); Woittiez et al. (2009);Dhariwal and Bragg (2018); Momenifar et al. (2019). The effect of gravity is quantified by the Froude number Fr which is defined as the ratio between the turbulent acceleration and the acceleration due to gravity. Bec et al. (2014) studied the interactions between (2) g 11 (r) = c 0r −c 1 , turbulence, gravity and particle sizes by combining DNS with theoretical results based on their asymptotic analysis. They observed that gravity acted in a non-uniform manner depending on the values of Fr and St. For Fr ≪ 1 , gravity suppresses clustering for St ≲ 1 and increases clustering for St ≳ 1 whereas for moderate values of Fr and St, clustering is enhanced.  analyzed the effect of gravity by using perturbation theory in the Kubo number (a dimensionless correlation time) and derived a theory to describe the clustering of particles falling through a turbulent flow. A model system was analyzed where particles are subject to gravity in a random two-dimensional velocity field. The inertial particle response to turbulence fluctuations and gravity are not found to the additive. They found a similar non-uniform effect of gravity behaviour and found that particles at small and intermediate values of St falling at finite Fr tend to cluster less due to the reduced correlation between particles and flow structures whereas for particles with higher values of St, clustering of settling particles is enhanced. Ayala et al. (2008a) extended the relations from Chun et al. (2005) for monodisperse sedimenting particles in a turbulent flow. The power law coefficients were assumed to be functions of a non-dimensional parameter for gravity g∕(u ∕t ) (this can also be interpreted as the inverse Froude number-Fr −1 , where g is the acceleration due to gravity, u and t are the Kolmogorov velocity and time scales respectively) in addition to St and the Taylor Reynolds number Re , and empirical expressions were generated by curvefitting to DNS data. A similar procedure was followed to obtain expressions for the bidisperse distribution and the cross over length scale r c . In addition to particle clustering, the effect of gravity on the relative velocities and preferential concentration of colliding particles was explored by Woittiez et al. (2009), Ireland et al. (2016a and Dhariwal and Bragg (2018) using DNS. The rate of collisions of sedimenting particles was found to be lower than non-sedimenting particles for a monodisperse distribution. For a bidisperse distribution, the higher relative velocities due to gravity enhanced particle accelerations resulted in a higher rate of collisions. In this study, particle clustering is analysed in decaying turbulence. In the unique setup adopted here, turbulent flow structures are evolving in both time and space and accordingly particles interact with evolving flow scales. Previous works mentioned in this section were carried out using deterministic forcing schemes which generated statistically stationary turbulence. The RDF which appears as a correction to the collision kernel is explored in significant detail in this work. This study marks the first time an attempt has been made to study the RDF across a range of scales in decaying turbulence to our knowledge. Additionally, we generate models for the kinematic collision kernel for both non-sedimenting and sedimenting particles. This is a significant contribution of this work. A complete analytical model is of course not created, however we have created a framework with a robust set of carefully chosen functions which can help determine the particle statistics. For this, we have relied on previously existing analytical and also generated new empirical relations using a series of non-linear regression analyses utilizing the significant DNS data generated specially for this work.
The paper is structured as follows. In Sect. 2, the case setup is presented along with the characteristics of the turbulence and the relevant droplet parameters. In Sect. 3, the results for both monodisperse and bidisperse RDFs are presented with comparisons to theoretical predictions. In Sect. 4, the effect of gravity on the geometric collision kernel is presented and models are created for the collision kernels. Section 5 includes a discussion of the results along with concluding remarks.

Case Setup
A novel decaying turbulence case is presented that closely resembles the flow in wind tunnel experiments. In wind tunnel experiments, grid-generated turbulence decays as the fluid moves downstream. The numerical case setup in this study aims to simulate the decaying turbulence by capturing the temporal evolution of a fluid section exiting the grid. Hence the simulations are temporal flow experiments and are similar to observing and following a fluid section through a camera which is being translated along an axis x at a velocity equal to the mean flow speed U.
The wind tunnel setup used by Bateson and Aliseda (2012) is taken as a reference for flow characteristics and grid dimensions. Simulating the motion of an active grid is difficult to achieve numerically. We circumvent this by simulating the turbulent flow just downstream of a passive grid with mesh spacing G = 0.1016 m . The average fluid velocity ũ just downstream of the grid is U g = Q∕(A − A obs ) between the grid, and zero at the grid. Here Q is the flow rate and A obs , A are the obstructed and total area of the grid respectively. The frame of reference for the numerical simulation is obtained by applying a Galilean transformation x → x − Ut where U = Q∕A . In the moving frame of reference, the initial condition u is given by The flow is then perturbed with random white noise and allowed to transition with time into isotropic, homogeneous turbulence.
The governing equations of the flow are given by the continuity equation and the incompressible Navier-Stokes equations where u = (u, v, w) is the velocity vector, P is the kinematic pressure and is the kinematic viscosity. Particle motion is governed by the simplified Maxey-Riley equations (Maxey and Riley 1983) where x(t) and v(t) are the particle position and velocity vector respectively, p = 2 p R 2 ∕9 is the particle inertial relaxation time with R the particle radius, p the particle density and the fluid density. Also, g = ge x where g is the acceleration due to gravity and e x is the unit vector along x . Here, x was chosen instead of the more conventional ẑ , as choosing ẑ caused heavy droplets to rain out in our smaller domain (compared to the wind tunnel cross-section) due to the non-periodicity in ẑ (more on this in the following subsection). Since the turbulence is homogeneous and isotropic, gravity can be defined in any direction.

Numerical Simulation Details
The DNS code SPARKLE is used, which solves the incompressible Navier-Stokes equations under the Boussineq approximation, and transport equations for scalars to fourth-order accuracy. Details of the numerical method used in SPARKLE and other details can be found in Craske and van Reeuwijk (2015). A tricubic Hermite interpolation scheme is used to interpolate the Eulerian flow velocities to the particle locations. A third order Adams-Bashforth timestepping scheme is used to advance the particle locations and velocities at each time step. The particles on processor boundaries are communicated across processors after each time integration step following the flow field (Perrin and Jonker 2015;Nair et al. 2021). The collision algorithm implemented in SPARKLE follows the general scheme in Allen and Tildesley (1987) where hard sphere molecular dynamics is solved. More details of the numerical implementation of the particle routine can be found in Nair (2021). Only particles that can collide within one time step are considered as prospective colliding pairs. Also, only binary collisions are considered and colliding particles are allowed to geometrically overlap post-collision (ghost collisions), i.e. no post-collision dynamics are calculated for the particles. Since we only look at geometric collisions, the effects of collision efficiencies are ignored. The code also calculates the kinematic collision kernel (using Eq. 16). The kinematic collision kernel has been validated against the dynamic collision kernel using DNS (Nair 2021, P. 67). A cubic domain is discretized using 1024 3 grid cells. Time-steps are dynamically calculated with a Courant-Friedrichs-Lewy (CFL) number of 0.4. The kinematic viscosity of air = 1 × 10 −5 m 2 s −1 . Periodic boundary conditions are imposed along the x and ŷ axes, and a free slip condition is imposed along the ẑ axis. The domain thus resembles a wind tunnel with frictionless side walls. The simulation was run for 1.1s. The computational domain is parallelized in the x and ŷ directions using 4096 processors.
The number of particles required to reach a target droplet volume fraction of O(10 −5 ) in a numerical domain whose sizes matches that of the wind-tunnel test section would be prohibitive (in terms of computational cost). Hence we calculate the two-point spatial correlation by determining the normalized spatial correlation tensor where u ′ i is the fluctuation from the mean for the instantaneous velocity field u. A length scale L 11 indicative of the size of the largest eddies in the flow can then be calculated by L 11 = ∫ L 0 R 11 dr , hence allowing us to determine domain sizes accordingly. We find that the value of L 11 = 0.06 m at t = 1s in a simulation without particles and a numerical domain of size 1.2 m × 1.2 m × 1.2 m . Accordingly, a numerical domain with a side length of 0.509 m is chosen which is around ten times the length scale L 11 and which is also almost half the size of the wind tunnel cross-section. This allows us to reach the target volume fraction with a computationally reasonable number of particles and also avoid any confinement effects due to the size of the domain. Table 1 shows characteristic values of the turbulent flow at x = Ut . The mean turbulent kinetic energy in the domain can be characterized as The velocity profiles averaged over the ŷ −ẑ plane and over a time interval Δt = 0.1s evolves and the effects of the grid slowly fade with time. After t = 0.6s , the flow can be considered to be isotropic. Also shown in Table 1 is the Taylor Reynolds number Re = u rms ∕ , with the Taylor micro-scale = √ 15 u 2 rms ∕ and the Froude number

Turbulence and Droplet Characterization
The evolution of the velocity variances u ′2 , v ′2 , and w ′2 with time are plotted in Fig. 1a. For isotropic homogeneous turbulence, these variances are equal to each other. The figure shows that the values are very close to each other after t = 0.5s thus attesting to the Table 1 Turbulence parameters at x = Ut for a simulation with an isotropic grid spacing Δx Here, is the dissipation rate of turbulent kinetic energy, , and u are the Kolmogorov length,time and velocity scales respectively, is the Taylor length scale, Re is the Taylor Reynolds number and Fr is the Froude number t = 0.5 s t = 0.6 s t = 0.7 s t = 0.8 s t = 0.9 s t = 1.0 s t = 1.  isotropy of the flow. The initial high values for u ′ compared to v ′ and w ′ are due to the effects of the initial conditions. Figure 1b shows normalized with the classical scaling u 3 e ∕L e . This value has been shown to be constant for decaying grid-turbulence in experiments (Sreenivasan 1984) and in forced as well as decaying DNS simulations of turbulence (Sreenivasan 1998). The magnitude of the constant C E = L e ∕u 3 e ≈ 1.8 shown in Fig. 1b is comparable with those from other DNS data where it falls in the range 1 − 3 (see Fig. 1 in the review by Vassilicos 2015). Figure 1c also shows the change of Re , Fr and St (for a droplet of radius 15μm ) as a function of time.
A polydisperse distribution of droplets with 36 different discrete radii between 10 and 55μm are randomly distributed across the whole computational domain at t = 0.4 s . This size range is selected since these are the most important sizes in the problem of the sizegap bottleneck in warm rain initiation. The number of droplets in each size category is decided by a volume-weighted distribution (such that the total volume of droplets in each size category is the same). Using such a distribution, a target volume fraction of 5 ×10 −5 is achieved by introducing a total of 62.5 million droplets. We assume that the system is well within the dilute limit and any modulation of fluid turbulence due to two-way coupling with the particles can be neglected. We also neglect particle-particle aerodynamic interactions.
The simplified form of the Maxey-Riley equations is valid only in the case where R∕ ≪ 1 . Table 2 shows the values of R∕ of six radius categories in the polydisperse system at = 215 cm 2 s −3 and t = 1 s . The values for all the categories are well within this limit. The added-mass, Basset history force and pressure forces are neglected under the assumption that p ∕ ≫ 1 . Since we are considering water droplets ( p = 1000 kg m −3 ) dispersed in air ( = 1.1436 kg m −3 ), p ∕ ∼ 1000 . Another important assumption is the Stokes-drag approximation for the viscous drag forces. This assumption is valid in the range where the particle Reynolds number Re p ≪ 1 . Ayala et al. (2008b) showed that for R > 30 μm , Re p is large enough to warrant the inclusion of a non-linear factor to model the departure from Stokes drag law. However, as was done in both Ayala et al. (2008b) and Woittiez et al. (2009), we employ the linear Stokes drag formulation to focus on the effects of gravity. Two simulations with the same flow characteristics and droplet parameters, but one with non-sedimenting particles and another with sedimenting particles, are performed.
Droplet-turbulence interactions are characterised using a series of non-dimensional parameters: the Stokes number St = p ∕ , the parameter p ∕ e (to characterize the interaction between particles and the larger turbulent flow scales), the settling parameter Sv = w t ∕u (to characterize the interaction between settling particles and the smaller scales of turbulence) where w t = p g is the droplet terminal velocity, and the parameter V T = w t ∕u e (to characterize the interaction between settling droplets and the larger scales Table 2 Parameters for particle turbulence interaction at = 215 cm 2 s −3 ( t = 1s)  Table 2 at t = 1s.

Radial Distribution Function
The three-dimensional RDF is calculated during the simulations using where 12 is the total number of droplet pairs with radii R 1 and R 2 with separation distance r detected in a spherical shell of inner radius equal to r − , outer radius equal to r + , and volume V s . Also, N 1 and N 2 are the total number of droplets of size R 1 and R 2 respectively. In this paper, = 0.02r . For an RDF of same-sized droplets, g 11 , N 2 is replaced with (N 1 − 1)∕2 . The RDF g 12 (r;t) is further averaged over time to obtain g 12 (r) . The calculation of the RDF and collision statistics is started at t = 0.5 s and averaged over every 0.1s thereafter. The very large number of particles ensures that converged particle statistics are obtained with this low averaging time. The maximum value of the inter-particle separation distance r over which the RDF is calculated is 4Δx . Considering particles at larger distances (beyond 4 cells in a direction) is restricted by the fact that the computational domain is parallelized in the ŷ and ẑ directions. Calculating particle separations at the processor boundaries requires exchange of data in the corresponding boundary cells of the neighbouring processor(s) which is a memory-and time-intensive procedure.

Monodisperse RDFs
As mentioned in the introduction, the RDF for same-sized particles can be expressed in terms of a power law. Figure 2a shows three dimensional RDFs in log-log coordinates for selected same-sized ( g 11 ), non-sedimenting (square markers) and sedimenting (plus markers) particles at t = 1 s ( = 215 cm 2 s −3 , Fr = 0.09 ). The RDFs for both non-sedimenting and sedimenting particles reach a saturation at r∕ < 1 at St = 0.62 (yellow squares and plus markers) after which the values start decreasing for higher St. However the magnitude of the maximum is higher in the absence of gravity. We observe a clear power law scaling for both cases for values of r∕ upto 4. The RDF at r∕ < 1 remains the same up to Sv ≈ 2 after which sedimenting particles show a lower RDF. At these intermediate St, gravity clearly suppresses particle clustering. At higher Stokes numbers, gravity increases the RDF for r ≪ as can be seen for St = 1.03 . To explore this in detail, the RDFs g 11 for droplets of size 40.8 μ m are shown as a function of time Fig. 2b. The figure shows how a droplet of a particular size interacts with the evolving flow scales of decaying turbulence. The power law is maintained throughout decaying turbulence with the slopes of the straight line changing as a function of time. As the particle is suspended in decaying turbulence the Stokes numbers vary from 1.97 at t = 0.6 s ( Fr = 0.34 ) to 0.81 at t = 1 s ( Fr = 0.09 ). The rate of change of the slopes is high corresponding to the change (of almost an order of magnitude) in from 0.1245m 2 s −3 (at t = 0.6 s ) to 0.0712 m 2 s −3 (t = 0.7 s ). From t = 0.8s to 1s , the change in is very small which corresponds to a small change in the slope. An interesting feature from these plots is the role gravity plays on the clustering or RDF as a function of both Fr and St.
For small values of Fr, which indicates a larger effect of gravity accelerations, the clustering is lesser (cyan and pink plus markers) as opposed to the non-sedimenting case ( Fr = ∞ , cyan and pink squares). As Fr increases, indicating weaker gravity accelerations compared to higher turbulent accelerations, particle clustering is higher for sedimenting particles (red, black and blue markers) for r∕ < 1 . Gravity reduces both preferential concentration and the path-history effects (Ireland et al. (2016a)). At low St, gravity weakens the preferential sampling drift effect leading to a reduction in the overall clustering. At higher Stokes numbers, gravity increases the RDF slightly by weakening the path-history effect whose effect is normally to decrease clustering. The coefficients c 0 and c 1 of the power law from Eq. (2) for the RDF values from t = 0.7s to t = 1.0s are shown in Fig. 3 for non-sedimenting (a, c) and sedimenting (b, d) particles. These are obtained by fitting Eq. (2) to the DNS data. The theory behind Eq. (2) was developed for steady-state realised by driven turbulence. Figure 2 clearly shows that a power law behaviour exists for decaying turbulence as well and Fig. 3 shows the coefficients c 0 and c 1 for the same. The dependency of the coefficients c 0 and c 1 on St is similar to that shown in Reade and Collins (2000). At low Stokes numbers, the coefficients are similar for both sedimenting and non-sedimenting particles. However it is interesting to note that for the coefficient c 0 , the peaks decrease with time for sedimenting particles (Fig. 3b), i.e. particles with the same Stokes number cluster to a lesser degree. The corresponding Fr decreases with time, indicating an increase in the gravitational accelerations compared to turbulent accelerations leading to a reduction in clustering in the intermediate St range of 0.3-1.
An empirical fit for the coefficients c 0 and c 1 was performed by Reade and Collins (2000) by curve fitting the following functions to their DNS data The analytical solution developed by Chun et al. (2005) was extended to sedimenting particles by Ayala et al. (2008a) who developed an empirical model using DNS results to quantify the RDF. In order to account for the effects of gravity, they assumed the coefficients c 0 and c 1 to be a function of St, Re and a dimensionless gravity parameter |g|∕(u ∕ ) . An empirical expression for c 1 for sedimenting particles was then obtained by curve-fitting to the DNS data with c 0 = 1 . However, Fig. 3 clearly shows that c 0 is a function of St. Even though the range of Re is small in our simulations, there is a clear change in the maxima for c 0 (Fig. 3b) as the simulation progresses. This behaviour can be well represented by the empirical relation where (7)

3
Curve-fitting this function with the data plotted in Fig. 3b gives the following parameters: b 0 = 11.7825, b 1 = −9.3844, b 2 = 0.8338 and b 3 = 0.9738 . The behaviour of coefficient c 1 is captured by (Ayala et al. 2008a) where This function is fitted to the data shown in Fig. 3d for c 1 and results in the following parameters: b 4 = −0.2817, b 5 = 1.6703, b 6 = −3.3334, b 7 = 2.3639, b 8 = −0.0632 and b 9 = −9.9914 . These values vary from those obtained by Ayala et al. (2008a), however the parameters for c 1 are of the same sign. The potential of using these empirical formulas for c 0 and c 1 to predict the RDF of colliding droplets will be explored in Sect. 4. Figure 3 also shows how the the coefficients change as the droplets interact with decaying turbulence. Collins and Keswani (2004) reported that with increasing Re (their study reached Re values up to 152), the RDF saturates and approaches a plateau which was in contrast to the linear growth with Re that was reported by Wang et al. (2000a). Falkovich and Pumir (2004) also reported an increase in the exponential coefficient c 1 with increasing Re and dissipation . This was attributed to an increase in intermittency in the gradients of the turbulent velocities with increasing Re . However, Saw et al. (2008Saw et al. ( , 2012b showed using wind tunnel experiments that there is a strong saturation in c 1 in the limit of large Re and the good agreement between experimental ( Re = 440 ) and DNS results ( Re = 140 , Saw et al. 2012a) is also attributed to this saturation and the weaker dependence on Re as compared to St. A Stokes number 'similarity' was also observed in the RDFs measured in Saw et al. (2008Saw et al. ( , 2012b where the plots of RDF for droplets with different sizes but same St (from measurements taken at different downstream locations of the wind tunnel) coincided even though these were calculated from different flow conditions. We see a similar St similarity for the exponential coefficient c 1 which is the slope of the RDF. It is possible that this reinforces the observation from Saw et al. (2012b) that particle clustering occurs at scales that are well separated from the inhomogeneity in the larger scales. But we are hesitant to present this as a conclusion as we do not attain the maximum or the range of Re in our numerical simulations ( ∼ 50-67) as in the experiments. Most of the previously stated DNS studies used different simulations of forced turbulence (at different turbulence intensities and hence Re ) to study the Re dependence. Future studies where results from our novel setup will be compared to experimental data has been planned where we will aim for the range of Re seen in wind tunnel flows. An attempt to provide a definitive answer for the Reynold's number dependence to clustering will be made in the future.
To conclude the analysis of the monodisperse RDFs, we plot the correlation dimension which is used extensively in dynamical systems theory (Bec et al. 2008;Gustavsson et al. 2015) in Fig. 4. The correlation dimension D 2 can be related to the exponent c 1 as D 2 = d − c 1 where d = 3 is the number of spatial dimensions. At low St, there is a decrease in the value of D 2 with increasing St which points to an increasing spatial correlation, i.e. a higher degree of clustering. This decrease is followed by the value of D 2 reaching a minimum between St = 0.4−0.6 indicating maximum particle clustering which is consistent with results from experiments (Monchaux et al. 2012). The effect of gravity is clearly seen at the intermediate St range where the value of D 2 is higher for sedimenting droplets (plus markers) indicating a lower spatial correlation and hence a lower degree of clustering. At high St, we see lower values of D 2 indicating a higher spatial correlation and hence a higher degree of clustering for non-sedimenting droplets (square markers) compared to sedimenting droplets.

Bidisperse RDFs
The analytical theory of Chun et al. (2005) predicts the RDF for inertial particles of different sizes using the following expression at r ≪ where c 0 and c 1 are interpreted similarly as in the monodipserse case (Eq. (2)) and r c is a cross-over length scale. They showed that shear-driven diffusion dominates for ≫ r ≫ r c with the resulting RDF exhibiting a power law behaviour, while an acceleration-driven diffusion dominates for r ≪ r c resulting in an RDF which is independent of r. Hence r c can be considered to be a length scale below which the power-law behaviour of the RDF crosses over to a plateau. A semi-empirical relation for r c (for Re = 47.1 ) is given by (Chun et al. 2005) and for c 1 as Eqs. (11) and (13) for different ΔSt . The bidisperse RDFs exhibit a power law at relatively larger scales which taper off to a plateau at smaller scales. It is also clear that this transition to a plateau where the RDF is independent of r clearly happens at much larger scales for sedimenting droplets (for the same ΔSt ). We observe a similar asymmetry in g 12 for non-sedimenting droplets when St 2 is increased or decreased relative to the fixed St 1 , as was shown by Saw et al. (2012a). Interestingly, a symmetry is observed in g 12 for sedimenting droplets. Any small deviation from a monodisperse distribution results in a sharp change in the RDF slope due to the difference in settling velocities for droplets with different sizes (or for small ΔSt = St 2 − St 1 ) irrespective of St 2 being greater or lesser than St 1 . This is explored in detail for colliding droplets in Sect. 4. Saw et al. (2012a) observed that for non-sedimenting droplets, the power-law coefficient c 1 can be recovered only at sufficiently small scales such that r ≪ , which makes the retrieval of c 1 difficult when the value of r c approaches for very small ΔSt . This happens in our study as well, as is evident by the dashed lines in Fig. 5. Saw et al. (2012a) extended the bidisperse theory and proposed a new power-law coefficient c ′ 1 which corresponds to the slope of the RDF in the range of scales between r c and the large-scale flattening of the RDF ( r∕ ≈ 1 − 6 in their study) and approximated this as the lesser of the two corresponding monodisperse values, i.e. This relation was shown to hold for ΔSt = |St 2 − St 1 | of the order of 0.01-0.1 and arises as a result of a saturation effect (with c ′ 1 increasing with St 2 to a maximum and thereafter remaining constant) which suggests that the value of c ′ 1 is always limited by the least clustered particles. We replace c 1 with c ′ 1 in Eq. (10) and calculate r c for non-sedimenting droplets by performing a nonlinear fit of Eq. (10) to the DNS data for bidisperse RDFs (for the entire range of r). The resulting values of r c are shown in Fig. 6 for different St 1 at t = 1s and are compared with the theoretical values (Eq. (11)). There is a very good agreement between the two values suggesting that the extension proposed by Saw et al. (2012a) to the theory of Chun et al. (2005) works well even for r∕ < 1 . This is important as Eq. (15) can be used to determine the collision kernel as will be done in Sect. 4. Unfortunately, extending this theory to sedimenting droplets is not straightforward and an attempt to do so failed since r c ∕ approaches 1 for really small ΔSt , and we do not see a power-law behaviour at larger scales (within 4 ) as is clear in Fig. 5b, d. We also highlight alternate results from Bec et al. (2005); Meibohm et al. (2017aMeibohm et al. ( , 2017b; Bhatnagar et al. (2018) which points towards r c scaling as |St 1 − St 2 |∕(St 1 + St 2 ) . We performed an alternate analysis to retrieve r c from the DNS data as before but with St 2 ) is the harmonic mean of the two Stokes numbers (again motivated by results from Bec et al. 2005;Meibohm et al. 2017a, b;Bhatnagar et al. 2018). The resulting values for r c is similar to that obtained above at small ΔSt and the same issue with r c approaching at small ΔSt prevented us from looking at the actual scaling of r c at larger range of ΔSt for non-sedimenting particles and for any ΔSt for sedimenting particles.

Geometric Collision Kernel
In this section, we analyse the effect of gravity on the collision kernel. The most common kinematic formulation for the geometric collision rate of a bidisperse distribution is given by (Sundaram and Collins 1997;Wang et al. 2000bWang et al. , 1998, where R c = R 1 + R 2 is the collision radius for two colliding droplets of radius R 1 and R 2 . The radial relative velocity w r is defined in terms of the relative velocity w between two droplets with separation vector r as w r = w ⋅ r∕|r| . The mean radial relative velocity ⟨�w r (r = R c )�⟩ is then the conditional average, i.e. the mean of |w r | conditional on r = R c . Figure 7 shows contour plots for the collision kernel Γ K (a, d), RDF g(R c ) (b, e) and the magnitude of the mean radial relative velocities ⟨�w r (R c )�⟩ (c, f) for all colliding droplets in the polydisperse distribution. For all the plots in this section, the statistics are shown for t = 1 s . It should be noted that for this figure, the DNS data for the RDF g(R c ) is calculated separately in the collision routine of SPARKLE and is not extracted from the RDF data used in Sect. 3 (details in Nair (2021)). All the figures show a striking symmetry which confirms the fact that the statistics for collisions between droplets of radii R 1 and R 2 are the same as those between R 2 and R 1 . Figure 7a, d show that sedimentation enhances the collision kernel except when collisions are between droplets of the same size. For collisions between same-sized particles (along the leading diagonal), the kernel is higher for non-sedimenting droplets (Fig. 7a) compared to sedimenting droplets (Fig. 7d), noticeably so for R > 25μm . For collisions between unequal-sized droplets (either side of the leading diagonal), the enhancement observed in Fig. 7d is obvious for R 1 > 25 μm, R 2 > 25μm . The higher velocity decorrelation due to the differential settling velocities for unequal-sized droplets becomes significant at sizes greater than 25μm as can be seen by comparing Fig. 7c, f in the same size range. The figures reveal a difference of at least an order of magnitude between ⟨�w r �⟩ for nonsedimenting and sedimenting colliding droplets. For both cases, bands are observed within which the collision kernels have a similar order of magnitude. Detailed analysis of the kernels are done later in the section using a cross-section from the contour plot.
A detailed analysis of the RDF at scales extending up to 4 was performed in Sect. 3. In Fig. 7b, e, we look at the effect of gravity on the RDF of colliding droplets, i.e. at r = R c . Figure 7b suggests that the gradient from the monodisperse peak is much more gradual for non-sedimenting colliding droplets when compared to sedimenting colliding droplets. In Fig. 7e, the RDF of sedimenting droplets falls rapidly from the monodisperse value for even a small ΔR = R 2 − R 1 . This is due to the rapid decorrelation of the droplet concentrations for different droplet sizes (Zhou et al. 2001). The decorrelation arises from the fact that particles with different inertia tend to cluster in different (slightly-offset) regions with respect to the flow. Differential sedimentation enhances this decorrelation leading to a very sharp diagonal with high values of RDF and the values dropping to one even for a small ΔR . In Fig. 7c, f, the plots are symmetric indicating that ⟨�w r �⟩ depends on ΔR (or ΔSt ), the difference in inertia, rather than R 1 or R 2 . The value of ⟨�w r �⟩ reaches a minimum for same-sized droplets and increases with ΔR for both cases. The effect of gravity is an overall increase in ⟨�w r �⟩ of almost an order of magnitude for unequal-sized droplets. The RDFs are now analysed separately in Fig. 8 by comparing the DNS data with values predicted from different analytical/empirical models for both non-sedimenting and sedimenting droplets. The DNS data in Fig. 8a for the monodisperse RDF g 11 (R c ) corresponds to the leading diagonal of Fig. 7b and that for the bidisperse RDF g 12 (R c ) in Fig. 8b corresponds to a slice through the contour plots at R 1 = 31.9μm(St = 0.5) . For g 11 (R c ) , the obvious model to predict the RDF would be using Eq. (2) where the coefficients c 0 and c 1 are calculated using the equations and fitted parameters from Sect. 3, i.e. Eq. (7) for nonsedimenting droplets and Eqs. (8) and (9) for sedimenting droplets. The fitted parameters Fig. 7 Comparisons between the polydisperse distribution of non-sedimenting droplets (top row) and sedimenting droplets (bottom row) for the (a, d) kinematic collision kernel Γ K , (b, e) RDF and (c, f) radial relative velocity ⟨�w r �⟩ 1 3 d 0 − d 7 and b 0 − b 9 for these equations remain the same as in Sect. 3. Similarly, the bidisperse RDF g 12 (R c ) can be predicted using Eq. (10) where r c is calculated using Eq. (11) for non-sedimenting droplets and (13) for sedimenting droplets. The coefficients c 0 and c 1 are selected from the corresponding monodisperse values as min(c 0 (St 1 ), c 0 (St 2 )) and min(c 1 (St 1 ), c 1 (St 2 )) for non-sedimenting droplets. For sedimenting droplets these are calculated using Eqs. (8) and (9) with St = max[St 1 , St 2 ] . This model is referred to as model M1. The corresponding values obtained using model M1 (dashed lines) are compared with the DNS data (circles) in Fig. 8.
For the monodisperse case, there is very good agreement between the DNS data and model M1 for colliding non-sedimenting (red) and sedimenting droplets (black) as shown in Fig. 8a. For sedimenting droplets (black dashed lines), a good agreement is observed for smaller St even though the peak seems to be shifted with the model slightly over predicting for St > 1 . For the bidisperse case, there is a good agreement for the non-sedimenting case for all sizes except for when R 2 > 35 which corresponds to ΔSt > 0.3 . Note that this is much better behaviour than expected since the relation given by Eq. (15) was shown to be valid for ΔSt up to 0.1 only. However, for sedimenting droplets, there is a marked disagreement for even small ΔSt . This is due to the values of r c ∕ rapidly approaching 1 even for small ΔSt which renders the retrieval of the power-law coefficient (and ultimately the RDF) using the model impossible.
We proceed to model the RDF using the procedure described in Ayala et al. (2008a) to investigate if this gives a better agreement with DNS results. This model will be named M2. We first curve-fit Eq. (2) to our DNS data for the monodisperse RDF g 11 (R c ) setting c 0 = 1 and c 1 following Eq. (7) for non-sedimenting droplets, and Eq. (9) for sedimenting droplets. This fit results in a new set of parameters for c 1 which are: . For bidisperse RDF g 12 (R c ) , this process is extended by assuming that r c follows Eq. (13) for both cases but with F(a 0g , Re ) replaced by F(Re ) = 11.5283 a 0 ∕Re 1∕2 and a 0 = (9.5686 + 0.1303Re )∕(205.29 + Re ) for nonsedimenting droplets. For sedimenting droplets, F(a 0g , Re ) has the same form and parameters as Eq. (14). The coefficient c 0 = 1 , and c 1 is calculated similar to model M1. Very good agreement is observed for the monodisperse case for non-sedimenting droplets (red dotted lines) and sedimenting droplets (black dotted lines). An excellent agreement is observed for the bidisperse case for both non-sedimenting and sedimenting droplets. Details of both models and the respective equations and parameters are summarized in Table 3.
To obtain a clearer picture of the droplet collision process, the kinematic collision kernel Γ K 11 , the RDF g 11 , and the radial relative velocities ⟨�w r �⟩ for colliding droplets of the same size are plotted in Fig. 9a-c. These plots correspond to the leading diagonals of Fig. 7a-c respectively. Figure 9d-f show the collision statistics Γ K 12 , the RDF g 12 and ⟨�w r �⟩ for unequal-sized droplets where St 1 = 0.5 and St 2 is varied. This corresponds to a slice through the contour plots at R 1 = 31.9 m in Fig. 7. Also plotted are the theoretical/model values for comparison. Model M2 is used for comparison with the RDF values obtained from DNS in Fig. 9b, e. The values of ⟨�w r �⟩ obtained from DNS are compared with models from Wang et al. (2000a) for non-sedimenting droplets and Ayala et al. (2008a) for sedimenting droplets. For non-sedimenting droplets, this is given by where where = max[ 2 ∕ 1 , 1 ∕ 2 ] , = 2.5 p ∕T e and w shear , w accel are the radial relative velocities due to the shear and acceleration mechanism respectively. For sedimenting droplets, the following expression for the formula from Ayala et al. (2008a) for ⟨�w r �⟩ is used (2) )⟩ and ⟨(v �(i) ) 2 ⟩ 1∕2 is the rms fluctuation velocity of the ith particle. The kinematic collision kernels Γ K 11 and Γ K 12 can then be predicted using model M2 to calculate the RDF, and Eqs. (17) and (18) for ⟨�w r �⟩ for nonsedimenting and sedimenting droplets respectively. Figure 9b shows that for same-sized colliding droplets, the RDF reaches a maximum at around St = 0.5 after which it reduces steadily. In the range St = 0.28 to St = 0.86 , gravity clearly suppresses particle clustering due to reasons discussed in Sect. 3.1. Figure 9e shows the RDFs for the bidisperse case for St 1 . The sharp fall of the RDF values from a maximum to one is clearly visible for sedimenting droplets (black circles) when compared to the more gradual fall for non-sedimenting droplets (red circles). Again, model M2 shows very good agreement with the DNS values.
Sedimenting droplets of the same size have a lower ⟨�w r �⟩ than their non-sedimenting counterparts as shown in Fig. 9c. Gravity effects the relative velocities of colliding particles through preferential sampling of the fluid velocity field and the temporally nonlocal path-history of the particle pair (Ireland et al. 2016a;Dhariwal and Bragg 2018). Due to particles preferentially sampling the fluid velocity field, they experience high differences in fluid velocities. When ΔSt = 0 , and for small values of St, gravity reduces the preferential sampling effect which leads to particles having dynamics closer to that of the fluid and a resulting reduction in the relative velocities. At higher St, non-sedimenting droplets experience a higher influence of the path-history effects. The particles remember the interactions with the different turbulent scales which are significantly different from that at the current separation. The contribution to the relative velocities for a particle pair will be dominated by the fluid velocities at larger separations along their path histories. This results in relative velocities that are higher than the local fluid velocity difference at the point of collision. At higher St, gravity reduces the path-history effect causing a lower relative velocity value for sedimenting particles compared to non-sedimenting particles. This can be seen by the higher difference between the two cases for St > 1 . When ΔSt > 0 , gravity can enhance inertial particle accelerations leading to values of relative velocities higher than in the non-sedimenting case (Dhariwal and Bragg 2018). This is clear from the order of magnitude difference between the values in Fig. 9f. The enhancements have been shown to be due to large settling velocities of particles with finite St and low Fr, which causes a rapid change in the fluid velocities experienced by the particles. Studies by Ireland et al. (2016a), Dhariwal and Bragg (2018), Momenifar and Bragg (2020) have explored these in detail.
A very good agreement between the model and DNS data is observed for same-sized droplets in Fig. 9c for both the non-sedimenting and sedimenting cases. An excellent agreement is observed for sedimenting droplets for unequal-sized droplets, whereas for the non-sedimenting case, the semi-empirical formula from Wang et al. (2000a) seems to over-predict the values of ⟨�w r �⟩ for small sizes. The kinematic collision kernel calculated using (16) is shown in Fig. 9a. Gravity clearly decreases the monodisperse collision kernel for St > 0.3 by influencing the accumulation effect(particle clustering) and the turbulence transport effect (particle relative velocities). The model for the monodisperse kinematic collisions kernel Γ K 11 works very well for both cases. For the bidisperse case, there is a slight over-prediction of Γ K 11 for St 2 < 0.5 as a result of the model for ⟨�w r �⟩ not exactly matching the DNS values. However, for the sedimenting case, there is an excellent agreement between the model and DNS.
It is important to highlight that even though the formulation of the kinematic collision kernel according to Eq. (16) suggests that particle clustering explicitly enhances the collision rates, an analysis by  suggests this may not be the case.  show that the collision rate is equal to m 1 (R c )∕2 where m p is the p-th moment of w r . They argue that since ⟨�w r �⟩ is defined as a conditional average (conditional to r = R c ), ⟨�w r �⟩ can be expressed as m 1 (R c )∕m 0 (R c ) and with g(R c ) ≈ m 0 (R c )∕R d−1 (where d is the number of spatial dimensions), by definition (Eq. 16), any direct enhancement to the collision kernel from clustering through g(R c ) will be cancelled by the probability that is present in the formulation of ⟨�w r (r + R c )�⟩ and the actual collision rate will be given by m 1 (R c ) . Also, any contribution from spatial clustering g(R) will be indirect; as the fractal evolves in phase space, the velocity moment m 1 may be influenced by the same degree of fractal clustering in velocity. This could explain why in Fig. 9d, the collision kernel appears to be smoother compared to g 12 and ⟨�w r �⟩ , i.e. the non-smooth parts are cancelled out. For sedimenting particles, we speculate that the contributions of differential settling likely makes m 1 (R c ) and hence the collision rates nonsmooth as seen in Fig. 9d. An analysis to confirm this is reserved for future work.
The enhancement of the geometric collision kernel due to air turbulence relative to the gravitational kernel is shown in Fig. 10a. The figure shows the ratio Γ K ∕Γ g where Γ g is the gravitational collision kernel. Also shown is the RDF (Fig. 10b) and the enhancement of the average radial relative velocity with respect to the purely gravitational case (Fig. 10c). Turbulence-induced enhancements of up to a factor of 3 are observed, especially for Fig. 10 a Turbulence enhancement of the geometric collision kernel Γ K , b RDF and c Turbulence enhancement of the radial relative velocity ⟨�w r �⟩ Table 3 Different models used for predicting RDFs for colliding droplets collisions between droplets with a small size difference. Devenish et al. (2012) showed enhancements between 1.5 and 4 at similar and size range. The value of Γ K ∕Γ g is not defined for particles of the same size (since Γ g = 0 ), however Fig. 10b shows that particle clustering can lead to enhancement in collisions up to 12 times that in the purely gravitational case. For larger droplets, as particle inertia increases, the effect of radial relative velocity is the significant contributor to the turbulence enhancement of the collision process. A difference in radial relative velocities between colliding particles can lead to enhancements of up to a factor of 2 for droplets of unequal size as shown in Fig. 10c.

Concluding Remarks
In this paper, we investigated the effect of gravity on the preferential concentration and geometric collisions of a polydisperse distribution of water droplets in decaying homogeneous isotropic turbulence. A polydisperse distribution of water droplets with radii ranging from R = 10 to R = 55μ m was randomly distributed in isotropic turbulent flow. The radial distribution function (RDF) was used to quantify particle clustering and simulations were performed with both non-sedimenting and sedimenting droplets.
Analyses of the RDFs for same-sized droplets reveal a power law behaviour for both non-sedimenting and sedimenting droplets in decaying turbulence. As the flow decays, droplets encounter and interact with flow scales that are evolving with time. As a result, the flow Froude number Fr and droplet Stokes numbers St change, and gravity plays a different role depending on the magnitudes of Fr and St. We see that the power law is maintained throughout the decay. For small values of Fr, clustering is less compared to the nonsedimenting case. At higher Fr, gravity enhances clustering. Empirical models have been generated for the exponential coefficients for the power law for both non-sedimenting and sedimenting droplets. We refrain from making conclusions on the dependence of particle clustering on the Reynolds number due to the small range of values of Re the flow attains during the time the particle statistics are calculated. Future studies will focus on attaining a higher maximum value as well as a wider range of Re when particle statistics are calculated.
For clustering in droplets of unequal sizes, the cross-over scale at which the power-law behaviour of the RDF tapers to a plateau is calculated from theoretical predictions for both non-sedimenting and sedimenting droplets. Excellent agreement is observed between the r c recovered from the DNS data and the theoretical formula of Chun et al. (2005) for nonsedimenting droplets. However, due to the difficulty in recovering the power-law coefficient c 1 when the value of r c approaches , comparison of r c with the empirical formula of Sundaram and Collins (1997) was not possible. Contrary to non-sedimenting droplets, we show that sedimenting droplets of unequal sizes exhibit a symmetric behaviour in RDF with droplets smaller and larger in size. The asymmetry in non-sedimenting droplets is due to the differential sedimenting velocities for particles of unequal sizes.
The effect of gravity on geometric collisions is analysed by looking at the kinematic collision kernel, the RDF and radial relative velocity at the point of collision ( r = R c ) of droplets. For same-sized droplets, gravity reduces the collision kernel for bigger droplets. This is mainly due to a reduction in relative velocity once the effects of droplet settling velocities dominate the particle inertial effects resulting in the droplet velocities being increasingly correlated. In the case of collisions between droplets of unequal sizes, gravity results in a significant increase in the kernel magnitude. A model to predict the kinematic collision kernel is developed by combining the models for the RDF and the already existing models for the conditional radial relative velocities in literature. Very good agreement is observed when compared with the DNS data. However, a smooth distribution of the kinematic collision kernel compared to the RDF and the conditional radial relative velocities suggest that any contribution due to particle clustering via the RDF is cancelled due to the definition of the conditional radial relative velocities (which includes a probability of the separation distance being at R c ) as was shown in . This causes the nonsmooth contribution to cancel out leading to the collisional kernel being determined by the first moment of the relative velocities. However confirmation of this has been reserved for future studies. We also explore the turbulence enhancement of the geometric collision kernel. Enhancement of up to factors of 3 is observed especially for monodisperse collisions.
Numerical studies involving a polydisperse distribution of droplets in decaying homogeneous isotropic turbulence are quite rare. The contour plots for the geometric collision kernel show how gravity plays an important role in enhancing collision between unequal sized droplets with radii greater than 30μ m. The results of this study can help in developing models of the collision kernel based on actual cloud conditions where gravity plays an important role. The value of typical mean turbulence dissipation rates in clouds is around 400 cm 2 s −3 . The simulations in this study consider a droplet size distribution which is comparable to cloud droplet size distributions. However, the simulations are at a much lower Reynolds number than in actual clouds. Establishing a saturation of the RDF with respect to growing Reynolds numbers would have allowed us to comment on the suitability for using these models for actual cloud conditions. However, past studies such as those by Collins and Keswani (2004), Rosa et al. (2013) have found that the RDF and the radial relative velocities reach a plateau with increasing Re . This suggests that these models can be applied to study the evolution of droplet size distribution in clouds. We also plan to compare the DNS data with data from wind tunnel experiments for similar turbulence and droplet parameters.

Ethics approval
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.