Lagrangian supersaturation fluctuations at the cloud edge

Evaporation of cloud droplets accelerates when turbulence mixes dry air into the cloud, affecting droplet-size distributions in atmospheric clouds, combustion sprays, and jets of exhaled droplets. The challenge is to model local correlations between droplet numbers, sizes, and supersaturation, which determine supersaturation fluctuations along droplet paths (Lagrangian fluctuations). We derived a statistical model that accounts for these correlations. Its predictions are in quantitative agreement with results of direct numerical simulations, and it explains the key mechanisms at play.

When dry air is mixed into a cloud, water droplets at the cloud edge evaporate.This causes the droplet-size distribution to broaden [1], a pre-requisite for rain formation [2].Similar processes occur in combustion sprays [3][4][5], and for respiratory droplets in exhaled jets of air [6][7][8][9][10].In these systems, an essential physical ingredient is that evaporating droplets saturate the surrounding air.But the subtle coupling between phase change and turbulent mixing at widely separated turbulent scales [11] makes it difficult to predict the local supersaturation and, as a consequence, droplet-size distributions.In this work, we study the interplay between these processes and their influence on droplet growth, focusing on the parameter regime relevant to the edge of a cloud.
The distribution of supersaturation at droplet positions -the Lagrangian distribution -can be strongly non-Gaussian.Direct numerical simulations (DNS) of transient mixing of a three-dimensional slab of cloudy air with the surrounding dry air [Fig.1(a)] show exponential tails [12][13][14][15].Non-Gaussian supersaturation fluctuations are also seen in cloud-chamber experiments [16][17][18].Without phase change, a passive scalar field mixed by turbulence can exhibit non-Gaussian concentration fluctuations in the presence of a mean scalar gradient [19][20][21].Without the mean gradient, however, the steady-state distribution is essentially Gaussian [13,14,22,23].Non-Gaussian tails may appear in transient mixing [24], but must eventually disappear in a homogeneous system.
Whether the tail of the Lagrangian supersaturation distribution is Gaussian or not makes a significant difference, because the tail determines how rapidly certain droplets evaporate, and thereby influences the sizes of the remaining droplets, and thus the droplet-size distibution in unkown ways.
The key question is thus how droplet phase change affects the Lagrangian supersaturation distribution.When phase change is frequent and rapid, no consideration of passive-scalar mixing can explain the non-Gaussian relaxation of the Lagrangian supersaturation distribution.
Large-eddy simulations of droplet growth by condensation in a cloud chamber [17] show an anticorrelation between supersaturation s(x, t) and the local dropletnumber density n(x, t) in the steady state: in regions with many droplets, the air is strongly subsaturated (very negative values of s), while it is less so in regions with few droplets (less negative s).It is plausible that this effect may change the tails of the Lagrangian supersaturation distribution, but existing stochastic models [25][26][27][28][29][30][31][32] cannot explain this, because they do not describe how s(x, t) is affected by phase change locally.
We derived a statistical model that describes transient Lagrangian supersaturation fluctuations developing from an initial inhomogeneity, such as the configuration shown in Fig. 1(a), used in earlier DNS studies [12][13][14][15]30].We show here that the model captures the key mechanisms that determine the shape of the Lagrangian supersaturation distribution.First, when mixing occurs on time scales much shorter than phase change, non-Gaussian tails form only during the initial transient, and large-time relaxation is characterised by a Gaussian distribution of s(x, t).Second, in the opposite limit of rapid phase change, the distribution is no longer Gaussian.Strong phase change drives the mean of the distribution close to its upper bound, while the variance decays more slowly.The distribution is squeezed and becomes non-Gaussian.This is reflected in a strong positive correlation between n and s, formed because phase change tends to establish saturation in regions with many droplets.This mechanism -the opposite of the effect described in Ref. [17] -explains the tails in the Lagrangian supersaturation distribution described in earlier studies [12][13][14][15].Beyond this qualitative explanation, the model predicts supersaturation distributions that are in excellent quantitative agreement with earlier DNS results.The key to success is that the model inherits its supersaturation dynamics from first principles, rather than imposing an external driving resulting in a Gaussian steady state [25,26,29].
We start from simplified microscopic equations [30] governing droplet evaporation in turbulent flow: Eq. (1a) is the Navier-Stokes equation for the incompressible fluid-velocity field u(x, t), where ϱ a is the mass density of air, and ν its kinematic viscosity.Eq. (1b) describes supersaturation s = q v /q vs −1, with water-vapour mixing ratio q v = ϱ v /ϱ a , the ratio of the mass densities of vapour and air, and κ is the diffusivity of supersaturation.Further, C d (x, t) = 4π 3 ρ w n(x, t) d dt r 3 is the local rate of change of droplet mass, averaged over all droplets in the vicinity of x with droplet radius r.Here ρ w is the liquid-water density, and n(x, t) is the droplet-number density.Eqs.(1c) state that the droplets follow the flow, and how the droplet radius r changes [33].
We emphasise that the dynamics (1) is transient, and tends towards a well-mixed steady state.The variance σ 2 s (t) of supersaturation fluctuations tends to zero due to dissipation (with scalar dissipation rate ε s ≡ 2κ⟨|∇s| 2 ⟩) and phase change: Here, the averages are over a given microscopic confiugration at time t.
How the steady state is approached depends on the non-dimensional parameters of the problem.We nondimensionalise Eqs.(1,2) as follows [30]: time, velocities, and positions with the large-eddy turnover time τ L = k/ε (with turbulent kinetic energy k and kinetic dissipation rate ε), and the turbulent r.m.s velocity u 0 = 2 k/3; supersaturation with |s e |, where s e is the initial subsaturation of the dry air outside the cloud; droplet radii with the initial average droplet radius r 0 .In the limit of large Reynolds number, the following non-dimensional parameter remain: the Damköhler numbers Da d = τ L /τ d and Da s = τ L /τ s (with supersaturation relaxation time τ s and droplet-evaporation time τ d defined as in [30]), and the volume fraction χ of cloudy air [Fig.1(a)].The Schmidt number Sc= ν/κ is of order unity [34].
Earlier attempts to analyse the process were based on statistical models of mixing and evaporation that describe droplet evaporating in direct response to a spatially inhomogeneous mean field given by an ensemble average, ⟨s(x, t)⟩ [31,32,35].A more sophisticated model [30] accounts for how droplet-phase change is affected by Lagrangian supersaturation fluctuations, but still assumes that they decay exponentially towards the mean.Both types of models explain how the extent of complete droplet evaporation depends on the Damköhler numbers, but fail to reproduce the far tail of cloud-droplet size distribution obtained in DNS [12][13][14][15].A likely reason is that these models underestimate the magnitude of supersaturation fluctuations.A further shortcoming is that these models assume exponential relaxation of supersaturation to ⟨s(x, t)⟩.As a consequence, they fail to reproduce the passive-scalar limit [36], namely Gaussian supersaturation fluctuations with exponentially decaying variance.
Model.The mapping-closure approximation [42,43] describes how the shape of a passive-scalar distribution changes as the scalar is mixed in turbulence.The approximation relies only on one-point statistics.Correlations are not needed.As a consequence, the approximation does not predict the speed of the mixing process, but yields accurate and robust predictions for the sequence of shapes of the distribution.Therefore it is ideally suited for our purposes.The mapping closure for Eulerian supersaturation fluctuations starts from where λ(t) a time-dependent length scale, ξ is a spatially smooth random Gaussian field with mean zero and unit variance, and X is the time-dependent mapping from ξ(x/λ(t)) to s(x, t).Inserting (3) into (1b), one obtains In comparison to the original model [42,43], Eq. (4) contains the phase-change term ⟨C d |s =X⟩ = X⟨rN/V |s = X⟩, with droplet number N and spatial volume V .We approximate this term by a mean-field decoupling of the conditional average (see SM [44] for details) The factor φ(t) = κ/[u 2 0 τ L λ 2 (t)] in Eq. ( 4) is the nondimensional relaxation rate of the Eulerian distribution.How φ(t) changes as a function of time is determined by processes at both small and large length scales, as the following argument shows.For passive-scalar mixing, the scalar variance decays exponentially in the selfsimilar regime [36] In this case, φ(t) approaches the steady-state value, φ * ∼ C ϕ /2.The steady state emerges as a balance between the scalar variance cascading towards large wave numbers and rapid dissipation at large wave numbers.In physical dimensions, the steady-state length scale λ * = [2κk/(C ϕ ε)] 1/2 equals λ T (5C ϕ Sc) −1/2 where λ T = (10νk/ε) 1/2 is the Taylor microscale.In other words, both large-scale mixing and small-scale diffusion matter.
shows how φ(t) evolves as a function of t.For passivescalar mixing, the predicted plateau at φ * is approached after two large-eddy turnover times, at t ≈ 2. With phase change, φ(t) is larger, corresponding to smaller λ(t).This is consistent with the notion that phase change generates supersaturation gradients by driving the air towards saturation where droplets exist, while subsaturated regions without droplets remain subsaturated.
To obtain the Lagrangian supersaturation fluctuations, Pope [43] suggested to use a Langevin equation for ξ(t) dξ = −R(t)ξdt+[2R(t)] 1 2 dv , where dv is the increment of a Gaussian random process, and to compute the supersaturation as s(t) = X[ξ(t), t].The Langevin equation ensures that the distribution of ξ(t) relaxes to a normalised Gaussian, and the function X(η, t) maps this Gaussian to the Eulerian supersaturation distribution.This ensures that the Lagrangian supersaturation distribution relaxes to the Eulerian one, as required.We set R(t) = Cφ(t), where C is a constant.This is motivated -at least for a passive scalar -by the fact that supersaturation fluctuations due to turbulent mixing experienced by a fluid element reflect the diffusive term κ∇ 2 s in Eq. (1b), and that the fluctuations of this term are proportional to φ(t) under the mapping closure.Comparison with DNS shows that R(t) = Cφ(t) works very well for the first two large-eddy turnover times, for Da s up to 8.0.We find that C decreases as Da s increases (see SI [44]), because phase change tends to maintain saturation in regions with droplets.For times much larger than the large-eddy turnover time, the precise form of R(t) does not matter because the Lagrangian distribution has almost relaxed to the Eulerian one.
lines), compared with earlier DNS results [13] (dashed lines).Shown are two cases: small and large Damköhler numbers.We see that the statistical model reproduces the DNS results quantitatively.
For small Da s , the effect of phase change is small at short times, the distribution is close to that of a passive scalar (blue dashed lines).Supersaturation behaves essentially like a passive scalar during the first largeeddy turnover time, t ∼ 1.At later times, droplet-phase change matters more, but its effect is straightforward, it causes the peak of the distribution to shift somewhat compared to the distribution for Da s = 0, to less negative values of s, while supersaturation fluctuations are still approximately Gaussian [Fig.2(a)].For larger values of Da s , by contrast, the evolution of the supersaturation distributions looks very different [Fig.2(b)].The distribution remains non-Gaussian at large times.
To pin down the precise mechanism, we followed the droplet-number density n(t) and the local supersaturation s(t) for different fluid parcels in DNS.The results are summarised in Fig. 3 which shows the conditional average of the local-droplet number density conditional on the surrounding supersaturation for the same parameters as in Fig. 2. For small Damköhler numbers, the average does not change much during the time shown, it is still strongly influenced by the initial condition.For large Da s , by contrast, the average changes rapidly, and the positive correlation between n(t) and s(t) increases significantly.The statistical model captures this very well.The mechanism is simply that a parcel containing many droplets cannot remain sub-or supersaturated for long, because phase change drives the air quickly towards saturation when Da s is large.As a consequence, parcels with few droplets tend to have much more negative values of s, compared with a parcel with small Da s .The strong suppression of the conditional average ⟨n(t)|s(t) = s⟩ at large Da s explains how the non-Gaussian tails evolve in Fig. 2(b): the left tail of P L (s; t) disappears quickly as time increases, because droplets saturate their surroundings, and therefore fewer of them experience very dry air.
For small Da s ‚ the mapping closure results in Gaussian relaxation.Phase change causes non-Gaussian tails.In order to describe these tails, it is necessary to condition C d on supersaturation, Eq. ( 5).The conditioning also ensures that the supersaturation fluctuations remain bounded, as they must because neither phase change nor mixing can turn subsaturated into supersaturated air.
Discussion.We begin by discussing in more detail, how our results relate passive-scalar mixing.Eswaran and Pope [37] described the shape change of the Eulerian passive-scalar distribution as a function of time.The initial condition [Fig.1(a)] dictates that the Eulerian distribution is, initially, the sum of two narrow peaks, located at s = s c and s e .It relaxes first to a U-shaped form.The left tail of the Lagrangian supersaturation distribution reflects how the Eulerian peak at s = s c broadens.At large times, the Eulerian U-shaped distribution relaxes to a Gaussian.
Our model predicts the same for small but not negligible Da s , with one important difference: phase change causes the mean of the Lagrangian supersaturation distribution to shift to the right [Fig.2(a)], while mixing causes the distribution to narrow, remaining approximately Gaussian.In this case, the mapping is approximately given by X(η, t) = σ s (t)η + µ s (t).Inserting this into Eq.( 4) and assuming passive-scalar relaxation of the width, dσ s /dt = −C ϕ σ s /2, yields dµ s /dt = −Da s µ s .So the standard deviation decays more rapidly than the mean for small Da s , consistent with Gaussian relaxation.At large Da s , the time evolution of the Lagrangian supersaturation distribution is strongly affected by phase change, resulting in persistent non-Gaussian tails [Fig.2(b)].Our model explains why the Lagrangian supersaturation distributions relax so differently for small and large Da s .Rapid phase change quickly drives the mean of the distribution towards the upper bound of the supersaturation distribution.As a result, the distribution is squeezed towards s = 0, thus preventing a Gaussian from forming.The distribution is bounded because subsaturated air can not obtain a positive supersaturation through droplet evaporation.Therefore saturation (s = 0) constitutes an upper bound for the Lagrangian supersaturation fluctuations.
We contrast our results with those of Prabhakaran et al. [17].They found a negative correlation between n(t) and s(t) in LES designed to model droplet condensation in a cloud chamber.Their system is statistically stationary, but the statistical model highlights the mechanism leading to their findings.In their case, the air is saturated or supersaturated, so droplets tend to grow by condensation.The resulting drive towards saturation gives rise to a positive correlation between n(t) and s(t).
The model describes not only the Lagrangian supersaturation fluctuations quantitatively, but also the dropletsize distribution (not shown, see SI [44]).Our earlier model [30] yielded qualitative but not quantitative agreement, highlighting the importance of correlations between n and s.
We recall that Damköhler numbers in atmospheric clouds tend to be large, simply because the relevant length scale L is large, causing large τ L .It is tempting to argue that the persistent left tail of the Lagrangian supersaturation distribution at large Damköhler numbers in Fig. 2(b) is more representative of atmospheric relaxation than the Gaussian relaxation at small Damköhler numbers.However, the local supersaturation field around individual droplets is hard to observe in situ.Observations resolving supersaturation at larger scales, of the order of one metre, indicate Gaussian distributions [45], but better resolved laboratory measurements reveal skewed distributions [46], as predicted by our model.
Here we analysed moist systems, with Da d /Da s ∝ (ρ w n 0 r 3 0 ) −1 ∼ 0.1.We expect the present model to apply equally well to dry clouds where complete droplet evaporation occurs frequently, but we have not yet explored this regime.
Villermaux et al. [5] measured the joint dynamics of vapour and droplets in a dense acetone spray.They analysed vapour concentrations for different flow configurations considering the limit of large droplet-number density n 0 and large Damköhler numbers, where the droplets in the spray prevent each other from evaporating, but evaporate instantenously in dry air.In this limit, correlations between n(x, t) and s(x, t) are extreme, and it remains to be seem whether they can be captured by our model.More generally, it is of interest to compare the accuracy of the present mapping-closure model with predictions of the linear-eddy model [47], where tur-bulent stretching and folding is represented by a onedimensional map [48][49][50].
Conclusions.We derived a statistical model for the transient supersaturation fluctuations around droplets near the cloud edge, where turbulence mixes dry with cloudy air, causing the droplets to evaporate.The model explains the key mechanisms determining Lagrangian supersaturation fluctuations, and its predictions are in quantitative agreement with earlier DNS studies of droplet evaporation at the cloud edge [12][13][14][15]30].This advance became possible because the model describes supersaturation dynamics and the local coupling due to phase change using a mapping closure, which is known to yield quantitative results for passive-scalar mixing.At the same time, the model is simple enough so that it can be used to resolve sub-grid scale effects in large-eddy simulations with high precision.
We stress that the present model, unlike earlier statistical models, accounts for local correlations between droplet numbers and supersaturation.This opens the possibility to model the dynamics of denser turbulent aerosols, such as industrial sprays.
JF was supported by grants from the Knut and Alice Wallenberg (KAW) Foundation (no.2014.0048) and Vetenskapsrådet (VR), no.2021-4452.G. Sardina acknowledges support from VR (grant no.2022-03939).The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at NSC and PDC partially funded by the Swedish Research Council through grant agreements no.2022-06725 and no.2018-05973.The collaboration of AP and BM during the Multiphase 22 program at the KITP was supported in part by the National Science Foundation under grant no.NSF PHY-1748958.

I. DESCRIPTION OF DIRECT NUMERICAL SIMULATIONS
In this Section, we describe the direct numerical simulations (DNS) used to obtain the data shown in Figs 1(b), 2, and 3.The governing equations (1) were solved in an Eulerian-Lagrangian approach: the Eulerian momentum and supersaturation equations [Eqs.(1a,b)] were solved using a pseudo-spectral DNS code coupled with the equation of motion for the droplets [Eq.(1c)].The numerical solver has been used in previous works on cloud droplet evaporations in turbulence [1,2].The Eulerian equations were solved in Fourier space with the nonlinear terms calculated in physical space using a standard 2/3 rule for de-aliasing.Direct and inverse fast Fourier transforms were employed.The equations were advanced in time with a low-storage third-order Runge-Kutta method with the diffusive terms analytically calculated, while an Adam-Bashforth scheme was employed for the nonlinear terms.Likewise, a lowstorage third-order Runge-Kutta temporal scheme was used to integrate the equations of motion for the droplets.
Air velocity and supersaturation at the droplet positions were calculated with a linear interpolation scheme.Linear extrapolation was used to evaluate the C d on the Eulerian grid.The simulation configurations and parameters were from Kumar et al. [3].More precisely, we used the parameters given in the middle line and the third line of their Table 2. Following Kumar et al. [3], the simulation domain was taken to be a cube with a length of L x = L y = L z = 25.6 cm, where homogeneous turbulence was forced to obtain a mean turbulent kinetic dissipation rate of ε = ν⟨|∇u| 2 ⟩ = 33.75cm 2 s −3 , or, equivalently, a Kolmogorov scale of η = (ν 3 /ε) 1/4 = 1 mm.The turbulent kinetic energy evaluated to k = 3u 2 0 /2 = 49 cm 2 s −2 , where u 0 = u x,rms cm s −1 is the turbulent r.m.s.velocity.We employed 512 3 collocation points to solve the equation on the Eulerian grid.The turbulence was forced to maintain a statistically steady state.In Fourier space, the forcing was is the Fourier transform of the velocity field, and the wave vectors are of the form as well as all permutations of the components.
All simulations had an initial droplet radius r 0 = 20µm.Initial supersaturation and droplet positions were chosen as illustrated in Fig. 1(a).In particular, the cloud slab had an initial droplet-number density of n 0 = 164 cm −3 , supersaturation s c = 0.02, and volume fraction χ = 0.4, while the dry region contained no particles initially, and had supersaturation s e = −0.2.Other relevant parameters of the two DNS runs were: supersaturation Schmidt number ν/κ = 0.7, A 2 = 1/(ϱ a q vs ) = 265 m 3 /kg (ϱ a is the mass density of dry air, and q vs is the saturation value of the vapor mixing ratio at 270 K), and A 3 = 50.7µm 2 /s −3 in Fig. 2a).Following Kumar et al. [3], we increased A 3 by a factor of ten for Fig. 2(b), in order to simulate inhomogeneous mixing conditions.The values of the Damköhler numbers are summarised in Table S1.

II. DETAILS OF MAPPING CLOSURE
In this Section, we derive the mapping closure with phase change, Eq. ( 4), and explain how we closed the phasechange term.The description of the mapping closure summarised below differs from Ref. [4] only in details.The main difference is that we assume that the physical supersaturation field is given by s(x, t) = X[ξ(x/λ(t)), t] [Eq.( 3)], rather than working with a surrogate field.
We start from the evolution equation for the Eulerian cumulative distribution function F (S; t) = Prob(s(x, t) < S).From Eq. (1b), one obtains an equation for F : Here and in the remainder of this Section, we use dimensional units, as in Eq. ( 1).In order to derive Eq. ( 4), we need to insert the ansatz s(x, t) = X[ξ(x/λ(t)), t], as described by Pope [4].We assume that X(η, t) is a non-decreasing function of η.Now recall that ξ is a standardised Gaussian field.This implies where G(η) is the probability density of a standardised Gaussian random variable.Differentiating this relation with respect to time yields Comparing the terms in Eqs.(S1) and (S3), one finds: < l a t e x i t s h a 1 _ b a s e 6 4 = " K + i 0 e d f 7 w  To arrive at Eq. ( 4), one needs to evaluate the average κ∇ 2 s|X .The gradients of supersaturation are computed using the ansatz s(x, t) = X[ξ(x/λ(t)), t].This yields: Here, we introduced the components z i of z = x/λ(t), the argument of the Gaussian field ξ.Furthermore, summation over repeated indices is implied, also in the following.Now we take the average conditional on ξ = η.Since ξ is a standardised Gaussian field, one has [4] Averaging Eq. (S4) with (S6) yields Eq. ( 4), if we set Λ = 1.This can be done without loss of generality, because for general Λ, the prefactor κ/λ 2 (t) in the first term on the right-hand side of Eq. ( 4) is replaced by κ/[Λλ(t)] 2 .Since λ(t) is an undetermined function (see next Section), this does not make any difference.We note that Pope [4] does not assume Λ = 1.As a consequence, a length scale λ θ appears in the evolution equations for the mapping X(η, t) in his formulation.Our length scale λ(t) equals the product of λ θ and their time-dependent non-dimensional function J(t), namely λ(t) = J(t)λ θ .
We now discuss how we approximated the conditional average ⟨C d |s = X⟩ in Eq. ( 4).It involves the average ⟨rn|s = X⟩ which is approximated as in terms of the average number N of droplets with supersaturation in the interval [X, X + dX], their sizes, and the volume of regions with that supersaturation (Fig. S1).The ratio ⟨N |s =X⟩/⟨V |s =X⟩ can be expressed as Here n 0 is the droplet-number density, and is the Eulerian probability density of supersaturation [5].The average radius ⟨r|s = X⟩ conditional on supersaturation s = X is given by ⟨r|s =X⟩ = dr r P L (r, s; t) P L (s; t) .(S10) Here P L (s; t) is the Lagrangian supersaturation distribution discussed in the main text, and P L (r, s; t) is the corresponding joint distribution of droplet radii and supersaturation.< l a t e x i t s h a 1 _ b a s e 6 4 = " x j P j 8 < l a t e x i t s h a 1 _ b a s e 6 4 = " n M L e m w Y P r 4 Z 6 z 6 < l a t e x i t s h a 1 _ b a s e 6 4 = " y 5 9 X X t c l P B X P r V i e g 8 K q L F m z W 3 J Z y s d 6 v l L 2 p 6 f 5 Z L z 4 X P l w y S n 9 b H 9 l C F T D y A G 4 6 I 1 9 d s 5 / D + 4 r r r e q V u / q l e a 1 Z m 9 q 2 g P H a A j 5 K E G a q J L 1 small Da Following [6], we use a Monte-Carlo method to simulate the statistical model, using an ensemble of Lagrangian fluid elements.With each such element there is an associated value of ξ(t) from which its supersaturation s(t) = X[ξ(t), t] is computed using the mapping X[η, t].There are two types of fluid elements, one that is used to sample the Eulerian distribution f (X; t), and the other to sample the Lagrangian distribution P L (s, r; t) and the corresponding marginal distribution P L (X; t).Each fluid element of the latter type is associated with a radius r(t), in addition to supersaturations s(t).In practice, f (X; t), P L (X; t) in Eq. (S8), and the average ⟨r|s = X⟩ in Eq. (S10) are evaluated by binning the continuous variable X.Note that the bin [X, X + dX] corresponds represent a rather complex shape in configuration space, as schematically illustrated in Fig. S1.

III. DETERMINATION OF φ(t) AND R(t)
In this Section, we describe how to determine the unknown rates φ(t) and R(t) (see main text).We first discuss the function λ(t).It describes how the length scale of supersaturation fluctuations decays.The mapping closure is a one-point approximation (Section II), so the rate at which the fluctuations relax must be obtained in another way.
Here we extract it from DNS (Section I).In order to make sure that the model predicts the correct scalar dissipation rate ε s (t) = 2κ ⟨∇s • ∇s⟩, we must first evaluate This and the following equations are written in the non-dimensional units introduced in the main text.The average on the r.h.s. is an unconditional average.But since the average conditioned on ξ = η in Eq. (S6) evaluates to unity, independently of η, both averages are the same, conditional and unconditional, and equal to unity.This implies: To compute φ(t) from Eq. (S12) in DNS is straightforward.First, the supersaturation dissipation rate is ε s (t) is readily computed from the gradient of s(x, t).Second, the mapping X(η, t) can be computed from the Eulerian cumulative distribution function of supersaturation, using Eq.(S2).We remark that the DNS results are subject to statistical variability.The evolution of φ(t) for three independent DNS at three different parameter combinations is shown Fig. (S2).Despite the variability, the trend for φ(t) to grow more rapidly at larger Damköhler numbers is evident.The functions φ(t) shown in Fig. 2(a) are averages of the independent functions φ(t) in Fig. S2.Each average in Fig. 1(b) is taken over the corresponding three independent functions in Fig. S2 with identical parameters.Now consider the rate R(t) at which the Lagrangian supersaturation distribution relaxes to the Eulerian one.We now show how to extract the rate R(t) from DNS, and that it obeys for not too large times, which allows us to extract the constant C. The rate R(t) is determined so that the Lagrangian supersaturation distribution in the model relaxes in the same way as in the DNS.We measure the overlap between the two distributions using the Bhattacharyya overlap [7].Computing this measure does not entail dividing with probability densities, whose magnitudes may be small.As a consequence, it is less sensitive to numerical fluctuations than the perhaps more familiar Kullback-Leibler divergence.More precisely, we extract the function from the DNS that minimises the overlap between the Lagrangian distributions in the statistical model and from the DNS.The Lagrangian supersaturation distribution is readily extracted from the DNS.The corresponding distribution in the statistical model is evaluated as follows.We note that it is given by the distribution of ξ(t), because random samples ξ are mapped to the Lagrangian supersaturation distribution by X(η, t).Solving the Langevin equation quoted in the main text yields a solution for the distribution of ξ(t) parameterised by t L (t).Mapping this with X(η, t) yields the model distribution, and we determine t L (t) so that it minimises the overlap between model and DNS distributions.The result is shown in Fig. S3.We plot the Lagrangian time as a function of the Eulerian time for the different parameter combinations used in the main text.We see that the two times are proportional to each other initially, for at least up to one large-eddy time (not shown).This implies that Eq. (S13) is a good approximation, and allows us to extract the constant C (Table S1).The values of C are subject to some statistical variability, as expected.Nevertheless, a tendency for C to decrease with increasing Damköhler number is evident.This reflects the fact that phase change suppresses Lagrangian supersaturation fluctuations.It is worth noting that it is important to model R(t) accurately at early times, but it is not very critical at late times, as explained in the main text.

IV. DETAILS REGARDING FIGS. 2 AND 3
The parameters of the simulations in Fig. 2(a) and (b) are taken from the DNS study of Ref. [3].The simulations in Fig. 2(a) uses parameters from the second row of Table 2 in Ref. [3], and the simulations shown in Fig. 2 (b) use parameters from the third row of that table.We extracted non-dimensional parameters from Ref. [3] as follows.
q M / 0 S j 6 H l 1 E P 6 5 K o 8 5 K 8 w q t R f T z E s N R 6 C U = < / l a t e x i t > '(t) < l a t e x i t s h a 1 _ b a s e 6 4 = " C P G 1 B N X M w y 0 b c g 6 6 X e O I y q L j q n M = " > Da s =0.82 < l a t e x i t s h a 1 _ b a s e 6 4 = " t B e u N E C J A F j o w K M T 7 z 6 J i C m P T w I 9 5 8 w 5 8 4 4 b M S q k Y X z k t I 3 8 5 t b 2 z m 6 h u L d / c F g q H z 2 I M O a Y 9 H H I Q j 5 w k S C M B q Q v q W R k E H G C f J e R R 3 d 6 k + U f n w k X N A x 6 c h Y R x 0 d P A f U o R l J J o 9 J p Y n M f 3 q J 0 N A e R 2 p V r u 2 L p j V G p a u h G 6 8 p s t a C h N 6 3 W 1 a W l w D A a b c u E p o I s q m A Z 3 V E 5 V 7 f H I Y 5 9 E k j M k B B D 0 4 i k k y A u K W Y k L d i x I B H C U / R E h g o D 5 B P h J P M n p L C m l D H 0 Q q 5 O I O F c X e 1 I k C / E z H f r P p I T L j y h W j I U v 4 s y 8 a / c M J a e 5 S Q 0 i G J J A r y Y 6 M U M y h B m x s A x 5 Q R L N l O A M K d q a Y g n i C M s l X 1 r U 3 5 W K N R W B / R M J 8 m W z 2 5 f y y T T K N N F X X 3 U l P B n x M Z p Q f n 7 b S L 8 H x 4 a u t n W m / f N a q e + d H o H n I F z c A F M c A k 6 4 A 5 0 Q R 9 g 8 A J e w R t 4 z 3 1 q e a 2 o 7 S 9 K t d y y 5 x i s h X b y B d l 4 t h k = < / l a t e x i t > Da s =8.2  < l a t e x i t s h a 1 _ b a s e 6 4 = " x j P j 8 < l a t e x i t s h a 1 _ b a s e 6 4 = " n M L e m w Y P r 4 Z 6 z 6 r V i e g 8 K q L F m z W 3 J Z y s d 6 v l L 2 p 6 f 5 Z L z 4 X P l w y S n 9 b H 9 l C F T D y A G 4 6 I 1 9 d s 5 / D + 4 r r r e q V u / q l e a 1 Z m 9 q 2 g P H a < l a t e x i t s h a 1 _ b a s e 6 4 = " X t 7 H z y R 8 2 g q G y H l g M g v U R A / g s e E = " > We extracted the turbulent dissipation rate ε, and used that to force the velocity field as described in Section I.
We then determined the stationary turbulent kinetic energy k in our DNS by averaging the instantaneous kinetic energy over the duration of the simulation.This yields the eddy-turnover time, τ L = k/ε.The values obtained are summarised in Table S1.To compute the Damköhler numbers Da d and Da s , we extracted from Ref. [3] the initial droplet radius r 0 = 20 µm, the droplet-number density n 0 = 164 cm −3 of the initially cloudy air, the density ϱ w = 1000 kg/m 3 of liquid water, the saturation water-vapour mixing ratio q vs = 0.00356, and the mass density ϱ a = 1.06 kg/m 3 of dry air.Comparing Eq. (1b) to Eq. ( 3) of Ref. [3], we inferred A 2 = 1/(ϱ a q vs ) = 265 m 3 /kg.Then we extracted the supersaturation s e = −0.2 of the initially dry air from Fig. 1 of Ref. [3].This allowed us to compute τ s = (4πA 2 A 3 ϱ w n 0 r 0 ) −1 = 1.806 s and 0.1806 s, as well as τ d = r 2 0 /(2A 3 |s e |) = 19.72 s and 1.972 s, for the second and third row of Table 2 of Ref. [3], using the values A 3 = 5.07 × 10 −11 m 2 /sand 5.07 × 10 −10 m 2 /s given there.The resulting Damköhler numbers are quoted in Table S1.Their values differ slightly from the corresponding values quoted in Ref. [3].This reflects that statistically independent DNS with identical values of τ s and τ d , and ε may have slightly different turbulent kinetic energies k, leading to different values of τ L , Da d and Da s .In order to study the effects of droplet phase change independently of the statistical variability of DNS, all DNS reported in Figs. 2 and 3 are shown for identically evolving velocity fields and droplet positions.
The initial mapping was computed using Eq.(S2) and the initial supersaturation profile in Fig. 1 of Ref. [3].From that figure, we extracted the initial volume fraction χ = 0.4 of cloudy air and used this volume fraction to initialise ξ(t)-values for the droplets.To this end, we sampled ξ(t = 0) from a standardised Gaussian distribution conditional on ξ(t) > Φ −1 (1 − χ), where Φ −1 was the inverse cumulative distribution function of a standardised normal distribution.Eq. (S2) implies that this initialisation corresponds to a uniform distribution of droplets in the fraction χ of air with the largest supersaturation [Fig.1(a)].
V. DETAILS REGARDING FIG. 3 The conditional average ⟨n(t)|s(t) = s⟩ in Figure 3 was obtained as follows from the DNS data: ⟨n(t)|s(t) = s⟩ = dn ′ n ′ P(n ′ |s) . (S16) The conditional probability was computed as P(n|s) = P(n, s)/P(s), where P(n, s) is the joint distribution of droplet number density n and supersaturation s, and P(s) is the distribution supersaturation, both obtained from the DNS.Since the average droplet-number density tended to be quite small, we found it necessary to use non-uniform bins in the droplet number n to represent P(n ′ , s), choosing smaller bins for smaller droplet number densities, to resolve the shape of the distribution at small values of n.To this end, we used half-interval Chebyshev nodes.

VI. MODEL PREDICTION FOR DROPLET-SIZE DISTRIBUTION
Fig. S4 shows droplet-size distributions for parameters values extracted from Kumar et al. [3].Shown are DNS results (thin black line), the results of the present model (thick black line), and results of the earlier statistical model of Fries et al. [8] (red line).Both models are in qualitative agreement with the DNS results, but the new model agrees better with the DNS in the tails, compared with the model from Ref. [8].The latter describes mixing of (b) which
g 5 o S 7 y 5 I 2 1 d q 0 c B v 0 7 m v W Q 6 3 e r l e b j c z e I r p E 1 + g W W e g R N d E z a q E u o g j Q C 3 p F b / g d f + I v / L 0 u L e C s 5 w L l A v / 8 A o S C r Q M = < / l a t e x i t > X < l a t e x i t s h a 1 _ b a s e 6 4 = " K 2 4 o q a D 1 U + x w + V B m D o d w W 4 j k A 7 g = " > A A A C Q X i c b V H L S s N A F J 3 4 r P X V 6 t L N a B E E p S R S t M u C G 5 c V + g i 2 o U w m k 3 b o Z B J m J o U S 8 h d u 9 W f 8 C j / B n b h 1 4 6 T N w r R e G D i c + z r 3 j B s x K p V p f h g b m 1 v b O 7 u l v f L + w e H R c a V 6 0 p N h L D D p 4 p C F w n a R J I x y 0 FIG. S1.Schematic.Shows contours of supersaturation s(x, t) = X and s(x, t) = X + dX.To evaluate the averages in Eq. (S7), one records the droplet number N and their sizes r in the hashed region of volume (here area) V .

<
l a t e x i t s h a 1 _ b a s e 6 4 = " g C C E 5 H K P V n S e e 9 O 7 K S m G t d I g W U s = " > A A A C W n i c b Z F b S 8 M w F M e z e t v F y 6 b i i y + Z Y + B D G a 0 M 9 U U Y 6 I O P E 3 a D t Y w 0 T b e w 9 E K S D k b e 9 m 4 Z 1 2 2 g + N W s t M 7 e 3 B M 7 A B a g D C 9 y B F n g E b d A F G H D w A l 7 B m / F u f B p f x v e y d M P I e 0 5 B I Y y f X 3 R A r + I = < / l a t e x i t > (b) FIG. S2.Time dependence of the non-dimensional rate φ(t) from DNS, for Damköhler numbers.Passive scalar (solid black lines), parameters corresponding to Fig. 2(a), green lines, and for the parameters corresponding to Fig. 2(b), violet lines.For all three parameter sets, results of three independent are shown.Also shown is the theoretical estimate φ * = 1 (see main text).
5 6 l s 2 I 0 Z P + Z 5 7 e 7 7 V T f r J M v B d k K 5 A F 6 3 i 6 H y 7 8 z c r F G t q k I 4 J Y u 0 4 T b T L W 2 I c Z w I W c d Z Y 0 I R V Z A J j D y W p w e b t 0 v I C 7 3 q m w K U y / k i H l + x N R U t q a + c 1 7 d X E T Y 0 t r Z c E a G 8 X B f K + 3 L h x 5 c e 8 5 V I 3 D i S 7 6 l g 2 A j u F w y B w w Q 0 w J + Y e E G a 4 N 4 3 Z l B j C n B / X W p d r C / H u z Q Y n a d 4 G 8 + H 1 t U x b 6 c D b n l 9 M B W Z G R L G 4 9 8 H e / 0 / 6 m 9 q 5 < l a t e x i t s h a 1 _ b a s e 6 4 = " g C C E 5 H K P V n S e e 9 O 7 K S m G t d I g W U s = " > A A A C W n i c b Z F b S 8 M w F M e z e t v F y 6 b i i y + Z Y + B D G a 0 M 9 U U Y 6 I O P E 3 a D t Y w 0 T b e w 9 E K S D k b t l / F V v 5 D g h z H d h u z i g c C P / z k n 5 + Q f J 2 J U S M P 4 z m l 7 + w e H R / l C s X R 8 c n p W r p z e 9 m 4 Z 1 2 2 g + N W s t M 7 e 3 B M 7 A B a g D C 9 y B F n g E b d A F G H D w A l 7 B m / F u f B p f x v e y d M P I e 0 5 B I Y y f X 3 R A r + I = < / l a t e x i t > (b) r o P e T c O 6 b T S f m r W W m d t b A m f g A t S B B e 5 A C z y C N u g C D D h 4 A a / g z X g 3 P o 0 v 4 3 t Z u m H k P a e g E M b P L 3 Y W r + M = < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " O h e q l P N k R m T B c c D / z + x z s B x H d 3 8 r o P e T c O 6 b T S f m r W W m d t b A m f g A t S B B e 5 A C z y C N u g C D D h 4 A a / g z X g 3 P o 0 v 4 3 t Z u m H k P a e g E M b P L 3 J q r + E = < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 W n y 4 m dO v j n D O 2 y Y n p f L F B Y a O A o = " > A A A C 2 H i c b V F N a 9 w w E N W 6 X 6 n 7 k a Q 9 9 i K 6 B F J Y F j u E t s c t p d B L I K H Z J H R t F k k e 7 w r L l p H G g c U Y e i u 9 9 t B r + y v 6 X / p v I m 2 2 I Z t k Q O j x Z p 7 m a Y b X S l q M o n + 9 4 N 7 9 B w 8 f b T w O n z x 9 9 n x z a / v F i d W N E T A W W m l z x p k F J S s Y o 0 Q F Z 7 U B V n I F p 7 z 4 6 P O n 5 2 C s 1 N U x L m p I S z a r Z C 4 F Q 0 c l O G 0 T U 9 J P 3 S 6 + m W 7 1 o 2 G 0 D H o b x C v Q J 6 s 4 n G 7 3 / i a Z F k 0 J F Q r F r J 3 E U Y 1 p y w x K o a A L k 8 Z C z U T B Z j B x s G I l 2 L R d m u 7 o j m M y m m v j T o V 0 y V 5 X t K y 0 d l H y Q c l w b m x u n c R D e 7 P I k 3 f l J g 3 m 7 9 N W V n W D U I n L j n m j K G r q R 0 E z a U C g W j j A h J H O N B V z Z p h A N 7 C 1 L l c W w p 3 r D Y 7 j t P X m / e tr m b a o P W 8 H b j U F m H O m s u 7 O B w f / P + l u b h e u q G Q F M L d G d C b C J I M 8 8 S n e e v X I y 2 FIG. S3.Comparison of the two time scales tL(t) and tE(t) extracted from DNS (thick lines).Also shown are linear fits in the region where tL(t) vs. tE(t) roughly proportional lines).Results of three independent simulations are shown.The resulting proportionality constants C (see text) are given in Table S1.(a) Das = 0. (b) Same (a), but at the low Damköhler numbers of Fig. 2(b).(c) Same as (a), but the large Damköhler numbers of Fig. 2(c).

TABLE S1 .
Parameter values for for the DNS runs, see Sections I and III in this supplemental information.The large-eddy time τL, the droplet evaporation time τ d , the supersaturation relaxation time τs, the Damköhler numbers Da d and Das, and the constant C extracted from four indepdent DNS of a passive scalar, four independent DNS at low Dast numbers, and three indepdent DNS at high Das numbers.The DNS in Fig.2are for identically evolving velocity fields and droplet positions (see Section IV).