Key parameters for droplet evaporation and mixing at the cloud edge

The distribution of liquid water in ice‐free clouds determines their radiative properties, a significant source of uncertainty in weather and climate models. Evaporation and turbulent mixing cause a cloud to display large variations in droplet number density, but quite small variations in droplet size (Beals et al., Science, 2015, vol. 350, pp. 87–90). However, direct numerical simulations of the joint effect of evaporation and mixing near the cloud edge predict quite different behaviours, and how to reconcile these results with the experimental findings remains an open question. To infer the history of mixing and evaporation from observational snapshots of droplets in clouds is challenging, because clouds are transient systems. We formulated a statistical model that provides a reliable description of the evaporation–mixing process as seen in direct numerical simulations and allows us to infer important aspects of the history of observed droplet populations, highlighting the key mechanisms at work and explaining the differences between observations and simulations.


INTRODUCTION
Clouds play a major role in regulating weather and climate on Earth, by modulating the incoming solar radiation. Despite substantial scientific advances in the last decades, clouds still represent the primary source of uncertainty in climate projections (Stocker et al., 2013;Pincus et al., 2018). A key challenge is to understand how entrainment of dry air at the edges of ice-free clouds affects the size distribution and number density of droplets (Blyth, 1993). This is important because the amount and distribution of liquid water determine cloud optical properties (Kokhanovsky, 2004) and precipitation efficiency (Burnet and Brenguier, 2007). As a consequence, weather and climate models are sensitive to how entrainment at the cloud edge is parameterised (Mauritsen et al., 2012).
The optical properties of clouds are of crucial importance for the radiation balance of the Earth's climate system (Dufresne and Bony, 2008;Caldwell et al., 2016). The size distribution and number density of droplets are key ingredients, because the light-extinction coefficient of the cloud is determined by the number of the droplets it contains times their average surface area (Kokhanovsky, 2004). F I G U R E 1 (a) Effect of mixing on local droplet populations (schematic). On large spatial scales ( 1 and 2 ), mixing and evaporation is not yet complete, but some smaller regions (of size 3 ) are in local steady states: saturated air with droplets, or subsaturated air without droplets. (b) Initial cloud configuration in our model. Before mixing, moist air and droplets reside in a w × L × L slab, contained in a cubic domain of side length L. Regions with dry air are dashed. The solid line is the initial profile of supersaturation s (see text).
Regarding precipitation efficiency, the mechanism of rain formation in ice-free clouds is a longstanding unresolved problem in atmospheric physics (Grabowski and Wang, 2013). A broad initial droplet size distribution is needed to activate the collisions and coalescences of droplets necessary for the rapid onset of rain formation observed empirically in warm clouds (Szumowski et al., 1997;1998;Devenish et al., 2012). Microscopic droplets grow by condensation of water vapour, or shrink by evaporation. However, when a droplet-containing parcel does not mix with its surroundings, condensation causes the droplet size distribution to narrow because the diffusional growth of a droplet is inversely proportional to its radius, so that small droplets grow faster than large ones (Rogers and Yau, 1989).
Turbulence has a strong influence upon droplet condensation and evaporation in clouds (Bodenschatz et al., 2010). Turbulent mixing causes water-vapour and liquid-water content to fluctuate on different length-and time-scales (Vaillancourt et al., 2001). As a consequence, nearby droplets may have experienced quite different growth histories. The droplets of a cloudy parcel that is mixed with dry air evaporate at different rates, so that the droplet size distribution broadens (Lanotte et al., 2009;Sardina et al., 2015;Li et al., 2020). Cloud-resolving simulations (Hoffmann and Feingold, 2019) show that this mechanism can have a strong effect on droplet number densities in turbulent clouds.
Entrainment of dry air at the cloud edge triggers rapid changes in the droplet size distribution (Perrin and Jonker, 2015;Abade et al., 2018). As turbulence mixes dry air into the cloud, it creates long-lasting regions of dry air where droplets can evaporate rapidly. Droplets in regions with higher water-vapour concentration, by contrast, may saturate the air and survive for a much longer time (Figure 1a). While droplets evaporate, turbulence mixes the cloud at many length-scales, ranging from the Kolomogorov length-of the order of millimetres (Devenish et al., 2012)-to a few kilometres (Rogers and Yau, 1989). Evaporation and mixing on a spatial scale depend on the turbulent mixing time at that scale, and upon the relevant thermodynamic time-scale . Their ratio forms a Damköhler number Da = ∕ (Dimotakis, 2005). The thermodynamic process parameterised by Da is limited by the rate of mixing if Da is large, and limited by thermodynamics if Da is small. The dynamics at large Damköhler numbers is referred to as inhomogeneous mixing (Baker et al., 1980), where some droplets evaporate completely while others do not evaporate at all. Small-Da mixing is called homogeneous, where droplets evaporate at approximately the same rate, so that the droplet size distribution remains narrow.
The notion of homogeneous and inhomogeneous mixing remains debated (Tölle and Krueger, 2014), but it can be given a precise meaning in terms of the fraction P e (t) of droplets that have evaporated completely at time t. However, it is not understood at present which mechanisms and parameters determine the transition from homogeneous to inhomogeneous mixing. Several authors have attempted to describe turbulent mixing in terms of one Damköhler number. Lehmann et al. (2009) used a combined microphysical response time r , a function of the two thermodynamic time-scales of the problem, d (droplet evaporation), and s (supersaturation relaxation). Lu et al. (2018), by contrast, suggest that d should be used to formulate a single-parameter criterion for inhomogeneous mixing. The direct numerical simulations (DNS) by Kumar et al. (2018) explored how the nature of mixing changes as L ∕ r increases with the linear size L of the simulated domain. However, no clear sign of inhomogeneous mixing was found. The authors mention that this may be a consequence of the thermodynamic setup used. Another possibility is that the simulated system was just not large enough. Jeffery (2007) emphasised that evaporation and mixing cannot be described by a single Da alone, because there are two thermodynamic time-scales, leading to three key nondimensional parameters: Da d , Da s , and the volume fraction of cloudy air. However, Jeffery only studied the case Da d = Da s and did not discuss the implications of varying the Damköhler numbers separately. Pinsky et al. (2016a) and Pinsky and Khain (2018a) modified an equation for advected liquid water (Rogers and Yau, 1989) into a diffusion-reaction equation for the droplet size distribution and emphasised the significance of a nondimensional parameter, their potential evaporation parameter.
Here we derive a statistical model for evaporation and turbulent mixing at the cloud edge from first principles. The model takes into account the multiscale turbulent dynamics, as turbulent clouds can have large Reynolds numbers; Re ∼ 10 7 is a conservative estimate for convective clouds (Devenish et al., 2012). The model quantitatively predicts the outcomes of the DNS by Kumar et al. (2012;. Furthermore, the model shows, in accordance with Jeffery (2007), that the evolution of the droplet size distribution is determined by Da d , Da s , and . We find that the potential evaporation parameter of Pinsky et al. (2016a) and Pinsky and Khain (2018a) is determined simply by the ratio of the Damköhler numbers. Finally, the model allows us to interpret the results of in situ measurements of clouds at the centimetre scale (Beals et al., 2015).

METHOD
We study mixing and evaporation of cloud droplets by mixing moist air, droplets, and dry air in a cubic domain of side length L with periodic boundary conditions. Initially, the saturated or slightly supersaturated moist air with supersaturation s c ≥ 0 is contained in a w × L × L slab, together with N 0 randomly distributed water droplets ( Figure 1b). The dry air, initially outside the slab, has negative supersaturation s e < 0. The mixing is driven by statistically stationary homogeneous isotropic turbulence, with turbulent kinetic energy TKE and mean dissipation rate per unit mass (Frisch, 1995). Essentially the same setup is used in the DNS of Kumar et al. (2012;, which allows us to understand their simulation results in terms of our model.

Microscopic equations
For the turbulent mixing, we start from the microscopic equations of Kumar et al. (2018) and earlier studies (Vaillancourt et al., 2002;Lanotte et al., 2009;Paoli and Shariff, 2009;Kumar et al., 2014;Perrin and Jonker, 2015). We neglect buoyancy, particle inertia and settling, temperature changes due to vertical motion, and temperature and pressure dependences of the thermodynamic coefficients, and subsume the joint effects of temperature and water vapour into a single supersaturation field. We denote fluid velocity and pressure by u(x, t) and p(x, t), and supersaturation by s(x, t). The spatial position of droplet is x (t), its radius equals r (t), and the index ranges from 1 to N 0 . We nondimensionalise as follows: u ′ = u∕U, x ′ = x∕(U L ), t ′ = t∕ L , p ′ = p/( a U 2 ), x ′ = x ∕(U L ), s ′ = s/|s e |, r ′ = r ∕r 0 , and s ′ c = s c ∕|s e |. Here U = √ 2 TKE∕3 is the turbulent root-mean-square velocity, L = TKE∕ is the large-eddy time [proportional to L/U if the size of the largest eddies is of the order L (Pope, 2000)], a is the reference mass density of air, is the initial volume-averaged droplet radius, and |s e | is the (positive) subsaturation of air outside the initial cloud slab. Dropping the primes, the microscopic equations take the nondimensional form: Equations 1-4 are the incompressible Navier-Stokes equations, with Lagrangian time derivative D Dt = t + (u ⋅ ∇). In DNS of Equations 1-4, a forcing is imposed to sustain stationary turbulence. This is not necessary in the model introduced below, and we therefore do not include a forcing term in Equation 1. Equation 2 is the equation for supersaturation. The first term on its right-hand side describes diffusion of the supersaturation s(x, t). The second term models the effect of condensation and evaporation through the average r (t)s(x , t), taken over all droplets in the vicinity of x. Droplets are advected by the turbulent flow (Equation 3), and Equation 4 models how the droplet radius r changes due to evaporation and condensation. When a droplet has evaporated completely, we impose that it must remain at r = 0. In the derivation of Equation 4, s(x (t), t) enters as the supersaturation at distances from the droplet much larger than the droplet radius (Rogers and Yau, 1989). In other words, Equation 4 relies on a scale separation between droplet sizes and the lengths that characterise supersaturation fluctuations generated by turbulent mixing (Vaillancourt et al., 2001). As a consequence, droplets interact locally with the supersaturation field over finite volumes through the average r (t)s(x , t) in Equation 2. Further details regarding Equations 1-4 are given in the Supporting Information, where we also show how to derive Equations 1-4 from the more detailed dynamical description of Vaillancourt et al. (2001;2002), Kumar et al. (2014;, and Perrin and Jonker (2015).
An advantage of writing the dynamics in nondimensional form is that this determines the independent nondimensional parameters. First, Re L = 2 3 TKE 2 ∕( ) is the turbulence Reynolds number (Pope, 2000), where is the kinematic viscosity of air. The Schmidt number is defined as Sc = ∕ , where is the diffusivity of s. The volume fraction of cloudy air is given by = w∕L, and V = [L∕(U L )] 3 is the dimensionless domain volume. The Damköhler number Da s is defined as Da s = L ∕ s , where s = (4 A 2 A 3 w n 0 r 0 ) −1 is the supersaturation relaxation time. This is the time-scale at which the supersaturation decays towards saturation, assuming that all droplets have the same radius r 0 , and for droplet number density n 0 = N 0 /(wL 2 ). Further, w is the density of pure liquid water, and A 2 and A 3 are thermodynamic coefficients, specified in the Supporting Information. The Damköhler number Da d is defined as is the droplet evaporation time, the time that it takes for a droplet of radius r 0 to evaporate completely in a constant ambient supersaturation s e < 0.
The Damköhler numbers determine the extent to which saturation and droplet evaporation are limited by the rate of mixing (Dimotakis, 2005). Saturation is mixing-limited at large Da s , since regions with evaporating droplets-created by mixing of cloudy and dry air at the time-scale L -saturate faster than L . When Da s is small, by contrast, evaporating droplets saturate the air more slowly than it is mixed. In this case, saturation is not limited by the rate of mixing. Droplet evaporation is mixing-limited at large Da d , since droplets then evaporate more rapidly than the exposure to subsaturated air changes. At small Da d , mixing is faster than droplet evaporation, and droplets tend to evaporate mainly after the system has been mixed. The droplets then experience roughly the same supersaturation as they evaporate.
In the limit Re L → ∞, three key nondimensional parameters remain in Equations 1-4: , Da d , and Da s . The system can be parameterised by , Da d , and In this way, the scale dependence of the mixing process is contained in Da d only. The Damköhler-number ratio ℛ is inversely proportional to the density of liquid water in the cloud slab; it regulates the moisture of the mixing process (details in the Supporting Information). The bifurcation between moist steady states, where droplets remain in saturated air, and dry steady states, where all droplets have evaporated completely (Jeffery, 2007;Kumar et al., 2013;Pinsky et al., 2016b), occurs at a critical value of ℛ, ℛ c . The critical ratio ℛ c can be computed from the conserved quantity = −⟨s(t)⟩ − 2 3ℛ [1 − P e (t)]⟨r 3 (t)⟩, which is analogous to the liquid-water potential temperature at fixed altitude (Gerber et al., 2008;Lamb and Verlinde, 2011;Kumar et al., 2014). Here, P e (t) is the fraction of completely evaporated droplets, the fraction of droplets for which r (t) = 0 at time t. Furthermore, is the mean cubed droplet radius, conditioned on r (t) > 0 by the factor [1 − P e (t)] −1 . The conservation of can be concluded by integrating Equation 2 for supersaturation over the domain volume: see the Supporting Information for details. Moist steady states have ⟨ r 3 (t) ⟩ > 0 and ⟨s(t)⟩ = 0, dry steady states have P e (t) = 1 and ⟨s(t)⟩ < 0, so the sign of determines whether the steady state is moist or dry. The value of is determined by the initial conditions, = −⟨s(0)⟩ − 2 ∕(3ℛ). Setting = 0, we find the critical Damköhler-number ratio ℛ c = − 2 3 ∕⟨s(0)⟩. DNS of Equations 1-4 for experimentally observed dissipation rates and droplet number densities are feasible only for quite small systems (Kumar et al., 2018). This restricts the range of scales that can be explored and makes it difficult to detect inhomogeneous mixing in DNS. We therefore pursue an alternative approach and adapt a PDF model (Pope, 2000)-commonly used to describe combustion processes (Haworth, 2010)-to the inhomogeneous cloud edge. As opposed to the kinematic statistical models reviewed by Gustavsson and Mehlig (2016), we must also take thermodynamic processes into account.

Statistical model
Statistical models have been used to describe droplets in a supersaturation field that fluctuates around zero, as in the cloud core (Paoli and Shariff, 2009;Sardina et al., 2015;Chandrakar et al., 2016;Siewert et al., 2017). At the cloud edge, there are large deviations from this equilibrium. Jeffery (2007), Pinsky et al. (2016a), and Pinsky and Khain (2018a) formulated models for the cloud edge where droplets evaporate in direct response to a mean supersaturation field. This does not take into account that mixing is local, and that small droplets are advected together with the supersaturation field. For an accurate description of mixing and evaporation, it is essential to describe how each droplet carries its own local supersaturation (Siewert et al., 2017). Our model does just that. It is derived from first principles using the established framework of PDF models (Pope, 2000). For the configuration shown in Figure 1b, we derive one-dimensional statistical-model equations from Equations 1-4 (details in the Supporting Information): Equation 6 describes the fluctuating acceleration of Lagrangian fluid elements in turbulence. Here, d are Brownian increments with zero mean and variance dt, and C 0 is an empirical constant (Pope, 2011). Each Lagrangian fluid element has a supersaturation s, and may contain a droplet (at position x, of size r). Equation 7 approximates the supersaturation dynamics as decay towards ⟨s(x, t)⟩, regulated by the empirical constant C (Pope, 2000). The second term on the right-hand side of Equation 7 represents the effect of condensation and evaporation through ⟨r(t)s(x, t)⟩. The position-dependent averages ⟨ … ⟩ in Equation 7 are taken over fluid elements located at x at time t (details in the Supporting Information). The statistical model in Equations 6-9 becomes independent of Re L at large Reynolds numbers, where C 0 approaches a definite limit (Pope, 2011). It is independent of Sc, in accordance with the known behaviour of advected scalars in fully developed turbulence (Shraiman and Siggia, 2000).
Equations 6-9 rest upon a probabilistic description of the dynamics of the two phases, droplets and air (Pope, 2000;Jenny et al., 2012). The corresponding evolution equations, dictated by Equations 1-4, contain unclosed terms that must be approximated (Pope, 1985;Haworth, 2010) to obtain a closed model such as in Equations 6-9. Following Pope (2000), we cast the model into the form of stochastic dynamical equations for Lagrangian fluid elements. Since the dynamics is statistically one-dimensional in our configuration, we can average over the y and z coordinates to obtain Equations 6-9 in one-dimensional form. For the closures, we rely on standard approximations, common and justified in PDF modelling of singlephase flows (Pope, 1985) and in models for turbulent combustion (Haworth, 2010;Jenny et al., 2012;Stöllinger et al., 2013). The explicit mathematical approximations for the closures provided in the Supporting Information render the interpretation of the statistical model definite, and indicate how to improve the model when necessary.
In the following, we briefly summarise the closures. Equation 6 contains the closure for fluid-element accelerations. It reproduces the empirically observed effect of turbulent diffusion of passive-scalar averages (Pope, 2000). Equation 7 contains two closure approximations. First, the decay towards ⟨s(x, t)⟩ approximates the diffusion of supersaturation. This closure ensures that the mean of a passive scalar is conserved, and that a passive scalar remains bounded between its minimal and maximal values (Pope, 2000). Furthermore, it describes the decay of passive-scalar variance in statistically homogeneous turbulent mixing of two scalar concentrations (Pope, 2000). However, this closure does not reproduce the relaxation of the single-point PDF of scalar concentration from a two-peaked distribution via a U-shaped distribution into a Gaussian (Eswaran and Pope, 1988;Pope, 1991). For our initial conditions (Figure 1b), the decay towards ⟨s(x, t)⟩ captures the supersaturation fluctuations experienced by a fluid element as it moves towards or away from the most cloudy region. Note that this closure does not account for large saturated cloud structures that tend to relax slowly towards ⟨s(x, t)⟩. Such events are most relevant during the initial stages of the mixing-evaporation process, and their effect is expected to diminish with time, as large cloud structures are mixed into smaller and smaller structures. Consequently, the statistical model may describe short initial transients only qualitatively, not quantitatively.
The second closure in Equation 7 approximates the effects of droplet phase change on the supersaturation field: the local average r (t)s(x , t) in Equation 2 is replaced by the ensemble average ⟨r(t)s(x, t)⟩. This is the simplest closure that preserves the conservation of the parameter , and it is therefore common in PDF models that describe combustion of particles in turbulence (Haworth, 2010;Jenny et al., 2012;Stöllinger et al., 2013). The average ⟨r(t)s(x, t)⟩ takes into account that droplet evaporation is delayed locally when nearby droplets saturate the surrounding air. Since we obtain closure by replacing the local average r (t)s(x , t) by an average over fluid elements with one-dimensional dynamics, variations in the rate of supersaturation relaxation in the y-and z-coordinates are not described.
It is expected that the large cloud structures mentioned above, and their three-dimensional forms, matter more at very large Damköhler numbers. Therefore it cannot be excluded that the statistical model is only qualitative in this extreme limit. Below we show that the model works very well even for the largest Damköhler numbers in DNS studies (Kumar et al., 2012;. Also, since the statistical model is derived using the established framework of PDF models (Pope, 2000), it can be improved straightforwardly by incorporating additional variables in the probabilistic description (Pope and Chen, 1990;Pope, 2000;Meyer and Jenny, 2008), or by using more refined approximations (Pope, 1991;Jenny et al., 2012).

Comparison with DNS
The statistical model can be used to understand the DNS results of Kumar et al. (2012;. Figure 2 shows good agreement for the time evolution of the fraction P e (t) of droplets that have completely evaporated, even though the statistical-model dynamics is slightly slower. Panels (a) and (b) in Figure 3 shows that the model reproduces the broadening of the droplet size distribution. The slightly slower dynamics in Figure 2 and the deviations in the tails in Figure 3 suggest that the statistical model does not reproduce the most rapid evaporation rates. This could be due to turbulent fluctuations in the supersaturation diffusion, neglected in Equation 7. Kumar et al. (2012; compute droplet size distributions with prominent exponential tails using DNS-some of them are seen in Figure 3 (black lines)-and connect these tails to corresponding exponential tails in the PDF of supersaturation at droplet positions. In our statistical-model simulations, we observe heavy tails in the PDF of supersaturation at droplet positions, but the tails are less pronounced than in the DNS (not shown). Heavy tails are consistent with the results of Eswaran and Pope (1988) mentioned above, who observed how an initially bimodal supersaturation relaxes. Despite these shortcomings, our model describes the time evolution of P e (t) very well ( Figure 2). It is also a significant improvement over models in which the droplets interact with a mean supersaturation field (Jeffery, 2007; Pinsky et al., 2016a;Pinsky and Khain, 2018a). In reality, the droplets react to the local supersaturation, as mentioned above, and this may be particularly important at large values of Da d , where locally saturated regions can persist for a long time. Figure 4 shows the steady-state value P * e of P e (t) computed from the statistical model as a function of Da d and ℛ∕ℛ c . We see how P * e increases with both Da d and ℛ∕ℛ c . The DNS results of Andrejczuk et al. (2006) and Kumar et al. (2012; form a pattern in Figure 4 that verifies these dependences: open symbols correspond to DNS with little or no complete evaporation in the steady state (P * e < 10%), and filled symbols to P * e > 10%. Figure 4 also explains why the DNS of Kumar et al. (2018) did not exhibit significant levels of inhomogeneous mixing: since their ℛ was quite small, small values of P * e require values of Da d much larger than unity (Da d ∼ 10 2 for P * e = 10%). Furthermore, the substantially different outcomes of the DNS of Kumar et al. (2012; are explained: their parameters lie on opposite sides of the bifurcation line. Figure 4 also explains, at least qualitatively, numerical results of DNS of transient turbulence with quite different initial conditions (Andrejczuk et al., 2006), namely how the amount of complete droplet evaporation increases with both ℛ∕ℛ c and Da d . There is no parameter corresponding directly to Da d in the simulations of Andrejczuk et al. (2006), because they are for different initial conditions and flows. We therefore place these simulations in Figure 4 by computing a time-scale ratio that, in a qualitative sense, incorporates the same physics as Da d (details in the Supporting Information). Key parameters of the DNS in Figure 4 are summarised in Table 2; a complete description is provided in the Supporting Information.

Mixing histories from observations
A common way of characterising the droplet content of a cloud is to plot the mean cubed radius r 3 and number density n of droplets for observed cloud-droplet populations in a mixing diagram. Figure 5(a) shows a mixing diagram with empirical data from Beals et al. (2015). Black crosses are values of r 3 and n extracted from snapshots (linear size 15 cm) of local droplet populations measured during an airplane flight through a convective cloud. Observational data in mixing diagrams are commonly discussed in relation to the homogeneous mixing line, a curve of the global steady states (r 3 * , n * ) that result from homogeneous mixing (no complete evaporation) between different proportions of undiluted cloudy and dry environmental air (Gerber et al., 2008; Table 1) for different times. The probability density of the nondimensional droplet radius r is shown for the statistical model (coloured lines) and DNS of Kumar et al. (2014) (black lines). The initial droplet size distribution is monodisperse and centred at r = 1. (c) Same as (b), but for a very moist case (parameters in Table 1) and DNS of Kumar et al. (2012). shown in Figure 5a. A fundamental problem, however, is that it is not clear how to interpret mixing diagrams such as Figure 5a, since it is not clear that the empirically observed droplet populations reflect global steady states (Pinsky and Khain, 2018a). It is nevertheless likely that most data points in Figure 5a sample local steady states, that is, locally well-mixed droplet populations that reside in saturated air. To describe such droplet populations, one must refer to the multiscale turbulent mixing process. We attempted this analysis using the statistical model, assuming that the statistical model with the initial condition shown Figure 1b describes how a cloud structure at the spatial scale L develops. Under this assumption, n * and r 3 * are given by the droplet number density (normalised by n 0 ) and the mean cubed droplet radius ⟨ r(t) 3 ⟩ in the steady state, and we can conjecture the mixing histories that formed the measured droplet populations.

TA B L E 1
We begin by noting that and P * e are completely determined for any steady-state point (r 3 * , n * ) in a mixing diagram. To show this, we write the volume-averaged initial supersaturation as ⟨s(0)⟩ = (1 + s c )( + 0 ) − 1, where 0 is a constant that depends on the initial supersaturation profile (details in the Supporting Information). Inserting Equations 10 and 11 determine how to map (r 3 * , n * ) to ( , P * e ). As a consistency check we note that one obtains the homogeneous mixing line (Pinsky et al., 2016a) from Equation 11 by setting P * e = 0. This allows us to infer that F I G U R E 4 Steady-state fraction P * e of completely evaporated droplets as a function of Da d and ℛ∕ℛ c ; details in the Supporting Information. The fraction P * e is colour-coded. The solid red line is the contour P * e = 10%, the black dashed line indicates the transition between moist and dry steady states, and symbols indicate DNS results from previous studies: ▵ (Kumar et al., 2012), ◽ (Kumar et al., 2013), ⬦ (Kumar et al., 2014), ▿ (Kumar et al., 2018), and • (Andrejczuk et al., 2006). Filled symbols indicate P * e > 10%. The DNS of Andrejczuk et al., (2006) should not be compared quantitatively with the statistical model results, since they are for different initial conditions and decaying turbulence (see text and Supporting Information). Red, blue, and light blue circles correspond to the dry, moist, and very moist simulations in Figures 2 and 3. ℛ = 0.17 and s c = 0 = 0 for the homogeneous mixing line of Beals et al. (2015).
Any point in the mixing diagram must correspond to a local steady state with certain values of P * e and . Each statistical-model simulation for given Da d , ℛ, and yields a certain value of P * e . This allows us to extract a value of Da d for each point in the mixing diagram from our statistical-model simulations. The result is shown in Figure 5a. We see that Da d increases rapidly above the homogeneous mixing line. Estimating s ∼ 1s from Beals et al. (2015) and conservatively estimating ∼ 1 cm 2 ⋅s −3 for a convective cloud, a value of Da d = 1000 implies that an observed droplet population was mixed at spatial scales of the order of L ∼ √ 3 L ∼ 5 km, larger than the size of the cloud. In other words, the rapid increase of Da d in Figure 5a suggests that most of the data in the mixing diagram cannot be in a global steady state of a mixing process parameterised by ℛ = 0.17.
We concluded above that most measurements of Beals et al. (2015) are likely to correspond to local steady states. As undiluted cloudy air is mixed with premixed air, such steady states are formed locally and temporarily as local mixing processes equilibrate at small spatial scales (Figure 1a). We now discuss how the analysis of local steady states may yield insight into possible local histories of the cloud. Air affected by earlier mixing events is not as dry as environmental air, so the mixing of undiluted cloud with premixed air is governed by smaller values of ℛ. We therefore ask: which values of ℛ are consistent with the assumption that the experimentally observed droplet population in the middle panel of figure 2 of Beals et al. (2015)-the red cross in Figure 5a-reflects a local steady state? Our model allows us to determine possible combinations of Da d , ℛ, and consistent with a local steady state. We know that ℛ must be smaller than 0.17, the upper limit dictated by the homogeneous mixing line of Beals et al. (2015). Furthermore, since the data cannot lie below the homogeneous mixing line of the global mixing process, a lower bound for ℛ is ℛ min = 2 3 (n * − r 3 * )n * ∕

[
(1 + s c )( 0 + n * ) − 1 ] = 0.0236. Figure 5b shows values of ℛ and Da d obtained from our statistical-model simulations that are consistent with these constraints. We see that the range of possible values of Da d covers several orders of magnitude. This means that local mixing processes consistent with the red cross in Figure 5a may have occurred over a large range of spatial scales. We also see that ℛ does not vary much in Figure 5b, only between 0.024 and 0.03. This allows us to conclude that some important aspects of the mixing dynamics are essentially independent of spatial scale. First, the fact that ℛ is substantially smaller than 0.17 indicates that the noncloudy air was premixed. Second, using Equation 11, we find that the reduction in droplet number density was caused primarily by dilution even at the largest scales, since P * e increases only up ∼1.4% for the largest value of Da d , at Da d ∼ 1000 and ℛ = 0.03. Put differently, ∼ n * for all values of Da d we considered.
How does the outcome of a local mixing process depend on its scale? Larger scales correspond to larger values of Da d , and Figure 5b shows that complete droplet evaporation begins to occur around Da d = 1, where ℛ starts to exceed ℛ min (blue circle). Estimating ∼ 10 cm 2 ⋅s −3 , a typical value for convective clouds (Devenish et al., 2012), we find that Da d = 1, s ∼ 1, and ℛ= ℛ min correspond to the spatial scale 9 m. Mixing processes leading to the red cross in Figure 5a that occurred at scales smaller than 9 m were therefore perfectly homogeneous: none of the droplets evaporated completely, as they were diluted by premixed air. At larger spatial scales, small but nonzero fractions of the droplets evaporated completely. Equation 11 gives ℛ = 0.028 for P * e = 1%, and from Figure 5b we read off Da d = 13 (green circle). For ∼ 10 cm 2 ⋅s −3 these values correspond to 300 m. This suggests that reductions in droplet number density are also  dominated by dilution, and not complete droplet evaporation, for mixing processes that range over hundreds of metres. Furthermore, since most data points in Figure 5a reside well above the region where equilibria are found for ℛ = 0.17, we conclude that they too resulted from mixing with premixed air.

DISCUSSION
A general conclusion from our analysis is that both Damköhler numbers are important for the transition to inhomogeneous mixing; Da d parameterises the mixing-limited nature of droplet evaporation, and the ratio ℛ = Da d ∕Da s regulates the self-limiting effect of droplet evaporation, namely that droplets cease to evaporate when they have saturated the surrounding air or evaporated completely. Analysing the parameters of our microscopic equations, Equations 1-4, we see that ℛ = − 3 2 R, where R is the potential evaporation parameter of Pinsky et al. (2016a) and Pinsky and Khain (2018a). Therefore R is in fact given by the ratio of Da d and Da s , consistent with our conclusion that both Damköhler numbers matter. Pinsky and Khain (2018b; concluded that the Damköhler-number ratio ℛ determines whether a cloud expands by dilution or shrinks by complete droplet evaporation. A mixing process that mixes equal proportions of saturated cloudy and subsaturated noncloudy air has ℛ c = 2 3 , so the symmetric configuration they adopted for the cloud edge implies that the cloud expands if ℛ < 2 3 and shrinks otherwise. We note that whether the cloud expands or shrinks depends on the position and scale at which one perceives it. A local mixing process with small ℛ∕ℛ c tends towards a moist steady state, so the cloud dilutes locally. A local mixing process that tends towards a dry steady state, by contrast, consumes the cloud. It is possible for a cloud to expand locally for some time, even if this local expansion is part of a mixing process that consumes the cloud at larger length-and time-scales. Although global mixing processes are transient, they contain local steady states. Diluted and saturated local droplet populations (such as the red cross in Figure 5a) can be a part of such transients. However, such local steady states must eventually be abandoned as mixing proceeds globally.
Adopting this multiscale picture of mixing, in which a large-scale mixing process consists of many mixing processes at smaller scales, it is natural to expect that large ranges of the parameters , ℛ, and Da d are relevant. If one moves the domain of a local mixing process from the interior of the cloud towards the cloud edge, the liquid water content decreases, so that ℛ∕ℛ c and Da d increase. Lehmann et al. (2009) point out that Da d increases as one perceives mixing processes at larger and larger spatial scales. This corresponds to moving to the right in Figure 4. We note that ℛ∕ℛ c and Da d tend to increase with distance from the interior of the cloud. Moving the sampling volume towards the cloud edge then corresponds to a motion upwards and to the right in Figure 4. The amount of complete droplet evaporation increases in this direction, consistent with the fact that complete droplet evaporation takes place at the cloud edge.
A number of assumptions may influence our interpretation of the empirical data in Section 3.2. First, the model configuration in Figure 1b is simplified compared with real clouds, which have irregular shapes that deform during the mixing-evaporation process. Second, the observational method may not detect droplets with radii smaller than 3 m (r 3 = 0.2), as stated by Beals et al. (2015). If many small droplets were not detected, the observations in Figure 5a are located too far from the homogeneous mixing line. Third, at the upper end of the Damköhler range in Figure 5, the statistical model may not be quantitatively accurate, as stated above. We nevertheless expect that the statistical model reproduces the evolution of droplet size distributions in DNS qualitatively in Figure 5. This expectation is corroborated by the robust tendency for P * e to increase with increasing values of Da d and ℛ in Figure 4, and follows directly from the roles of the Damköhler numbers in mixing-evaporation dynamics.
The deviations in the tails of droplet size distributions at moderate Damköhler numbers in Figure 3 suggest that the next step in improving the statistical model should aim at reproducing the fastest evaporation rates in the transient mixing-evaporation process. A better agreement in the tails could be achieved by refining the closure for supersaturation diffusion by using a dynamic C in Equation 7 (Jenny et al., 2012), or by introducing additional fluctuations (Pope, 1991). Another possibility is to refine the description of the spatial structure of the supersaturation field by improved closures (Pope, 1991;Vedula et al., 2001;Meyer and Jenny, 2008;Jenny et al., 2012).

CONCLUSIONS
We derived a statistical model for evaporation and turbulent mixing at the cloud edge from first principles. The model explains results of earlier DNS studies of mixing (Andrejczuk et al., 2006;Kumar et al., 2012; and shows that two thermodynamic time-scales are important for a mixing process, the droplet evaporation time and the supersaturation relaxation time. This means that one must consider two Damköhler numbers in order to quantify the mixing-evaporation dynamics. We concluded that the simulations of Kumar et al. (2018) did not exhibit a transition to inhomogeneous mixing with increasing spatial scale, because the supersaturation relaxation time was too small compared to the droplet evaporation time.
Our analysis supports general conclusions regarding in-situ observation of droplets in turbulent clouds. First, most of the local and instantaneous snapshots of droplet configurations observed by Beals et al. (2015) cannot be in the steady states of a global mixing process that mixed undiluted cloud with dry environmental air. However, a local droplet population may still be in a local steady state, established as the droplets saturated the air locally. Such local steady states belong to the transient of a global mixing process. In order to understand the nature of this transient, it was necessary to consider the whole range of possible steady states at different length-scales ( Figure 1a). In short, clouds are not equilibrated at large scales, yet local steady states occur at small scales.
Our analysis also indicates that most of the droplet populations observed by Beals et al. (2015) are likely to have resulted from mixing with premixed air, and we concluded that the corresponding local steady states arose by dilution rather than complete evaporation. Our model indicates that only a very few droplets evaporated completely.
We found that the statistical-model dynamics is somewhat slower than the DNS of Kumar et al. (2012; and that the tails of our droplet size distributions are somewhat lighter. We speculated that this may be due to the supersaturation dynamics being oversimplified. Since our model belongs to the family of established PDF models (Pope, 2000), it is clear how to address this question in the future (Pope, 1991;Vedula et al., 2001;Jenny et al., 2012).
Last but not least, our analysis highlights which additional observational data are needed for a more quantitative statistical-model analysis of in situ cloud-droplet measurements. To determine the three key parameters, the volume fraction of cloudy air as well as the two Damköhler numbers, one needs joint measurements of local droplet populations, supersaturation levels in their vicinity, and the sizes of the local cloud structures. This will allow us to characterise and understand the mechanisms underlying local mixing processes observed on different length-and time-scales and at different distances from the cloud edge. A challenge for the future is to understand the global picture, how evaporation distributes in the cloud, and where complete droplet evaporation takes place. This is necessary to improve the parameterisation of mixing and evaporation at the cloud edge in subgrid-scale models, in order to better represent the radiative effects of clouds.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.