Extreme responses of a coupled scalar–particle system during turbulent mixing

One signature of turbulent flow is its efficiency in stirring and mixing [1], with the simplest scenario being the advection and diffusion of a continuous, passive scalar [2, 3]. In this work we consider advection–diffusion with a discrete scalar constituent, particles in a turbulent flow, which in turn are coupled to a continuous 'vapor' field. Active feedback to the turbulence itself is not considered, so that we may concentrate solely on the influence of turbulent mixing on the particle–continuous scalar coupling. The problem is of general interest as discrete advection–diffusion, but has direct relevance to a variety of applications; we cast the problem in terms of the mixing of cloudy air and clear air in the atmosphere.

Even the simplest models of cloudy convection depend sensitively on the extent and nature of entrainment of ambient, non-cloudy air, with implications for cloud depth, cloud optical properties, precipitation time scales and cloud lifetime [4]. These cloud properties play an important role in climate dynamics. The entrainment of clear air and its subsequent mixing with cloudy air occurs during the entire life of a cloud and introduces strong inhomogeneities at spatial scales ranging from 10² m down to 1 mm and at time scales from hours to seconds [5]. The entrainment and mixing process changes the water droplet size distribution through dilution, evaporation and perhaps enhanced collision. This distribution is of direct consequence to rain formation and cloud radiative properties and will eventually determine the cloud lifetime itself [6].

The mixing of cloudy and clear air is a complex process, not only because of the inherent role of turbulence with its broad range of coupled scales, but also because the scalar fields (e.g. droplet number density, water vapor concentration and temperature) are no longer 'passive': latent heating effects couple the energy and condensed phase equations, mass conservation couples the condensed and vapor phases and the condensed phase itself can respond through various pathways even for a single bulk thermodynamic outcome. In this work we neglect the first aspect, the active thermal feedback, in order to isolate the influence of the mixing process on the evolution of the droplet population. For example, upon mixing of cloudy air and clear air, droplets will evaporate until the mixture becomes saturated (assuming the initial presence of sufficient condensed water), but this could occur by all droplets evaporating by the same amount, or by a subset of droplets evaporating completely, leaving the remaining droplets unchanged. It has been argued that the microphysical response to mixing can be characterized by the Damköhler number, the ratio of a fluid time scale to a characteristic thermodynamic time scale associated with the evaporation process (phase relaxation time), i.e. $Da=\frac {\tau _{\mathrm {fluid}}}{\tau _{\mathrm {phase}}}$. The two limits, Da ≪ 1 and Da ≫ 1, describe homogeneous and inhomogeneous mixing, respectively [7–9]. Homogeneous mixing occurs when the condensational growth or evaporation of cloud water droplets is slow compared to the mixing and therefore takes place in a well-mixed environment. Inhomogeneous mixing occurs when the evaporation proceeds much faster than the flow structures evolve, with the result that droplets near the clear air–cloud interface experience evaporation while others do not. Both processes can coexist in a turbulent cloud because of the broad spectrum of fluid time scales that are present, with inhomogeneous mixing dominating at large scales, and homogeneous mixing occurring at fine scales [10]. Concomitantly, Damköhler numbers can be defined using the large-eddy time (Da_L) and the dissipation time (Da_η), and the length scale at which the transition from inhomogeneous to homogeneous mixing is expected to occur can be estimated as l_c = Da^−3/2_ηη [10, 11], where η is the Kolmogorov dissipation length scale.

We pose two questions for this work: firstly, does the Damköhler number capture all aspects of the mixing or, said another way, is there Damköhler number similarity? Secondly, how does the droplet size distribution respond differently under conditions of strongly homogeneous and inhomogeneous mixing? These two questions are directly motivated by recent work in which we investigated the response of the cloud droplet size distribution to mixing using direct numerical simulation (DNS) of turbulence containing a population of cloud droplets [11]. The approach, which we now extend to the questions posed here, was a combined Eulerian representation of the fluid turbulence and water vapor concentration fields and Lagrangian representation of the discrete droplets. Previous attempts at these questions which extend classical ideas of inhomogeneous mixing [9] were conducted by Krueger et al in a one-dimensional (1D) explicit mixing parcel model [12] and by Jeffery in an idealized eddy-diffusivity model [13].

Thermodynamic and microphysical conditions in our DNS model were representative of those existing in a cumulus cloud, but of course because of computational limitations the turbulence Reynolds number was much lower than in a real cloud, with the implication that the range of spatial scales was also limited. This limitation on the largest simulated scales constrained the large-eddy Damköhler number to values of Da_L ∼ 1 (see table 2 in [11]), and therefore limited our ability to investigate the full range of Damköhler numbers, from strongly homogeneous to strongly inhomogeneous. The inhomogeneous limit Da_L ≫ 1 is of particular interest because of the intriguing expectation that the size distribution will maintain a nearly constant mean diameter but with a steadily decreasing droplet number density (e.g. through complete evaporation of a subset of the droplets, with the remaining droplets not influenced). In this work, we reach the extremes of strongly homogeneous and inhomogeneous mixing, even with limited computation size, by directly varying the rate coefficient that couples the discrete, condensed phase and the continuous, vapor phase. While admittedly not a direct simulation of atmospheric clouds, we will argue that the study captures essential aspects of the response of the particle size distribution to mixing in the extreme Da limits.

While the mixing process has been the focus of previous computational investigations, those studies have primarily considered the problem solely from the continuum microphysics perspective [9, 14, 15]. We explicitly treat the Lagrangian nature of the discrete droplet field. Lagrangian treatment of cloud droplets has been carried out elsewhere [16, 17], with an emphasis on the evolution of the supersaturation field and the droplet population during steady growth, whereas this study is focused on the microphysical response to a transient mixing event. This paper is organized as follows. Section 2 reviews the governing equations and the simulation parameters relevant to this work. Section 3 considers the first question regarding Damköhler similarity, and section 4 describes the results of the investigation of mixing with extreme Damköhler numbers. The paper closes with a discussion and summary.

2.1. The Euler–Lagrangian equations of motion

The governing equations, simulation parameters and simulation geometry, except where noted below, are directly taken from our previous study and are described in detail there [11]. As an overview, the simulation is initialized with an idealized slab-cloud geometry that could be considered analogous to a cloud edge or cloud filament boundary formed through inhomogeneous mixing at larger scales. The initial slab-like filament of supersaturated vapor and cloud droplets fills approximately the middle one-third of the cubic simulation box, which has periodic boundary conditions. The incompressible turbulent flow is described by the velocity field u(x,t) and the pressure field p(x,t), and within this flow the vapor mixing ratio field q_v(x,t) is transported and diffuses. The vapor mixing ratio is defined as $q_{\mathrm {v}}({\bf x},t)=\frac {\rho _{\mathrm {v}}}{\rho _{\mathrm {d}}}$, where ρ_v and ρ_d are the mass densities of vapor and dry air, respectively. As discussed already, the advection–diffusion equation for temperature is not considered in this study. The Eulerian equations for the turbulent fields are

$$\begin{equation} \boldsymbol{\nabla}\cdot{\bf u}=0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \partial_t {\bf u} + ({\bf u}\cdot{\bf\nabla}) {\bf u} =-\frac{1}{\rho_0}{\bf\nabla}p+ \nu{\bf\nabla}^2{\bf u}+{\bf f}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \partial_t q_{\mathrm{v}} + {\bf u}\cdot{\bf\nabla} q_{\mathrm{v}} = D {\bf\nabla}^2 q_{\mathrm{v}} -C_{\mathrm{d}}, \end{equation} \tag{ 3 }$$

where f(x,t) is a bulk forcing that sustains the statistically stationary turbulence and C_d is the condensation rate. The Lagrangian evolution of the discrete scalar field, consisting of N droplets in the volume V , is described simply as

$$\begin{eqnarray} \frac{\mbox{d}{\bf X}(t)}{\mbox{d}t}={\bf V}(t) & = & {\bf u}(x,t).\, \end{eqnarray} \tag{ 4 }$$

Here, X is the droplet position and V its velocity, and we neglect droplet inertia and gravitational settling again, in order to isolate the essential aspects of the particle–scalar coupling during the mixing process. In our previous study [11], we investigated the effect of inertia and gravitational settling on the mixing process. For the parameters accessible, we concluded that particle inertia and gravitational settling lead to a more rapid initial evaporation, but ultimately lead only to a slight depletion of both tails of the droplet size distribution.

As droplets are advected by the fluid, they can grow or evaporate in response to the local vapor field (recalling that in this study the temperature is constant). Direct droplet interactions through collision are neglected in order to focus solely on the initial stage of the entrainment and mixing process. The condensation rate field C_d(x,t) is directly coupled to the discrete droplets, as in [16] and [11]. In detail, the condensation rate is given by

$$\begin{equation} C_{\mathrm{d}}({\bf x},t)=\frac{1}{m_{\mathrm{a}}} \frac{\mbox{d}m_{\mathrm{l}}({\bf x},t)}{\mbox{d}t} =\frac{4\pi\rho_{\mathrm{l}} K}{\rho_0 a^3} \sum^{\triangle}_{\beta=1}S({\bf X}_{\beta},t) r(t), \end{equation} \tag{ 5 }$$

where m_a is the mass of air per grid cell and the sum collects the droplets inside each of the grid cells of size a³ that surround the (grid) point x. This relation closes the system of Eulerian–Lagrangian equations (1)–(4). Switching between the Eulerian field values at grid positions and the enclosed droplet position is done by trilinear interpolation. The inverse procedure is required for the calculation of the condensation rate which is evaluated first at the droplet position and then redistributed to the nearest eight grid vertices.

For the derivation of the final expression in (5), we used the equation for diffusional droplet growth by condensation (detailed derivations are found in [4, 11]). Its derivation is based on the assumption of steady fluxes of vapor mass to the droplet surface and of released latent heat due to condensation away from the droplet. Energy conservation in a steady flux regime requires

$$\begin{equation} F_{\mathrm{v}} + L^{-1} F_Q=0, \end{equation} \tag{ 6 }$$

with L being the latent heat of vaporization. It is further assumed that the water vapor pressure at the droplet surface is at the saturation value e_s(T_r) and thus ρ_v,r = ρ_vs. Saturation pressure and density are connected via the ideal gas law e_s = R_vρ_vsT, where R_v is the vapor gas constant. The vapor flux across a sphere of radius R has to be equal to the change in liquid water mass inside the sphere

$$\begin{equation} 4\pi\!R^2 D \frac{\mbox{d}\rho_{\mathrm{v}}}{\mbox{d}R}=\frac{\mbox{d}M_{\mathrm{l}}}{\mbox{d}t}=4\pi\rho_{\mathrm{l}} r^2\frac{\mbox{d}r} {\mbox{d}t}. \end{equation} \tag{ 7 }$$

Integration from R = r to ∞ results in combination with an expression for saturation pressure (and thus for the vapor mass density at saturation) which follows from the Clausius–Clapeyron equation to

$$\begin{equation} r \frac{\mbox{d}r}{\mbox{d}t} \simeq \frac{\mathscr{D}}{\rho_{\mathrm{l}}}\, (\rho_{\mathrm{v},0}-\rho_{\mathrm{vs}}(T_0)). \end{equation} \tag{ 8 }$$

$\mathscr{D}$ is a modified diffusion coefficient for water vapor in air, accounting for the coupled transport of energy and vapor [11]. Quantities ρ_l and ρ_vs(T₀) denote the liquid water mass density and the vapor mass density at saturation at reference temperature T₀, respectively. With the definition of the supersaturation

$$\begin{equation} S({\bf x},t)=\frac{q_{\mathrm{v}}({\bf x},t)}{q_{\mathrm{vs}}(T_0)}-1,\, \end{equation} \tag{ 9 }$$

in which q_vs(T₀) = ρ_vs(T₀)/ρ_d is the vapor mixing ratio at saturation at the reference temperature T₀, the growth equation can be written finally as

$$\begin{eqnarray} &&r(t)\frac{\mbox{d}r(t)}{\mbox{d}t}= K\!S({\bf X},t), \end{eqnarray} \tag{ 10 }$$

in which the constant K is given as [4]

$$\begin{eqnarray} &&K= \left[ \rho_{\mathrm{l}} \left(\frac{R_{\mathrm{v}} T_0}{D e_{\mathrm{s}}(T_0)}+\frac{L^2 }{k R_{\mathrm{v}} T_0^2}\right)\right]^{-1}=\frac{\rho_{\mathrm{vs}}(T_0)}{\rho_{\mathrm{l}}}{\mathscr{D}}\, . \end{eqnarray}\noindent \tag{ 11 }$$

Parameter k denotes the thermal conductivity.

Thermal effects are accounted for at the individual droplet scale, but are neglected in their influence on the turbulent flow itself. Instead, the temperature is held fixed at a reference value of T₀ in the simulation domain, thus leaving the saturation vapor mixing ratio constant. This is clearly a simplification which decouples the temperature field completely from the dynamics, but allows us to focus on the mixing aspects in the present model. As a consequence, the turbulence is not driven by buoyancy effects as in similar studies [14–16], but by a volume forcing that mimics a cascade of kinetic energy from larger scales [18]. Finally, we wish to also mention that volume considered in this study will be fairly small, such that strong temperature variations over a time span of 1 min or less are not very probable.

2.2. Initial conditions

The initial condition of the turbulent velocity field was generated by running a separate flow simulation with the same volume forcing as stated in equation (2). Such a forcing term f injects a fixed amount of turbulent kinetic energy into the flow at the largest scales and ensures a statistically stationary homogeneous isotropic turbulence. The preparation run is conducted for a time lag of 15 s. The turbulent kinetic energy and the volume average of the kinetic energy dissipation rate fluctuate then about their temporal means. This flow field served as an initial state for all parameter studies reported here. In order to focus on mixing from a state possessing a single, clearly defined initial length scale, we have chosen to begin with an idealized slab-cloud geometry. This slab configuration can be considered as an approximation for one small filament at the cloud edge which can be formed through a larger-scale entrainment. More precisely, the initial condition is a slab-like filament of supersaturated vapor that fills about one third of the cubic simulation box. The vapor mixing ratio profile across the slab is given by

$$\begin{equation} q_{\mathrm{v}}(x,y,z,t=0)=(q^{\max}_{\mathrm{v}}-q^{\mathrm{e}}_{\mathrm{v}}) \exp\left[-A(x-x_0)^6\right] + q^{\mathrm{e}}_{\mathrm{v}}. \end{equation} \tag{ 12 }$$

Here, $q^{\max }_{\mathrm {v}}$ is the maximum amplitude of q_v, which exceeds the saturation value q_vs(T₀) by 2%. The variable q^e_v stands for the environmental vapor mixing ratio, representing the subsaturated clear air outside the supersaturated cloudy air filament. Initially, droplets are seeded randomly in the supersaturated part of the profile, which is shown in figure 1, as a monodisperse particle ensemble. The simulation domain is a cube with a side length of L_x = 25.6 cm. The Kolmogorov dissipation length η will be resolved with one grid cell. All parameter studies started with the same initial conditions.

Zoom In Zoom Out Reset image size

The vapor mixing ratio profile is chosen such that the saturation value of q_vs = 0.003 56 corresponds to the typical conditions at the reference temperature T₀ = 270 K. The other turbulence parameters have been chosen corresponding to typical quantities in cumulus clouds. The mean energy dissipation is 〈ε〉 = 33.75 cm² s⁻³. Root mean square values of the velocity field are about 10 cm s⁻¹. The Schmidt number Sc = ν/D = 0.7 (ν is the kinematic viscosity of air and D the diffusivity of water vapor in air).

2.3. Numerical implementation

Kumar et al [11] investigated the phase relaxation process during a mixing event under a variety of realistic microphysical conditions for a cumulus cloud. The characteristic time scale associated with this process is the phase relaxation time, defined by

$$\begin{equation} \tau_{\mathrm{phase}} \approx \frac{1}{4\pi n_{\mathrm{d}} {\mathscr{D}} r}. \end{equation} \tag{ 13 }$$

The mixing of the clear air with the cloudy air causes a decay of the initial q_v profile. In figure 2, we plot 1D cuts of the vapor mixing ratio fluctuations at a fixed position (y₀,z₀). They are defined by

$$\begin{equation} q_{\mathrm{v}}^{\prime}({\bf x},t)=q_{\mathrm{v}}({\bf x},t)-\overline{q}_{\mathrm{v}}(t)\quad \mbox{with}\;\,\overline{q}_{\mathrm{v}}(t)=\frac{1}{V}\int_V q_{\mathrm{v}}({\bf x},t)\,\mathrm{d}V. \end{equation} \tag{ 14 }$$

Entrainment and mixing are observable by the reduction of the minimum and maximum amplitudes and the deformation of the curves. After about 2 s the mixing process is coming to the end and the fluctuations have completely decayed in the simulation domain. Although there is a source term in the advection–diffusion equation given by the condensation rate C_d, the scalar field q_v behaves practically as a conserved scalar. The fluctuations decay while the volume mean $\overline {q}_{\mathrm {v}}(t)$ approaches the saturation value q_vs for cases with a sufficiently large amount of liquid water. We always found that the impact of the droplet evaporation as a source of vapor is small within the present mixing model.

Zoom In Zoom Out Reset image size

In the studies by Kumar et al [11], the Damköhler numbers describing the mixing events changed according to the specified initial droplet number density n_d, initial droplet radius r₀ and Taylor microscale Reynolds number R_λ. Within the investigated range there were several distinct combinations of microphysical and turbulence conditions that resulted in the same value of Da, so we compare those results here in order to determine whether the Da alone is sufficient to characterize the mixing process. We refer to this as Da-similarity.

Eight sets of simulation conditions that correspond, two each, to Da_L = 0.20, 0.55, 0.81 and 1.09 are summarized in table 1. Because the Taylor microscale Reynolds number R_λ is the same in all eight simulations, the Da_L variation results solely from changes in the phase relaxation time through microphysical initial conditions. In the pair of simulations with Da_L = 0.55, for example, the initial radius and number density are varied by a factor of 2 in opposite directions, such that the phase relaxation time remains constant. Results are shown in figure 3, where we plot the size pdfs in terms of the fluctuation of the square radius

$$\begin{equation} \delta r^2 = \frac{r^2-\overline{r^2}}{\sigma_{r^2}},\, \end{equation} \tag{ 15 }$$

where σ_r² is the standard deviation of the square radius. In this way, the comparison accounts properly for the droplet-diameter-dependent growth or evaporation rates. Specifically, for a given supersaturation the single-droplet growth law equation (10) suggests that dr²/dt∝K, i.e. a constant, so when a size distribution is plotted in r² coordinates the shape is preserved under simple condensation and evaporation. By removing the distribution mean and scaling by the distribution width, we emphasize the comparison of the dispersion of the r² distribution. The distributions are plotted at t = 30 s, which is large compared to the mixing time scales so that the mixing can be taken as fully relaxed. As seen in the figure, the microphysical responses are essentially identical when Da_L is held fixed, albeit over a limited range of investigation. Nevertheless, this supports the concept that the Damköhler number captures the essential physics of the evolution of the droplet size distribution. This Damköhler similarity suggests, in turn, that we can investigate the microphysical response to mixing just by varying Da, regardless of the specific microphysical or turbulence conditions.

Zoom In Zoom Out Reset image size

Ideally, extremes of mixing would be studied directly by going to extremely large computational domains, i.e. very large R_λ, but this is greatly limited with current computational resources. In the simulations reported here, for example, the simulation domain is 25.6 cm in spatial extent. Large-eddy simulation (LES) allows larger scales, for which inhomogeneous mixing likely dominates, to be resolved, but then the discrete, Lagrangian approach to resolving the microphysics is no longer possible. It is therefore impossible to fully resolve the mixing process over a range of spatial scales that allows the transition from strongly inhomogeneous mixing (large scales) to homogeneous mixing (at the smallest scales) to be observed. As pointed out in the introduction, in our previous work [11] we achieved large-eddy Damköhler numbers of order Da_L ∼ 1. Even if the transition from inhomogeneous to homogeneous mixing is challenging to simulate with current computational abilities, can extremely homogeneous or inhomogeneous scenarios be achieved without going to greatly unrealistic microphysical conditions?

Here we investigate further extremes in Da_L within our spatial constraints by varying the transport coefficient that couples the continuous and discrete scalar fields, i.e. the thermally modified vapor diffusion coefficient $\mathscr{D}$ [11, 22] or, equivalently, the single-droplet growth constant K (see equations (10) and (11)). This directly influences Da_L through the phase relaxation time. The Damköhler number Da_L = T_L/τ_phase is therefore directly proportional to $\mathscr{D}$ or K. Da_L can also be varied by suitably changing n_d or r₀ in the simulation, as in the approach taken by Kumar et al [11], but we have chosen instead to maintain constant initial liquid water content and microphysical properties between the simulations to isolate the Damköhler effect alone. Lastly, it is important to point out that in spite of changing $\mathscr{D}$ and the associated microphysical response between the simulations, we hold the water-vapor Schmidt number (Sc = ν/D) to be constant. Thus, the fine-scale properties of mixing, including both viscous and scalar dissipation, are fixed throughout the simulations, except insofar as they are modified through direct coupling to the discrete droplet field.

The three simulations have an initial droplet radius r₀ = 20 μm, initial number density with respect to the slab cloud n^(c,0)_d = 164 cm⁻³, initial number density with respect to the whole volume n^(g,0)_d = 66 cm⁻³ and a global liquid water content w^(g,0) = 2.2 g m⁻³. To achieve stronger homogeneous and inhomogeneous mixing we have arbitrarily decreased and increased by a factor of 10 the original values of $\mathscr{D}$ and K that were taken for typical cumulus cloud conditions [11]. The values are listed in table 2, along with other relevant turbulence and cloud parameters (the middle row corresponds to the realistic cloud conditions from [11]). Simulations with identical microphysical initial conditions (e.g. the same initial droplet radius, droplet number density and therefore the same liquid water content) were then repeated using the values on the first and third rows of the table. One could consider the scenario with low $\mathscr{D}$ to be analogous to thermodynamically 'sluggish' droplets, for example due to low temperature or low mass accommodation coefficient (e.g. coated with a thin organic film), and the scenario with high $\mathscr{D}$ to be 'highly responsive' droplets, for example due to very high temperature or the presence of an unusually volatile substance. We once again emphasize that achieving high and low Da_L in this manner is still not the same as extending the range of spatial scales by simulating larger computational volumes, mainly because it still does not allow for the transition from inhomogeneous to homogeneous mixing to be observed. Nevertheless, it allows us to investigate the influence of inhomogeneous mixing on the droplet population in a Lagrangian-resolved simulation—beyond what is capable in LES, which generally can only resolve the strongly inhomogeneous portion of the mixing process [21].

Table 2. Relevant parameters for the study of Damköhler number, including the modified diffusion coefficient $\mathscr{D}$ in m² s⁻¹, the single-droplet growth constant K in μm² s⁻¹, the single-droplet evaporation time τ_r = R²₀/2K, the phase relaxation times (13) based on the number densities n^(c,0)_d and n^(g,0)_d and the large-scale eddy turnover time T_L. Four possible Damköhler numbers, based on dissipation and large-eddy time scales and on local and global phase relaxation times, are also listed. All times are given in seconds. The Kolmogorov time scale in the three DNS runs is τ_η = 0.067 s and the Taylor microscale Reynolds number is R_λ = 59.

$\mathscr{D}$	K	τr	τ(c,0)phase	τ(g,0)phase	TL	Da(c,0)η	Da(c,0)L	Da(g,0)η	Da(g,0)L
1.31×10−6	5.07	364	19	46	2.5	0.0036	0.14	0.0014	0.055
1.31×10−5	50.7	36.4	1.9	4.6	2.5	0.036	1.4	0.014	0.55
1.31×10−4	507	3.64	0.19	0.46	2.5	0.36	14	0.14	5.5

The temporal evolution of the droplet size distributions for the three mixing scenarios is shown in the left columns of figure 4, with the slowest microphysical response (homogeneous mixing) in the top left panel and fastest (inhomogeneous mixing) in the bottom left panel. The initial size distribution of the droplets is delta-distributed at 20 μm and in all three scenarios the same mass of condensed water evaporates until the full volume of air is brought to saturation. The microphysical response is starkly different, however, in the two extreme mixing scenarios, with the size distributions both rapidly achieving a distinctive shape: an extremely narrow distribution (width of approximately 0.1 μm) with a steadily decreasing mean for the strongly homogeneous mixing. For times up to 200 s, this distribution relaxes to a mean radius of approximately 19.1 μm as seen in the upper left panel of figure 4. A negatively skewed distribution for the strongly inhomogeneous mixing is found. The negative tail of the droplet diameter pdf for Da_L ≫ 1 is approximately exponential, suggesting that it could be related to the growth of the interface for a passive scalar in a turbulent flow [2, 23]. Because the droplet size distribution is inherently coupled to the supersaturation field sampled along Lagrangian droplet paths, we show the corresponding supersaturation pdfs conditioned on droplet position in the right column of figure 4. The early development of the supersaturation pdf, e.g. at times t ≲ T_L, is remarkably similar over the three simulations, always showing the development of an approximately negative exponential tail. As required for phase relaxation, the supersaturation pdf relaxes to delta-distributed at S = 0 at t ≫ τ_phase for all scenarios. The intermediate behavior, however, is quite distinct: for strongly homogeneous mixing the negative tail quickly collapses and P(S) becomes narrow and symmetric as it steadily approaches a mean of S = 0, whereas for strongly inhomogeneous mixing the negatively skewed distribution approaches a mean of S = 0 through the gradual collapse of the negative tail of P(S). It is reasonable to surmise that the distinctive, negative-exponential shape of the droplet size distribution for strongly inhomogeneous mixing is a direct reflection of the similarly skewed supersaturation pdf, and this idea is further explored in section 5.

Zoom In Zoom Out Reset image size

The size distributions in figure 4 suggest that significant changes in the shape of the distribution occur early in the mixing process, after which the shape becomes 'frozen in'. We investigate this early development by plotting the size pdfs in terms of δr², as shown in figure 5, recalling that in r²-coordinates the size distribution maintains a constant shape during uniform condensation or evaporation. The size distributions are plotted for 0.2, 0.4, 0.6, 1.0, 2.0 and 4.0 times the large-eddy time T_L (=2.5 s), and it is immediately seen that essentially all the change in the pdf shape occurs for t ⩽ T_L. The response of droplets during this transient period of mixing therefore dominates the final shape of the size distribution. This provides a relatively simple explanation for the distinctive shapes resulting from strongly homogeneous and inhomogeneous mixing: during homogeneous mixing the negative supersaturation tail collapses within time T_L and the size distribution therefore does not have sufficient time to respond significantly, whereas during inhomogeneous mixing the droplets respond rapidly to the negatively skewed supersaturation pdf and then maintain the resulting shape as the vapor field is thoroughly mixed.

Zoom In Zoom Out Reset image size

To understand further the appearance of the negative-exponential tail in the droplet size distribution, we identify three droplet groups, as shown in the lower left panel of figure 4. The groups are chosen from the left tail, center and right tail of a size distribution near the final stages of mixing. There are 3000 droplets in each of the groups, taken from among the total 1.1 million droplets. The positions of these droplets are displayed in figure 6 for times t = 0 and t = 0.2T_L, with droplets in groups 1, 2 and 3 being dark blue, light blue and green, respectively. As might have been expected, droplets from group 3 lie mostly away from the interface, within the 'core' of the unmixed cloud, while droplets from groups 1 and 2 are mostly positioned near the clear air–cloud interface. The right panel of figure 6 clearly shows that droplets in group 1, those that experience the most extreme evaporation, tend to lie at the very edge of eddies that are ejected from the initial slab cloud into the clear air during early stages of the mixing process.

Zoom In Zoom Out Reset image size

The Lagrangian behavior of individual droplets from each group is illustrated in figures 7 and 8. Five droplets are selected from each group defined in figure 4 and this is further extended to the other two simulations (droplets similarly sampled from the left tail, center and right tail of each distribution). In figure 7, each row displays radii versus time for a particular value of $\mathscr{D}$, homogeneous mixing for the top row and inhomogeneous mixing for the bottom row. Each column corresponds to a particle group, with group 1 on the left, group 2 in the center, and group 3 on the right. Figure 8 shows supersaturation sampled by the individual droplets versus time, with panels similarly arranged as in figure 7. In all radius versus time plots the phase relaxation is apparent, and as expected, the strongly inhomogeneous mixing scenario shows the most rapid response of droplet size. In that case (bottom row), full relaxation is reached by t ≈ T_L, with a nearly linear in time decrease of the droplet radius. The final radius is lowest for group 1 and highest for group 3, as can be understood from the corresponding supersaturation evolution plots, which show a progressively shallower minimum from group 1 to group 3 (figure 8). The supersaturation histories for all scenarios and particle groups are qualitatively similar, with a rapid decrease followed by a somewhat slower, asymptotic increase towards S = 0. The main dip occurs within t ≈ 2T_L and is therefore clearly associated with the transient mixing event during the first large-eddy turnover times.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The most intriguing regime in this study is that corresponding to strongly inhomogeneous mixing. As we have seen in the last sections, the droplet size distribution then exhibits a pronounced exponential tail. In the following, we give a short explanation of how the left tail of the pdf of r² follows as a consequence of the pdf of the supersaturation S along the droplet tracks. Equation (5) can be rewritten as

$$\begin{equation} \frac{\mathrm{d} r^2(t)}{\mathrm{d}t}=2\,K S(t), \end{equation} \tag{ 16 }$$

and we then obtain the general solution

$$\begin{equation} r^2(t)=r_0^2+ 2\,K \int_0^t S(t^{\prime})\,\mathrm{d}t^{\prime}. \end{equation} \tag{ 17 }$$

This requires a knowledge of the full time history of the supersaturation along the droplet track from the beginning of the mixing to time t. The solution can be simplified to

$$\begin{equation} r^2(t)\approx r^2(\overline t)+ 2\,K S (t-\overline t) , \end{equation} \tag{ 18 }$$

if the supersaturation S is roughly constant with an average $\overline {S}$ along the droplet trajectories for a certain time lag starting at $\overline t$. Figure 6 (bottom row) shows that this indeed seems to be the case for all three groups between t ≲ 1 and approximately 2.5 s, the latter of which is T_L. Thus we take $\overline t\approx 0.6\,\mathrm {s}$. The pdfs of the radius square and supersaturation can be related by the following expression:

$$\begin{equation} p(\tilde r^2,t)=\int P(S) \,\delta(\tilde r^2-r^2(S,t)) \,\mathrm{d}S , \end{equation} \tag{ 19 }$$

The substitution rule for the delta distribution states that

$$\begin{equation} \delta(\tilde r^2-r^2(S,t))=\frac{1}{|r^{2\,\prime}(\skew5\tilde{S})|}\delta(S-\skew5\tilde{S}), \end{equation} \tag{ 20 }$$

where the prime denotes the derivative of r² with respect to S taken in (18) at the value $\skew5\tilde{S}$ at which $r^2(S,t)=\tilde {r}^2$. This gives $2\!K(t-\overline {t})$ and $\skew5\tilde{S}=(\tilde r^2-r^2(\overline t))/2\,K (t-\overline t)$. For simplicity, we substitute, furthermore, $r^2(\overline t)$ by $\overline {r^2}(t)$ that we have determined for figure 5 and thus

$$\begin{equation} p(\tilde r^2,t)=\frac{1}{2\!K (t-\overline t)}P\left(\frac{\tilde r^2-\overline{r^2}(t))}{2\,K (t-\overline t)}\right). \end{equation} \tag{ 21 }$$

Supposing that to a good approximation an exponential tail of the supersaturation pdf is established, then for S < 0 we can set

$$\begin{equation} \fl P(S,t)\sim \exp\left(-\frac{|S|}{\alpha}\right) \Rightarrow p(r^2,t)\sim \frac{1}{2\!K (t-\overline t)}\exp\left(-\frac{|r^2-\overline{r^2}(t)|}{2\!K\alpha (t-\overline t)}\right). \end{equation} \tag{ 22 }$$

The parameter α will be a function of t and K. For the strongly inhomogeneous mixing regime, the tail of the supersaturation pdf can be matched with α = 0.02, 0.023 and 0.018 for times t = 1, 1.5 and 2 s. Figure 9 demonstrates that the supersaturation pdf with the exponential tail indeed results in a size distribution with a shape that can be matched reasonably well to that prescribed by (22).

Zoom In Zoom Out Reset image size

It remains to be shown why the entrainment process generates exponential tails of the supersaturation pdf. As was seen in figure 4 (right panels), however, its appearance during the transient mixing for t ≲ T_L is largely independent of details of the particle response. We draw attention, however, to work done in the context of the so-called Kraichnan model, a passive scalar advection model in a synthetic turbulent flow without time correlations [24] (see [2, 25] for comprehensive reviews). For both spatially rough and smooth advecting flows, it can be shown within this model that exponential tails of the passive scalar evolve. Numerical and experimental studies on Navier–Stokes turbulence can give rise to super-Gaussian or exponential tails as well as to sub-Gaussian tails as discussed in [26]. This difference might be attributed to the finite-time correlations that every real turbulent flow obeys, a point that is still open.

The idealized numerical experiments presented here aimed at exploring extreme cases of mixing and entrainment at the edge of a cloud within a simplified numerical model that disentangles mixing dynamics from the strong temperature feedback via the latent heat release as would be present during the evolution in a cloud. The regimes of strongly homogeneous and inhomogeneous mixing were established here by a manipulation of the droplet response to local super- or subsaturation which is established by a variation of the constant K (or ${\mathscr{D}}$) in growth equation (6) for the droplet radius. The strongly inhomogeneous mixing case is characterized by pronounced exponential tails of the pdfs of supersaturation and radius which can be quantitatively related to each other.

As we stated at the beginning of this work, our present study should be considered only as a first step of a longer-term analysis on the complex subject of turbulent mixing in clouds. There is no doubt that it remains to be seen which aspects will be observable when the feedback of the temperature field is included and the numerical experiments are scaled up to larger Reynolds numbers and box sizes which automatically lead to larger Damköhler numbers Da_L (e.g. fully resolved cloud volumes of the order of 10 m³). Studies that incorporate the active role of temperature on the flow and the latent heat release in the mixing dynamics are currently under way and will be reported in the near future. However, the present analysis already provides what in our view are interesting implications for models in which groups of 10⁵–10⁶ individual cloud droplets are represented by one quasi-particle, commonly referred to as a super droplet [21, 27]. In those models that are usually coupled to LES the present studies can provide new strategies to model the fluctuations of S on the subgrid scales. Regardless of the mode of implementation, ultimately it will be important to couple the effects of local mixing and dilution, such as those considered here, with the fluctuations in humidity that result from larger-scale motions within the cloud. Recent efforts to perform detailed microphysical calculations for ensembles of Lagrangian parcel trajectories (e.g., [28, 29]) either neglect mixing between parcels or represent the parcel-scale mixing in simple ways that do not fully capture the Damköhler-number dependence. Similar statements could be made regarding recent developments on stochastic condensation theory (e.g. [30, 31]). These approaches give very important insight into the crucial role of large-scale mixing and dynamics in the evolution of droplet size distributions, and it remains to be seen how an accurate representation of the microphysical response to entrainment and dilution will alter the emerging picture.

The authors are grateful to Alan Kerstein, Steven Krueger, Alexander Kostinski and Holger Siebert for helpful discussions. The authors acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) within the Research Focus Program Metström (SPP 1276). JS acknowledges additional support from the Heisenberg Program under grant no. SCHU 1410/5-1. Furthermore, support from COST Action MP0806 is kindly acknowledged. The numerical simulations have been carried out at the Jülich Supercomputing Centre (Germany) on the computer cluster JUROPA under grant no. HIL03. RAS acknowledges support from the National Science Foundation under grant no. AGS1026123.