Key parameters for droplet evaporation and mixing at the cloud edge
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b RESEARCH ARTICLE
Key parameters for droplet evaporation andmixing at the cloud edge
J. Fries | G. Sardina | G. Svensson | B. Mehlig Department of Physics, GothenburgUniversity, SE-41296 Gothenburg, Sweden Department of Mechanics and MaritimeSciences, Division of Fluid Dynamics,Chalmers University of Technology,SE-41296 Gothenburg, Sweden Department of Meteorology, StockholmUniversity and Swedish e-science ResearchCentre, Stockholm, Sweden
Correspondence
Bernhard Mehlig, Department of Physics,Gothenburg University, SE-41296Gothenburg, SwedenEmail: [email protected]
Funding information
Vetenskapsrådet, grant number:2017-3865; Formas, grant number:2014-585; Knut and Alice WallenbergFoundation, Dnr. KAW 2014.0048
The distribution of liquid water in ice-free clouds determinestheir radiative properties, a significant source of uncertaintyin weather and climate models. Evaporation and turbulentmixing cause a cloud to display large variations in droplet-number density, but quite small variations in droplet size (Bealset al., 2015). Yet direct numerical simulations of the joint ef-fect of evaporation and mixing near the cloud edge predictquite different behaviors, and it remains an open questionhow to reconcile these results with the experimental find-ings. To infer the history of mixing and evaporation fromobservational snapshots of droplets in clouds is challengingbecause clouds are transient systems. We formulated a sta-tistical model that provides a reliable description of the eva-poration-mixing process as seen in direct numerical simula-tions, and allows to infer important aspects of the historyof observed droplet populations, highlighting the key mech-anisms at work, and explaining the differences between ob-servations and simulations.
KEYWORDS cloud micro-physics, turbulent mixing, droplet evaporation,Lagrangian droplet dynamics | 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 uncertaintyin climate projections (Pincus et al., 2018; Stocker et al., 2013). A key challenge is to understand how entrainmentof dry air at the edges of ice-free clouds affects the size distribution and number density of droplets (Blyth, 1993).This is important because 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 aresensitive 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(Caldwell et al., 2016; Dufresne and Bony, 2008). The size distribution and number density of droplets are key ingre-dients, 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).Regarding precipitation efficiency, the mechanism of rain formation in ice-free clouds is a longstanding unresolvedproblem in atmospheric physics (Grabowski and Wang, 2013). A broad initial droplet-size distribution is needed toactivate the collisions and coalescences of droplets that are necessary for the rapid onset of rain formation observedempirically in warm clouds (Devenish et al., 2012; Szumowski et al., 1997, 1998). Microscopic droplets grow bycondensation of water vapour, or shrink by evaporation. Yet when a droplet-containing parcel does not mix with itssurroundings, condensation causes the droplet-size distribution to narrow because the diffusional growth of a dropletis 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 (Vail-lancourt et al., 2001). As a consequence, nearby droplets may have experienced quite different growth histories. Thedroplets of a cloudy parcel that is mixed with dry air evaporate at different rates, so that the droplet-size distributionbroadens (Lanotte et al., 2009; Sardina et al., 2015, 2018; Li et al., 2020). Cloud-resolving simulations (Hoffmann andFeingold, 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 wheredroplets can rapidly evaporate. Droplets in regions with higher water-vapour concentration, by contrast, may saturatethe air and survive for a much longer time, Figure 1( a ). While droplets evaporate turbulence mixes the cloud at manylength scales, ranging from the Kolomogorov length – of the order of millimetres (Devenish et al., 2012) – to a fewkilometres (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, andlimited by thermodynamics if Da is small. The dynamics at large Damköhler numbers is referred to as inhomogeneousmixing (Baker et al., 1980), where some droplets evaporate completely while others do not evaporate at all. Small-Damixing is called homogeneous, where droplets evaporate at approximately the same rate, so that the droplet-sizedistribution remains narrow.The notion of homogeneous and inhomogeneous mixing remains debated (Tölle and Krueger, 2014), but it can begiven a precise meaning in terms of the fraction P e ( t ) of droplets that have completely evaporated at time t . However,it is at present not understood which mechanisms and parameters that determine the transition from homogeneousto inhomogeneous mixing. Several authors have attempted to describe turbulent mixing in terms of one Damköhlernumber. Lehmann et al. (2009) used a combined microphysical response time τ r , a function of the two thermodynamic . Fries et al. 3 (a) (b) ℓ ℓ ℓ ss c s e − w w − L L x FIGURE 1 ( a ) Effect of mixing on local droplet populations (schematic). On large spatial scales ( ℓ and ℓ ) mixingand evaporation is not yet complete, but some smaller regions (of size ℓ ) are in local steady states: saturated airwith droplets, or subsaturated air without droplets. ( b ) Initial cloud configuration in our model. Before mixing, moistair and droplets reside in a w × L × L slab, contained in a cubic domain of side length L . Regions with dry air aredashed. The solid line is the initial profile of supersaturation s (see text).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 nu-merical simulations (DNS) by Kumar et al. (2018) explored how the nature of mixing changes as τ L / τ r increases withthe linear size L of the simulated domain. However, no clear sign of inhomogeneous mixing was found. The authorsmention that this may be a consequence of the thermodynamic setup used. Another possibility is that the simulatedsystem was just not large enough. Jeffery (2007) emphasised that evaporation and mixing cannot be described by asingle Da alone, because there are two thermodynamic time scales, leading to three key non-dimensional parameters,Da d , Da s , and the volume fraction χ of cloudy air. However, Jeffery only studied the case Da d = Da s and did notdiscuss 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 forthe droplet-size distribution, and emphasised the significance of a non-dimensional parameter, their potential evapo-ration parameter.Here we derive a statistical model for evaporation and turbulent mixing at the cloud edge from first principles. Themodel takes into account the multi-scale turbulent dynamics, as turbulent clouds can have large Reynolds numbers;Re ∼ is a conservative estimate for convective clouds (Devenish et al., 2012). The model quantitatively predictsthe outcomes of the DNS by Kumar et al. (2012, 2013, 2014, 2018). Furthermore, the model shows in accordancewith Jeffery (2007) that the evolution of the droplet-size distribution is determined by Da d , Da s , and χ . We find thatthe potential evaporation parameter of Pinsky et al. (2016a) and Pinsky and Khain (2018a) is simply determined by theratio of the Damköhler numbers. Finally, the model allows to interpret the results of in-situ measurements of cloudsat 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 ofside length L with periodic boundary conditions. Initially, the saturated or slightly supersaturated moist air withsupersaturation s c ≥ is contained in a w × L × L slab together with N randomly distributed water droplets, Figure 1( b ). J. Fries et al.
The dry air, initially outside the slab, has negative supersaturation s e < . The mixing is driven by statistically stationaryhomogeneous 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, 2013, 2014, 2018), which allows us tounderstand 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 (Vaillan-court et al., 2002; Paoli and Shariff, 2009; Lanotte et al., 2009; Kumar et al., 2014; Perrin and Jonker, 2015). Weneglect buoyancy, particle inertia and settling, temperature changes due to vertical motion, temperature and pressuredependencies of the thermodynamic coefficients, and subsume the joint effects of temperature and water vapourinto 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 positionof droplet α is x α ( t ) , its radius equals r α ( t ) , and the index α ranges from to N . We non-dimensionalise as follows: u ′ = u / U , x ′ = x /( Uτ L ) , t ′ = t / τ L , p ′ = p /( ̺ a U ) , x ′ α = x α /( Uτ L ) , s ′ = s /| s e | , r ′ α = r α / r , and s ′ c = s c /| s e | . Here U = p TKE / is the turbulent r.m.s. velocity, τ L = TKE / ε is the large-eddy time [ ∝ L / U if the size of the largest eddiesis of the order L (Pope, 2000)], ̺ a is the reference mass density of air, r = [ N − Í N α =1 r α ( ) ] / is the initial volume-averaged droplet radius, and | s e | is the (positive) subsaturation of the air outside the initial cloud slab. Dropping theprimes, the microscopic equations take the non-dimensional form: DD t u = − + p + Re − L + u with + · u =0 , (1a) DD t s = ( Re L Sc ) − + s − Da s χV r α ( t ) s ( x α , t ) , (1b) dd t x α = u ( x α , t ) , (1c) dd t r α = Da d s ( x α ( t ) , t )/( r α ) if r α > . (1d)Equations 1a are the incompressible Navier-Stokes equations, with Lagrangian time derivative DD t = ∂ t + ( u · + ) . In DNSof Equations 1, a forcing is imposed to sustain stationary turbulence. This is not necessary in the model introducedbelow, and we therefore do not include a forcing term in Equation 1a. Equation 1b is the equation for supersaturation.The first term on its r.h.s. describes diffusion of the supersaturation s ( x , t ) . The second term models the effect ofcondensation and evaporation through the average r α ( t ) s ( x α , t ) , taken over all droplets in the vicinity of x . Dropletsare advected by the turbulent flow (Equation 1c), and Equation 1d models how the droplet radius r α changes due toevaporation and condensation. When a droplet has evaporated completely, we impose that it must remain at r α = 0 .In the derivation of Equation 1d, s ( x α ( t ) , t ) enters as the supersaturation at distances from the droplet that are muchlarger than the droplets radius (Rogers and Yau, 1989). In other words, Equation 1d relies on a scale separation betweendroplet sizes and lengths that characterise supersaturation fluctuations generated by turbulent mixing (Vaillancourtet al., 2001). As a consequence, droplets interact locally with the supersaturation field over finite volumes through theaverage r α ( t ) s ( x α , t ) in Equation 1b. Further details regarding Equations 1 are given in the Supporting Information(SI), where we also show how to derive Equations 1 from the more detailed dynamical description of Vaillancourt et al.(2001, 2002), Kumar et al. (2014, 2018), and Perrin and Jonker (2015).An advantage of writing the dynamics in non-dimensional form is that this determines the independent non-dimensional parameters. First, Re L = TKE /( εν ) is the turbulence Reynolds number (Pope, 2000), ν is the kinematicviscosity of air. The Schmidt number is defined as Sc = ν / κ , where κ is the diffusivity of s . The volume fraction of . Fries et al. 5 cloudy air is given by χ = w / L , and V = [ L /( Uτ L ) ] is the dimensionless domain volume. The Damköhler numberDa s is defined as Da s = τ L / τ s , where τ s = ( πA A ̺ w n r ) − is the supersaturation relaxation time. This is the timescale at which the supersaturation decays towards saturation, assuming that all droplets have the same radius r ,and for droplet number density n = N /( w L ) . Further, ̺ w is the density of pure liquid water, and A and A arethermodynamic coefficients, specified in the SI. The Damköhler number Da d is defined as Da d = τ L / τ d , where τ d = r /( A | s e |) is the droplet evaporation time, the time that it takes for a droplet of radius r to evaporate completelyin a constant ambient supersaturation s e < .The Damköhler numbers determine the extent to which saturation and droplet evaporation are limited by therate 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 rateof mixing. Droplet evaporation is mixing limited at large Da d , since droplets then evaporate more rapidly than theexposure to subsaturated air changes. At small Da d , mixing is faster than droplet evaporation, and droplets tend toevaporate mainly after the system has been mixed. The droplets then experience roughly the same supersaturationas they evaporate.In the limit Re L → ∞ , three key non-dimensional parameters remain in Equations 1: χ , Da d , and Da s . The systemcan be parameterised by χ , Da d , and R = Da d / Da s . (2)In this way, the scale dependence of the mixing process is contained in Da d only. The Damköhler-number ratio R is inversely proportional to the density of liquid water in the cloud slab; it regulates the moisture of the mix-ing process (details in SI). The bifurcation between moist steady states, where droplets remain in saturated air, anddry steady states, where all droplets have evaporated completely (Jeffery, 2007; Kumar et al., 2013; Pinsky et al.,2016b), occurs at a critical value of R , R c . The critical ratio R c can be computed from the conserved quantity θ = −h s ( t ) i − χ R [ − P e ( t ) ] h r ( t ) i , 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 evapo-rated droplets, the fraction of droplets for which r α ( t ) = 0 at time t . Furthermore, h s ( t ) i = V − ∫ V s ( x , t ) d x is thevolume average of supersaturation, and (cid:10) r ( t ) (cid:11) = {[ − P e ( t ) ] N } − Í N α =1 r α ( t ) is the mean cubed droplet radiusconditioned on r α ( t ) > by the factor [ − P e ( t ) ] − . The conservation of θ can be concluded by integrating Equa-tion 1b for supersaturation over the domain volume, see the SI for details. Moist steady states have (cid:10) r ( t ) (cid:11) > and h s ( t ) i = 0 , dry steady states have P e ( t ) = 1 and h s ( t ) i < , so the sign of θ determines whether the steady-state ismoist or dry. The value of θ is determined by the initial conditions, θ = −h s ( ) i − χ /( R ) . Setting θ = 0 , we find thecritical Damköhler-number ratio R c = − χ /h s ( ) i .DNS of Equations 1 for experimentally observed dissipation rates and droplet-number densities are feasible onlyfor quite small systems (Kumar et al., 2018). This restricts the range of scales that can be explored, and makes itdifficult to detect inhomogeneous mixing in DNS. We therefore pursue an alternative approach and adapt a PDFmodel (Pope, 2000) – commonly used to describe combustion processes (Haworth, 2010) – to the inhomogeneouscloud edge. As opposed to the kinematic statistical models reviewed by Gustavsson and Mehlig (2016), we must alsotake thermodynamic processes into account. J. Fries et al. | Statistical model
Statistical models have been used to describe droplets in a supersaturation field that fluctuates around zero, as in thecloud core (Paoli and Shariff, 2009; Sardina et al., 2015; Chandrakar et al., 2016; Siewert et al., 2017). At the cloudedge, 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 supersaturationfield. This does not take into account that mixing is local, and that small droplets are advected together with thesupersaturation field. For an accurate description of mixing and evaporation, it is essential to describe how eachdroplet carries its own local supersaturation (Siewert et al., 2017). Our model does just that. It is derived from firstprinciples using the established framework of PDF models (Pope, 2000).For the configuration shown in Figure 1( b ) we derive one-dimensional statistical-model equations from Equa-tions 1 (details in the SI): d u = − C u d t + (cid:0) C (cid:1) / d η , (3a) dd t s = − C φ [ s − h s ( x , t ) i] − Da s χV h r ( t ) s ( x , t ) i , (3b) dd t x = u , (3c) dd t r = Da d s /( r ) if r > . (3d)Equation 3a describes the fluctuating acceleration of Lagrangian fluid elements in turbulence. Here, d η are Brownianincrements with zero mean and variance d t , and C is an empirical constant (Pope, 2011). Each Lagrangian fluidelement has a supersaturation s , and may contain a droplet (at position x , of size r ). Equation 3b approximates thesupersaturation dynamics as decay towards h s ( x , t ) i , regulated by the empirical constant C φ (Pope, 2000). The secondterm on the r.h.s. of Equation 3b represents the effect of condensation and evaporation through h r ( t ) s ( x , t ) i . Theposition-dependent averages h· · · i in Equation 3b are taken over fluid elements located at x at time t (details in theSI). The statistical-model Equations 3 becomes independent of Re L at large Reynolds numbers, where C approachesa definite limit (Pope, 2011). It is independent of Sc, in accordance with the known behaviour of advected scalars infully developed turbulence (Shraiman and Siggia, 2000).Equations 3 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, contain unclosed terms thatmust be approximated (Pope, 1985; Haworth, 2010) to obtain a closed model such as Equations 3. Following Pope(2000), we cast the model into the form of stochastic dynamical equations for Lagrangian fluid elements. Since thedynamics is statistically one-dimensional in our configuration, we can average over the y and z coordinates to obtainEquations 3 in one-dimensional form. For the closures, we rely on standard approximations, common and justifiedin PDF modeling of single-phase flows (Pope, 1985), and in models for turbulent combustion (Haworth, 2010; Jennyet al., 2012; Stöllinger et al., 2013). The explicit mathematical approximations for the closures provided in the SI renderthe interpretation of the statistical model definite, and indicate how to improve the model when necessary.In the following, we briefly summarise the closures. Equation 3a contains the closure for fluid-element acceler-ations. It reproduces the empirically observed effect of turbulent diffusion of passive-scalar averages (Pope, 2000).Equation 3b contains two closure approximations. First, the decay towards h s ( x , t ) i approximates the diffusion ofsupersaturation. This closure ensures that the mean of a passive scalar is conserved, and that a passive scalar remainsbounded between its minimal and maximal values (Pope, 2000). Furthermore, it describes the decay of passive-scalarvariance in statistically homogeneous turbulent mixing of two scalar concentrations (Pope, 2000). However, thisclosure does not reproduce the relaxation of the single-point PDF of scalar concentration – from a two-peaked distri- . Fries et al. 7 TABLE 1 Summary of statistical-model simulations, details in SI. Damköhler numbers Da d and Da s ,Damköhler-number ratio R , critical ratio R c , and volume fraction χ of cloudy air. Simulation Da d Da s R R c χ Figure 2 and 3( a ) [dry] 2.44 0.968 2.52 0.859 0.428Figure 2 [moist] 1.09 1.43 0.76 0.859 0.428Figure 3( b ) [very moist] 0.754 8.20 0.092 0.683 0.4Figure 4 5E-3-4E2 1E-3-9E3 4E-2-4E0 0.913 0.429Figure 5( a ) 1E-2-1E3 6E-2-6E3 0.17 0.18-2.7 0.2-0.8Figure 5( b ) 1E-2-8E2 3E-3-4E4 2.4E-2-2.9E-2 0.38-0.41 0.369-0.374bution via a U-shaped distribution into a Gaussian (Eswaran and Pope, 1988; Pope, 1991). For our initial conditions[Figure 1 (b) ], the decay towards h s ( x , t ) i captures the supersaturation fluctuations experienced by a fluid element asit moves towards or away from the most cloudy region. Note that this closure does not account for large saturatedcloud structures that tend to relax slowly towards h s ( x , t ) i . Such events are most relevant during the initial stagesof the mixing-evaporation process, and their effect is expected to diminish with time, as large cloud structures aremixed into smaller and smaller structures. Consequently, the statistical model may describe short initial transientsonly qualitatively, not quantitatively.The second closure in Equation 3b approximates the effects of droplet phase change on the supersaturationfield: the local average r α ( t ) s ( x α , t ) in Equation 1b is replaced by the ensemble average h r ( t ) s ( x , t ) i . This is thesimplest closure that preserves the conservation of the parameter θ , and it is therefore common in PDF modelsthat describe combustion of particles in turbulence (Haworth, 2010; Jenny et al., 2012; Stöllinger et al., 2013). Theaverage h r ( t ) s ( x , t ) i takes into account that droplet evaporation is delayed locally when nearby droplets saturate thesurrounding air. Since we obtain closure by replacing the local average r α ( t ) s ( x α , t ) by an average over fluid elementswith one-dimensional dynamics, variations in the rate of supersaturation relaxation in the y and z -coordinates are notdescribed.It is expected that the large cloud structures mentioned above, and their three-dimensional forms, matter moreat very large Damköhler numbers. Therefore it cannot be excluded that the statistical model is only qualitative inthis extreme limit. Below we show that the model works very well even for the largest Damköhler numbers in DNSstudies Kumar et al. (2012, 2014). Also, since the statistical model is derived using the established framework of PDFmodels (Pope, 2000), it can be straightforwardly improved by incorporating additional variables in the probabilisticdescription (Pope and Chen, 1990; Pope, 2000; Meyer and Jenny, 2008), or by using more refined approximations(Pope, 1991; Jenny et al., 2012). | RESULTS3.1 | Comparison with DNS
The statistical model can be used to understand DNS results of Kumar et al. (2012, 2013, 2014, 2018). Figure 2 showsgood agreement for the time evolution of the fraction P e ( t ) of droplets that have completely evaporated, even thoughthe statistical-model dynamics is slightly slower. Panels ( a ) and ( b ) in Figure 3 show that the model reproduces the J. Fries et al.
TABLE 2 Parameters of DNS shown in Figure 4: Damköhler number Da d , Damköhler-number ratio R , critical ratio R c , and volume fraction χ of cloudy air. Some dimensional parameters are also shown: domain size L , meandissipation rate ε , and droplet-number density n of the initially cloudy air. Non-dimensional parameters Dimensional parameters
L ε n Reference Da d Da s R R c χ [cm] [cm /s ] [cm − ] ◦ Andrejczuk et al. (2006) 8E − -1E -3E − -9E -1E △ Kumar et al. (2012) 8E − -8E − − -8E − (cid:3) Kumar et al. (2013) 0.14, 0.31 0.62, 0.41 0.22, 0.73 0.68 0.4 26 34 164 ⋄ Kumar et al. (2014) 0.61-2.4 0.97-1.9 0.31-2.5 0.84 0.42 51 34 153 ▽ Kumar et al. (2018) 0.12-0.91 0.51-4.0 0.23 0.90-0.95 0.42-0.45 1E -2E P e ( t ) very well (Figure 2). It is also asignificant 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 mentionedabove, and this may be particularly important at large values of Da d , where locally saturated regions can persist for along 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 R / R c . We see how P ∗ e increases with both Da d and R / R c . The DNS results of Andrejczuk et al. (2006) and Kumar et al.(2012, 2013, 2014, 2018) form a pattern in Figure 4 that verifies these dependencies: open symbols correspond toDNS with little or no complete evaporation in the steady state ( P ∗ e < %), and filled symbols to P ∗ e > %. Figure 4 alsoexplains why the DNS of Kumar et al. (2018) did not exhibit significant levels of inhomogeneous mixing: since their R was quite small, small values of P ∗ e require values of Da d much larger than unity (Da d ∼ for P ∗ e = 10 %). Furthermore,the substantially different outcomes of the DNS of Kumar et al. (2012) and Kumar et al. (2014) are explained, theirparameters lie on opposite sides of the bifurcation line. Figure 4 also explains, at least qualitatively, numerical resultsof DNS of transient turbulence with quite different initial conditions (Andrejczuk et al., 2006), namely how the amountof complete droplet evaporation increases with both R / R c and Da d . There is no parameter corresponding directlyto Da d in the simulations of Andrejczuk et al. (2006), because they are for different initial conditions and flows. Wetherefore place these simulations in Figure 4 by computing a time-scale ratio that, in a qualitative sense, incorporatesthe same physics as Da d (details in the SI). Key parameters of the DNS in Figure 4 are summarised in Table 2, a completedescription is provided in the SI. . Fries et al. 9 moistdry tP e FIGURE 2 Fraction P e ( t ) of droplets that have completely evaporated as a function of non-dimensional time t ,parameters in Table 1. Coloured lines are statistical-model simulations, black lines are DNS of Kumar et al. (2014). -4 -2 (a) (b)dry very moist r r P D F t = 2.11 t = 0.84 steady state t = 1 . t = 0 . FIGURE 3 ( a ) Evolution of the droplet-size distribution (parameters in Table 1) for different times. The probabilitydensity of the non-dimensional droplet radius r is shown for the statistical-model (coloured lines) and DNS of Kumaret 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). | Mixing histories from observations
A common way of characterising the droplet content of a cloud is to plot the mean cubed radius r and numberdensity n of droplets for observed cloud-droplet populations in a mixing diagram. Figure 5( a ) shows a mixing diagramwith empirical data from Beals et al. (2015). Black crosses are values of r and n extracted from snapshots (linear size 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, acurve of global steady states ( r ∗ , n ∗ ) that result from homogeneous mixing (no complete evaporation) between differ-ent proportions of undiluted cloudy and dry environmental air (Gerber et al., 2008; Kumar et al., 2014; Pinsky et al.,2016a). Beals et al. (2015) calculated this line, it is also shown in Figure 5( a ). A fundamental problem is however thatit is not clear how to interpret mixing diagrams such as Figure 5( a ), since it is not clear that the empirically observeddroplet populations reflect global steady states (Pinsky and Khain, 2018a).It is nevertheless likely that most data points in Figure 5( a ) sample local steady states, i.e. locally well-mixeddroplet populations that reside in saturated air. To describe such droplet populations, one must refer to the multi- -2 -4 -3 -2 -1 drymoist R / R c Da d P ∗ e FIGURE 4 Steady-state fraction P ∗ e of completely evaporated droplets, as a function of Da d and R / R c , details inSI. The fraction P ∗ e is colour coded. The solid red line is the contour P ∗ e = 10 %, the black dashed line indicates thetransition between moist and dry steady states, and symbols indicate DNS results from previous studies: △ (Kumaret al., 2012), (cid:3) (Kumar et al., 2013), ⋄ (Kumar et al., 2014), ▽ (Kumar et al., 2018), and ◦ (Andrejczuk et al., 2006).Filled symbols indicate P ∗ e > %. The DNS of Andrejczuk et al. (2006) should not be quantitatively compared to thestatistical model results, since they are for different initial conditions and decaying turbulence (see text and SI). Red,blue and light blue circles correspond to the dry, moist, and very moist simulations in Figures 2 and 3.scale turbulent mixing process. We attempted this analysis using the statistical model, assuming that the statisticalmodel with the initial condition shown Figure 1( b ) describes how a cloud structure at the spatial scale L develops.Under this assumption, n ∗ and r ∗ are given by the droplet-number density (normalised by n ) and the mean cubeddroplet radius (cid:10) r ( t ) (cid:11) in the steady state, and we can conjecture the mixing histories that formed the measured dropletpopulations.We begin by noting that χ and P ∗ e are completely determined for any steady-state point ( r ∗ , n ∗ ) in a mixing diagram.To show this, we write the volume-averaged initial supersaturation as h s ( ) i = ( s c ) ( χ + χ )− , where χ is a constantthat depends on the initial supersaturation profile (details in the SI). Inserting χ = n ∗ /( − P ∗ e ) (4)into θ = −h s ( t ) i − χ R [ − P e ( t ) ] h r ( t ) i we find: P ∗ e = 1 − n ∗ [ R ( s c ) ] n ∗ r ∗ + R [ − χ ( s c ) ] . (5)Equations 4 and 5 determine how to map ( r ∗ , n ∗ ) to ( χ , P ∗ e ). As a consistency check we note that one obtains thehomogeneous mixing line (Pinsky et al., 2016a) from Equation 5 by setting P ∗ e = 0 . This allows us to infer that R = 0 . and s c = χ = 0 from 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 χ . Eachstatistical-model simulation for given Da d , R , and χ yields a certain value of P ∗ e . This allows us to extract a value ofDa d for each point in the mixing diagram from our statistical-model simulations. The result is shown in Figure 5( a ).We see that Da d increases rapidly above the homogeneous mixing line. Estimating τ s ∼ s from Beals et al. (2015)and conservatively estimating ε ∼ cm s − for a convective cloud, a value of Da d = 1000 implies that an observeddroplet population was mixed at spatial scales of the order of L ∼ q ετ L ∼ km, larger than the size of the cloud. In . Fries et al. 11 (a) n ∗ r ∗ Da d .
024 0 .
027 0 . -2 -1 (b) R Da d P ∗ e =0 ( L ∼ m ) P ∗ e =1 % ( L ∼ m ) R m i n R m i n FIGURE 5 ( a ) Mixing diagram. Empirical data from Beals et al. (2015) (black crosses). The homogeneous mixing line(black) from top panel of Figure 3 of Beals et al. (2015) corresponds to R = 0 . . The coloured region shows wheresteady states are found in the mixing diagram for R = 0 . . The corresponding values of Da d are colour coded(legend). The red cross is the measurement shown in the middle panel of Figure 2 of Beals et al. (2015). ( b ) Values of R and Da d consistent with a steady state at the red cross in panel ( a ) (red). We estimate L ∼ m for Da d = 1 (bluecircle). These local mixing processes have R = R min (black dashed line), so no droplets evaporated completely( P ∗ e = 0 ). The green circle corresponds to Da d = 13 and R = 0 . with P ∗ e = 1 % and L ∼ m.other words, the rapid increase of Da d in Figure 5( a ) suggests that most of the data in the mixing diagram cannot bein a global steady state of a mixing process parameterised by R = 0 . .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 localmixing processes equilibrate at small spatial scales [Figure 1( a )]. We now discuss how the analysis of local steadystates may yield insight into possible local histories of the cloud. Air affected by earlier mixing events is not as dry asenvironmental air, so the mixing of undiluted cloud with premixed air is governed by smaller values of R . We thereforeask: which values of R are consistent with the assumption that the experimentally observed droplet population in themiddle panel of Figure 2 of Beals et al. (2015) – red cross in Figure 5( a ) – reflects a local steady state? Our model allowsus to determine possible combinations of Da d , R , and χ consistent with a local steady state. We know that R mustbe 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 R is R min = ( n ∗ − r ∗ ) n ∗ /[( s c ) ( χ + n ∗ ) − ] = 0 . .Figure 5( b ) shows values of R and Da d obtained from our statistical-model simulations that are consistent withthese constraints. We see that the range of possible values of Da d covers several orders of magnitude. This meansthat local mixing processes consistent with the red cross in Figure 5( a ) may have occurred over a large range of spatialscales. We also see that R does not vary much in Figure 5( b ), only between 0.024 and 0.03. This allows us to concludethat some important aspects of the mixing dynamics are essentially independent of spatial scale. First, the fact that R is substantially smaller than 0.17 indicates that the non-cloudy air was premixed. Second, using Equation 5, wefind that the reduction in droplet-number density was primarily caused by dilution even at the largest scales, since P ∗ e increases only up ∼ d , at Da d ∼ and R = 0 . . Put differently, χ ∼ n ∗ for all valuesof Da d we considered.How does the outcome of a local mixing process depend on its scale? Larger scales correspond to larger valuesof Da d , and Figure 5( b ) shows that complete droplet evaporation begins to occur around Da d = 1 , where R starts toexceed R min (blue circle). Estimating ε ∼ cm s − , a typical value for convective clouds (Devenish et al., 2012), wefind that Da d = 1 , τ s ∼ , and R = R min correspond to the spatial scale m. Mixing processes leading to the red cross in Figure 5( a ) that occurred at scales smaller than m were therefore perfectly homogeneous: none of the dropletsevaporated completely as they were diluted by premixed air. At larger spatial scales, small but non-zero fractionsof the droplets evaporated completely. Equation 5 gives R = 0 . for P ∗ e = 1 %, and from Figure 5( b ) we read offDa d = 13 (green circle). For ε ∼ cm s − these values correspond to m. This suggests that reductions in dropletnumber concentration are dominated by dilution, and not complete droplet evaporation, also for mixing processesthat range over hundreds of meters. Furthermore, since most data points in Figure 5( a ) reside well above the regionwhere equilibria are found for R = 0 . , 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 inhomo-geneous mixing; Da d parameterises the mixing-limited nature of droplet evaporation, and the ratio R = Da d / Da s regulates the self-limiting effect of droplet evaporation, namely that droplets cease to evaporate when they have sat-urated the surrounding air or evaporated completely. Analysing the parameters of our microscopic equations 1 we seethat R = − R where R is the potential evaporation parameter of Pinsky et al. (2016a) and Pinsky and Khain (2018a).So 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, 2019) concluded that the Damköhler number ratio R determines whether a cloudexpands by dilution or shrinks by complete droplet evaporation. A mixing process that mixes equal proportions ofsaturated cloudy and subsaturated non-cloudy air has R c = , so the symmetric configuration they adopted for thecloud edge implies that the cloud expands if R < , and shrinks otherwise. We note that whether the cloud expandsor shrinks depends on the position and scale at which one perceives it. A local mixing process with small R / R c tendstowards 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 expansionis part of a mixing process that consumes the cloud at larger length and time scales. Although global mixing processesare transient, they contain local steady states. Diluted and saturated local droplet populations [such as the red crossin Figure 5( a )] can be a part of such transients. However, such local steady states must eventually be abandoned asmixing proceeds globally.Adopting this multi-scale picture of mixing, in which a large-scale mixing process consists of many mixing pro-cesses at smaller scales, it is natural to expect that large ranges of the parameters χ , R and Da d are relevant. If onemoves the domain of a local mixing process from the interior of the cloud towards the cloud edge, the liquid watercontent decreases, so that R / R c and Da d increase. Lehmann et al. (2009) point out that Da d increases as one per-ceives mixing processes at larger and larger spatial scales. This corresponds to moving to the right in Figure 4. Wenote that R / R c and Da d tend to increase with distance from the interior of the cloud. Moving the sampling volumetowards the cloud edge then corresponds to a motion upwards and to the right in Figure 4. The amount of completedroplet evaporation increases in this direction, consistent with the fact that complete droplet evaporation takes placeat the cloud edge.A number of assumptions may influence our interpretation of the empirical data in Section 3.2. First, the modelconfiguration in Figure 1( b ) is simplified compared to real clouds, which have irregular shapes that deform during themixing-evaporation process. Second, the observational method may not detect droplets with radii smaller than 3 µ m( r α = 0 . ), as stated by Beals et al. (2015). If many small droplets where not detected, the observations in Figure 5( a )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 . Fries et al. 13 model reproduces the evolution of droplet-size distributions in DNS qualitatively in Figure 5. This expectation iscorroborated by the robust tendency for P ∗ e to increase with increasing values of Da d and R in Figure 4, and followsdirectly 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 thatthe next step in improving the statistical model should aim at reproducing the fastest evaporation rates in the transientmixing-evaporation process. A better agreement in the tails could be achieved by refining the closure for supersatu-ration diffusion by using a dynamic C φ in Equation 3b (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 byimproved closures (Vedula et al., 2001; Pope, 1991; 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 modelexplains results of earlier DNS studies of mixing (Andrejczuk et al., 2006; Kumar et al., 2012, 2013, 2014, 2018), andshows that two thermodynamic time scales are important for a mixing process, the droplet evaporation time and thesupersaturation relaxation time. This means that one must consider two Damköhler numbers in order to quantify themixing-evaporation dynamics. We concluded that the simulations of Kumar et al. (2018) did not exhibit a transitionto inhomogeneous mixing with increasing spatial scale because the supersaturation relaxation time was too smallcompared 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 inthe steady states of a global mixing process that mixed undiluted cloud with dry environmental air. However, a localdroplet population may still be in a local steady state, established as the droplets saturated the air locally. Such localsteady states belong to the transient of a global mixing process. In order to understand the nature of this transient, itwas necessary to consider the whole range of possible steady states at different length scales, Figure 1 ( a ). 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 haveresulted from mixing with premixed air, and we concluded that the corresponding local steady states arose by dilutionrather than complete evaporation. Our model indicates that only very few droplets evaporated completely.We found that the statistical-model dynamics is somewhat slower than the DNS of Kumar et al. (2012, 2014), andthat the tails of our droplet-size distributions are somewhat lighter. We speculated that this may be due to that thesupersaturation dynamics is 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 (Vedula et al., 2001; Pope, 1991; Jenny et al., 2012).Last but not least our analysis highlights which additional observational data is needed for a more quantitativestatistical-model analysis of in-situ cloud-droplet measurements. To determine the three key parameters, the volumefraction of cloudy air as well as the two Damköhler numbers, one needs joint measurements of local droplet popula-tions, supersaturation levels in their vicinity, and the sizes of the local cloud structures. This will allow to characteriseand understand the mechanisms underlying local mixing processes observed on different length and time scales, andat different distances from the cloud edge. A challenge for the future is to understand the global picture, how evapo-ration distributes in the cloud, and where complete droplet evaporation takes place. This is necessary to improve theparameterisation of mixing and evaporation at the cloud edge in sub-grid scale models, in order to better representthe radiative effects of clouds. | SUPPORTING INFORMATION
The SI contains Tables S1-S3 and Appendix S1. Table S1 lists all parameters of our statistical model simulations.Tables S2 and S3 specify the DNS results from Kumar et al. (2012, 2013, 2014, 2018), and Andrejczuk et al. (2006)in Figure 4. Appendix S1 describes the relation between our microscopic Equations 1, and those used in the DNS ofothers (Vaillancourt et al., 2002; Paoli and Shariff, 2009; Lanotte et al., 2009; Perrin and Jonker, 2015; Kumar et al.,2014, 2018). We discuss the approximations made in deriving Equations 1, and quantify their accuracy. We alsodescribe how to derive the statistical model, Equations 3, from Equations 1, using the framework of PDF methods(Pope, 2000). We detail our computer simulations of the statistical model, and address their numerical convergence.Finally, Appendix S1 contains details concerning our interpretation of the data of Beals et al. (2015), and of directnumerical simulations conducted by other authors.
ACKNOWLEDGEMENTS
JF and BM thank B. Kumar and J. Schumacher for providing details of their simulations. We acknowledge supportby Vetenskapsrådet (grant number 2017-368 3865), Formas (grant number 2014-585), and by the grant ‘Bottlenecksfor particle growth in turbulent aerosols’ from the Knut and Alice Wallenberg Foundation, Dnr. KAW 2014.0048.Simulations were performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE)provided by the Swedish National Infrastructure for Computing (SNIC). . Fries et al. 15 | LIST OF SYMBOLS t time x = ( x , y , z ) spatial position N number of droplets at initialisation P e fraction of completely evaporated droplets u fluid velocity p pressure L side length of cubic simulation domain w width of initially cloudy region in simulation domain s c supersaturation within initial cloud slab s e supersaturation outside initial cloud slab TKE turbulent kinetic energy U root-mean-square of fluid velocity ε turbulent dissipation rate per unit mass ν kinematic viscosity κ diffusivity of supersaturation s supersaturation r droplet radius r initial volume radius of droplets n droplet-number density of intially cloudy region r α ( t ) s ( x α , t ) average of r α ( t ) s ( x α , t ) for droplets in the vicinity of x at time t DD t Lagrangian time derivative τ d droplet evaporation time τ s supersaturation relaxation time τ ℓ time scale for mixing at the length scale ℓτ L large-eddy turnover time in simulation domainRe L turbulence Reynolds numberSc Schmidt number V non-dimensional volume of simulation domainDa d Damköhler number based on droplet evaporation timeDa s Damköhler number based on supersaturation relaxation time R Damköhler-number ratio R c Critical Damköhler-number ratio R min Lower bound for Damköhler-number ratio related to mixing diagrams χ volume fraction of cloudy air χ contribution to the initial volume average of supersaturation θ conserved quantity that reflects the conservation of water and energy C , C φ empirical constants h . . . i volume average or ensemble average in statistical model REFERENCES
Abade, G. C., Grabowski, W. W. and Pawlowska, H. (2018) Broadening of cloud droplet spectra through eddy hopping: Tur-bulent entraining parcel simulations.
Journal of the Atmospheric Sciences , , 3365–3379.Andrejczuk, M., Grabowski, W. W., Malinowski, S. P. and Smolarkiewicz, P. K. (2006) Numerical simulation of cloud–clear airinterfacial mixing: Effects on cloud microphysics. Journal of the Atmospheric Sciences , , 3204–3225.Baker, M. B., Corbin, R. G. and Latham, J. (1980) The influence of entrainment on the evolution of cloud droplet spectra: I. Amodel of inhomogeneous mixing. Quarterly Journal of the Royal Meteorological Society , , 581–598.Beals, M. J., Fugal, J. P., Shaw, R. A., Lu, J., Spuler, S. M. and Stith, J. L. (2015) Holographic measurements of inhomogeneouscloud mixing at the centimeter scale. Science , , 87–90.Blyth, A. M. (1993) Entrainment in cumulus clouds. Journal of applied meteorology , , 626–641.Bodenschatz, E., Malinowski, S. P., Shaw, R. A. and Stratmann, F. (2010) Can we understand clouds without turbulence? Science , , 970–971.Burnet, F. and Brenguier, J.-L. (2007) Observational study of the entrainment-mixing process in warm convective clouds. Journal of the Atmospheric Sciences , , 1995–2011.Caldwell, P. M., Zelinka, M. D., Taylor, K. E. and Marvel, K. (2016) Quantifying the sources of intermodel spread in equilibriumclimate sensitivity. Journal of Climate , , 513–524.Chandrakar, K. K., Cantrell, W., Chang, K., Ciochetto, D., Niedermeier, D., Ovchinnikov, M., Shaw, R. A. and Yang, F. (2016)Aerosol indirect effect from turbulence-induced broadening of cloud-droplet size distributions. Proceedings of the NationalAcademy of Sciences , , 14243–14248.Devenish, B. J., Bartello, P., Brenguier, J.-L., Collins, L. R., Grabowski, W. W., IJzermans, R. H. A., Malinowski, S. P., Reeks,M. W., Vassilicos, J. C., Wang, L.-P. and Warhaft, Z. (2012) Droplet growth in warm turbulent clouds. Quarterly Journal ofthe Royal Meteorological Society , , 1401–1429.Dimotakis, P. E. (2005) Turbulent Mixing. Annual Review of Fluid Mechanics , , 329–356.Dufresne, J.-L. and Bony, S. (2008) An assessment of the primary sources of spread of global warming estimates from coupledatmosphere–ocean models. Journal of Climate , , 5135–5144.Eswaran, V. and Pope, S. (1988) Direct numerical simulations of the turbulent mixing of a passive scalar. Physics of Fluids , ,506–520.Frisch, U. (1995) Turbulence: The Legacy of A. N. Kolmogorov . Cambridge University Press.Gerber, H. E., Frick, G. M., Jensen, J. B. and Hudson, J. G. (2008) Entrainment, mixing, and microphysics in trade-wind cumulus.
J. Meteor. Soc. Japan , , 87–106.Grabowski, W. W. and Wang, L.-P. (2013) Growth of cloud droplets in a turbulent environment. Annual Review of Fluid Me-chanics , , 293–324.Gustavsson, K. and Mehlig, B. (2016) Statistical models for spatial patterns of heavy particles in turbulence. Advances inPhysics , , 1–57.Haworth, D. C. (2010) Progress in probability density function methods for turbulent reacting flows. Progress in Energy andCombustion Science , , 168–259.Hoffmann, F. and Feingold, G. (2019) Entrainment and mixing in stratocumulus: Effects of a new explicit subgrid-scale schemefor large-eddy simulations with particle-based microphysics. Journal of the Atmospheric Sciences , , 1955–1973. . Fries et al. 17 Jeffery, C. A. (2007) Inhomogeneous cloud evaporation, invariance, and Damköhler number.
J. Geophys. Res. , , D24S21.Jenny, P., Roekaerts, D. and Beishuizen, N. (2012) Modeling of turbulent dilute spray combustion. Progress in Energy andCombustion Science , , 846–887.Kokhanovsky, A. (2004) Optical properties of terrestrial clouds. Earth-Science Reviews , , 189–241.Kumar, B., Götzfried, P., Suresh, N., Schumacher, J. and Shaw, R. A. (2018) Scale dependence of cloud microphysical responseto turbulent entrainment and mixing. Journal of Advances in Modeling Earth Systems , , 2777–2785.Kumar, B., Janetzko, F., Schumacher, J. and Shaw, R. A. (2012) Extreme responses of a coupled scalar–particle system duringturbulent mixing. New J. Phys. , , 115020.Kumar, B., Schumacher, J. and Shaw, R. A. (2013) Cloud microphysical effects of turbulent mixing and entrainment. Theor.Comp. Fluid Dyn. , , 361–376.— (2014) Lagrangian mixing dynamics at the cloudy–clear air interface. Journal of the Atmospheric Sciences , , 2564–2580.Lamb, D. and Verlinde, J. (2011) Physics and Chemistry of Clouds . Cambridge University Press.Lanotte, A. S., Seminara, A. and Toschi, F. (2009) Cloud droplet growth by condensation in homogeneous isotropic turbulence.
Journal of the Atmospheric Sciences , , 1685–1697.Lehmann, K., Siebert, H. and Shaw, R. A. (2009) Homogeneous and inhomogeneous mixing in cumulus clouds: Dependenceon local turbulence structure. Journal of the Atmospheric Sciences , , 3641–3659.Li, X.-Y., Brandenburg, A., Svensson, G., Haugen, N. E. L., Mehlig, B. and Rogachevskii, I. (2020) Condensational and collisionalgrowth of cloud droplets in a turbulent environment. Journal of the Atmospheric Sciences , , 337–353.Lu, C., Liu, Y., Zhu, B., Yum, S. S., Krueger, S. K., Qiu, Y., Niu, S. and Luo, S. (2018) On which microphysical time scales to usein studies of entrainment-mixing mechanisms in clouds. J. Geophys. Res. Atmos. , , 3740–3756.Mauritsen, T., Stevens, B., Roeckner, E., Crueger, T., Esch, M., Giorgetta, M., Haak, H., Jungclaus, J., Klocke, D., Matei, D.,Mikolajewicz, U., Notz, D., Pincus, R., Schmidt, H. and Tomassini, L. (2012) Tuning the climate of a global model. Journalof Advances in Modeling Earth Systems , , 1–18.Meyer, D. W. and Jenny, P. (2008) An improved mixing model providing joint statistics of scalar and scalar dissipation. Com-bustion and Flame , , 490–508.Paoli, R. and Shariff, K. (2009) Turbulent condensation of droplets: Direct simulation and a stochastic model. Journal of theAtmospheric Sciences , , 723–740.Perrin, V. E. and Jonker, H. J. J. (2015) Lagrangian droplet dynamics in the subsiding shell of a cloud using direct numericalsimulations. Journal of the Atmospheric Sciences , , 4015–4028.Pincus, R., Winker, D., Bony, S. and Stevens, B. (2018) Shallow Clouds, Water Vapor, Circulation, and Climate Sensitivity . Springer.Pinsky, M. and Khain, A. (2018a) Theoretical analysis of mixing in liquid clouds – Part IV: DSD evolution and mixing diagrams.
Atmos. Chem. Phys. , , 3659–3676.— (2018b) Theoretical analysis of the entrainment–mixing process at cloud boundaries. Part I: Droplet size distributions andhumidity within the interface zone. Journal of the Atmospheric Sciences , , 2049–2064.— (2019) Theoretical analysis of the entrainment–mixing process at cloud boundaries. Part II: Motion of cloud interface. Journal of the Atmospheric Sciences , , 2599–2616.Pinsky, M., Khain, A. and Korolev, A. (2016a) Theoretical analysis of mixing in liquid clouds – Part 3: Inhomogeneous mixing. Atmos. Chem. Phys. , , 9273–9297. Pinsky, M., Khain, A., Korolev, A. and Magaritz-Ronen, L. (2016b) Theoretical investigation of mixing in warm clouds – Part 2:Homogeneous mixing.
Atmos. Chem. Phys. , , 9255–9272.Pope, S. (1991) Mapping closures for turbulent mixing and reaction. Theoretical and Computational Fluid Dynamics , , 255–270.Pope, S. B. (1985) PDF methods for turbulent reactive flows. Progress in Energy and Combustion Science , , 119–192.— (2000) Turbulent Flows . Cambridge University Press.— (2011) Simple models of turbulent flows.
Physics of Fluids , , 1–20.Pope, S. B. and Chen, Y. L. (1990) The velocity-dissipation probability density function model for turbulent flows. Physics ofFluids , , 1437–1449.Rogers, R. R. and Yau, M. K. (1989) A Short Course in Cloud Physics . Pergamon Press.Sardina, G., Picano, F., Brandt, L. and Caballero, R. (2015) Continuous growth of droplet size variance due to condensation inturbulent clouds.
Physical Review Letters , , 184501.Sardina, G., Poulain, S., Brandt, L. and Caballero, R. (2018) Broadening of cloud droplet size spectra by stochastic condensation:Effects of mean updraft velocity and CCN activation. Journal of the Atmospheric Sciences , , 451–467.Shraiman, B. I. and Siggia, E. D. (2000) Scalar turbulence. Nature , , 639–646.Siewert, C., Bec, J. and Krstulovic, G. (2017) Statistical steady state in turbulent droplet condensation. Journal of Fluid Me-chanics , , 254–280.Stocker, T., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V. and Midgley, P. (eds.) (2013) Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Inter-governmental Panel on Climate Change . Cambridge University Press.Stöllinger, M., Naud, B., Roekaerts, D., Beishuizen, N. and Heinz, S. (2013) PDF modeling and simulations of pulverized coalcombustion–Part 1: Theory and modeling.
Combustion and Flame , , 384–395.Szumowski, M. J., Rauber, R. M., Ochs III, H. T. and Beard, K. V. (1998) The microphysical structure and evolution of Hawaiianrainband clouds. Part II: Aircraft measurements within rainbands containing high reflectivity cores. Journal of the Atmo-spheric Sciences , , 208–226.Szumowski, M. J., Rauber, R. M., Ochs III, H. T. and Miller, L. (1997) The microphysical structure and evolution of Hawaiianrainband clouds. Part I: Radar observations of rainbands containing high reflectivity cores. Journal of the AtmosphericSciences , , 369–385.Tölle, M. H. and Krueger, S. K. (2014) Effects of entrainment and mixing on droplet size distributions in warm cumulus clouds. Journal of Advances in Modeling Earth Systems , , 281–299.Vaillancourt, P. A., Yau, M. K., Bartello, P. and Grabowski, W. W. (2002) Microscopic approach to cloud droplet growth bycondensation. Part II: Turbulence, clustering, and condensational growth. Journal of the Atmospheric Sciences , , 3421–3435.Vaillancourt, P. A., Yau, M. K. and Grabowski, W. W. (2001) Microscopic approach to cloud droplet growth by condensation.Part I: Model description and results without turbulence. Journal of the Atmospheric Sciences , , 1945–1964.Vedula, P., Yeung, P. and Fox, R. O. (2001) Dynamics of scalar dissipation in isotropic turbulence: a numerical and modellingstudy. Journal of Fluid Mechanics , , 29. r X i v : . [ phy s i c s . f l u - dyn ] F e b SUPPORTING INFORMATION
Supporting Information for "Key parameters fordroplet evaporation and mixing at the cloud edge"
J. Fries | G. Sardina | G. Svensson | B. Mehlig Department of Physics, GothenburgUniversity, SE-41296 Gothenburg, Sweden Department of Mechanics and MaritimeSciences, Division of Fluid Dynamics,Chalmers University of Technology,SE-41296 Gothenburg, Sweden Department of Meteorology, StockholmUniversity and Swedish e-science ResearchCentre, Stockholm, Sweden
Correspondence
Bernhard Mehlig, Department of Physics,Gothenburg University, SE-41296Gothenburg, SwedenEmail: [email protected]
This Supporting Information (SI) contains:
Appendix S1Tables S1 to S3SI References
Appendix S1
S1.1 | MICROSCOPIC EQUATIONS
We denote fluid velocity and pressure by u ( x , t ) and p ( x , t ) . Additional fields are the supersaturation s ( x , t ) , and thecondensation rate C d ( x , t ) . Initially, at t = 0 , there are N spherical droplets, labeled by α = 1 , . . . , N , with positions x α ( t ) , velocities v α ( t ) , and radii r α ( t ) . Our microscopic equations read: ∂ u ∂t + ( u · + ) u = − ̺ a + p + ν + u , + · u = 0 , (S1.1a) ∂s∂t + ( u · + ) s = κ + s − A C d ( x , t ) (S1.1b)with C d ( x , t ) = N Õ α =1 G ( | x − x α ( t ) |) π ̺ w d r α ( t ) d t (S1.1c)d x α d t = u ( x α ( t ) , t ) , (S1.1d)d r α d t = 2 A s ( x α ( t ) , t ) . (S1.1e)Here, ̺ a is the density of air and ν is its kinematic viscosity, ̺ w is the density of pure liquid water, κ is the diffusivity ofsupersaturation, and A and A are thermodynamic coefficients. In Equation S1.1c, G ( | x − x α ( t ) |) is a spatial kernel,normalised to unity. It makes it possible to construct the continuous condensation-rate field from the disperseddroplets, and ensures that the condensation rate is computed locally (Jenny et al., 2012). The supersaturation isdefined as s = e v / e vs ( T ) − , (S1.2)where e v is the partial pressure of water vapour and e vs ( T ) is the equilibrium water-vapour pressure at the tempera-ture T (Rogers and Yau, 1989). The partial pressure of water vapour depends on the density ̺ v and gas constant R v of water vapor, and is given by the ideal-gas law (Rogers and Yau, 1989): e v = ̺ v R v T . (S1.3)Equations S1.1 describe an incompressible flow with advected droplets that condense or evaporate in response tothe supersaturation in their vicinity. Evaporation or condensation causes the droplet radii r α to decrease or increase, ata rate regulated by the thermodynamic coefficient A . The supersaturation s ( x , t ) evolves according to an advection-diffusion equation, with a source term that describes the response of s ( x , t ) to droplet phase change. The strength ofthis response is regulated by the thermodynamic coefficient A . The parameters ̺ a , ν , κ , ̺ w , A , and A are assumedconstant in time. In addition, we impose that a droplet that has evaporated completely must remain at r α = 0 .Equation S1.1e relies on a scale separation: spatial scales of supersaturation fluctuations induced by turbulentmixing are much larger than the radii of droplets (Vaillancourt et al., 2001). This means that droplets do not imposeboundary conditions on the supersaturation s ( x , t ) in Equation S1.1b. They do interact locally with the supersatura-tion field, but over a finite volume determined by the spatial kernel in Equation S1.1c. It also means that the super-saturation s ( x , t ) is not defined on length scales comparable to the radii of droplets, and therefore no a microscopicfield in the most strict sense. . Fries et al. 3 The microscopic equations S1.1 can be derived from more fundamental descriptions (Vaillancourt et al., 2001,2002; Kumar et al., 2014, 2018; Perrin and Jonker, 2015). In the following we discuss the assumptions underlyingEquations S1.1.
S1.1.1 | Droplet inertia and settling
Equation S1.1d assumes that the droplets are small enough so that droplet inertia and settling are negligible (Gustavs-son et al., 2014; Gustavsson and Mehlig, 2016). To neglect droplet inertia is a reasonable approximation when theStokes number St = τ p / τ η is small enough. Here τ p = 2 ̺ w r α /( ̺ a ν ) is the Stokes time for a small water droplet inair with ̺ w ≫ ̺ a , and τ η = ( ν / ε ) is the Kolmogorov time (of the order of the turnover time of the smallest eddies).Settling is negligible at small settling numbers Sv = gτ p / u η , where u η = ( νε ) is the Kolmogorov velocity and g > isthe gravitational acceleration. For typical cloud conditions (dissipation rate per unit mass ε ∼ − m /s and viscosity ν ∼ . × − m /s), the Stokes number is of the order St ∼ − ( r α / µ m ) , and Sv ∼ − ( r α / µ m ) . In Figure 5 weanalyse observational data from Beals et al. (2015). Droplet-size distributions in Beals et al. (2015) suggest that veryfew of the measured droplets exceed r α = 6 µ m, corresponding to St = 1 . × − and Sv = 2 . × − . Droplet inertiaand settling can therefore be neglected for most of these droplets. Out of the five DNS studies referred to in the maintext, Kumar et al. (2013, 2014, 2018) incorporate droplet inertia and settling. Droplet inertia and settling can haveminute effects on the droplet-size distribution, despite quite large settling numbers, Sv ∼ . (Kumar et al., 2013). S1.1.2 | Buoyancy
Buoyancy is neglected in the momentum equation, Equation S1.1a. The neglected buoyancy terms read (Bannon,1996): gB ( x , t ) = g (cid:20) T − T s ( z ) T s ( z ) + 0 . q v − q ℓ (cid:21) . (S1.4)Here, B = B ( x , t ) is buoyancy, T = T ( x , t ) is temperature, and T s ( z ) is a static base profile of temperature as afunction of altitude z . Furthermore, q v = q v ( x , t ) is the water-vapour mixing ratio, and q ℓ = q ℓ ( x , t ) is the liquid-water mixing ratio. In a dry atmosphere ( q v = q ℓ = 0 ), buoyancy variations are caused only by vertical displacements.If T s ( z ) depends linearly on z , then a parcel that is neutrally buoyant at z = 0 and is displaced adiabatically to analtitude z has the buoyancy B = 1 T s ( ) (cid:18) Γ − d T s d z (cid:19) z + O( z ) . (S1.5)Here, Γ is the dry adiabatic lapse rate, the vertical temperature gradient d T / d z in a dry adiabatic atmosphere athydrostatic equilibrium (Rogers and Yau, 1989). One finds an upper bound for the buoyancy variations within a systemof spatial scale ℓ by substituting z = ℓ in Equation S1.5. An analogous constraint is derived in Dougherty (Dougherty,1961). Comparing buoyancy acceleration gB at a spatial scale ℓ to the inertial acceleration ( ε / ℓ ) at that scale, oneobtains the Dougherty-Ozmidov length scale (Grachev et al., 2015): ℓ B = ε (cid:16) g T s ( ) (cid:12)(cid:12)(cid:12)(cid:12) Γ − d T s d z (cid:12)(cid:12)(cid:12)(cid:12) (cid:17) − . (S1.6) J. Fries et al.
For typical atmospheric conditions [ ε ∼ − m /s , Γ ∼ − − K/m, d T s / d z ∼ − × − K/m (Dougherty, 1961)], onefinds ℓ B ∼
45 m. Under these conditions, buoyancy effects may matter on spatial scales larger than 45 m.In a moist atmosphere, buoyancy varies not only as a consequence of vertical displacements, but also becausedroplets condense and evaporate. The largest impact of droplet phase change occurs in regions where mixing be-tween cloudy and subsaturated air takes place (Vaillancourt and Yau, 2000). In such regions, buoyancy reduces asevaporating droplets absorb latent heat from the air. Grabowski (1993) show that sedimentation of droplets from7 cm wide stationary filaments of cloudy air can amplify buoyancy reductions due to droplet evaporation by a factor ∼ within 5 s. We argue however that such amplifications can not be expected in turbulent clouds, because buoy-ancy fluctuations at spatial scales ℓ ∼ cm are smoothened by turbulent dissipation at time scales τ ℓ ∼ ( ℓ / ε ) / that are smaller than 5 s, already at quite low dissipation rates. We therefore analyse buoyancy reductions causedby evaporation using the mixing curve of buoyancy (Siems et al., 1990), which does not take into account dropletsedimentation. This mixing curve gives an upper bound on the buoyancy reduction ∆ B that phase change induceswhen a volume fraction χ of cloudy air is mixed with a volume fraction − χ of non-cloudy air. This upper boundis achieved when the mixture equilibrates at a moist or dry steady state, and it is a definite function of the volumefraction χ , the temperatures T c and T e , the water-vapour mixing ratios q vc and q ve , and the liquid-water mixing ratio q ℓ c (Grabowski, 1993). The subscripts c and e indicate values for the cloudy and non-cloudy mixing substrates. Weestimate T c ∼ . K, q vc ∼ . · − , q ℓ c ∼ . · − , T e ∼ . , and q ve ∼ . · − for the cloud and cloudenvironment observed by Beals et al. (2015) (see Section S1.8). From these estimates we compute the upper bound ∆ B ∼ − for buoyancy reductions caused by droplet evaporation, using the mixing curve of buoyancy in Grabowski(1993). Comparing g ∆ B with the a typical inertial acceleration ( ε / ℓ ) at the spatial scale ℓ , we find that buoyancyaccelerations due to evaporation are smaller than inertial accelerations for ℓ < m.One caveat is that buoyancy may cause large-scale anisotropies in the fluid motion. This can introduce additionalcomplexity to droplet evaporation and mixing (Perrin and Jonker, 2015). Nevertheless, buoyancy is not a necessaryingredient in a model for these processes. Buoyancy is neglected in two of the DNS studies that we compare to(Kumar et al., 2012, 2013), as well as in the studies of Pinsky and Khain (2018) and Jeffery (2007). S1.1.3 | Temperature changes due to vertical air motion
Neglecting temperature changes due to vertical air motion is justified for systems that remain at a constant altitudeand that are not too large or too moist, as we now explain. Given the spatial scale ℓ of a system that remains at a fixedaltitude, the maximal temperature changes caused by vertical air motion are comparable to ℓ Γ . For the temperatureto change by 1 K due to vertical air motion, Γ ∼ − K/m dictates that we must have ℓ ∼ m.Temperature regulates droplet evaporation through supersaturation. Supersaturation is determined by the tem-perature T and the water-vapour density ̺ v according to Equations S1.2 and S1.3. For the values of ̺ v that implysaturation ( s = 0 ) at temperatures above 273.15 K, we find that a temperature change of 1 K causes the supersat-uration to change by 0.07 or less in magnitude. In a quite moist system, where the supersaturation is smaller thanthis in magnitude everywhere, temperature changes due to vertical air motion can be more important than mixingfor droplet evaporation. Since temperature changes due to vertical air motion increase with the spatial scale of thesystem, results obtained using Equations S1.1 may be inaccurate for systems at constant altitudes that are large andmoist. Temperature changes due to vertical air motion are often negelcted in studies of droplet evaporation (Kumaret al., 2012, 2013, 2012, 2014, 2018; Andrejczuk et al., 2006; Pinsky and Khain, 2018; Siewert et al., 2017; Jeffery,2007). . Fries et al. 5 S1.1.4 | Temperature and pressure dependence of the coefficient A We use a constant thermodynamic coefficient A . This parameter is given by A = (cid:20) (cid:18) Λ R v T − (cid:19) Λ ̺ w K a T + ̺ w R v TD v e s ( T ) (cid:21) − , (S1.7)where Λ is the latent heat of water vapor, K a is the thermal conductivity of air, and D v is the diffusivity of water vapour(Rogers and Yau, 1989). The parameter A decreases with pressure through D v , and increases with the temperature T (Rogers and Yau, 1989). The values that we use for A are computed from temperatures and pressures that arerepresentative of the systems that we consider. To use a constant value for A is usually a good approximation (Lanotteet al., 2009; Sardina et al., 2015; Lehmann et al., 2009; Kumar et al., 2012; Andrejczuk et al., 2009; Siewert et al.,2017; Jeffery, 2007; Pinsky and Khain, 2018; Devenish et al., 2016). The reason is that the temperature and pressuredependencies are quite weak. From Figure 7.1 of Rogers and Yau (1989) we see, for example, that A increases by 26%as the temperature increases from 10 ◦ C to 20 ◦ C and the pressure decreases from 100 kPa to 70 kPa. For observationsof a convective cloud at 4000 m altitude (Beals et al., 2015) we estimate the temperature variations to be smaller than ∼ ◦ C (see Section S1.8). These temperature variations cause A to change by less than 10%, which suggests that itis a good approximation to keep the coefficient constant in analyses of convective clouds at fixed altitude. S1.1.5 | Linearisation of the supersaturation profile
We treat the joint effect of temperature and water vapour at the level of a single field, the supersaturation field. Thisone-field treatment is obtained by linearising the supersaturations dependence on water-vapour mixing ratio q v andtemperature T , and by setting the diffusivities of water vapour and temperature equal to the same value κ . Withthe same diffusivity for temperature and water vapor, and with s given by a linear combination of q v and T , the twoseparate advection-diffusion equations for temperature and water-vapour mixing ratio (Vaillancourt et al., 2001) canbe concatenated into a single equation for s .We denote the heat capacity of dry air at constant pressure by c p and estimate ̺ a c p ∼ J/(m K) (Rogers andYau, 1989). With this estimate, Table 7.1 of (Rogers and Yau, 1989) implies that the diffusivity of water vapour roughlyequals the diffusivity of temperature for atmospheric conditions, so that it is sufficient to consider just one value forthe diffusivity. The supersaturation lineqrised around a reference state with temperature T r and water-vapour mixingratio q vr reads: s ( q v , T ) = s ( q vr , T r ) + ∂s∂q v (cid:12)(cid:12)(cid:12) q vr T r ( q v − q vr ) + ∂s∂T (cid:12)(cid:12)(cid:12) q vr T r ( T − T r ) . (S1.8)Using d q v = − d q ℓ and d T = ( Λ / c p ) d q ℓ , we find that the coefficient A in Equation S1.1b is given by A = 1 ̺ a (cid:16) ∂s∂q v (cid:12)(cid:12)(cid:12) q vr T r − Λ c p ∂s∂T (cid:12)(cid:12)(cid:12) q vr T r (cid:17) . (S1.9)In Figures 2 to 4 we compare our statistical-model results with DNS results (Kumar et al., 2014, 2018; Andrejczuket al., 2009) that are obtained using models where temperature and water vapour are treated as two separate fields. Inorder to match their non-linear supersaturation to our linear supersaturation given by Equation S1.8, we linearise theirsupersaturation around a saturated reference state ( T r , q vr ) . We use T r = 270 . K (Kumar et al., 2014), T r = 282 . K(Kumar et al., 2018) , and T r = 293 K (Andrejczuk et al., 2009), determining our values of q vr through Equations S1.2, J. Fries et al.
S1.3, and q v = ̺ v / ̺ a . For a convective cloud observed by Beals et al. (2015) we estimate the temperatures T c ∼ . and T e ∼ . , and the water-vapour mixing ratios q v c ∼ . × − and q v e ∼ . × − , wihtin and outside the cloud(see Section S1.8). Over these ranges, the absolue error of the supersaturation linearised at ( T r , q vr ) = ( T c , q v c ) is lessthan 1%. The relative error at ( T e , q v e ) is around 2%. S1.2 | INITIAL CONDITIONS
To allow for quantitative comparisons with the results of Kumar et al. (2013, 2012, 2014, 2018), we use the same(or almost identical) initial conditions: a slab of cloudy air inside a cubic domain of size L with periodic boundaryconditions. The cloudy air occupies a fraction χ of the domain, and it is surrounded by subsaturated air [Figure 1( b )].The initial droplet-size distribution is either monodisperse or Gaussian, with mean r µ and standard deviation σ . Whenwe compare to DNS from Kumar et al. (2013, 2012, 2014, 2018) in Figures 2 to 4, we use their initial droplet-sizedistribution, which is monodisperse ( σ = 0 ). When we analyse measurements of Beals et al. (2015), we use theGaussian droplet-size distribution of the undiluted cloud that corresponds to the measurements ( σ , , see SectionS1.5). The simulations of Andrejczuk et al. (2006) are for decaying turbulence with quite strong buoyancy effects,and a different initial geometry. Nevertheless, some of our results can be compared qualitatively to the results ofAndrejczuk et al. (2006), as explained in Section S1.7. Our analysis of the measurements of Beals et al. (2015) is basedon our initial condition. This analysis rests upon the assumption that our system models a structure of size L formedat the cloud edge, although its detailed composition and geometry can deviate from our slab-like initial conditions.Our initial conditions are homogeneous in two of the three coordinate axes. We take x to denote the coordinatein the spatially inhomogeneous direction. Initially the droplets are distributed uniformly and randomly within a three-dimensional cloud slab, the volume for which − χL / < x < χL / [Figure 1( b )]. The initial droplet-number densityof the slab is n = N /( χV ) , where V = L is the domain volume. To obtain the statistical-model results shown inFigures 2 to 4 we use the initial supersaturation profiles of Kumar et al. (2013, 2012, 2014, 2018): s ( x , ) = ( s c − s e ) exp (cid:20) − ζ (cid:16) xL (cid:17) ζ (cid:21) + s e . (S1.10)Here, s c > and s e < denote the initial supersaturation at the center of the slab and of the dry air [Figure 1( b )]. Theshape of the initial supersaturation profile is contained in the parameters ζ and ζ . The homogeneous mixing line ofBeals et al. (2015) passes through the top-right corner of the mixing diagram in Figure 5( a ). To be able to compare withthis homogeneous mixing line, we use a sharp initial supersaturation profile with s c = 0 when analysing the measureddata: s ( x , ) = for | x / L | < χ / , s e otherwise . (S1.11) S1.3 | NON-DIMENSIONAL EQUATIONS AND PARAMETERS
We non-dimensionalise Equations S1.1 using U = p TKE / and the large-eddy time τ L = TKE / ε , which is proportionalto L / U if the size of the largest eddies scale with L (Pope, 2000). In addition, we use the following quantities to non-dimensionalise: the droplet-number density of the slab [ n = N /( χV ) ], the (positive) supersaturation | s e | , and the . Fries et al. 7 initial volume radius r = N N Õ α =1 r α ( ) . (S1.12)We non-dimensionalise as follows: u ′ = u / U , x ′ = x /( Uτ L ) , t ′ = t / τ L , p ′ = p /( ̺ a U ) , x ′ α = x α /( Uτ L ) , s ′ = s /| s e | , r ′ α = r α / r , G ′ = G ( Uτ L ) , n ′ = n / n , s ′ c = s c /| s e | , σ ′ = σ / r , L ′ = L /( Uτ L ) and V ′ = V /( Uτ L ) . Dropping the primes,Equations S1.1 take the non-dimensional form: ∂ u ∂t + ( u · + ) u = − + p + 1 Re L + u , (S1.13a) + · u = 0 , (S1.13b) ∂s∂t + ( u · + ) s = 1 Re L Sc + s − Da s χV r α ( t ) s ( x α , t ) , (S1.13c)d x α d t = u ( x α ( t ) , t ) , (S1.13d)d r α d t = Da d s ( x α ( t ) , t ) . (S1.13e)Here r α ( t ) s ( x α , t ) denotes the average of G ( | x − x α ( t ) |) r α ( t ) s ( x α , t ) over the droplets: r α ( t ) s ( x α , t ) = 1 N N Õ α =1 G ( | x − x α ( t ) |) r α ( t ) s ( x α , t ) . (S1.14)The average r α ( t ) s ( x α , t ) is position dependent, it depends on number and sizes of the local droplets, and upon thelocal supersaturation s ( x α , t ) tied to these droplets. The dynamics in Equations S1.13 is identical to the dynamics inEquations S1.1, but non-dimensional.Equations S1.13 have the following non-dimensional parameters: the domain volume V , the initial cloud fraction χ , the turbulence Reynolds number Re L = TKE /( εν ) (Pope, 2000), the Schmidt number Sc = ν / κ [of order unitysince the diffusivities of both temperature and water vapour are roughly the same as the kinematic viscosity of air(Perrin and Jonker, 2015)], and the two Damköhler numbersDa s = τ L / τ s and Da d = τ L / τ d , (S1.15)where τ s and τ d denote the supersaturation relaxation time and the droplet evaporation time. These time scales aregiven by τ s = ( πA A ρ w n r ) − and τ d = r A | s e | . (S1.16)The time scales τ s and τ d can be very different. This can be understood by analysing their ratio R = τ s τ d = Da d Da s = | s e | πA ρ w n r = 2 | s e | A ̺ ℓ . (S1.17)Here, we introduce ̺ ℓ = 4 πr n ̺ w / as a scale for liquid water density (mass per volume), the liquid water density in J. Fries et al. a cloud with droplet-number density n and mean volume radius r . The scale for liquid water density determines R together with the supersaturation s e of the dry mixing substrate, and the thermodynamic coefficient A . For typicalcloud conditions, we find A ∼ m /kg using Equation S1.9.Consider, for example, a cumulus cloud with a typical (Rogers and Yau, 1989) liquid water density . g/m inthe bulk that contains premixed and slightly subsaturated air with supersaturation − . . Substituting s e = − . and ̺ ℓ = 0 . g/m into Equation S1.17, we find that R ∼ . for the local mixing processes in the bulk, so τ d is morethan ten times larger than τ s . At the cloud edge, by contrast, the liquid water density is lower, and mixing with dry airoccurs. Assuming, say, s e = − . and ̺ ℓ = 0 . g/m for a local mixing process yields R = 2 . , so τ d is smaller than τ s . As we continue to move away from the cloud core, ̺ ℓ tends to zero, and τ s / τ d tends to infinity. In conclusion,the supersaturation relaxation time and the droplet evaporation time for local mixing processes can differ by ordersof magnitude. S1.4 | STATISTICAL MODEL
The statistical model, Equations 3, rests upon on a probabilistic description of the microscopic Equations 1, in termsof two probability-density functions (PDF:s). The first PDF, denoted by F , describes the droplets, and the secondone, denoted by f , describes the air. The dynamics of Lagrangian fluid elements in Equations 3 constitutes a closedset of evolution equations for these PDF:s. The presence of droplets makes it necessary to consider two types offluid elements, one for each PDF. Similar extensions of single-phase PDF models have already been made to describecombustion problems, where a gas interacts with droplets or other dispersed particles (Jenny et al., 2012; Stöllingeret al., 2013).In the following, we derive the statistical model in terms of F and f , since this allows us to highlight similaritiesand differences to PDF models for single-phase flows (Pope, 1985). The first step is to define the two PDF:s. Theirexact evolution equations are not closed, as they contain conditional averages that are not known in terms of thePDF:s. We employ standard closures, designed to approximate the effects of acceleration and diffusive scalar fluxin single-phase flows (Pope, 1985) and multiphase combustion (Jenny et al., 2012; Stöllinger et al., 2013; Haworth,2010). The resulting closed set of equations constitutes a statistical model that describes the dynamics dictated bythe microscopic Equations 1. As a final step, the model is cast into the form of Langevin equations (Pope, 2000),invoking the concept of Lagrangian fluid elements. S1.4.1 | Probabilistic description
We describe the droplets by the joint PDF of droplets and air, F ( U , S , R , x ; t ) , which is the probability density of theevent E d : { x α ( t ) = x , r α ( t ) = R > , u ( x α ( t ) , t ) = U , and s ( x α ( t ) , t ) = S } (S1.18)for a droplet α . Our definition of probability densities of events is that of Pope (1985). The PDF F is Lagrangian , sinceit is a PDF evaluated along Lagrangian trajectories (Pope, 2000), the trajectories of the advected droplets. By defining F as a density in squared droplet radius r α ( t ) , instead of droplet radius r α ( t ) , we avoid complications that followfrom that d r α / d t diverges in the limit r α ( t ) → when we formulate the evolution equation of F below. Droplets that . Fries et al. 9 evaporate according to Equation S1.13e may adopt a zero radius, r α ( t ) = 0 . When r α ( t ) becomes zero the droplethas evaporated completely and no longer plays any role for the dynamics. Therefore only droplets with R > areconsidered in F . In this setup F is normalised to − P e ( t ) , reflecting that the total probability of F decreases asdroplets completely evaporate. The observables that we solve for are contained in F . We compute the droplet-sizedistribution as F r ( R ; t ) = 2 R ∫ F ( U , S , R , x ; t ) d U d S d x , (S1.19)and we compute the fraction of completely evaporated droplets as P e ( t ) = 1 − ∫ F ( U , S , R , x ; t ) d U d S d R d x . (S1.20)The air is described by the joint PDF of velocity and supersaturation, f ( U , S ; x , t ) , which is the probability densityof the event E a : { u ( x , t ) = U and s ( x , t ) = S } . (S1.21)As opposed to F , the PDF f is normalised to unity at all times. Note also that f is Eulerian , since it is a PDF of fluidproperties at a fixed position (Pope, 2000). PDF models for single-phase flows describe turbulent flows with one orseveral advected scalar fields (Pope, 2000). They solve for an Eulerian joint PDF of velocity and scalars. This PDF isanalogous to f , but may describe more than one advected scalar field. As in some combustion problems (Stöllingeret al., 2013; Jenny et al., 2012; Haworth, 2010), we must consider a dispersed phase (the droplets) that interacts witha fluid phase (the air). As a consequence, we need two PDF:s to describe the dynamics.Exact evolution equations of Lagrangian and Eulerian PDF:s are given in Pope (1985). To obtain the evolutionequation for F , we integrate a corresponding transition PDF in Pope (1985) over F ( U , S , R , x ; t = 0 ) , which isknown from the intial conditions. The microscopic equations S1.13 imply that the dynamics of F and f read: ∂ F ∂t + U i ∂ F ∂x i + Da d S ∂ F ∂R = − ∂∂U i (cid:28) Re L + u i − ∂p∂x i E d (cid:29) F − ∂∂S (cid:28) Re L Sc + s E d (cid:29) F − Da s χV ∂∂S D r α ( t ) s ( x α , t ) E d E F , (S1.22a) ∂f∂t + U i ∂f∂x i = − ∂∂U i (cid:28) Re L + u i − ∂p∂x i E a (cid:29) f − ∂∂S (cid:28) Re L Sc + s E a (cid:29) f − Da s χV ∂∂S D r α ( t ) s ( x α , t ) E a E f . (S1.22b)Here, we follow Pope (1985) and write the components of positions and velocities by ( x , x , x ) , ( u , u , u ) and ( U , U , U ) , and sum over repeated indices. Averages conditioned upon the events in Equations S1.18 and S1.21 aredenoted by h . . . | E d i and h . . . | E a i .We have imposed that a droplet that has evaporated completely must remain at r α ( t ) = 0 . Since F describesdroplets with positive radii, this dynamics implies a boundary condition for F at R = 0 , described by the probability current (Gardiner, 2009) J e ( U , S , x ; t ) = − Da d S lim R → + F ( U , S , R , x ; t ) if S < otherwise (S1.23)in the negative R direction. Here, the limit R → + is required, since R > in F . The probability current inEquation S1.23 is non-negative, J e ( U , S , x ; t ) ≥ , which reflects that droplets can reach r α ( t ) = 0 from above, butnot from below. The rate of complete droplet evaporation can be written in terms of the probability current asd P e d t = ∫ J e ( U , S , x ; t ) d U d S d x . (S1.24)In addition to the boundary condition at R = 0 , F and f inherit the spatial periodicity of the slab configuration[Figure (1 b )]. The variables U and S are not bounded, and require no boundary conditions. The initial conditions for F and f are determined by the initial supersaturation profile (Equation S1.10 or S1.11), the initial spatial distributionof droplets, the initial droplet-size distribution, and the PDF of velocity ϕ ( U ; x , t ) , which is the probability density ofthe event u ( x , t ) = U . We discuss the PDF ϕ ( U ; x , t ) below.The evolution equations S1.22, together with their initial and boundary conditions, constitute an exact descriptionof the dynamics of F and f . The terms on the left-hand sides of Equations S1.22a and S1.22b are in closed form, asthey can be computed in terms of the PDF:s. The terms on the right-hand sides are not in closed form, and must beapproximated in order to solve for the joint evolution of F and f . They correspond to the fluctuating accelerationand diffusive flux of supersaturation that a fluid element experiences. S1.4.2 | Closure
We follow (Pope, 2000) and model the fluctuating acceleration of fluid elements with velocities u ( t ) in statisticallystationary, homogeneous and isotropic turbulence as a three-dimensional Ornstein-Uhlenbeck process:d u = − C u d t + (cid:18) C (cid:19) d η . (S1.25)In this Langevin equation, C is a constant that depends on the Taylor-scale Reynolds number Re λ = p Re L (Pope,2011), and d η is the increment of a three-dimensional isotropic Wiener process (Pope, 2000). We model the Re λ -dependence of C according to the empirical fit of Pope (2011), C = 6 . /( Re − / λ ) . Equation S1.25 dictatesthat each velocity component of a fluid element evolves stochastically according to an independent one-dimensionalOrnstein-Uhlenbeck process with dimensional variance TKE / . Another consequence of Equation S1.25 is that theauto-correlation function ρ ( s ) = h u ( t ) u ( t + s ) i (dimensional t and s ) of the velocity component u ( t ) of a fluid elementdecays exponentially at the time-scale τ L /( C ) . In our non-dimensional units, the Ornstein-Uhlenbeck processeshave variance unity and auto-correlation time /( C ) . With Equation S1.25, the PDF of fluid-element velocity relaxesto a joint normal distribution at this time scale (Pope, 2000).The joint normal PDF of fluid-element velocity components and the exponential decay of ρ ( s ) dictated by Equa-tion S1.25 is observed empirically (Pope, 2000, 2011). In the limit Re λ → ∞ , C approaches 6.5, a measured valuefor the Kolmogorov constant that specifies how the second order Lagrangian structure function depends on the La-grangian auto-correlation time and the mean dissipation rate in high-Reynolds number turbulence (Pope, 2011). Con- . Fries et al. 11 sistently with the Kolmogorov theory of turbulence (Kolmogorov, 1941), the statistics pertaining to the motion of afluid element subject to Equation S1.25 is Reynolds-number independent at high Reynolds numbers. Equation S1.25does not take into account that the dissipation rate experienced by a fluid element fluctuates intermittently aroundthe mean dissipation rate ε , as predicted by Kolmogorov’s refined theory of turbulence (Kolmogorov, 1962). Thisintermittency can be taken into account using the refined Langevin model of Pope and Chen (1990), which is anextension of Equation S1.25. Since Equation S1.25 reproduces the observed variance and autocorrelation time offluid-element velocities, Equation S1.25 ensures that the empirically observed turbulent diffusion of passive scalarsis correctly described (Pope, 2000).Equation S1.25 for the velocity of fluid elements provides closure for the first terms on the right-hand sides ofEquations S1.22a and S1.22b. The joint PDF of droplets and air F is a probability density in the velocity of a fluidelement, a fluid element that moves together with one of the advected droplets. The closure for Equation S1.22atherefore follows directly from Equation S1.25. The joint PDF of velocity and supersaturation f is not a probabilitydensity in the velocity of a fluid element. The closure implied for Equation S1.22b is therefore somewhat more intricate.As explained by Pope (2000), the crucial step in the derivation of this closure is to consider the relation between(the Eulerian) f and a corresponding Lagrangian PDF that is conditioned upon the initial position of a fluid element.In an incompressible flow, f is obtained by integrating this Lagrangian PDF over initial positions of fluid elements.Equation S1.25 prescribes a Fokker-Planck equation for the Lagrangian PDF. Since this Fokker-Planck equation doesnot depend upon the initial position of the fluid element, a Fokker-Planck equation of the same form follows for f . Following Pope (2000), we conclude that Equation S1.25 implies that the first terms on the right-hand sides ofEquations S1.22a and S1.22b are given by two Fokker-Planck equations: − ∂∂U i (cid:28) Re L + u i − ∂p∂x i E d (cid:29) F = 34 C ∂∂U i U i F + 34 C ∂ F ∂U i ∂U i , and (S1.26a) − ∂∂U i (cid:28) Re L + u i − ∂p∂x i E a (cid:29) f = 34 C ∂∂U i U i f + 34 C ∂ f∂U i ∂U i . (S1.26b)In addition to providing closure to the first terms on the right-hand sides of Equations S1.22a and S1.22b, Equa-tion S1.25 sets the initial conditions of F and f . Statistical homogeneity and stationarity dictates that the PDF ofvelocity ϕ ( U ; x , t ) discussed above is independent of x and t . To obtain statistically stationary modeled turbulence itfollows from Equation S1.25 that ϕ ( U ; x , t ) must be joint normal in the velocity components.For the second terms on the right-hand sides of Equations S1.22a and S1.22b, we make the following closure: (cid:28) Re L Sc + s E d (cid:29) = (cid:28) Re L Sc + s E a (cid:29) = − C φ ( S − h s ( x , t ) i) . (S1.27)Here, C φ is an empirical constant, and h s ( x , t ) i is defined as the position-dependent average h s ( x , t ) i = ∫ Sf ( U , S ; x , t ) d U d S . (S1.28)Equation S1.27 ensures that the combined effect of molecular diffusion and turbulent mixing results in the correctdecay of the variance Var [ s ( x , t ) ] = (cid:10) s ( x , t ) (cid:11) − h s ( x , t ) i of supersaturation, at least under homogeneous conditions( ∂f / ∂x i = 0 ). Here, (cid:10) s ( x , t ) (cid:11) is the mean-squared supersaturation, given by Equation S1.28 with S replaced by S inthe integrand. Experiments and DNS support that, under homogeneous conditions, the variance of supersaturation in fully developed turbulence decays as dd t Var [ s ( x , t ) ] = − C φ Var [ s ( x , t ) ] , (S1.29)with C φ ∼ , independently of the molecular diffusivity (independently of Sc) (Pope, 2000). This decay rate is re-produced by Equation S1.27. Furthermore, Equation S1.27 respects the boundedness condition of a passive scalar(Pope, 2000). This means that, in the absence of evaporation or condensation, the supersaturation remains boundedbetween its maximal and minimal values. Other aspects of the decay of supersaturation fluctuations may not becorrectly reproduced by Equation S1.27 (Pope, 2000).For the third terms on the right-hand sides of Equations S1.22a and S1.22b, we make the following closure: D r α ( t ) s ( x α , t ) E d E = D r α ( t ) s ( x α , t ) E a E = h r ( t ) s ( x , t ) i , (S1.30)where h r ( t ) s ( x , t ) i is the position-dependent average h r ( t ) s ( x , t ) i = ∫ R S F ( U , S , R , x ; t ) d U d S d R . (S1.31)Equation S1.30 ensures that the liquid-water potential temperature analogue θ is conserved under the statistical-model. With this closure we replace the local average r α ( t ) s ( x α , t ) , which describes local neighborhoods of droplets,by h r ( t ) s ( x , t ) i , which is an expectation for a single droplet. As a consequence, the self-limiting nature of dropletevaporation may not be accurately described by the statistical model. To model local conditional averages as in Equa-tion S1.30 is a standard way to obtain closure for PDF models that describe combustion of particles in turbulence(Jenny et al., 2012; Stöllinger et al., 2013; Haworth, 2010).Inserting Equations S1.26, S1.27 and S1.30 into Equations S1.22, we obtain a model for the joint evolution of F and f : ∂ F ∂t + U i ∂ F ∂x i + 2 Da d S ∂ F ∂R = 34 C ∂∂U i U i F + 34 C ∂ F ∂U i ∂U i + 12 C φ ∂∂S ( S − h s ( x , t ) i) F − Da s χV ∂∂S h r ( t ) s ( x , t ) i F , (S1.32a) ∂f∂t + U i ∂f∂x i = 34 C ∂∂U i U i f + 34 C ∂ f∂U i ∂U i + 12 C φ ∂∂S ( S − h s ( x , t ) i) f − Da s χV ∂∂S h r ( t ) s ( x , t ) i f . (S1.32b)These equations are the statistical model, described as a closed set of evolution equations for F and f . We notethat F and f couple to each other through the averages h s ( x , t ) i and h r ( t ) s ( x , t ) i . This coupling is a consequenceof that droplets are affected by the supersaturation of the air, and that the supersaturation of the air is affected byevaporating droplets. S1.4.3 | Lagrangian dynamics
The PDF dynamics in Equations S1.32 can be cast into a set of equations that describe the dynamics of Lagrangianfluid elements (Pope, 2011). The dynamics of Lagrangian fluid elements with positions x ( t ) , velocities u ( t ) , and . Fries et al. 13 supersaturations s ( t ) that corresponds to Equation S1.32b is given by Equation S1.25, together with:d x d t = u , (S1.33a)d s d t = − C φ ( s − h s ( x , t ) i) − Da s χV h r ( t ) s ( x , t ) i . (S1.33b)Since droplets are advected, they move together with Lagrangian fluid elements. Therefore, some Lagrangian fluidelements coincide with droplets. We describe such fluid elements by the radii r ( t ) of the droplets that they coincidewith, in addition to their positions, velocities and supersaturations. Their dynamics corresponds to Equation S1.32a,and is given by Equations S1.25 and S1.33, together withd r d t = Da d s if r ( t ) > if r ( t ) = 0 . (S1.34)The relations between fluid-element dynamics and evolution equations of Eulerian and Lagrangian PDF:s arederived by Pope (1985). To conclude why Equations S1.25, S1.33, and S1.34 correspond to Equations S1.32, wegive a brief summary of the derivations for our statistical-model dynamics. First, one recognises that fluid elementswith droplets sample F , and that fluid elements without droplets sample f . Second, one formulates the Fokker-Planck equations for the two types of fluid elements. The Fokker-Planck equation for fluid elements with dropletsis Equation S1.32a, subject to the boundary condition S1.23. To obtain Equation S1.32b from the fluid elementswithout droplets, one integrates their Fokker-Planck equation over initial positions. Equation S1.32b then follows, asa consequence of the relation between Lagrangian and Eulerian PDF:s mentioned above. S1.4.4 | Computation of observables and statistical one-dimensionality
Since the PDF dynamics of our statistical model is implied by the dynamics of Lagrangian fluid elements, we solveEquations S1.32 for the PDF:s by evolving fluid elements according to Equations S1.25, S1.33 and S1.34. The aver-ages h s ( x , t ) i and h r ( t ) s ( x , t ) i are known in terms of the fluid elements, since the fluid elements sample F and f .In practice, these averages are computed using kernel estimates, as described by Pope (2000). Since h s ( x , t ) i and h r ( t ) s ( x , t ) i are known, Equations S1.25, S1.33 and S1.34 form a closed set of equations that govern the simultane-ous evolution of an ensemble of fluid elements. Since fluid elements with droplets sample F , we can compute thedroplet-size distribution and the fraction of completely evaporated droplets from them according to Equations S1.19and S1.20.To compute our results, we use the one-dimensional form of the statistical model shown in the main text. This ispossible, since the dynamics that we model is statisticallyone-dimensional. Statistical one-dimensionality follows fromthat the turbulence is statistically homogeneous and isotropic, and from that our initial conditions are one-dimensional.The direction of statistical inhomogeneity is x . The PDF:s F ( U , S , R , x ; t ) and f ( U , S ; x , t ) , as well as the averages h r ( t ) s ( x , t ) i and h s ( x , t ) i that they form, depend on the three-dimensional position x only through x . As explainedby Pope (1985), it is therefore sufficient to evolve only the x -components of fluid element positions and velocities. Inthe main text, we have taken the statistical one-dimensionality into account and replaced h r ( t ) s ( x , t ) i and h s ( x , t ) i by h r ( t ) s ( x , t ) i and h s ( x , t ) i . S1.4.5 | Averaging
The mean cubed radius of droplets and the volume average of supersaturation are required to compute the liquid-water potential temperature analogue θ . In the statistical model the mean cubed radius of droplets is given by: (cid:10) r ( t ) (cid:11) = 11 − P e ( t ) ∫ R F ( U , S , R , x ; t ) d U d S d R d x . (S1.35)The volume average of supersaturation is given by: h s ( t ) i = 1 V ∫ Sf ( U , S ; x , t ) d U d S d x . (S1.36) S1.5 | SIMULATION PARAMETERS
The parameters of our simulations are listed in Table S1. The parameters for Figures 2 to 4 are taken from the DNSof Kumar et al. (2012, 2014, 2018). The results of (Kumar et al., 2014, 2018) and Andrejczuk et al. (2006) shown inthese figures are not obtained using Equations 1, the microscopic equations that we model. The interpretation of thesimulation setups of Kumar et al. (2014, 2018) and Andrejczuk et al. (2006) in terms of our parameters is discussed inSections S1.1 and S1.7. Kumar et al. (2018) did not give their values of ζ , ζ , and s c , but we obtained them throughprivate communication (Kumar, 2019). For Figure 4, we use the values of C and L of Run 5 of Kumar et al. (2018),since we wanted to compare with DNS and are interested in large-Reynolds number regimes. All simulations use C φ = 2 , a common estimate for this empirical parameter (Pope, 2000).In Figure 5 we analyse measurements from the top panel of Figure 3 of Beals et al. (2015). In the supplementalmaterial of that paper, the authors show droplet-size distributions and droplet-number densities of the undilutedcloud air. The distribution of droplet radii that we use to analyse the top panel in Figure 3 of Beals et al. (2015) isvery close to Gaussian, as explained in Section S1.8. For our Figure 5, we fit a Gaussian to this distribution, and thefit is parameterised by the non-dimensional standard deviation σ = 0 . . We also use C = 6 . and L = 2 . forFigure 5. This value of L is an estimate of the constant of proportionality between the dimensional domain size and Uτ L at large Reynolds numbers. We compute this estimate by in two steps. First, we estimate the turbulent kineticenergy associated to a spatial scale ℓ by integrating the energy-spectrum function (Pope, 2000) from the wave number k = 2 π / ℓ : TKE = ∫ ∞ π / ℓ C ε k − d k = C (cid:18) εℓ π (cid:19) . (S1.37)Here, C is a Kolmogorov constant that has been measured to C = 1 . (Pope, 2000). Second, we use this value of theKolmogorov constant and find ℓUτ L = 4 π C − = 2 . . (S1.38)Since 2.28 is the constant of proportionality between an arbitrary spatial scale and Uτ L , it is the constant of propor-tionality between the dimensional domain size and Uτ L . We therefore conclude that our estimate in Equation S1.37gives L = 2 . and V = L = 11 . .Table S1 provides a set of independent parameters that, together with C φ = 2 , specify our simulations. Also . Fries et al. 15 shown are a few dependent parameters. These include the Damköhler number ratio R = Da d / Da s and the criticalDamköhler number ratio R c = − χ h s ( ) i , (S1.39)where h s ( ) i is the initial volume average of the supersaturation (defined below). Also included is the coefficient χ that determine h s ( ) i together with χ and s c . The initial volume average of the supersaturation can be written h s ( ) i = ( s c ) ( χ + χ ) − , (S1.40)a replacement that is done in the main text. For Figure 5, we have χ = s c = 0 , so h s ( ) i = χ − . For Figures 2 to 4, wehave a non-zero χ , since the initial supersaturation profile is smooth. Here, the parameters ζ and ζ that determinethe initial supersaturation profile are tied to specific values of χ . It is only for these values of χ that the initial super-saturation is approximately zero at x = ± χL / , the edges of the cloud slab [Figure 1( b )]. For a general χ , the volumeaverage of the initial supersaturation is given by Equation S1.40, with a value of χ that can be computed from ζ , ζ ,and s c . S1.6 | NUMERICS
We solve the statistical model using computer simulations that evolve an ensemble of Lagrangian fluid elements ac-cording to Equations 3. Fluid elements with droplets are initialised uniformly over the interval − χL / < x < χL / . Toensure that h s ( x , t ) i can be computed in the whole simulation domain, we initialise fluid elements without dropletsuniformly over − L / < x < L / (Pope, 2000). We solve for the supersaturation s ( t ) and squared radius r ( t ) of fluidelements using Euler’s method (Rade and Westergren, 2013). We use a dynamical time step that ensures that theabsolute errors in s ( t ) and r ( t ) at each time step are smaller than a given threshold. At each time step, the position x ( t ) and velocity u ( t ) of a fluid element are drawn from their exact joint distribution, which is conditioned on theposition and velocity before the time step (Gillespie, 1996). The averages h s ( x , t ) i and h r ( t ) s ( x , t ) i are computedfrom the fluid elements on a regular mesh. We obtain the mean-field values at the positions of fluid elements by linearinterpolation. The spatial variations of these averages reduce with time due to mixing, and we increase the computa-tional efficiency of our simulations by neglecting these spatial variations when they become very small. That h s ( x , t ) i and h r ( t ) s ( x , t ) i become essentially independent of x does not mean that the simulated system becomes spatiallyuniform. In particular, individual realisations of the simulated system may still exhibit supersaturation fluctuations thataffect the evaporation of droplets for some time.Figures 2 and 3 show droplet-size distributions and time series of P e ( t ) for three simulations. To conclude conver-gence for these simulations, we simply check that the results do not change as we vary the parameters and thresholdsthat control the numerics. Figures 4 and 5 show results from hundreds of simulations. To address the convergence ofthese simulations, we test the accuracy of our simulations for representative initial conditions and representative com-binations of Da d , R and χ . In each test, we vary the numerical parameters separately to exclude systematical errors.To estimate the statistical errors, we compute several independent realisations for each combination of the numericalparameters. These tests allow us to conclude that no systematical errors pertain to the results in Figures 4 and 5,and that the values of P ∗ e that we compute have relative errors that are less than 5 %, and/or absolute errors that aresmaller than − . The simulation results in Figures 4 and 5 are either P ∗ e or functions of P ∗ e , and we conclude that none of our conclusions in the main text are affected by numerical errors. S1.7 | PARAMETERS OF DNS IN FIGURE 4
To place DNS of Andrejczuk et al. (2006); Kumar et al. (2012, 2013, 2014, 2018) in our phase diagram (Figure 4), weextract their dimensional parameters that determine Da d and R / R c . The computed values are listed in Table S2 forKumar et al. (2012, 2013, 2014, 2018), and in Table S3 for Andrejczuk et al. (2006). S1.7.1 | DNS of Kumar et al. (2012, 2013, 2014, 2018)
It is straightforward to compute Da d , R , and R c for the simulations in Kumar et al. (2012, 2013, 2018) and threesimulations in Kumar et al. (2014), since they are for stationary turbulence with Lagrangian droplets, just as Equations 1.First, one casts the dynamics into the dynamics of Equations S1.1, using the simplifications described in Section S1.1.After that one non-dimensionalise as described in Section S1.3. The simulations of Kumar et al. (2012, 2013, 2014,2018) shown in Figure 4 are the simulations for which one can conclude if P ∗ e > % or if P ∗ e < %. S1.7.2 | DNS of Andrejczuk et al. (2006)
Figure 4 shows 42 simulations of Andrejczuk et al. (2006). Andrejczuk et al. (2006) reports of 58 simulations, but wecan only read off whether P ∗ e > % or not for 50 of them. Out of these 50 simulations, eight simulations fall outsidethe plot range of Figure 4. We therefore show 42 simulations of Andrejczuk et al. (2006) in Figure 4.The simulations of Andrejczuk et al. (2006) are for decaying turbulence and initial conditions in the form of ran-domly distributed cloudy filaments in dry air. The parameter Da d in Equation S1.15 is defined for the stationaryturbulence and the initial conditions of our simulations, so it has no direct counterpart in Andrejczuk et al. (2006). Tonevertheless discuss the simulations of Andrejczuk et al. (2006) we compute a time-scale ratio that incorporates thesame physics as Da d . We then place a simulation in Figure 4 by assuming that this time-scale ratio is equal to Da d .The time-scale ratio is computed as Da d = τ L / τ d , but with τ L taken as TKE / ε at the time when TKE decayed fastest inAndrejczuk et al. (2006). We estimate these values of
TKE and ε from plots in Andrejczuk et al. (2006).The coefficient A enters the dynamics multiplied by an unnamed factor γ A in Andrejczuk et al. (2006), and weread off A / γ A = 10 − m /s. We extract the the density ̺ a of air and the density ̺ w of pure liquid water Andrejczuket al. (2006, 2004). The expression for saturation pressure of water vapour is not given in Andrejczuk et al. (2006),so we can not compute their value of A using Equation S1.9. We therefore use the value A = 413 m /kg thatwe compute for Kumar et al. (2014). To nevertheless analyse DNS results of Andrejczuk et al. (2006) qualitativelyis justified since the temperatures and water-vapour densities in Andrejczuk et al. (2006) and Kumar et al. (2014)imply that thair values of A are similar. We compute s e as one minus the relative humidity of the non-cloudy initialfilaments, which we extract from tables in Andrejczuk et al. (2006). We also extract the volume fraction χ of cloudyair, and the liquid-water mixing ratio q ℓ of the initially cloudy filaments from these tables. We interpret the initiallycloudy filaments as saturated and possesing sharp edges, even though this is not stated explicitly in Andrejczuk et al.(2006), so that R c = − χ /[ ( χ − ) ] .Andrejczuk et al. (2006) do not use Lagrangian droplets. Instead, they employ a field description in which dropletsare binned into 16 size categories. The three initially populated bins are centered at droplet radii , . , and . µ m,and they contained 25%, 50%, and 25% of the liquid-water mixing ratio q ℓ . We compute the number density of . Fries et al. 17 droplets within a bin centered at the droplet radius r that contains a fraction ξ r of liquid-water mixing ratio as n r =3 ξ r q ℓ ̺ a /( πr ̺ w ) . By summing up the number densities for r = 8 , . , and . µ m, we find the droplet-number den-sity n of the initially cloudy filaments. We compute the initial volume radii of droplets as r = [ q ℓ ̺ a /( πn ̺ w ) ] / S1.8 | ESTIMATES USED TO ANALYSE DATA FROM Beals et al. (2015)
The black crosses in Figure 5 (a) are from the top panel of Figure 3 of Beals et al. (2015), and represent droplets mea-sured during a research flight through a convective cloud. Droplet-size distributions and thermodynamic conditionsare not given for this flight, so we estimate them using data from a typical flight detailed in the supplementary materialsof Beals et al. (2015), namely pass 2 of the research flight RF05.To estimate the parameters of our model, we need actual values for physical coefficients. We use
Λ = 2 . · J/kg, c p = 1005 J/(kg · K), and R v = 461 . J/(kg · K) from Rogers and Yau (1989). We also set the density of pure liquidwater to ̺ w = 1000 kg/m . Figure S3 of Beals et al. (2015) shows that pass 2 of RF05 traverses a cloud at a constantaltitude of 4000 m. We therefore assume the pressure of the International Standard Atmosphere (Cavcar, 2000) atthis altitude, p ∼ Pa. The precise value of the pressure is not important, since conclusions that rely on theestimates in this Section are unchanged if we assume a pressure that is 10 % lower or higher. We extract K a and D v atour assumed pressure from Table 7.1 of Rogers and Yau (1989). Furthermore, we assume that the saturation pressureof water vapour is given by Equation (2.12) in Rogers and Yau (1989). We also use the relation between potentialtemperature Θ , temperature T , and pressure p of Rogers and Yau (1989), Θ = T ( Pa / p ) . .The droplet-size distribution of undiluted cloud air traversed during pass 2 of RF05 is shown in Figure S7 of Bealset al. (2015). The distribution of droplet radii is very close Gaussian, normalised to the number-density of droplets. Wemake a Gaussian fit with mean . µ m, and standard deviation . µ m. The corresponding volume radius is r = 4 . µ m, which gives the non-dimensional standard deviation σ = 0 . . We use this value of σ to parameterise theinitial droplet-size distribution in the statistial-model simulations shown in Figure 5. The normalisation of the Gaussianfit gives the droplet-number density n = 764 cm − . We compute the liquid-water density of the undiluted cloud as ̺ ℓ c = ( πr / ) n ̺ w = 2 . · − kg/m . This value of ̺ ℓ c is consistent with the liquid-water density plot in Figure S4of Beals et al. (2015).We estimate temperatures from the liquid-water potential temperature shown in Figure S4 of Beals et al. (2015).Here, the liquid-water potential temperature varies between Θ ℓ c = 294 . K within the cloud and Θ ℓ e = 295 . Koutside the cloud. The liquid-water potential temperature is the same as the potential temperature Θ e outside thecloud, Θ e = 295 . K. We compute the temperature T e ∼ . outside the cloud, using Θ e and the pressure of theInternational Standard Atmosphere (Cavcar, 2000). We estimate the density of dry air as ̺ a = p /( R v T e ) ∼ . kg/m . This density gives the liquid-water mixing ratio q ℓ c = ̺ ℓ c / ̺ a ∼ . · − within the cloud. The potentialtemperature within the cloud can now be estimated as Θ c = Θ ℓ c + ( Λ / c p ) q ℓ c ∼ . K (Lamb and Verlinde, 2011),and the corresponding temperature estimate is T c ∼ . K. It is seen in Figure S4 of Beals et al. (2015) that thesupersaturation within the cloud is smaller than ∼
2% in magnitude. Assuming saturation within the cloud, the estimate T c ∼ . K gives the water-vapour mixing ratio q vc = e vs ( T c )/( ̺ a R v T c ) ∼ . · − . We extract the (negative)supersaturaion s e ∼ − . from Figure S4 of Beals et al. (2015), and estimate the water-vapour mixing ratio outsidethe cloud, q ve = e vs ( T e ) ( s e + 1 )/( ̺ a R v T e ) ∼ . · − .In our analysis of empirical data from Beals et al. (2015) we estimate τ s ∼ s based on our estimates of T c , q vc , and p . Equation S1.9 gives A ∼ m /kg when linearising the supersaturation around the temperature T c ∼ . Kand the water-vapour mixing ratio q vc ∼ . · − . We compute A ∼ · − m /s at the pressure p = 62000 Pa and temperature T c ∼ . K using Equation S1.7. Using the volume radius r and droplet-number density n estimatedabove, we find τ s ∼ s using Equation S1.16. . F r i e s e t a l. Table S1
TABLE S1 Simulation parameters. Our statistical-model simulations are completely specified by the two Damköhler numbersDa d and Da s , the constant C that regulates the auto-correlation time of fluid elements, the domain size L , the volume fraction χ of cloudy air, the constants ζ and ζ that determine the shape of the initial supersaturation profile in Figures 2 to 4, and thesupersaturation s c at the center of the initial cloud slab, together with the empirical constant C φ = 2 . From these parameters, wecompute the Damköhler number ratio R , the critical Damköhler number ratio R c , and the contribution χ to the initial volumeaverage of the supersaturation from the shape of the initial supersaturation profile. Simulation Independent parameters Dependent parametersDa d Da s C L χ σ ζ ζ s c R R c χ Figures 2 and 3, dry 2.44 0.968 5.22 2.96 0.428 0 4722 8 0.021 2.52 0.859 0.226Figure 2, moist 1.09 1.433 5.22 2.96 0.428 0 4722 8 0.021 0.760 0.859 0.226Figure 3, very moist 0.754 8.2 4.50 2.99 0.4 0 1410 6 0.1 0.092 0.683 0.154Figure 4 5E -3 -4E -3 -9E -2 -4E a ) 1E -2 -1E -2 -6E b ) 1E -2 -1E -1 -4E -2 -3E -2 Table S2
TABLE S2 Same as Table 2, but separately for each DNS from Kumar et al. (2012, 2013, 2014, 2018) in Figure 4.Non-dimensional parameters: Damköhler number Da d , Damköhler-number ratio R , critical ratio R c , and volumefraction χ of cloudy air. Dimensional parameters: domain size L , mean dissipation rate ε , and droplet-numberdensity n of the initially cloudy air. Non-dimensional parameters Dimensional parameters
L ε n Reference Da d Da s R R c χ [cm] [cm /s ] [cm − ]Table 2 in Kumar et al. (2012), Row 1 0.0075 0.0820 0.0916 0.6829 25.6 33.8 164Table 2 in Kumar et al. (2012), Row 2 0.0751 0.8204 0.0916 0.6829 25.6 33.8 164Table 2 in Kumar et al. (2012), Row 3 0.7511 8.2041 0.0916 0.6829 25.6 33.8 164Table 2 in Kumar et al. (2013), Row 3 0.3065 0.4185 0.7324 0.6829 25.6 33.8 164Table 2 in Kumar et al. (2013), Row 7 0.1362 0.6277 0.2170 0.6829 25.6 33.8 164Simulation S1 in Kumar et al. (2014) 2.4320 0.9660 2.5177 0.8377 51.2 33.8 153Simulation S2 in Kumar et al. (2014) 1.0809 1.4490 0.7460 0.8377 51.2 33.8 153Simulation S3 in Kumar et al. (2014) 0.6080 1.9319 0.3147 0.8377 51.2 33.8 153Run 1 in Kumar et al. (2018) 0.1191 0.5185 0.2296 0.9150 12.8 31.9 118Run 2 in Kumar et al. (2018) 0.1914 0.8333 0.2297 0.9155 25.6 34.6 118Run 3 in Kumar et al. (2018) 0.3227 1.4042 0.2298 0.9163 51.2 34.7 118Run 4 in Kumar et al. (2018) 0.5842 2.4509 0.2384 0.9503 102.4 32.1 113Run 5 in Kumar et al. (2018) 0.9122 4.0429 0.2256 0.9000 204.8 33.6 120 . Fries et al. 21 Table S3
TABLE S3 Same as Table 2, but separately for each DNS from Andrejczuk et al. (2006) in Figure 4.Non-dimensional parameters: Damköhler number Da d , Damköhler-number ratio R , critical ratio R c , and volumefraction χ of cloudy air. Dimensional parameters: domain size L , mean dissipation rate ε , and droplet-numberdensity n of the initially cloudy air. Non-dimensional parameters Dimensional parameters
L ε n Da d Da s R R c χ [cm] [cm /s ] [cm − ]S2a SI References
Andrejczuk, M., Grabowski, W. W., Malinowski, S. P. and Smolarkiewicz, P. K. (2004) Numerical simulation of cloud–clear airinterfacial mixing.
Journal of the atmospheric sciences , , 1726–1739.— (2006) Numerical simulation of cloud–clear air interfacial mixing: Effects on cloud microphysics. Journal of the AtmosphericSciences , , 3204–3225.— (2009) Numerical simulation of cloud–clear air interfacial mixing: Homogeneous versus inhomogeneous mixing. Journal ofthe Atmospheric Sciences , , 2493–2500.Bannon, P. R. (1996) On the anelastic approximation for a compressible atmosphere. Journal of the Atmospheric Sciences , ,3618–3628.Beals, M. J., Fugal, J. P., Shaw, R. A., Lu, J., Spuler, S. M. and Stith, J. L. (2015) Holographic measurements of inhomogeneouscloud mixing at the centimeter scale. Science , , 87–90.Cavcar, M. (2000) The international standard atmosphere (ISA). Anadolu University, Turkey , , 1–6.Devenish, B. J., Furtado, K. and Thomson, D. J. (2016) Analytical solutions of the supersaturation equation for a warm cloud. Journal of the Atmospheric Sciences , , 3453–3465.Dougherty, J. P. (1961) The anisotropy of turbulence at the meteor level. Journal of Atmospheric and Terrestrial Physics , ,210–213.Gardiner, C. (2009) Stochastic Methods , vol. 4. Springer Berlin.Gillespie, D. T. (1996) Exact numerical simulation of the Ornstein-Uhlenbeck process and its integral.
Physical Review E , ,2084–2091.Grabowski, W. W. (1993) Cumulus entrainment, fine-scale mixing, and buoyancy reversal. Quarterly Journal of the Royal Mete-orological Society , , 935–956.Grachev, A. A., Andreas, E. L., Fairall, C. W., Guest, P. S. and Persson, P. O. G. (2015) Similarity theory based on the Dougherty–Ozmidov length scale. Quarterly Journal of the Royal Meteorological Society , , 1845–1856.Gustavsson, K. and Mehlig, B. (2016) Statistical models for spatial patterns of heavy particles in turbulence. Advances inPhysics , , 1–57.Gustavsson, K., Vajedi, S. and Mehlig, B. (2014) Clustering of particles falling in a turbulent flow. Physical Review Letters , ,214501.Haworth, D. C. (2010) Progress in probability density function methods for turbulent reacting flows. Progress in Energy andCombustion Science , , 168–259.Jeffery, C. A. (2007) Inhomogeneous cloud evaporation, invariance, and Damköhler number. J. Geophys. Res. , , D24S21.Jenny, P., Roekaerts, D. and Beishuizen, N. (2012) Modeling of turbulent dilute spray combustion. Progress in Energy andCombustion Science , , 846–887.Kolmogorov, A. N. (1941) The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Doklady Akademii Nauk SSSR , , 299–303.— (1962) A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluidat high Reynolds number. Journal of Fluid Mechanics , , 82–85.Kumar, B. (2019) Private Communication. . Fries et al. 23 Kumar, B., Götzfried, P., Suresh, N., Schumacher, J. and Shaw, R. A. (2018) Scale dependence of cloud microphysical responseto turbulent entrainment and mixing.
Journal of Advances in Modeling Earth Systems , , 2777–2785.Kumar, B., Janetzko, F., Schumacher, J. and Shaw, R. A. (2012) Extreme responses of a coupled scalar–particle system duringturbulent mixing. New J. Phys. , , 115020.Kumar, B., Schumacher, J. and Shaw, R. A. (2013) Cloud microphysical effects of turbulent mixing and entrainment. Theor.Comp. Fluid Dyn. , , 361–376.— (2014) Lagrangian mixing dynamics at the cloudy–clear air interface. Journal of the Atmospheric Sciences , , 2564–2580.Lamb, D. and Verlinde, J. (2011) Physics and Chemistry of Clouds . Cambridge University Press.Lanotte, A. S., Seminara, A. and Toschi, F. (2009) Cloud droplet growth by condensation in homogeneous isotropic turbulence.
Journal of the Atmospheric Sciences , , 1685–1697.Lehmann, K., Siebert, H. and Shaw, R. A. (2009) Homogeneous and inhomogeneous mixing in cumulus clouds: Dependenceon local turbulence structure. Journal of the Atmospheric Sciences , , 3641–3659.Perrin, V. E. and Jonker, H. J. J. (2015) Lagrangian droplet dynamics in the subsiding shell of a cloud using direct numericalsimulations. Journal of the Atmospheric Sciences , , 4015–4028.Pinsky, M. and Khain, A. (2018) Theoretical analysis of mixing in liquid clouds – Part IV: DSD evolution and mixing diagrams. Atmos. Chem. Phys. , , 3659–3676.Pope, S. B. (1985) PDF methods for turbulent reactive flows. Progress in Energy and Combustion Science , , 119–192.— (2000) Turbulent Flows . Cambridge University Press.— (2011) Simple models of turbulent flows.
Physics of Fluids , , 1–20.Pope, S. B. and Chen, Y. L. (1990) The velocity-dissipation probability density function model for turbulent flows. Physics ofFluids , , 1437–1449.Rade, L. and Westergren, B. (2013) Mathematics Handbook for Science and Engineering . Springer Science & Business Media.Rogers, R. R. and Yau, M. K. (1989)
A Short Course in Cloud Physics . Pergamon Press.Sardina, G., Picano, F., Brandt, L. and Caballero, R. (2015) Continuous growth of droplet size variance due to condensation inturbulent clouds.
Physical Review Letters , , 184501.Siems, S. T., Bretherton, C. S., Baker, M. B., Shy, S. and Breidenthal, R. E. (1990) Buoyancy reversal and cloud-top entrainmentinstability. Quarterly Journal of the Royal Meteorological Society , , 705–739.Siewert, C., Bec, J. and Krstulovic, G. (2017) Statistical steady state in turbulent droplet condensation. Journal of Fluid Me-chanics , , 254–280.Stöllinger, M., Naud, B., Roekaerts, D., Beishuizen, N. and Heinz, S. (2013) PDF modeling and simulations of pulverized coalcombustion–Part 1: Theory and modeling. Combustion and Flame , , 384–395.Vaillancourt, P. A. and Yau, M. K. (2000) Review of particle-turbulenceinteractions and consequences for cloud physics. Bulletinof the American Meteorological Society , , 285–298.Vaillancourt, P. A., Yau, M. K., Bartello, P. and Grabowski, W. W. (2002) Microscopic approach to cloud droplet growth bycondensation. Part II: Turbulence, clustering, and condensational growth. Journal of the Atmospheric Sciences , , 3421–3435.Vaillancourt, P. A., Yau, M. K. and Grabowski, W. W. (2001) Microscopic approach to cloud droplet growth by condensation.Part I: Model description and results without turbulence. Journal of the Atmospheric Sciences ,58