Cloud droplet diffusional growth in homogeneous isotropic turbulence: bin microphysics versus Lagrangian superdroplet simulations

. Increase of the spectral width of initially monodisperse population of cloud droplets in homogeneous isotropic turbulence is investigated applying a finite-difference fluid flow model combined with either Eulerian bin microphysics or Lagrangian particle-based scheme. The turbulence is forced applying a variant of the so-called linear forcing method that maintains the 15 mean turbulent kinetic energy (TKE) and the TKE partitioning between velocity components. The latter is important for maintaining the quasi-steady forcing of the supersaturation fluctuations that drive the increase of the spectral width. We apply a large computational domain, 64 3 m 3 , one of the domains considered in Thomas et al. (2020). The simulations apply 1 m grid length and are in the spirit of the implicit large eddy simulation (ILES), that is, with explicit 20 small-scale dissipation provided by the model numerics. This is in contrast to the scaled-up direct numerical simulation (DNS) applied in Thomas et al. (2020). Two TKE intensities and three different droplet concentrations are considered. Analytic solutions derived in Sardina et al. (2015), valid for the case when the turbulence time scale is much larger than the droplet phase relaxation time scale, are used to guide the comparison between the two microphysics simulation 25 techniques. The Lagrangian approach reproduces the scalings relatively well. Representing the spectral width increase in time is more challenging for the bin microphysics because appropriately high resolution in the bin space is needed. The bin width of 0.5 μm is only sufficient for the lowest droplet concentration, 26 cm -3 . For the highest droplet concentration, 650 cm -3 microphysics schemes represent similar departures. Finally, because the fluid flow is the same for all simulations featuring either low or high TKE, one can compare point-by-point simulation results. Such a comparison shows very close temperature and water vapor point-by-point values across the computational domain, and larger differences between simulated mean droplet radii 35 and spectral width. The latter are explained by fundamental differences in the two simulation methodologies, numerical diffusion in the Eulerian bin approach and relatively small number of Lagrangian particles that are used in the particle-based microphysics. model numerics provide the required small-scale dissipation of the turbulent kinetic energy (TKE) and scalar variance. Details of the fluid flow model are presented in the next section with the emphasis on the forcing to maintain the quasi- 120 steady turbulence, the key element of the homogeneous isotropic turbulence DNS. Two turbulence cases are considered, the low TKE case (following Lanotte et al., 2009 and Thomas et al., 2020) and the high TKE case, the latter featuring hundred times larger TKE than the former. TKE, simulations completed,

spectral broadening comes from the coupling between vertical advection in a stratified environment (that provides the supersaturation source) and advection in the bin space that represents response of the droplet population to the supersaturation forcing. Cumulus congestus 95 case from G20b is based in a modeling study by Lasher-Trapp et al. (2005) that considered a cloud observed by a radar and an instrumented aircraft during the Small Cumulus Microphysics Study (SCMS) near Cape Canaveral, Florida, during July-August of 1995. A unique aspect of the G20b study is application of the piggybacking methodology. Piggybacking refers to using two microphysics schemes in a single cloud simulation, one scheme driving the dynamics and 100 the other one piggybacking the simulated flow, see Grabowski (2019) for a review. Operating the two schemes in the same cloud-scale flow allows point by point comparison of droplet spectra predicted by the two schemes. A significantly larger mean spectral width simulated by the bin scheme across the entire cloud depth is the largest difference between the two schemes in G20b simulations. 105 In this paper, we discuss differences between the two methodologies for representing droplet spectral evolution in numerical homogeneous isotropic turbulence. Li et al. (2017) present similar comparisons applying dynamic and kinematic simulations, and including collisioncoalescence. Here, we consider diffusional growth of cloud droplets only. The direct numerical 110 simulation (DNS) methodology (e.g., Vaillancourt et al., 2001Vaillancourt et al., , 2002Lanotte et al., 2009;Li et al., 2019) and scaled-up DNS technique (Thomas et al., 2020) allow representation of turbulence impact on the droplet spectral width with an unprecedented fidelity. Sardina et al. (2015) and Grabowski and Abade (2017) provide stochastic model reference for such studies (cf. Fig. 10 in Thomas et al., 2020). In contrast to Thomas et al. (2020) who used a traditional spectral DNS 115 code, we apply a finite-difference fluid flow model that does not require small-scale dissipation to maintain computational stability. It follows that the simulations are in the spirit of the implicit large eddy simulation (ILES) where the model numerics provide the required small-scale dissipation of the turbulent kinetic energy (TKE) and scalar variance. Details of the fluid flow model are presented in the next section with the emphasis on the forcing to maintain the quasi-120 steady turbulence, the key element of the homogeneous isotropic turbulence DNS. Two turbulence cases are considered, the low TKE case (following Lanotte et al., 2009 andThomas et al., 2020) and the high TKE case, the latter featuring hundred times larger TKE than the former. https://doi.org/10.5194/acp-2020-1106 Preprint. Discussion started: 6 November 2020 c Author(s) 2020. CC BY 4.0 License.
Section 3 introduces the temperature and water vapor equations for moist simulations and presents results from simulations without droplets. Section 4 introduces numerical representation 125 of cloud droplets applying either Eulerian bin microphysics or Lagrangian superdroplets. Results of simulations with droplets are presented in section 5 focusing on the ability of either scheme to represent theoretical scalings derived by Sardina et al. (2015). Section 6 shows grid-volume by grid-volume comparison of model results facilitated by the simulation methodology, exposing additional limitations of the two microphysics simulation approaches. A brief summary in 130 section 7 concludes the paper.

The model and model forcing
The EULerian--semi-LAGrangian (EULAG) anelastic finite-difference fluid flow model (http://www.mmm.ucar.edu/eulag/) is used in this study in the ILES mode (Margolin and Rider, 2002;Andrejczuk et al., 2004;Margolin et al., 2006;Grinstein et al., 2007). ILES implies that 140 the model uses no explicit dissipation and removes small-scale velocity and scalar fluctuations through numerical diffusion provided by the monotone advection scheme. The fluid flow equations for homogeneous isotropic turbulence simulations are (e.g., Lanotte et al., 2009, Li et al., 2017: where u is the fluid flow velocity, p is pressure, = 1 kg m -3 is the air density, and f is the 150 turbulence forcing term. The forcing term ensures that the turbulence is maintained throughout the simulation with TKE flowing from large-scales towards the small-scale dissipation. The traditional technique to force the quasi-equilibrium homogeneous isotropic turbulence, convenient for spectral models, is to consider the forcing only for a few low-wavenumber modes. However, such an approach is not practical for the finite-difference model used here. Instead, we 155 apply a method in the spirit of the so-called linear forcing of Rosales and Meneveau (2005) and Onishi et al. (2011). In the homogeneous isotropic turbulence, TKE increases with the eddy size L as L 2/3 , that is, TKE is dominated by contributions from the largest eddies. Hence, one can force the turbulence by simply ensuring that the TKE does not change from one model time step to the next one because such forcing affects mostly large eddies. This implies that u (n+1) = u (n) 160 with = (Et/E (n) ) 1/2 , where n and n+1 represent time levels, E (n) is TKE at the n time level, and Et is the target TKE (see Eq. 3 in Onishi et al., 2011). The finite difference representation of such a forcing is given by 165 that is, as in the case of the linear forcing of Rosales and Meneveau (2005). The TKE dissipation rate can be derived assuming that the TKE does not change in time [see (5)  Eq. A4 is particularly useful in ILES because it allows diagnosing the TKE dissipation rate that is otherwise not known.

175
The forcing described above was initially applied in dry ILES simulations, that is, the EULAG model solving (1) and (2). For those tests (and for other simulations described in this paper), the initial flow field (scaled-up to approximately match the required TKE) was taken as the initial flow pattern in decaying turbulence simulations in Andrejczuk et al. (2004), see Fig. 1 therein.
Other parameters of the test simulations correspond to low TKE setup as described in the next 180 section. Those initial tests forced as in (3) revealed the need for additional forcing modifications as discussed below.
As in other studies of forced homogeneous isotropic turbulence (e.g., Lanotte et al., 2009), the model applies computational domain with triply-periodic lateral boundary conditions. Such boundary conditions together with (2) imply that the mean flow across the domain has to be uniform. For instance, d<uz>xy/dz (where uz is the vertical velocity and <.>xy is the horizontal average) has to vanish, and the same is true for the other two spatial directions. However, the vertically-uniform <uz>xy can evolve in time. In the initial tests, a gradual development of the mean flow across the domain was noticed. In other words, in addition to driving the turbulence 190 inside the computational domain, the forcing (3) resulted in a gradual development of the mean flow across the domain. To eliminate this undesirable behavior, an additional forcing term is included in the model equations that controls the mean flow across the domain. The additional forcing term is a simple relaxation towards the vanishing mean flow, that is, where u = (ux, uy, uz), <.> is the 3D average, and is the relaxation time scale taken as 10 model time steps. Note that (5) does not dump the flow perturbations, but only prevents the mean flow development. In simulations presented here, the mean flow across the domain after applying (5)  200 was limited to about 10 -16 m s -1 .
Second, although (3) maintains the mean TKE, the TKE partitioning between the three velocity components is allowed to evolve. As a result, the magnitude of the root mean square (rms) vertical velocity can vary in time and affect supersaturation fluctuations and droplet growth. 205 Thus, (3) needs to be modified to maintain the uniform TKE partitioning between the three components. The idea is to apply the forcing for each component separately and assuming equipartition of the TKE between all three components. For instance, for the x velocity component, ux, the modified forcing should be ux where, as before, Et is the target TKE. The forcing term for ux is then given by: 210 and similar for the two other velocity components.
To document the necessity of using the more elaborate approach (6), two simulations applying either (3) combined with (5) or (6) combined with (5) were run (details of the simulations are provided in the next section). Figure 1 compares TKE and rms vertical velocity in the two simulations. The figure shows that the modified forcing, that is, applying (6) in place of (3), maintains not only the TKE, but also a uniform in time rms vertical velocity. The latter implies 220 even partitioning of the TKE between the velocity components in agreement with the forcing formulation. The fluctuating in time rms vertical velocity in simulations applying forcing (3) leads to evolving supersaturation standard deviation in moist simulations and thus more complex mean droplet size evolution (not shown). In summary, the forcing term driving the isotropic homogeneous turbulence applied in this study is the sum of (5) and (6).

The setup of dynamic simulations
The triply-periodic computational domain is 64 3 m 3 with the model grid length of 1 m. This is one of the domains considered in Thomas et al. (2020) and close to the turbulent rising parcel extent of 50 m considered in Grabowski and Abade (2017). Such a domain size is also similar to 240 the grid volume of LES simulations of natural clouds (e.g., Siebesma et al., 2003, Stevens et al., 2005, VanZanten et al., 2011. For the fluid flow, we consider two turbulence intensities as expressed by the prescribed TKE. The "low TKE" simulations assume the TKE of 5.2 x 10 -2 m 2 s -2 . Such a TKE corresponds to the TKE dissipation rate of 10 -4 m 2 s -3 in 64 3 m 3 scaled-up DNS simulations in Thomas et al. (2020)  Simulations with droplets discussed in section 5 were run for 20 min or about 6 integral time scales starting from minute 28 of 40-min simulations presented in the next section. The "high TKE" simulations consider hundred times larger TKE, i.e., 5.2 m 2 s -2 . High TKE simulations feature ten times larger rms and maximum velocities (i.e., around 2 m s -1 and 5 to 8 m s -1 , 255 respectively) together with ten times smaller integral time scale of about 19 sec. Model time step in high TKE simulations was proportionally reduced to 0.025 sec. The high TKE simulations were run for the same number of time steps as the low TKE simulations, that is, for the total time of either 4 or 2 minutes, that is, either 12 or 6 integral time scales. Model data were saved every 15/1.5 sec for low/high TKE simulations and they are used in the analysis presented here. 260 Eq. 4 allows estimation of the TKE dissipation rate. The parameter in the forcing term is monitored from time step to time step during the simulations. The typical value in both low and high TKE is -1 ≈ 2 x 10 -4 . With the target TKE Et = 5.2 x 10 -2 m 2 s -2 and Δ = 0.25 sec in low TKE simulations, the TKE dissipation diagnosed from (4) is ≈ 4 x 10 -4 m 2 s -3 . This value 265 approximately agrees with the assumed low TKE turbulence simulation setup. For the high TKE, the target TKE and the model time step imply ≈ 4 x 10 -1 m 2 s -3 . This is a rather extreme TKE dissipation rate for small cumulus dynamics (e.g., Siebert et al., 2006), but perhaps not unusual for deep convection as simulated, for example, by Benmoshe et al. (2012).  Figure 2 shows the TKE spectra for low and high TKE at the minute 28/2.8 that are used as initial conditions for moist 275 simulations with droplets. The spectra have a classical shape characteristic of a relatively low-Reynolds-number homogeneous isotropic numerical turbulence (e.g., Fig. 2 in Rosales and Meneveau, 2005). Spectra at different times are similar to those in Fig. 2 (not shown). The turbulence dynamics in moist simulations described in the next section is exactly as 285 described above, that is, the impact of cloud droplets and of the latent heating on the flow is neglected. This is because the air density is assumed constant and the flow equations exclude the buoyancy term as typical in the homogeneous isotropic DNS simulations (see, for instance, eq. 1 in Lanotte et al., 2009). Because the flow is exactly the same in all simulations, one can compare model results grid volume by grid volume as in the piggybacking methodology (Grabowski,290 2019 and references therein). This allows a comprehensive comparison of simulation results as illustrated in section 6.

295
In addition to the momentum, the moist ILES of homogeneous isotropic turbulence solves the temperature T and water vapor mixing ratio qv equations in the form: where Lv = 2.5 x 10 6 J kg -1 is the latent heat of condensation, cp = 1015 J kg -1 K -1 is the specific heat of air at constant pressure, g = 9.81 m s -2 is the gravitational acceleration, and Cd is the condensation rate, the rate of change of the cloud water mixing ratio resulting from the 305 diffusional growth of cloud droplets. Calculation of the condensation rate depends on the microphysics scheme as explained below.
In the spirit of DNS studies of homogeneous isotropic turbulence, the initial temperature and water vapor mixing ratio in moist simulations are assumed spatially uniform. The actual values 310 are taken as in Thomas et al. (2020), that is, T = 283 K and qv at saturation assuming environmental pressure of 1000 hPa. The spatially-uniform initial conditions justify the triplyperiodic computational domain. The last term in (6) represents the temperature change due to adiabatic air expansion resulting from the vertical motion in the stratified environment. This term drives small-scale supersaturation fluctuations in the otherwise uniform environment; see 315 discussion in section 3 in Vaillancourt et al. (2001). Such a modeling framework is a simplification of a truely stratified environment where the environmental temperature and pressure are functions of height and typically the potential temperature (an invariant for dry adiabatic vertical displacements) rather than the temperature is applied as the model variable. In some DNS studies, the temperature and moisture equations are combined into the supersaturation 320 equation that includes the source due to the vertical motion (as the last term in Eq. 6) and the sink due to droplet growth, Eq. (2) in Lanotte et al. (2009)

330
The initial test of the moist ILES framework considers a 40-min long low TKE and 4-min long high TKE simulations without droplets (i.e., both up to about 12 integral time scales) and initiated in the same way as the dry simulation illustrated in Figs. 1 and 2. The simulations apply the fluid flow as described above and solving (7) and (8) without the condensation term Cd. In such simulations, the largest temperature change is possible when the air parcel rises across the entire computational domain depth, that is, 64 m, with the corresponding temperature change of about 0.64 K as given by (7). Such a maximum temperature change leads to the supersaturation change from the initial zero to about 4%. In the numerical simulation, the maximum temperature deviations from the uniform initial 283 K are typically smaller than 0.5 K. The evolutions of 340 supersaturation statistics are shown in Fig. 3. Despite the dramatic difference in the TKE levels, the statistics are similar regardless of the TKE level, in agreement with the parcel argument.
Only after including the source due to condensation, the evolutions become different depending on the droplet characteristics and the TKE level. Small differences in the evolutions in Fig. 3 come from different flow realizations between low and high TKE cases. 345

ILES moist simulations with droplets
The general microphysical setup considers initially monodisperse population of cloud droplets with the radius of 13 µm present in three different concentrations: 26, 130, and 650 cm -3 . The 350 concentration 130 cm -3 was considered in Lanotte et al. (2009) and Thomas et al. (2020) and it corresponds to the mean cloud water content of about 1.2 g m -3 . In addition, the five-times smaller and five-times larger droplet concentrations are considered to document how the two microphysics schemes represent the expected scalings. The mean cloud water content in simulations with the increased (decreased) droplet concentration is about 6.0 (0.24) g m -3 . Moist 355 simulations with droplets start at minute 28 of the moist simulation without droplets, that is, applying the flow field analyzed in Fig. 2, together with the temperature and moisture field at minute 28 in Fig. 3. The simulations are run for additional 20 minutes for the low TKE setup, saving data every 15 sec. The high TKE simulations are run for 2 minutes and the data are saved every 1.5 sec. The data are used in the analysis of both macrophysical (e.g., the supersaturation 360 characteristics) and microphysical (e.g., spectral width) simulated by the two microphysics schemes. Sardina et al. (2015) derived scalings for the case when the droplet phase relaxation time is much shorter than the turbulence integral time scale. This is the case for the selected domain size and the low TKE for all droplet concentrations. For the high TKE and low droplet concentration (i.e., 26 cm -3 ), the phase relaxation time (about 10 seconds) is the closest to the turbulence integral time scale (about 19 seconds), so some deviations from the theoretical scaling should be expected as shown below.

370
The particle-based Lagrangian scheme considers on average 40 Lagrangian particles (superdroplets) per grid volume, each featuring the initial radius of 13 µm. The number of superdroplets per grid volume together with the assumed droplet concentration dictates the multiplicity that is assumed the same for all superdroplets. Although the average number of superdroplets per grid volume is small when compared to millions real droplets within a 1 cubic 375 meter grid volume, G20a and G20b document that the number as small as 10 per grid box provides physically-meaningful results; see also Li et al. (2017). At the simulation onset, each superdroplet is placed at a random position within a grid volume. To be consistent with the bin microphysics, superdroplets grow in response to the mean supersaturation predicted inside a grid volume it occupies. Superdroplets are advected applying a model flow field interpolated to the 380 droplet position as in Arabas et al. (2015). The interpolation scheme maintains the incompressibility of the flow at subgrid scales; see discussion of this aspect in section 2.4 in Grabowski et al. (2018). Droplet inertia and droplet sedimentation are not considered.
Condensation rate Cd is calculated by summing up the mass change of all superdroplets present within a given grid volume. The Eulerian bin microphysics considers spectral density function represented by 40 equallyspaced bins with the bin size modified in different simulations as described below. The reason 395 for modifications of the bin resolution is to match the results from the Lagrangian microphysics as shown in the results section. In each bin setup, there is a bin centered at 13 µm that is filled with droplets at the simulation onset. It should be stressed that the monodisperse initial droplet size distribution is impossible to be accurately represented using the spectral density function because the monodisperse distribution corresponds to the delta function. However, even with a 400 finite width of the initial distribution, the broadening of the distribution as time progresses (Sardina et al., 2015;Li et al., 2019;Thomas et al., 2020) can be appropriately represented as illustrated below provided that the bin width is appropriately small. In the bin microphysics, each bin is independently advected in the physical space using the same advection scheme that is applied to the momentum, temperature, and water vapor mixing ratio. Neither droplet 405 sedimentation nor droplet inertia are considered as in the Lagrangian scheme. All bins are combined at each grid volume to calculate evolution of the droplet spectrum due to the local sub-  Table 1. However, for the high TKE, the scaling breaks down 450 because the phase relaxation time is no longer much smaller than the turbulence integral time scale. For a given droplet concentration, shift from low to high TKE should result to ten-fold increase of σS because of the <w 2 > 1/2 increase. This is approximately valid for 650 cm -3 droplet concentration, but reduces to a factor of only about 5 for 26 cm -3 concentration.
The supersaturation standard deviation shown in Fig. 4 can be compared to the standard deviation resulting from the quasi-equilibrium supersaturation fluctuations. Evolution of the supersaturation S = qv/qvs -1 (where qvs is the saturated water vapor mixing ratio) can be derived by combining Eqs. 7, 8 and 9 as where relax is the phase relaxation time that depends of the mean droplet radius and concentration: 465 where w = 10 3 kg m -3 is the water density, Rv = 461 J kg -1 K -1 is the water vapor gas constant, and <.> in (11) depicts averaging over all droplets within a given grid volume. The quasiequilibrium supersaturation is obtained by setting the left-hand-side of (10) to zero that leads to 470 Seq = a1 w relax . For the mean temperature and humidity of the simulations and specific numerical values of the relevant constants, a1 = 6.54 x 10 -4 m -4 and the phase relaxation time for the 13 µm droplets and their concentration of 130 cm -3 is relax = 1.98 sec. The phase relaxation time is five times smaller/larger for the droplet concentration of 650/26 cm -3 . Since the quasiequilibrium supersaturation is proportional to the vertical velocity, its standard deviation is 475 proportional to the rms vertical velocity. For the low TKE, the quasi-equilibrium supersaturation standard deviation is 0.120, 2.39e-2, 4.79e-3% for 26, 130 and 650 cm -3 droplet concentrations.
These are in a relatively good agreement with simulated standard deviation shown in Table 1.    Figure 6 shows evolutions of the radius squared standard deviation. Sardina et al. (2015) show that the standard deviation of the radius squared distribution should increase in time as square 520 root of time as long as the phase relaxation time is much smaller than the turbulence integral time. The rate of increase is proportional to the supersaturation standard deviation (see Eq. 13 therein). The former has been shown in other numerical simulations, such as in Li et al. (2019) and Thomas et al. (2020). Fig. 6 shows that the Lagrangian microphysics reproduces the t 1/2 scaling and that the differences between various simulations for the low TKE can be explained 525 by the differences in the supersaturation standard deviation shown in Fig. 4 (note that on the loglog plot the rate of increase change corresponds to a vertical shift as shown in Fig. 6). For the https://doi.org/10.5194/acp-2020-1106 Preprint. Discussion started: 6 November 2020 c Author(s) 2020. CC BY 4.0 License.
high TKE, small deviations for the expected scaling can be explained by the phase relaxation time being no longer much smaller than the turbulence integral time. This is especially evident for the high-TKE SDS.26 case in the right panel of Fig. 6. Left panels of Figs. 4 and 6 show that 530 increasing the number of superdroplets from 40 to 150 per grid volume has virtually no impact on the results. Overall, the Lagrangian microphysics seems to represent expected scalings (or departures from them) without much difficulty.   vertical velocity (i.e., Fig. 5) for the bin microphysics is similar to that shown for the Lagrangian scheme and is not shown. In summary, Eulerian bin microphysics is capable in appropriately representing turbulent temperature and moisture fluctuations, but fails to simulate their impact on droplet spectra unless 580 appropriately high bin resolution is used. But even with sufficient resolution, the spatial variability in the droplet spectra differs between Lagrangian and Eulerian microphysics as documented in the next section.

Grid volume by grid volume analysis of marco-and microphysical properties 585
The analysis presented in the previous section concerned domain-averaged characteristics.
Because all simulations with either low or high TKE feature exactly the same evolving flow field, simulated thermodynamic variables (i.e., the temperature, water vapor, and cloud droplet characteristics) can be compared grid volume by grid volume and thus provide a comprehensive 590 comparison of simulated local conditions. Such a comparison is in the spirit of the piggybacking methodology applied in G20b. σ(R 2 ) (!m 2 ) https://doi.org/10.5194/acp-2020-1106 Preprint. Discussion started: 6 November 2020 c Author(s) 2020. CC BY 4.0 License. Figure 9 compares the temperature, water vapor, and cloud water mixing ratios (the latter derived from the predicted droplet spectra within each grid volume) between SDS.130 and BIN.130 low 595 TKE simulations at time of 20 minutes. For the temperature and moisture, plots at earlier times are similar to upper panels in Fig. 9 except for smaller ranges between minima and maxima.
Cloud water plots at earlier times are also similar to those shown in Fig. 9 except for the initial couple minutes. The temperature and water vapor values are extremely close between the two simulations: the root mean square difference between temperatures are 2.6 x 10 -4 K. For the 600 water vapor mixing ratios, the root mean square difference is 1.1 x 10 -4 g kg -1 . However, the cloud water mixing ratio can differ significantly between the bin and the superdroplet simulations. This is because of statistical fluctuations in the number of superdroplets per grid volume. With on average 40 superdroplets per grid volume, the standard deviation of the droplet number is around 6, or about 15% of the mean. Assuming that one can find grid volumes with 605 three times the standard deviation, the range of the cloud water mixing ratio can be as high as close to 50%. This can explain the spread seen in the lower left panel of   microphysics is the one that shows the correct scaling in Fig. 8. The mean radius comparison features some scatter that is reduced when increased resolutions simulations are compared.
However, the scatter is asymmetric with respect to 1:1 line and similar to the cloud water scatter 640 in the lower right panel of Fig. 9. The asymmetry shows that the bin microphysics tends to simulate larger droplets than the Lagrangian microphysics for droplets smaller than the mean, and the reverse is true for droplets larger than the mean. The spectral width panels show that the bin simulations feature smaller spread of the spectral width across the computational domain than the Lagrangian scheme. This seems independent of the bin resolution. In other words, spectral width simulated by the bin scheme varies less across the computational domain. In  In summary, the grid volume by grid volume comparison between Eulerian and Lagrangian results shows that the simulated spatial variability is smaller in the bin microphysics when compared to the superdroplets. Arguably, this comes from a combination of the numerical diffusion in the bin microphysics (i.e., smoothing bin results similar to other Eulerian fields) and small-scale fluctuations of the Lagrangian microphysics due to a relatively small mean number of superdroplets per grid volume.

685
Because we apply a finite-difference fluid flow model, we had to develop a turbulence forcing scheme that led to the quasi-steady homogeneous isotropic turbulence similar to that simulated by a spectral model. The forcing scheme applied here is in the spirit of the linear forcing of Rosales and Meneveau (2005) and Onishi et al. (2011). The idea is to ensure that the mean turbulence kinetic energy (TKE) remains constant in time and that the partitioning between the three TKE components does not change. The latter is important for the simulations of the turbulence impact on the droplet spectra because the driving mechanism in the homogeneous environment comes from the vertical velocity fluctuations affecting the local supersaturation.
The simulations are in the spirit of the implicit large eddy simulation (ILES), that is, without modeling the small-scale TKE and scalar variance dissipation (Margolin and Rider, 2002;695 Andrejczuk et al., 2004;Margolin et al., 2006;Grinstein et al., 2007). ILES applying a finite difference model relies on the model numerics to provide the required dissipation, in contrast to the scaled-up spectral model DNS in Thomas et al. (2020) where appropriately scaled-up molecular dissipation coefficients have to be used to ensure stable simulations. The TKE dissipation rate diagnosed from the forcing algorithm (see Eq. 6) is approximately correct for the 700 considered turbulence intensities.
There are two drastically different simulation techniques that can be applied to investigate the impact of cloud turbulence on the droplet spectra: the Eulerian bin microphysics and the Lagrangian particle-based approach, the latter often referred to as the superdroplet method (see 705 review in the introduction). We apply the ILES homogeneous isotropic turbulence setup to compare the two techniques following Li et al. (2017), , and Grabowski (2021) for the diffusional growth of cloud droplets only. The computational domain is 64 3 m 3 , one of the domain sizes considered in Thomas et al. (2020) and similar to grid volumes of a typical LES simulation of natural clouds. Two TKE intensities are considered, low and high, different by a 710 factor of one hundred. The latter implies that the velocity fluctuation differ by a factor of ten and the TKE dissipation rates differ by a factor of one thousand. The TKE dissipation spans the range observed in natural clouds. The initially monodisperse droplet size distribution poses a challenge for the Eulerian bin microphysics that needs to be kept in mind when comparing the two simulations methodologies. 715 The Lagrangian approach reproduces the expected scalings derived in Sardina et al. (2015) for the case when the turbulence integral time scale is much longer that the phase relaxation time of cloud droplets. Representing the scalings is more challenging for the bin microphysics, not only because of the ill-posedness of the initial conditions, but also because appropriately high 720 resolution in the bin space is needed. In fact, the standard bin resolution, with the bin width of 0.5 μm and covering the range up to 20 μm, similar to Grabowski ( , 2021, is only sufficient for the lowest droplet concentration (26 cm -3 ). For the highest droplet concentration, 650 cm -3 , even an order of magnitude smaller bin size is not sufficient to reproduce the expected scaling.
Such a bin resolution is impossible to use when collisional growth is also considered as in Li et 725 al. (2017). For the lowest droplet concentration (26 cm -3 ) and the high TKE case, the phase relaxation time is about 10 sec and the turbulence integral time is around 19 sec, so some departures from the expected scaling are expected. This is indeed the case, and the two simulation methodologies represent similar supersaturation and spectral width departures.

730
Because the fluid flow is the same for all simulations featuring either low or high TKE, one can compare model results point-by-point as in the piggybacking technique of Grabowski (2019).
Such a comparison shows miniscule differences between temperature and water vapor fields across the computational domain, and larger differences between simulated mean droplet radii and spectral width. These are consistent with fundamental differences in the two simulation 735 methodologies, numerical diffusion in the Eulerian approach and relatively small number of Lagrangian particles (superdroplets) that can be afforded in the particle-based microphysics.
Either one can be limited by either increasing the model resolution or increasing the number of Lagrangian particles, both significantly increasing computational cost. But there are additional options for the particle-based microphysics, for instance, assuming that a particle within a given 740 grid volume represents a cloud of particles spread over a prescribed halo. We are pursuing those ideas in ongoing research.
Author contributions. WWG developed the idea and ran the simulations. WWG and LT performed data analysis and prepared the manuscript.