A DNS study of aerosol and small-scale cloud turbulence interaction

The purpose of this study is to investigate the interaction between small-scale turbulence and aerosol and cloud microphysical properties using direct numerical simulations (DNS). We consider the domain located at the height of about 2000 m from the sea level, experiencing transient high supersaturation due to atmospheric fluctuations of temperature and humidity. To study the effect of total number of particles (Ntot) on air temperature, activation and supersaturation, we vary Ntot. To investigate the effect of aerosol dynamics on small-scale turbulence and vertical air motion, we vary the intensity of turbulent fluctuations and the buoyant force. We find that even a small number of aerosol particles (55.5 cm−3), and therefore a small droplet number concentration, strongly affects the air temperature due to release of latent heat. The system comes to an equilibrium faster and the relative number of activated particles appears to be smaller for largerNtot. We conclude that aerosol particles strongly affect the air motion. In a case of updraught coursed by buoyant force, the presence of aerosol particles results in acceleration of air motion in vertical direction and increase of turbulent fluctuations.

periencing transient high supersaturation due to atmospheric fluctuations of temperature and humidity. To study the effect of aerosol dynamics on the turbulence we vary the total number of particles (N tot ). In turn, to investigate the effect of small-scale turbulence on evolution of aerosol particles we 10 vary the intensity of turbulent fluctuations and the buoyant force. We find that even small amount of aerosol particles (55.5 cm −3 ) increases the air temperature by 1 K under supersaturated conditions due to release of latent heat. The system comes to an equilibrium faster and the relative number 15 of activated particles appears to be smaller for larger N tot . We conclude that the presence of aerosol particles results in deceleration of air motion in vertical direction and damping of turbulent fluctuations. 20

Introduction
Interaction of atmospheric turbulence with aerosol and cloud formation processes has been studied extensively. Due to non-linearity of particle formation and other aerosol dynamical processes, the fluctuations of temperature and rel-25 ative humidity can have strong effect on aerosol formation. Large scale fluctuations of atmospheric properties, which occur for example in the atmospheric boundary layer, can be the drivers for initiation of particle formation (Easter et al., 1994). Mixing of air with different properties, including tem-30 perature and relative humidity, have been shown to enhance atmospheric nucleation significantly (Nilsson and Kulmala, 1998). More specifically, the atmospheric waves can increase the nucleation rate several orders of magnitude and affect also the size spectum of the particles (Nilsson et al., 2000). 35 Also the activation of atmospheric particles in cloud areas is affected by the fluctuation of supersaturation, where some droplets were shown to grow also in undersaturated conditions due to fluctuations and the bi-modal particle size distribution could be observed after initially unimodal particle 40 population experienced fluctuating supersaturation (Kulmala et al., 1997).
In-cloud turbulence has been shown to intensify cloudmicrophysical processes determining cloud properties (Benmoshe and Khain, 2014). Similarly, aerosol loadings also in-45 fluence the numbers and sizes of cloud-particles, thereby influencing precipitation of formation, cloud extent and hence the climate (Forster, 2007). Both turbulence and aerosol particles can influence the chemical reactions in the atmosphere providing surfaces for aqueous phase chemistry and promot-50 ing uptake of gaseous species. This can have repercussions for chemical reactions in the air. Yet both environmental factors, namely aerosol composition and turbulence, have been usually considered as separate influences.
The large-scale atmospheric turbulence is well known to 55 affect aerosol processes significantly. In turn, aerosol properties may influence the buoyancy and turbulence. A motivation for a potential aerosol-turbulence interaction is that the most advanced cloud-microphysical models can now resolve the cloud edges where mixing occurs and where spatial 60 gradients are great in turbulence and concentrations of hydrometeors, including aerosols. Consequently, by including a more complete set of such interactions in the models, a better simulation of aerosols, turbulence and microphysics may be possible.

65
The main goal of this direct numerical simulation (DNS) study is to investigate interaction between small-scale turbulence and aerosol and cloud microphysical properties. Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2015-913, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 22 January 2016 c Author(s) 2016. CC-BY 3.0 License. case 1 N tot = 55.5 cm −3 case 2 N tot = 555 cm −3 case 3 N tot = 5550 cm −3 case 4 no particles  Table 2. The summary of key parameters for studying the effect of aerosol on turbulence. All four sets of parameters are used to compare cases 1 and 4 (see Table 1).
The choosen DNS domain is realistic for a small volume at cloud edge where turbulent mixing is a dominat fea-70 ture (Katzwinkel et al., 2014) and can create fluctuations in temperature as well as humidity. Such conditions however appear transient because presence of atmospheric aerosols leads to depletion of supersaturation via condensation and droplet activation. 75 We use the high order public domain finite difference PENCIL Code for compressible hydrodynamic flows. The code is highly modular and comes with a large selection of physics modules. It is widely documented in the literature (Dobler et al., 2006;Pencil Code, 2001, and references 80 therein). The chemistry module is responsible for the detailed description of the necessary quantities in a case of complicated chemical composition, such as diffusion coefficients, thermal conductivity, reaction rates etc. (Babkovskaia et al., 2011). The detailed description of aerosol module is 85 in Babkovskaia et al. (2015). The paper is constructed as follows. Section 2 is devoted to the description of the model. Results are presented in Section 3. Section 4 is the summary.
2 Description of the model As a starting point for our simulations we consider values typical for observations made in trade wind cumuli. During the CARRIBA project (Siebert et al. 2013, hereafter called SI13) total aerosol number concentrations in the marine subcloud layer up to a few 100 cm −3 have been observed ( Fig   95   6f in SI13). Although it was argued that the highest values are due to local biomass burning we consider a concentration of 550 cm −3 as typical and 55 and 5550 cm −3 as two extreme values for sensitivity test. The initial normalized aerosol number size distribution as shown as a dotted line 100 in Fig. 3 in our manuscript compares well with the shape of the observed distribution shown as red line in Fig. 8 We assume a soluble aerosol (NaCl) which will dilute inside droplet. We take 50 size bins logarithmically distributed in the range [10 nm, 1000 µm]. As an initial distribution 105 of particles we take the observational data at the sea level and assume that it is the same everywhere in the domain (see the distribution in Fig. 3). To analyse the effect of total number of aerosol particles, N tot , on the structure and properties of turbulent motion we consider the following cases: 110 N tot = 55.5 cm −3 (case 1), N tot = 555 cm −3 (observed data, case 2), N tot = 5550 cm −3 ( case 3), and no particles (case 4), see Table 1.

Air composition
We assume the following air composition O 2 + H 2 O + N 2 , 115 where nitrogen mass fraction is taken to be Y N2 = 70%. The observations provide us the absolute humidity to set the initial value for mass fraction of water vapour, Y H2O . Oxygen mass fraction is recalculated from the normalization conditions, i.e. Y O2 + Y H2O + Y N2 = 1, where Y O2 is oxygen mass 120 fraction.

Initial conditions
Based on data of CARRIBA observations typical for the upper parts of clouds / cloud edges in a height of 2000 m, we set the initial conditions for air temperature (T 0 = 285.4 K) and 125 water vapour mixing ratio (q 0 = 0.0124). The small vertical gradients of temperature (∆T = 0.001 K) and water content (∆q = 4 × 10 −5 ) are also based on the CARRIBA measurements.
Direct observations of supersaturation fluctuations are not 130 possible so far and one has to rely on estimates based on measurements of air temperature and water vapor mixing ratio. The measured temperature and absolute humidity led us to a range of supersaturation (S ) from 10.17 % to 10.54 %. Such values of S are extremely high in a case of an adia-135 batic lifting of a cloudy air parcel. The "quasi-steady state" supersaturation depends mainly on vertical updraft velocity and appears to be on the order of a few tenths of percent for vertical updrafts of about 1 m s −1 . However, it is well known that atmospheric clouds are mainly non-adiabatic due 140 to turbulent mixing, and a few percent supersaturation are realistic for higher updraft velocities (e.g., Korolev and Mazin, 2003). This is particular true at cloud edges where entrainment of unsaturated air into the cloud results in strong mixing and fluctuations of the water vapor and temperature field. It  provided some convincing arguments that at the cloud base high fluctuations of supersaturation on the order of up to several percent can exist. It is also well known that turbulence in high Reynolds number flows (typical in convective clouds) is 155 highly intermittent (Siebert et al., 2010). Shaw (2000) argued that under such conditions long-living vortex tubes could produce small areas with decreased droplet number concentrations and, therefore, high supersaturation resulting in secondary activation. Summarizing, there are good arguments 160 that supersaturation of a few percent can be generated without strong updrafts. On the other hand, the supersaturation excess would be eliminated by condensation onto droplets and the usual equilibrium supersaturation would be restored. Therefore, the key 165 question about the temporal time scale is under consideration now. The ratio of two time scales is important here: the phase relaxation time which describes how fast the supersaturation can react on the new thermodynamic condition and the turbulent mixing time scale which describes how fast turbulence 170 can mix a certain volume (eddy with the length scale l). If the phase relaxation time is smaller compared to the turbulent mixing time than the quasi-steady state solution would be applicable. However, for scales where the turbulent mixing is faster we expect strong supersaturation fluctuations to 175 "survive".
Let us now assume a small eddy of size l = 1 m and a local energy dissipation rate of ǫ = 0.1m 2 s −3 . These dissipation values are typical peak values for cumulus clouds on that small scales. The highest dissipation can be found at Kol-180 mogorov size (see discussion in Siebert et al., 2006Siebert et al., , 2010. The eddy turn-over time is τ eddy = (l 2 /ǫ) 1/3 ≈ 2 s. If we now take a phase relaxation time of the order of one second which is typical for cumulus clouds (see again Korolev and Mazin, 2003) we see that these two time scales are of the same or- der. Therefore, we conclude that on scales below one meter, strong supersaturation fluctuations can exist and the quasisteady state solution is not applicable anymore. This argumentation partly follows the discussion in Section 8e by Korolev and Mazin (2003).
190 Kulmala et al. (1997) estimated standard deviations of supersaturation of up to 5 % based on aircraft observations at cloud level but outside the clouds. Ditas et al. (2012) observed supersaturation fluctuations in a field of stratocumulus clouds and estimated from highly collocated temperature 195 and water vapor observations peak-to-peak values of up to 1.5 %, which is much higher compared to the quasi-steadystate solution. It is straightforward to assume that higher S values can be expected for parts of (shallow) cumulus clouds. Thus, we argue that for our small modeled volume a transient 200 supersaturation of 10 % can be realistic for specific mixing event.
The pressure, p, is assumed to be constant everywhere in the domain and it is also based on the measurements. The air density, ρ, is calculated from the equation of state p = 205 ρRT/m, where R is gas universal constant, m is air molar mass. The initial velocity is taken to be zero.
To generate the initial turbulent field we make first 100 iterations without evaporation/activation of aerosol particles (further "aerosol dynamics"), but with included randomly di-210 rected external forces (see next section). After that the external forces are set to zero, whereas the particles start to evolve. Thus, the turbulence is decaying for the analysed time. The maximal timestep allowed by the Courant condition for convergence is ∆t c = 10 −6 s (Courant time step). Since ∆t c is 215 much larger than the time step needed to move the smallest particle to the neighbour size bin (∆t a = 2 × 10 −7 s), at every Courant time step we make 5 substeps and calculate the particle evolution equation only.

Boundary conditions 220
In all three directions we set periodic boundary conditions for all variables, including the number density function. It means, that at every time step ∆t c the number of particles appearing on the bottom/left boundaries is equal to the number of particles disappearing through the top/right boundaries.

225
While the periodic boundary conditions modify the initial temperature stratification, they allow us to consider this domain as an isolated volume, i.e. the total mass, energy and number of particles in the domain does not change with time. In turn, it makes possible to compare the results of simula-230 tions, varying the key parameters of the model and to carry out the detailed quantitative analysis of the interaction between aerosol and turbulence.

Basic equations
The detailed description of the main equations is presented 235 by Babkovskaia et al. (2015). We consider the standard Navier-Stokes system including equation for conservation of mass, momentum, energy and chemical species. In the present study we modify the model, describing the evolution of aerosol number density function: water vapour pressure 240 over a droplet of radius r is defined as (Seinfeld and Pandis, 2006;Sorjamaa and Laaksonen, 2007) where p 0 is water vapor pressure over a flat surface at the same temperature, A = 0.66/T (in µm), where r 0 = 10 nm is 245 the radius of the droplet core, d w = 0.3 nm is the size of water molecule.
For generation of the initial turbulence the external forces f is used in a form 250 where x is the position vector. The wave vector k(t) and the random phase −π < φ(t) ≤ π change at every time step; N = f 0 c s (|k|c s /∆t c ) 1/2 is the normalization factor, c s is the sound speed, f 0 is a non-dimensional forcing amplitude; f k = (k × e)/ k 2 − (k · e) 2 ), where e is an arbitrary unit vector that is 255 real and not aligned with k (see detail in Pencil Code (2001)).
To study the effect of aerosol on the turbulence we consider two initial turbulent fields, taking different nondimensional forcing amplitudes as f 0 = 10 (low intensive turbulence) and f 0 = 100 (high intensive turbulence). Also, we 260 compare the effect of the aerosol in the case when the domain is in equilibrium with environment (equilibrium case) and when the temperatures inside and outside the domain are different (unequilibrium case), taking the environment temperature T 0 = 285.4 K for equilibrium case and T 0 = 293 K 265 for unequilibrium case . The buoyancy force is expressed as follows where g = 9.81 m s −2 is the acceleration of gravity, q 0 is the reference water vapor mixing ratio, ǫ + 1 = R v /R d is the ra-270 tio of the gas constant for water vapor and dry air, and q c is the cloud water mixing ratio (Andrejczuk et al., 2004;Babkovskaia et al., 2015). The summary of four considered cases is in Table 2.

275
Effect of total number of particles on air temperature and supersaturation distributions In Fig. 1 we present the air temperature averaged in ydirection. In Fig. 1 (a, b, c) we present the temperature distribution in the cases where the aerosol dynamics is included. 280 We find that without aerosol (case 4) after 3 s the difference   case 1 case 2 case 3 N tot ,cm −3 55.5 555 5550 N act , cm −3 49 493 4700 LWC, g/cm 3 2.3 10 −7 3.69 10 −7 3.76 10 −7 ∆ T, K 0.66 1.04 1.04 ∆ S max , % -6.77 -10.5 -10.5 ∆ S min , % -6.48 -10.13 -10.12 Table 3. Total number of particles (N tot ), number of activated particles (N act ), liquid water content (LWC), change of the temperature (∆ T), maximal values of supersaturation (∆ S max ), minimal values of supersaturation (∆ S min ) at t = 3 s for three considered cases. Unequilibrium case is considered. between the absolute value of maximal (285.381 K) and minimal (285.380 K) temperatures in the domain is still about 0.001 K, while in case 3 this difference achieves 0.02 K (see Fig. 1, d). Moreover, the temperature increases by about 1 K achieves maximum. It happens because of the different effect of aerosol dynamics in dependence of N tot . The small amount of aerosol particles (case 1) does not make any substantial 300 effect on the temperature distribution, i.e. both temperature and supersaturation are identically shifting with time in vertical direction. In turn, in cases 2 and 3 aerosol dynamics crucially changes the temperature distribution. One can see in Figs. 1, 2 that layers with larger temperature correspond 305 to layers with smaller supersaturation. It can be interpreted as follows. More intensive condensation occurs in initially warmer layers because supersaturation is larger there and respectively the temperature grows faster in these layers. In turn, since supersaturation exponentially depends on temper-310 ature, S ∝ exp (−T ), at some moment S appears to be smaller in warmer than cooler layers. Thus, equilibrium supersaturation is higher in the layers with temperature minimum (and vice versa). Also, we find that in case 1 at t = 3 s the supersaturation is about 5 % and the aerosol is still activating, 315 whereas in case 3 for the first 3 s the supersaturation almost drops to zero and the system is coming to an equilibrium.
Additionally, we investigate how fast the system with initially high value of supersaturation (S = 10 %) comes to an equilibrium (S = 0 %) in dependence on different total num-320 bers of droplets to answer the question how N tot affects the phase relaxation time. Analysing Fig. 2 we find that in a case of N tot = 55 cm −3 it takes more than 3 s for the system to come to equilibrium, i.e. phase relaxation time appears to be larger than theoretically estimated turbulent mixing time 325 (τ eddy = 2 s) (see discussion in section 2). Thus, strong supersaturation fluctuations can "survive" longer if the total number of droplets is smaller.

Effect of total number of particles on activation
In this study we consider that all particles with radius larger 330 than r cr = 1.75 µm are activated. This value is somewhat arbitrary but the results of our study were not sensitive on the choice of r cr provided that r cr ≤ 1.75 µm. In Fig. 3 we present the normalized particle distribution at t = 0 s and t = 3 s. We find that the smaller is the total number, the larger particles 335 are produced for the similar time period. Comparing the supersaturation at t = 3 s in these three cases (Fig. 2) one can see that in the case of the largest N tot the supersaturation appears to be close to zero and particles stop to grow, while for the smaller N tot the supersaturation is about 3 % and par-340 ticles continue growing. Therefore, the system is coming to an equilibrium (i.e. S = 0) faster and aerosol particles cease to grow earlier for the larger N tot . Thus, our DNS results for a case without vertical displacement and resulting convective cooling confirm the Twomey effect and show that aerosol 345 particles can grow and possibly achieve the size of the rain drop only in the case of small total number of particles. This is consistent with the fact that in natural clouds, deep ascent is needed to drive the prolonged condensational growth of drops, in order for rain to form (e.g. Rogers and Yau 1989, 350 Fig. 6. Averaged vertical velocity as a function of time with aerosol for N tot = 5500 cm −3 (dotted curve); and without aerosol (solid curve). Unequilibrium and low intensive turbulence case is considered (see Table 3). Fig. 7. The dependence of the average kinetic energy on time with aerosol for N tot = 5500 cm −3 (dotted curve); and without aerosol (solid curve). The other parameters correspond to unequilibrium case with low intensive turbulence (see Table 2). Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2015-913, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 22 January 2016 c Author(s) 2016. CC-BY 3.0 License. Fig. 8. The difference between turbulent kinetic energy averaged over time period t in a case without aerosol (E noa T ) and with aerosols (E a T ), where E T (t) = 1/t t 0 T KE(t ′ )dt ′ . All other parameters are the same as in Fig.7. a short course in cloud physics). In real clouds, the ascent is the source of larger-scale supersaturation. One should also mention, that in the scope of this model we neglect collision and coalescence of aerosol particles, which can be crutial in creation of rain drops.

355
In addition, in Fig. 4 we analyse the amount of liquid water, LW(r) = 4/3πr 3 ρ w N(r), on particles with radius r (where N(r) is the number of particles with radius r, ρ w is liquid water density) and in Fig. 5 we show the number of activated particles averaged over the domain as a function of 360 time in cases 1, 2, 3. In Table 3 we collect the most important quantities, such as number of activated particles, liquid water content LWC = LW(r)/∆r dr (where ∆r is a size of corresponding bin), and the changes of temperature and supersaturation averaged over domain at t = 3 s for cases 1, 2, 365

3.
We find that while the total number in case 2 is ten times smaller than in case 3, liquid water content is similar in these two cases. On the other hand, LWC appears to be smaller in case 1 than in cases 2, 3 (see Table 3 and Fig. 4). It happens 370 because the probability of water molecules to catch a particle is much smaller in a case of the smallest particle concentration (case 1) than in cases 2, 3. In Fig. 5 one can see that in case 1 and case 2 the number of activated particles does not grow after 1.2 s and 0.6 s, correspondingly, and in case 375 3 it happens after 0.4 s. Since in case 3 the equilibrium is achieved earlier than in cases 1 and 2, the maximum of final particle distribution in case 3 is shifted (see Fig. 3) and the final relative number of activated particles appears to be smaller than in cases 1 and 2. We find that the number of 380 activated particles approximately linearly depends on the total number, whereas the change of N tot in 100 times leads to increase of LWC just by 40 % (Table 3).

Effect of aerosol on the turbulent motion
We analyse the effect of aerosol on the turbulent motion.

385
First, we consider the equilibrium case, taking T 0 = 285.4 K and q 0 = 0.0124 in Eq. 3. In that case the turbulent field appears to be isotropic. Next, we increase T 0 = 293 K and study the developing of turbulence in unequilibrium case. Then the intensive vertical motion is generated because of buoyancy 390 force. Note, that the model domain is not vertically displaced during the simulation time but the vertical motion is generated within domain. Also, we vary parameter f 0 in Eq. 2 to compare the developing of low intensive ( f 0 = 10) and high intensive ( f 0 = 100) turbulence in both equilibrium and un-395 equilibrium cases. The summary of parameters is presented in Table 2. Finally, we compare the results of simulations with particles (taking N tot = 5500 cm −3 ) and without them in all four cases described above.
We find that the vertical motion of air is decelerated be- without them. We conclude that turbulence is damped because of presence of aerosols. We interpret these results as 405 follows. The air temperature increases because of latent heating caused by condensation onto drops, and the temperatures inside and outside the domain is reduced. It results in decreased buoyancy force and deceleration in vertical direction. Thus, we conclude that aerosol affects the turbulence 410 via buoyancy.
Also, we find that deceleration does not depend on intensity of turbulent fluctuations, i.e. deceleration in vertical direction is the same in low intensive turbulence and high intensive turbulence cases. Moreover, turbulent fluctuations are 415 damped because of presence of aerosol particles in all four considered cases. However, the damping effect is stronger in unequilibrium with low intensive turbulence than in other three cases. The dependences of vertical velocity and turbulent kinetic energy (TKE) averaged over the domain as a 420 function of time for unequilibrium case with low intensive turbulence are presented in Figs. 6 and 7, correspondingly.
Finally, we find that there is a strong variation of TKE with time (see Fig. 7). To interpret this fact we plot x-,y-and z-components of TKE in Fig. 9 and find strong time vari-425 ations only in z-component of TKE (x-and y-components are smoothly decreasing with time). We conclude that it results from fluctuations of temperature caused by aerosol dynamics, and therefore, changes of buoyancy force with time, which results in perturbation of vertical motion.

Summary
Turbulence, aerosol growth and microphysics of hydrometeors in clouds are intimately coupled. In the present study a new modelling approach was applied so as to quantify this linkage. This model represents the 3D fluid flow on the mi-435 croscale inside a volume of 10 cm x 10 cm x 10 cm, just inside the cloud in the mid-troposphere. We study the interaction in the cloud area under transient, high supesaturation conditions, using direct numerical simulations. As the initial conditions we take the observational data obtained by the 440 ACTOS platform. To analyse the effect of aerosol on the turbulence we choose the part of the cloud where the supersaturation is about 10 % and condensation is the dominant process. As an initial distribution of particles we take the data of measurements at the sea level and analyse the drops activated 445 by the aerosols in simulations.
We compare the results of simulations with particles and without them. We find that the total number of particles in the domain is crucial for distribution of temperature and for developing of turbulence. Even small amount of aerosol par-450 ticles (55.5 cm −3 ) increases the air temperature by 1 K because of latent heating caused by condensation onto drops. The system comes to an equilibrium faster for the larger total number of particles. Also, we find that the larger is N tot the smaller is the relative number of activated particles (averaged 455 over the domain). Analysing the effect of aerosol dynamics on the turbulent kinetic energy and on vertical velocity we conclude that the presence of aerosol has an effect on vertical motion and in our case tends to reduce downward velocity. We conclude that aerosols quite strongly influences the 460 dynamics in the cloud area. Such effect of aerosols can be crucial also for large scales usually studied with Large Eddy Simulation (LES) and the LES parametrization can be improved with direct numerical simulations.