Diffusional growth of cloud droplets in homogeneous isotropic turbulence: DNS, scaled-up DNS, and stochastic model

This paper presents a novel methodology to use direct numerical simulation (DNS) to study the impact of isotropic homogeneous turbulence on the condensational growth of cloud droplets. As shown by previous DNS studies, the impact of turbulence increases with the computational domain size, that is, with the Reynolds number, because larger eddies generate higher and longer-lasting supersaturation fluctuations that affect growth of individual cloud droplets. The traditional DNS can only simulate a limited range of scales because of the excessive computational cost that comes from resolving all scales involved, that is, from large scales at which the turbulent kinetic energy (TKE) is introduced down to the Kolmogorov microscale, and from following every single droplet. The novel approach is referred to as the “scaled-up DNS”. The scaling up is done in two parts, first by increasing both the computational domain and the Kolmogorov microscale and second by using super-droplets instead of real droplets. To ensure proper dissipation of TKE and scalar variance at small scales, molecular transport coefficients are appropriately scaled up with the grid length. For the scaled-up domains, say, meters and tens of meters, one needs to follow billions of real droplets. This is not computationally feasible, and so-called super-droplets are applied in scaled-up DNS simulations. Each super-droplet represents an ensemble of identical real droplets, and the number of real droplets represented by a super-droplet is referred to as the multiplicity attribute. After simple tests showing the validity of the methodology, scaled-up DNS simulations are conducted for five domains, the largest of 64³ m³ volume using a DNS of 256³ grid points and various multiplicities. All simulations are carried out with vanishing mean vertical velocity and with no mean supersaturation, similarly to past DNS studies. As expected, the supersaturation fluctuations as well as the spread in droplet size distribution increase with the domain size, with the droplet radius variance increasing in time t as t^1∕2 as identified in previous DNS studies. Scaled-up simulations with different multiplicities document numerical convergence of the scaled-up solutions. Finally, we compare the scaled-up DNS results with a simple stochastic model that calculates supersaturation fluctuations based on the vertical velocity fluctuations updated using the Langevin equation. Overall, the results document similar scaling to previous small-domain DNS simulations and support the notion that the stochastic subgrid-scale model is a valuable tool for the multi-scale simulation of droplet spectral evolution applying a large-eddy simulation model.

1 Introduction

The impact of turbulence on the growth of cloud droplets is an important and still poorly understood aspect of cloud physics. This is because of the wide range of spatial scales that affect droplet growth, from the Kolmogorov microscale (about a millimeter for typical atmospheric turbulence levels) to the scale of the entire cloud or cloud system. Cloud droplets grow by the diffusion of water vapor and by gravitational collision/coalescence, with the former dominating growth until droplets are large enough so that the collisional growth can be initiated and eventually lead to drizzle and rain formation. For the gravitational collision/coalescence, the frequency of droplet collisions depends on the droplet spectrum width. It follows that understanding processes leading to the observed droplet spectra is important for the understanding of the rain onset. Observations of natural droplet spectra go back to the early days of aircraft cloud studies (e.g., Warner, 1969) and continue in numerous subsequent investigations (e.g., Jensen et al., 1985; Brenguier and Chaumat, 2001; Pawlowska et al., 2006; Prabha et al., 2012, among many others; see also references in Grabowski and Wang, 2013). Those observations typically show that observed droplet spectra are wider than predicted by simple models of cloud dynamics and microphysics. In many instances, such a discrepancy can be explained by cloud entrainment (e.g., Warner, 1973; Paluch and Knight, 1984; Su et al., 1998; Lasher-Trapp et al., 2005, among many others). However, the presence of a significant spectral broadening in undiluted and weakly diluted cloudy volumes is more difficult to explain. One can wonder whether the presence of small-scale turbulence can lead to appreciable widening of the droplet spectra during diffusional growth within otherwise uniform cloudy volumes.

Vaillancourt et al. (2001, 2002) were the first to apply the direct numerical simulation (DNS) approach to study diffusional growth of cloud droplets in homogeneous isotropic turbulence by applying a rising adiabatic parcel setup. DNS was initially applied to turbulent particle-laden flows to study the so-called preferential concentration (or clustering) of inertial particles in turbulence (e.g., Eaton and Fessler, 1994, see also references in Shaw, 2003). For the DNS of cloud droplets growing by the diffusion of water vapor, droplets respond to the supersaturation fluctuations in their immediate environment as suggested by Srivastava (1989). Limited by the computational resources, Vaillancourt et al. (2002) were only able to consider small volumes of a turbulent cloud, around 1 L. Three sets of simulations were performed with turbulent kinetic energy dissipation rates relevant to cloud conditions and including droplet sedimentation. Earlier studies of particle-laden turbulent flows typically exclude sedimentation (see Eaton and Fessler, 1994), but this is inappropriate for weak to moderate turbulence intensities typical for natural clouds (Grabowski and Vaillancourt, 1999). Vaillancourt et al. (2002) simulations show a small impact of local supersaturation fluctuations on the droplet spectra: the standard deviation of the initially monodisperse droplet distribution increases very slowly with time, on the order of 0.01 µm min⁻¹. Similar simulations reported in Lanotte et al. (2009) applying larger domains and no mean ascent clearly show that the impact, although still relatively small (a few tenths of 1 µm), does increase with the domain size (see Fig. 3 therein). In similar DNS simulations, Li et al. (2019) demonstrate the increase in spectral broadening with the increase in the domain size (i.e., the Reynolds number) and the increase in the length of the simulations; see Figs. 3 and 4 therein. For the largest domain of 512³ and about a minute of the simulation time, the initially monodisperse 10 µm droplets evolve into a spectrum with about 1 µm width.

Homogeneous isotropic turbulence simulations of Vaillancourt et al. (2002), Lanotte et al. (2009), and Li et al. (2019) are limited by the computational domain size. As a result, simulations featuring domains larger than a fraction of a cubic meter are simply not yet possible. At the same time, as argued in Grabowski and Wang (2013) and documented in Grabowski and Abade (2017; see Fig. 4 therein) and Li et al. (2019; see Fig. 4 therein), the impact of supersaturation fluctuations in homogeneous isotropic turbulence on the spectral width increases with the domain size. A simple argument is that this is because the largest turbulent eddies feature the largest vertical velocity perturbations that result in the largest and longest-lasting supersaturation fluctuations and thus have the largest impact on the spread of droplet growth histories.

To this end, we propose to use what we refer to as the “scaled-up DNS” approach. Since the largest eddies are the key for the condensational growth, one would like to apply the DNS technique in simulations with domains much larger than currently possible. For instance, taking a 128³ DNS simulation with 10 cm grid length gives a computational domain of 12.8³ m³, that is, comparable to the grid volume of an LES cloud simulation. To ensure a proper dissipation of the turbulent kinetic energy (TKE) at the smallest scales, one needs to scale up the molecular viscosity with the increase in the model grid length. The increase in the small-scale dissipation is critical for traditional DNS models applying spectral techniques to simulate homogeneous isotropic turbulence as applied in this study. This is different from past turbulence-related studies applying finite-difference models with large domains and spatially uniform diffusion coefficients (e.g., Grabowski and Clark, 1993; Mellado et al., 2018; Rotunno and Bryan, 2018) or no explicit diffusion at all as in the so-called implicit large-eddy simulation (ILES; e.g., Margolin et al., 2006; Grinstein et al., 2007). The increased number of droplets in the large domain can be accounted for by the so-called super-droplet technique in which each super-droplet represents an appropriately scaled-up number of natural droplets (referred to as the multiplicity factor, Shima et al., 2009) as already applied in the appendix of Lanotte et al. (2009) and in Li et al. (2019).

The paper is organized as follows. The next section presents the model and modeling setup. Section 3 presents a general methodology of the scaled-up DNS and discusses numerical tests of this approach. Cloud droplets are added to scaled-up DNS simulations in Sect. 4 applying the super-droplet method. Section 5 compares DNS and scaled-up DNS supersaturation fluctuations with those obtained from a simple stochastic model. Concluding discussion is the focus of Sect. 6.

2 The model and modeling setup

The numerical code used here is that of Kumar et al. (2012, 2014). It solves evolution equations for the three velocity components (u, v, w), the temperature T, and the water vapor mixing ratio q_v. Cloud droplets are represented as point particles followed in space and they grow or evaporate as dictated by their local environment. Droplet collisions are not considered. The Navier–Stokes equations are solved by a pseudo-spectral method over a cubic volume with periodic lateral boundary conditions in three directions using the fast Fourier transforms. Time stepping is performed using a second-order predictor–corrector method. The code is parallelized in two dimensions and the cubic domain is decomposed into so-called pencils. The same procedures as in Kumar et al. (2012) are followed for the initial turbulent state preparation and the turbulence maintenance. See Kumar et al. (2012, 2014) for more details.

Two modifications have been made to the code to carry out the present study. First, we included an additional source/sink term in the temperature equation that was missing in the original code. The term describes evolution of temperature fluctuations affected by the vertical velocity. This effect is incorporated in the DNS through the source/sink term $-gw/C_{\mathrm{p}}$, where g is the gravitational acceleration, w is the local vertical velocity, and C_p is the specific heat capacity of air at constant pressure. The complete equation for the evolution of temperature fluctuations is

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial T{}^{\prime }}{\partial t}}+\mathbf{u}\cdot \mathrm{\nabla }T{}^{\prime }=K{\mathrm{\nabla }}^{\mathrm{2}}T{}^{\prime }+{\displaystyle \frac{L_{\mathrm{v}}}{C_{\mathrm{p}}}}{C}_{\mathrm{d}}-{\displaystyle \frac{g}{C_{\mathrm{p}}}}w{}^{\prime },\end{array}$$

where K is the thermal diffusion coefficient, L_v is the latent heat of vaporization, and C_d is the condensation rate. Without the last term, the vertical velocity simulated by the DNS has no impact on the supersaturation fluctuations. Since the emphasis in Kumar et al. (2014, 2018) was on the mixing between cloudy and clear air, this omission has a negligible impact on results presented there. However, this term is critical for the current study.

Second, we modified the way condensation rate is calculated for a single droplet. The analytic formulation applied originally has the form

$$\begin{array}{}\text{(2)}& C_{\mathrm{d}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }{K}_{\mathrm{r}}{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\Delta }V}}SR,\end{array}$$

where S is the supersaturation, R is the droplet radius, ρ₀ is the air density, ρ_w=10³ kg m⁻³ is the liquid water density, $K_{\mathrm{r}}=\mathrm{5.00}\times {\mathrm{10}}^{-\mathrm{11}}$ m² s⁻¹ is the condensational growth constant (i.e., $\mathrm{d}R/\mathrm{d}t=K_{\mathrm{r}}S/R$), and $\mathrm{\Delta }V=\mathrm{d}x\times \mathrm{d}y\times \mathrm{d}z$ is the grid box volume. To ensure mass conservation, Eq. (2) is modified to

$$\begin{array}{}\text{(3)}& C_{\mathrm{d}}\mathrm{\Delta }t={\displaystyle \frac{\mathrm{4}\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{3}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\Delta }V}}\left[R^{\mathrm{3}}\left(t\right)-R^{\mathrm{3}}\left(t-\mathrm{\Delta }t\right)\right],\end{array}$$

where Δt is the time step and R(t) and R(t−Δt) are droplet radii at time t and t−Δt, respectively. This means that droplet growth is calculated first, and then Eq. (3) is used to derive the condensation rate.

The coupling of the Eulerian fields and the droplets is done using trilinear interpolation. The condensation rate is calculated for each droplet by interpolating the values of T and q_v from the grid points enclosing the droplet. The condensation rate is calculated at the droplet position and then redistributed to the nearest eight grid points through a reverse procedure. The condensation rate provides a feedback on the temperature and water vapor evolutions. Inertial effects and gravitational settling are included in the droplet motion. More details can be found in Kumar et al. (2012).

The modeling setup follows one of the simulations discussed in Lanotte et al. (2009). We consider an initial mono-disperse droplet distribution of 13 µm radius and a concentration of 130 cm⁻³. The liquid water content (LWC) is 1.19 gm⁻³. The initial conditions include a uniform temperature of 283 K and zero supersaturation. The latter gives a water vapor mixing ratio of 7.65 g kg⁻¹. Since the mean velocity inside the DNS domain is zero, the total cloud water does not change with time, but the initial monodisperse droplet size distribution broadens because the supersaturation fluctuates in time and space, affecting the distribution (cf. Li et al., 2019; Saito et al., 2019). The two specific aspects are discussed in the next two sections that allow one to extend the DNS methodology to apply large spatial domains.

3 Scaled-up DNS simulation

The intensity of turbulence is typically expressed by the turbulent kinetic energy (TKE) dissipation rate ε. Increasing the domain size L for the same ε increases kinetic energy of turbulent motions. The TKE determines velocity fluctuations and controls the supersaturation variations that play the key role in the condensational growth of cloud droplets. The TKE dependence on ε and L is typically expressed as (e.g., Pope, 2000)

$$\begin{array}{}\text{(4)}& E\sim (L\mathit{\epsilon }{)}^{\mathrm{2}/\mathrm{3}}.\end{array}$$

In the classical DNS, the grid length has to be close to the Kolmogorov microscale η to allow proper TKE dissipation at the smallest scales. Increasing the domain size L without changing the number of grid points implies that the grid length increases as well. We will refer to such simulations as “scaled-up DNS”. With the increased grid length, one needs to increase molecular transport coefficients to maintain proper TKE dissipation as well as the removal of scalar fluctuations. Note that this is different from LES and ILES methodologies mentioned in the introduction. Assuming that the domain size L represents the appropriate scale of energy-containing eddies, the L∕η ratio represents the flow Reynolds number Re:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathit{\eta }}{L}}\sim {\mathit{Re}}^{-\mathrm{3}/\mathrm{4}}\end{array}$$

(e.g., Pope, 2000; Grabowski and Clark, 1993). Keeping the Reynolds number the same for the actual and scaled-up DNS implies that

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{L_{\mathrm{1}}}{{\mathit{\eta }}_{\mathrm{1}}}}={\displaystyle \frac{L_{\mathrm{2}}}{{\mathit{\eta }}_{\mathrm{2}}}},\end{array}$$

where subscripts 1 and 2 refer to different domain sizes and corresponding Kolmogorov microscales. The Kolmogorov microscale is given by $\mathit{\eta }={\left(\frac{{\mathit{\nu }}^{\mathrm{3}}}{\mathit{\epsilon }}\right)}^{\mathrm{1}/\mathrm{4}}$, where ν is the viscosity. Applying η with the same TKE dissipation rate ε for both ν₁ and ν₂ leads to

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{L_{\mathrm{1}}}{{\mathit{\nu }}_{\mathrm{1}}^{\mathrm{3}/\mathrm{4}}}}={\displaystyle \frac{L_{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{2}}^{\mathrm{3}/\mathrm{4}}}}.\end{array}$$

Equation (7) implies that when the TKE dissipation rate is supposed to be the same in DNS and scaled-up DNS, the molecular viscosity needs to be scaled up as

$$\begin{array}{}\text{(8)}& {\mathit{\nu }}_{\mathrm{2}}={\mathit{\nu }}_{\mathrm{1}}{\left({\displaystyle \frac{L_{\mathrm{2}}}{L_{\mathrm{1}}}}\right)}^{\mathrm{4}/\mathrm{3}},\end{array}$$

where ν₁ is the viscosity in the real DNS, ν₂ is the scaled-up viscosity, L₁ is the real-DNS domain size, and L₂ is the scaled-up domain size. Thermal and water vapor diffusivities are also changed in the same way as the viscosity.

We used a DNS with 256³ grid points to study scaling-up simulations without droplets. Real DNS was run for L=0.256 m and scaled-up DNS was run with domains of sizes L=2.56 m, L=25.6 m, and L=256 m. According to Eq. (8), the viscosity (taken as ν=0.15 cm² s⁻¹ for the real DNS) has to be scaled up by 21.54, 464.16, and 10 000 times for L=2.56 m, L=25.6 m, and L=256 m, respectively. All simulations are forced as described in Kumar et al. (2012) by applying a TKE dissipation rate ε of 10 cm² s⁻³ as in Lanotte et al. (2009).

Figure 1 shows energy spectra for the real 25.6 cm DNS and the three scaled-up DNS. The black dashed lines represent the $-\mathrm{5}/\mathrm{3}$ slope expected in the inertial range. The spectral peak shifts to the left and its value increases as the domain size increases (i.e., the wavenumber k decreases). The slope remains approximately similar for the four simulations.

Figure 1 Comparison of energy spectra for real and scaled-up DNS.

TKE in the scaled-up simulations for the same TKE dissipation rate ε should increase following the scaling originating from Eq. (4), that is,

$$\begin{array}{}\text{(9)}& E_{\mathrm{1}}=E_{\mathrm{2}}{\left(L_{\mathrm{1}}/L_{\mathrm{2}}\right)}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where E₁ and L₁ are for the scaled-up DNS and E₂ and L₂ are for real DNS (e.g., L₂=0.256 m and E₂=20 cm² s⁻²). Figure 2 shows the evolution of TKE and TKE dissipation rate for the four simulations in Fig. 1. For the TKE evolution, dashed lines show the expected scaling based on Eq. (9). TKEs from the scaled-up DNS simulations agree with the theoretical scaled-up TKE values. To show that the DNS and scaled-up DNS feature the same TKE dissipation rate, we also show the dissipation rate calculated as twice the product of the scaled-up viscosity and the mean simulated enstrophy as typically done in DNS studies. The plots show that the forcing is approximately correct in the scaled-up simulations. The scaled-up simulations need to be run for longer times, with the timescale following the $L/E^{\mathrm{1}/\mathrm{2}}$ scaling of the large-eddy turnover time. The simulations show that the scaled-up DNS with viscosity modified according to Eq. (8) produces the expected TKE.

Figure 2 Evolution of TKE (upper panels) and TKE dissipation rate (lower panels) for four simulations mentioned in the text. The dashed lines are theoretical values. Bottom panels show nondimensional time using eddy turnover time; see Eq. (13) in Sect. 5.

The simulations shown in Figs. 1 and 2 feature the same dynamic range, that is, the same Reynolds number and the L∕η ratio. However, one may also consider scaled-up DNS simulations where the dynamic range is changed. For instance, one may compare simulations with the same ε and L and different numbers of grid points N covering L. For such simulations, the change in the Kolmogorov microscale $\mathit{\eta }=L/N$ suggests the required rescaling of the dissipation coefficients. Since $\mathit{\eta }={\left(\frac{{\mathit{\nu }}^{\mathrm{3}}}{\mathit{\epsilon }}\right)}^{\mathrm{1}/\mathrm{4}}$, assuming ε=const gives the scaling similar to Eq. (8):

$$\begin{array}{}\text{(10)}& {\mathit{\nu }}_{\mathrm{2}}={\mathit{\nu }}_{\mathrm{1}}{\left({\displaystyle \frac{N_{\mathrm{1}}}{N_{\mathrm{2}}}}\right)}^{\mathrm{4}/\mathrm{3}},\end{array}$$

that is, with the number of grid points rather than the domain size providing the scaling.

The scaling Eq. (10) is illustrated in Fig. 3 that shows the spectra in simulations with the domain size of either 0.512 or 1.024 m and applying either DNS or scaled-up DNS. The spectra are obtained at final simulation times. The red lines represent spectra for the real DNS, and green and blue lines show spectra for scaled-up DNS. Scaling up accurately predicts the energy at the largest scales, but some energy at smaller scales, still far from the dissipation, is lost. This means that the total TKE for a scaled-up DNS is slightly lower than the real DNS within the same volume. For the simulations shown in Fig. 3, TKE for L of 0.512 m is 34.2, 32.0, and 26.0 cm² s⁻² for real DNS (N=512) and scaled-up DNS with N=256 and 128, respectively. For L of 1.024 m, TKE is 55.0, 50.0, and 41.0 cm² s⁻² for real DNS (N=1024) and scaled-up DNS with N=256 and 128, respectively. Because for the condensational growth the interest is in the largest scales as discussed in the introduction, the energy loss at smaller scales can be considered less important. However, this aspect is relevant for the comparison between scaled-up DNS and the stochastic model as discussed in Sect. 5.

Figure 3 Energy spectrum comparison of real DNS and scaled-up DNS. Panels (a, b) correspond to 512³∕1024³ real DNS.

4 Applying super-droplets for the scaled-up DNS

For a scaled-up DNS, one needs to follow a significantly larger number of droplets when compared to DNS. For instance, for the droplet concentration of 130 cm⁻³ one needs to follow 1.3×10¹¹ droplets for a domain of L=10 m. This is not computationally feasible. To overcome this problem, one can use the so-called super-droplets (Shima et al., 2009) instead of real droplets, where each super-droplet represents an ensemble of real droplets with the same radius. Position and velocity of each super-droplet are predicted in the same way as for the real droplet. The number of real droplets represented by a single super-droplet is referred to as the multiplicity attribute μ (Shima et al., 2009).

At the onset of simulations, super-droplets are inserted into the computational domain in the same way as regular droplets; that is, they are randomly positioned inside the domain and subsequently followed in space and time as regular droplets. The condensation rate for a super-droplet is calculated as in Eq. (3) except for an additional multiplicity factor μ. The evolution of the temperature and water vapor mixing ratio fluctuations is affected by the condensation rate of super-droplets within a grid box in the same way as regular droplets.

The super-droplet approach is first tested in the real DNS. Figure 4 shows evolutions of the standard deviation of the supersaturation spatial fluctuations¹ for real DNS of L=0.064 m (64³ grid points) and L=0.128 m (128³ grid points) with different multiplicity parameters. Number of super-droplets used and their multiplicities are listed in Table 1. For randomly distributed real droplets, the two simulations have droplets in one out of about eight grid boxes. For super-droplets in 64³ simulations, this number changes to 1 out of about 15 grid boxes for μ=2 and about 38 grid boxes for μ=5. For the 128³ simulations, super-droplets are on average in 1 out of about 43 and 87 grid boxes for μ=5 and 10, respectively. The mean supersaturation is close to zero, as expected (not shown). Supersaturation standard deviations fluctuate similarly in all simulations, with the mean values close among all multiplicities. The mean value of the standard deviation is larger for the larger domain, in agreement with simulations discussed in Lanotte et al. (2009).

Figure 4 Standard deviation of supersaturation fluctuations for 6.4³ cm³ real DNS (a) and 12.8³ cm³ real DNS (b). Colors represent different multiplicity as marked inside each panel.

Table 1 Number of super-droplets and their multiplicity for real DNS domains of volume 6.4³ and 12.8³ cm³.

In general, the multiplicity value should be decided on carefully because too large multiplicity results in too many grid boxes without droplets when compared to real droplets, and this may cause undesirable effects in the mean supersaturation and its spatial variability. In the two DNS cases, slight deviations in the mean supersaturation are present, although the simulations are not long enough to document the impact with confidence. For the scaled-up DNS, the number of droplets is in billions, and we have to select higher multiplicity values to make computations feasible.

The evolution of the radius squared (R²) standard deviation (${\mathit{\sigma }}_{R^{\mathrm{2}}}$) from the above simulations with droplets and super-droplets is shown in Fig. 5. Initially (i.e., at t=0), the distributions are monodisperse (i.e., ${\mathit{\sigma }}_{R^{\mathrm{2}}}=\mathrm{0}$). Supersaturation fluctuations in response to local vertical velocity fluctuations lead to the increase in ${\mathit{\sigma }}_{R^{\mathrm{2}}}$ in time. After some initial adjustment, the increase approximately follows the t^1∕2 scaling, with t being the time from the start of the simulation. This agrees with the study by Sardina et al. (2015), who applied a stochastic model and DNS. A similar result is also shown in Li et al. (2019) and Saito et al. (2019). As expected, the ${\mathit{\sigma }}_{R^{\mathrm{2}}}$ values are larger for the larger domain, in agreement with Fig. 3 in Lanotte et al. (2009) and Li et al. (2019).

Figure 5 Evolutions of the radius squared standard deviation (${\mathit{\sigma }}_{R^{\mathrm{2}}}$) for real DNS with different multiplicity parameters. Panels (a, b) correspond to 64³∕128³ real DNS ($L=\mathrm{6.4}/\mathrm{12.8}$ cm). Colors represent different multiplicities as marked inside each panel.

After the super-droplet technique is tested in DNS, the same method is used in scaled-up DNS. In general, one may expect that if the multiplicity is increased beyond a certain value, the results will start to deviate from those with a low multiplicity featuring a larger number of super-droplets. However, high multiplicity is desirable to reduce the number of super-droplets that need to be followed. For the DNS, the number of super-droplets was shown to be relatively low to maintain similarity between real droplet and super-droplet solutions (see Figs. 4 and 5; as low as one super-droplet in a few dozen grid volumes). With scaled-up DNS, one might expect a different requirement because of a stronger local forcing of the supersaturation due to higher TKE and thus larger vertical velocities.

For the scaled-up DNS study, we apply a 256³ domain to represent volumes with characteristic lengths of several meters. The TKE is scaled as explained in Sect. 3 with relevant parameters listed in Table 2. As the table shows, scaled-up DNS simulations typically have a relatively small number of super-droplets per grid box, similarly to DNS. This is because of computational efficiency considerations. However, one may question such an approach because scaled-up DNS include a large number of real droplets (e.g., ∼10⁹ for scaled-up DNS with 1 cm grid length and much larger numbers for scaled-up DNS with larger grid lengths). To show that the standard deviation of the supersaturation spatial distribution is not affected by the small number of super-droplets considered in the scaled-up DNS simulations, we included additional simulations (shown in bold in Table 2) that include about 10 super-droplets per grid volume and follow about 160 million super-droplets. Although arguably still a small number, 10 super-droplets per grid box is the number considered in one of the sensitivity simulations in Grabowski (2020). By comparing results of simulations with various numbers of super-droplet per grid volume, Grabowski (2020) shows that the number as small as 10 is sufficient to reasonably represent condensational growth of natural droplets in idealized simulations of laboratory chamber experiments.

Table 2 Details of DNS and scaled-up DNS. From left to right: domain length L, grid length l, viscosity ν, turbulent kinetic energy E, number of super-droplets in the domain N_s, multiplicity μ, and number of super-droplets per grid volume N_s∕N³. Numbers in bold correspond to the simulations with the highest number of super-droplets per grid volume.

Figures 6 and 7 present evolutions of the mean supersaturation and standard deviation of its spatial distribution for the scaled-up simulations from Table 2. The five scaled-up domains shown in the table and figures correspond to the domain size L of 2.56, 6.4, 12.8, 25.6, and 64 m. Note that the simulations extend to times of several minutes, that is, a significant fraction of a small convective cloud life cycle. All scaled-up simulations are run with three different multiplicities for super-droplets as listed in Table 2. As Fig. 6 shows, the mean supersaturation for all five scaled-up cases is close to zero after the initial spike. The spike magnitude, about 100 times smaller than the standard deviations shown in Fig. 7, increases as the domain size increases, and it is slightly larger for the higher multiplicity. Higher multiplicity also causes larger fluctuations after the initial spike, but the mean does not seem to be significantly affected. The standard deviation shown in Fig. 7 increases with the domain size, as expected. For all domains, standard deviations are similar for various multiplicities. In particular, based on 12.8 and 64 m simulations, the low number of super-droplets per grid volume (desirable for computational efficiency) seems to insignificantly impact the supersaturation statistics.

Figure 6 Evolution of the mean supersaturation for various scaled-up domains. Colors represent different domain sizes; different line styles correspond to different multiplicities. The additional simulation of 10 super-droplets per grid volume is only shown for 12.8³ and 64.0³ m³ volumes.

Figure 7 Evolution of standard deviation of supersaturation fluctuations for different domain sizes. Colors represent different domain sizes; different line styles correspond to different multiplicities. The additional simulation of 10 super-droplets per grid volume is only shown for 12.8³ and 64.0³ m³ volumes.

To further study the impact of the multiplicity, additional scaled-up DNS simulations are run with 256³ grid points for a domain of size 12.8³ m³. All simulations are listed in Table 3, with some already considered in Table 2 and Figs. 6 and 7. The total number of real droplets for the 12.8³ m³ domain with a droplet concentration of 130 cm⁻³ is about 2.7×10¹¹. The grid volume of the 256³ grid points and 12.83 m³ simulation is 125 cm³. When the multiplicity is 1625, the number of super-droplets is close to 170 million, and there are on average 10 super-droplets per grid volume. When μ is 7×10³, the number of super-droplets is close to 40 million, and there are on average about 2.3 super-droplets per grid volume. When μ is further increased to 1×10⁴, the number of super-droplets per grid volume decreases to about 1.6. For $\mathit{\mu }=\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$, the number further decreases to about 0.32 (i.e., a super-droplet in about three grid volumes). Finally, for $\mathit{\mu }=\mathrm{8.192}\times {\mathrm{10}}^{\mathrm{5}}$, a super-droplet is approximately in 1 out of 50 grid volumes.

Table 3 Number of super-droplets and multiplicity for different 12.8³ m³ scaled-up domain simulations.

Results obtained from these simulations are shown in Fig. 8, with some results already shown in Fig. 7. As the figure shows, only the largest multiplicity with a super-droplet in 1 out of 50 grid volumes differs significantly from other simulations. The highest multiplicity simulation also results in the non-zero mean supersaturation (not shown). Note that for real DNS (Table 2 and Fig. 4), having a droplet in one of several dozens of grid volumes still results in supersaturation fluctuations in agreement with real droplets. This suggests that the maximum multiplicity that can be used in scaled-up DNS depends on the domain size. This perhaps should not be surprising because the magnitude of the vertical velocity perturbation and thus the supersaturation forcing increases with the domain size. Results for the largest domain considered in the current study (64³ m³) suggest that the multiplicities selected for the scaled-up DNS provide robust (i.e., independent of the multiplicity) outcomes.

Figure 8 Standard deviation of supersaturation fluctuations for the 12.8³ m³ scaled-up domain. Colors indicate different multiplicities.

As shown in Fig. 9, evolutions of the radius squared standard deviation ${\mathit{\sigma }}_{R^{\mathrm{2}}}$ for scaled-up DNS domains follow the same trend as in the real DNS; that is, the standard deviations increase in time t as t^1∕2. The results are shown for the five scaled-up cases mentioned above. Scaled-up DNS for each domain was run for three different multiplicity values, one of them being 10 super-droplets per grid volume. The error bars correspond to the standard deviation among realizations with different multiplicities. Overall, the scatter resulting from different multiplicities is relatively small. The key result in Fig. 9 is that the spectral width increases with the domain size. For domain sizes of a few tens of meters, the spectral width after a few minutes reaches values of 1–2 µm, which is comparable to those observed in near-adiabatic cores of small cumuli (e.g., Jensen et al., 1985) or subtropical stratocumulus (e.g., Pawlowska et al., 2006).

Figure 9 Evolutions of the radius squared standard deviation ${\mathit{\sigma }}_{R^{\mathrm{2}}}$ for different domain sizes in the scaled-up DNS simulations. Horizontal bars along each line show variability resulting from the different multiplicity used for each domain size.

5 Stochastic model

We apply the stochastic model similar to that in Grabowski and Abade (2017) to simulate fluctuating supersaturation and compare results to the real and scaled-up DNS. The fluctuating in space supersaturation in the dynamic simulations (i.e., real DNS or scaled-up DNS) is modeled in the stochastic model as independent realizations of the fluctuation in time supersaturation as described below. For each realization, the supersaturation fluctuations are driven by the vertical velocity fluctuations as given by the Ornstein–Uhlenbeck process (e.g., Pope, 1994). In its finite difference implementation, the velocity perturbations are updated as in Grabowski and Abade (2017):

$$\begin{array}{}\text{(11)}& w{}^{\prime }(t+\mathit{\delta }t)=w^{\prime }\left(t\right)e^{-\mathit{\delta }t/\mathit{\tau }}+\sqrt{\mathrm{1}-e^{-\frac{\mathrm{2}\mathit{\delta }t}{\mathit{\tau }}}}{\mathit{\sigma }}_{w^{\prime }}\mathit{\psi },\end{array}$$

where δt is the model time step and ${\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}$ is the vertical velocity variance obtained from TKE as

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_{w^{\prime }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}E.\end{array}$$

ψ is a Gaussian random number with zero mean and unit variance generated at every time step; τ is the eddy turnover time calculated as

$$\begin{array}{}\text{(13)}& \mathit{\tau }={\displaystyle \frac{L}{(\mathrm{2}\mathit{\pi }{)}^{\mathrm{1}/\mathrm{3}}}}{\left({\displaystyle \frac{C_{\mathit{\tau }}}{E}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where C_τ is a constant equal to 1.5 as in Lasher-Trapp et al. (2005). Supersaturation fluctuations evolve according to the equation

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}S{}^{\prime }}{\mathrm{d}t}}=a_{\mathrm{1}}w{}^{\prime }-{\displaystyle \frac{S{}^{\prime }}{{\mathit{\tau }}_{\mathrm{relax}}}},\end{array}$$

where w^′ is the vertical velocity perturbation fluctuating in time as in Eq. (11), a₁ is a temperature-dependent numerical coefficient, and τ_relax is the phase relaxation time that depends on the temperature and pressure as well as on the droplet concentration and mean radius. For the conditions considered in this study, $a_{\mathrm{1}}=\mathrm{4.753}\times {\mathrm{10}}^{-\mathrm{4}}$ m⁻¹ and τ_relax=3.513 s.

The stochastic model used here applies 1000 realizations, each starting from a random velocity perturbation (i.e., ${\mathit{\sigma }}_{w{}^{\prime }}\mathit{\psi }$ as in Eq. 11) and zero supersaturation and run for six eddy turnover times. The time step in Eqs. (11) and (14) is taken as one-thousandth of the eddy turnover time. The number of realizations is sufficient to give results that change insignificantly when the number is further increased. Subsequently, the standard deviation of the supersaturation temporal evolution for each realization is derived. Its mean value averaged over all realizations together with the standard deviation among realizations is used in the comparison with the DNS and scaled-up DNS simulations.

Figure 10 shows the standard deviation of supersaturation fluctuations (σ_S) derived from the stochastic model as explained above for different domain sizes together with similar results from the DNS and scaled-up DNS dynamic simulations. The first five points (L=0.064, 0.128, 0.256, 0.512, 1.024 m) correspond to the real DNS, whereas the last five points (L=2.56, 6.4, 12.8, 25.6, 64.0 m) correspond to scaled-up DNS. The stochastic model uses TKE simulated by either DNS or scaled-up DNS, and for L=0.512 and 1.024 m by both. The vertical lines for the stochastic model show twice the standard deviation among the realizations. The red circles in the left part of the figure are from DNS simulations. Standard deviations from different multiplicities as shown in Fig. 4 are smaller than the circle radius, and thus they are not shown. The blue symbols are for scaled-up 256³ DNS simulations; standard deviations from different multiplicities are again smaller than the radius. The red circles for L=64 m show σ_S in scaled-up DNS simulations with grids of 128³ and 512³. The red line shows results from the stochastic model without droplets, that is, with no last term on the right-hand side of Eq. (14). Finally, the blue line is for the stochastic model assuming quasi-equilibrium supersaturation $S{{}^{\prime }}_{\mathrm{eq}}$, that is, with the left-hand side of Eq. (14) equal to zero and thus $S{}^{\prime }=S{{}^{\prime }}_{\mathrm{eq}}=a_{\mathrm{1}}w{}^{\prime }{\mathit{\tau }}_{\mathrm{relax}}$.

Figure 10 Standard deviation of supersaturation fluctuations in DNS and scaled-up DNS (color circles) and the stochastic model (black stars). Vertical lines for the stochastic model represent variability among individual realizations. For DNS and scaled-up DNS, the variability comes from different multiplicity for super-droplets; it is not shown as it is smaller than the symbol size. Red circles in the left half are for DNS with three different multiplicities as in Fig. 4. Blue circles are for 256³ scaled-up DNS simulations. Red circles for L=64 m are for 128³ (lower symbol) and 512³ (upper symbol) scaled-up DNS simulations. The two data points for the stochastic model with L=0.512 m and L=1.024 m come from applying TKE from either DNS or scaled-up DNS. The red and blue lines show results from the stochastic model without droplets and with the quasi-equilibrium supersaturation, respectively. See text for details.

Overall, the stochastic model seems to reasonably represent the scale dependence of the supersaturation fluctuations. At small scales (i.e., L=0.064 and 0.128 m), presence of droplets has a small impact on supersaturation fluctuations and the scale dependence is approximately as for the case without droplets as shown by the red line. At those scales, DNS seems to underestimate supersaturation fluctuations. Arguably, this is because of the small Reynolds number and thus a poor separation between forcing and dissipation scales. For scaled-up DNS, the stochastic model underestimates supersaturation fluctuations and the spread between the scaled-up DNS and stochastic model increases with the increase in the spatial scale. At the largest scales considered (i.e., L=12.8, 25.6, and 64 m), the quasi-equilibrium supersaturation provides a good estimate of the supersaturation fluctuations as shown by the blue line. This agrees with the argument put forward in Grabowski and Wang (2013; see the discussion in the last paragraph of Sect. 3.4).

There are a few reasons for the discrepancy between the stochastic model and scaled-up DNS. First, the stochastic model uses TKE obtained from the scaled-up DNS. However, scaled-up DNS features reduced TKE when compared to the real DNS as documented in Sect. 3. Allowing more TKE on input for the stochastic model would shift the stochastic model results upwards, that is, closer to the scaled-up DNS. But increasing the Reynolds number in the scaled-up DNS increases σ_S as well. This is illustrated by three data points for L=64 m scaled-up DNS with 128³, 256³, and 512³ simulations. Second, scaled-up DNS excludes scales of motion that are smaller than the scaled-up Kolmogorov microscale. For instance, for L=25.6 m and 256³ simulation, the scaled-up Kolmogorov microscale is 0.1 m. Hence, scales of motion between 10 cm and 1 mm are excluded when compared to the real DNS. Arguably, these small-scale motions in real DNS can affect supersaturation fluctuations and reduce σ_S. Such an argument seems to be contradicted by the results with L=64 m because σ_S increases, not decreases, between 128³, 256³, and 512³ scaled-up DNS. However, it is unclear whether the increase in σ_S with the further increase in the Reynolds number (i.e., the number of grid points) continues once real DNS is approached with further increase in the simulation towards the 64 000³ real DNS limit. Finally, one might argue that assuming a Gaussian vertical velocity distribution in Eq. (11) is an increasingly poor assumption with the increase in the domain size. Higher frequency of large vertical velocity perturbations (i.e., above the Gaussian distribution) should result in larger supersaturation fluctuations.

6 Discussion and conclusions

This study presents a novel modeling methodology that extends the traditional technique to simulate homogeneous isotropic turbulence, the direct numerical simulation (DNS). DNS is typically used for small-scale simulations applying grid lengths on the order of the Kolmogorov microscale, that is, about a millimeter for typical levels of atmospheric turbulence. Such a choice allows proper dissipation of the turbulent kinetic energy (TKE) that cascades through the inertial range from large scales where TKE is introduced. To reach domain sizes of about 1 m³ and beyond with a grid length of about 1 mm requires tremendous computation resources, with simulations featuring spatial scales of tens of meters and beyond (i.e., volumes of 1000s of cubic meters and larger) impossible for the foreseeable future. At the same time, one should expect the largest turbulent eddies to affect the diffusional growth of cloud droplets most significantly because such eddies feature the largest and longest-lasting vertical velocity and supersaturation perturbations. It is thus desirable to have a modeling approach similar to the traditional DNS but capable of reaching significantly larger spatial scales, say, tens and hundreds of meters.

This paper presents such an approach. The key idea is simple: rather than assuming that the dynamic model grid length is the Kolmogorov microscale η, we start with the DNS domain size L and adjust the Kolmogorov microscale given the computational resources. For instance, for L=100 m and 512³ simulation, η≃0.2 m. However, to use a traditional DNS code, one needs to allow proper TKE dissipation (as well as the scalar variance removal) at the smallest scales. It follows that the molecular transport coefficients from the traditional DNS need to be properly increased. The Reynolds number similarity is applied to develop a proper scaling; see Eqs. (8) and (10). The modified modeling approach is referred to as the scaled-up DNS. Section 2 presents numerical simulations applying the spectral DNS code that document robustness of the scaled-up DNS technique. We show that DNS and scaled-up DNS simulations with the same Reynolds number (i.e., the same L∕η ratio) forced to maintain the same TKE dissipation rate feature the expected TKE scaling (Eq. 4). However, when real DNSs are replaced by scaled-up DNSs with a reduced Reynolds number (i.e., keeping L the same and increasing η), a small fraction of the TKE is lost. As one might expect, the closer scaled-up DNS L∕η ratio is to the real DNS ratio, the closer TKEs are between the two simulations (cf. Fig. 3 and its discussion).

For simulations targeting growth of cloud droplets in homogeneous isotropic turbulence, the scaled-up DNS faces the problem of a large number of droplets that need to be followed inside the computational domain. For instance, a cube volume with L=100 m and droplet concentration of 100 cm⁻³ contains about 10¹⁴ droplets. Following all of them is computationally not possible. We apply a method already used in Lanotte et al. (2009) and in Li et al. (2019) and referred to as the super-droplet method in Shima et al. (2009). A super-droplet represents an ensemble of real droplets with the same radius; position and velocity of each super-droplet are predicted in the same way as for the real droplet. The number of real droplets represented by a single super-droplet is referred to as the multiplicity attribute (Shima et al., 2009). The multiplicity attribute is included in the condensation rate calculations. The super-droplet method is first tested in the real DNS and then implemented in the scaled-up DNS. Real DNS with L=0.064 m (64³ grid points) and L=0.128 m (128³ grid points) and with different multiplicity parameters gives consistent results even if the multiplicity parameter results in a super-droplet present in one out of a few dozen grid boxes. For scaled-up DNS (and likely for the real DNS as well), there is an upper limit for the multiplicity parameter before supersaturation fluctuations start deviating from the value obtained with lower multiplicities. Scaled-up DNSs presented here suggest that there should be at least a few super-droplets per grid box for approximately converged solutions. Such an estimate agrees with the result of idealized laboratory cloud chamber simulations reported in Sect. 4 of Grabowski (2020).

The scaled-up DNSs starting from unimodal droplet distribution with no mean ascent (i.e., as in Lanotte et al., 2009; Li et al., 2019) extend the validity of the scaling relationship obtained previously in either DNS simulations (e.g., Li et al., 2019; Saito et al., 2019) or in stochastic model simulations (Sardina et al., 2015). The scaling implies that the standard deviation of the droplet radius squared increases in time t as t^1∕2. DNS results of Li et al. (2019) show that the evolution of the droplet distribution spread depends on the Reynolds number (i.e., the DNS domain) and is insensitive to the TKE dissipation rate. The Reynolds number dependence is consistent with the eddy-hopping argument (see Sect. 3.5 in Grabowski and Wang, 2013) and the dominating impact of the largest eddies for the spread of the droplet size distribution. The standard deviation of the droplet radius squared increases in our simulations as t^1∕2 as well, with systematically larger values for larger scaled-up DNS domains as shown in Fig. 9.

Finally, we also consider supersaturation fluctuations in a simple stochastic model of a droplet ensemble (Grabowski and Abade, 2017) and compare the fluctuations to those simulated in DNS and scaled-up DNS. The key advantage of the stochastic model is that its computational cost is just a tiny fraction of a DNS simulation. The simulation time of the stochastic model is typically a mere few seconds on a laptop computer compared to hours of wall clock time of a high-performance massively parallel computer applied in DNS and scaled-up DNS simulations. As argued in Grabowski and Abade (2017), the stochastic model provides a simple and physically appealing approach to multiscale large-eddy simulation of a cloud applying Lagrangian particle-based microphysics (see Grabowski et al., 2019, for a discussion of the Lagrangian microphysics).

The scaled-up DNS methodology presented here was developed with diffusional growth of cloud droplets in mind. The next step can be to apply this approach in a rising parcel simulations as in Grabowski and Abade (2017) to understand the impact of turbulence on the cloud condensation nuclei activation/de-activation near the cloud base (see discussion in Abade et al., 2018). The rising parcel framework includes a prescribed ascent rate that affects the mean temperature inside the DNS domain and a reduction of the prescribed environmental pressure as the parcel rises; see Eqs. (1) and (3) in Grabowski and Abade (2017). One can argue that scale-dependent supersaturation fluctuations can induce significant droplet concentration heterogeneities at the cloud base that arguably affect droplet growth aloft. One may also consider applying the scaled-up DNS to the problem of droplet collisions. However, since collisions between cloud droplets take place at sub-Kolmogorov scales, applying scaled-up DNS for turbulent collisions is not straightforward. Finally, one can also consider applying scaled-up DNS in simulations of the turbulent entrainment and mixing similar to those discussed in Kumar et al. (2018) and Paoli and Shariff (2009). Such simulations would extend the still relatively small-domain DNS simulations into domain sizes comparable to the large entraining eddies in natural cumuli as discussed in Grabowski and Clark (1993). We hope to explore some of these research directions in the future.