# Diffuselet method for three-dimensional turbulent mixing of a cloudy air filament

Vladyslav Pushenko, Simone Scollo, Patrice Meunier, Emmanuel Villermaux,

Jörg Schumacher

► To cite this version:

Vladyslav Pushenko, Simone Scollo, Patrice Meunier, Emmanuel Villermaux, Jörg Schumacher. Diffuselet method for three-dimensional turbulent mixing of a cloudy air filament. Physical Review Fluids, 2025, 10, 10.1103/1bt3-mlpm. hal-05479000

HAL Id: hal-05479000

https://hal.science/hal-05479000v1

Submitted on 27 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Diffuselet method for three-dimensional turbulent mixing
of a cloudy air filament

Vladyslav Pushenko,¹ Simone Scollo,² Patrice Meunier,²

Emmanuel Villermaux,²,³ and Jörg Schumacher¹,⁴

¹Institute of Thermodynamics and Fluid Mechanics, P.O.Box 100565, Technische Universität Ilmenau,
D-98684 Ilmenau, Germany

²Institut de Recherche sur les Phénomènes Hors Equilibre (IRPHE), CNRS, Aix-Marseille Université,
13384 Marseille Cedex 13, France

³Institut Universitaire de France, 75005 Paris, France

⁴Tandon School of Engineering, New York University, New York, New York 11201, USA

<figure><img src="image_2.png" /></figure>

(Received 5 March 2025; accepted 18 June 2025; published 16 July 2025)

The mixing properties of vapor content, temperature, and particle fields are of paramount importance in cloud turbulence as they pertain to essential processes, such as cloud water droplet evaporation and entrainment. Our study examines the mixing of a single cloudy air (which implies droplet-laden) filament with its clear air environment, a characteristic process at the cloud edge, in two ways. The first consists of three-dimensional combined Euler-Lagrangian direct numerical simulations which describe the scalar supersaturation as an Eulerian field and the individual cloud water droplets as an ensemble of Lagrangian tracers. The second way builds on the recently developed diffuselet method, a kinematic Lagrangian framework that decomposes a scalar filament into a collection of small sections subject to deformation by local stirring and cross-sheet diffusion. The Schmidt number is $Sc = 0.7$. The entrainment process causes a deformation of the supersaturated cloud filament in combination with diffusion until the system reaches a well-mixed equilibrium state, which implies for the present configuration that all droplets are evaporated. We compare the time dependence of the mean square and probability density function of the supersaturation field. For the initial period of the mixing process, they agree very well; at later stages, deviations caused by nonzero mean of the conserved scalar are observed. For the cases including cloud water droplets, we also investigate the impact of droplet number density and condensation growth response. Turbulence causes deviations from the $d^2$ law similar to recent experiments in sprays. A simulation at a Schmidt number that is by a factor of 100 larger than in clouds improves the agreement between simulation and diffuselet method significantly. The latter result promotes the diffuselet framework as an efficient parametrization for turbulent high-Sc mixing which can reduce the resolution efforts of the viscous-convective range of scalar turbulence.

DOI: 10.1103/1bt3-mlpm

# I. INTRODUCTION

Turbulent mixing and the coupling to droplet microphysics in atmospheric clouds remains one major source of uncertainty for more reliable predictions of the ongoing climate change and the

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
connected global warming [1–4]. It is expected that a better understanding of these small-scale processes will improve their convective parametrization starting from large-eddy simulations of individual clouds all the way up to global atmospheric circulation models. Water exists in warm clouds in the form of water vapor and cloud water droplets [5]. The entrainment of dry clear into moist cloudy air has a significant impact on the statistical properties inside the clouds as well as their lifetime [6]. This also determines number density and size distribution of the droplets [7–14] and thus eventually the appearance of precipitation [15].

Turbulent mixing [16–18] in general combines two fundamental processes, the stirring of a scalar substance due to advection by the (turbulent) flow and molecular diffusion [19]. A turbulent flow at its smallest scales is characterized by spatially smooth velocity fields; they stretch and fold scalar filaments in a range of decreasing scales that starts around the viscous Kolmogorov length $\eta_K$ [20,21]. In this way, fluid motion leads to an ongoing aggregation of scalar filaments if the scalar diffusion is sufficiently small (which is quantified in Sec. II by the dimensionless Schmidt number $Sc$). A specific form of the probability density function of the scalar concentration field in the form of a gamma distribution follows from this picture [22–24]. At the smallest spatial scales in the mixing problem, molecular diffusion of the scalar will always dominate the dynamics. This is the beginning of the inertial-diffusive and viscous-diffusive scale range for Schmidt (or Prandtl) numbers $Sc \ll 1$ and $Sc \gg 1$, respectively.

The primary goal of the present work is to investigate whether a kinematic Lagrangian framework that incorporates these elementary building blocks of the turbulent mixing process—the diffuselet model [25]—is applicable to the entrainment at the edge of a turbulent cloud. Our model system is a thin initially slab-like filament of cloudy supersaturated air in a clear air environment at the same turbulent kinetic energy level, similar to previous numerical studies [10,26–29]. To this end, we will proceed in two subsequent steps. First, we consider the passive scalar mixing model of the supersaturation field [11,30,31] and investigate how the diffuselet model compares to fully resolved direct numerical simulations (DNS) of three-dimensional scalar turbulence in a periodic box, see, e.g., Refs. [31,32]. Two cases are studied here, one at a Schmidt number of $Sc = 0.7$, which corresponds to atmospheric conditions, and one at a hundred times larger value of $Sc = 70$, which develops a viscous-convective range for scales $r \lesssim \eta_K$. The latter case is added even though this high Schmidt number cannot be related to cloud turbulence. The early period of the mixing process is well described by the diffuselet model; at later times we will observe growing deviations and discuss the reasons.

Secondly, we add an initially monodisperse ensemble of cloud water droplets to the scalar filament and numerically solve the cloud model [31]. We investigate how strongly the deviations from the diffuselet model grow with time in comparison to the pure scalar field case. Therefore, we enhance the droplet number density beyond realistic values of $n_d \approx 10^2 \text{ cm}^{-3}$ [1,5]. Furthermore, we increase the cloud water droplet response to super- or subsaturation by an increase of the diffusion constant $K$ in the droplet growth equation [26]. Our results are also compared with dense spray mixing experiments [33]. In summary, we find that the additional impact of the droplets on the presently chosen configuration remains small, such that the diffuselet method in its present form is even applicable (in the same limits as for the passive scalar case) to a two-phase mixing process. A major scope of this work is to explore the potential but also the limitations of the diffuselet model in the absence of a viscous-convective scale range of the scalar field for the transport in a three-dimensional complex Navier-Stokes flow.

The outline of the manuscript is as follows: In Sec. II, we describe the equations, parameters, and physically relevant scales and present a brief summary of the central ideas of the diffuselet model. Section III presents all results, starting from the pure scalar case to the filament mixing with droplets. The paper ends with our conclusion and an outlook in Sec. IV. Technical details on the simulations are summarized in the Appendix.
## II. METHODS

### A. Eulerian fields and direct numerical simulations of turbulent mixing

The present cloud mixing model is based on three fields: the vapor mixing ratio $q_v(x, t)$, the liquid water mixing ratio $q_l(x, t)$, and the temperature field $T(x, t)$, see Ref. [31] for further details. The fallout of rain is not included, such that $q_v + q_l = \text{const}$. The modeling of the liquid water mixing ratio $q_l$ in the Lagrangian framework will be detailed in Sec. II B. The Eulerian fields $q_v$ and $T$ are combined in the scalar supersaturation field $s(x, t)$, which is approximated as a passive scalar that determines the diffusion growth of the cloud water droplets. It is given by

$$s(x, t) = \frac{q_v(x, t)}{q_{vs}[T(x, t)]} - 1, \qquad (1)$$

where $q_{vs}$ is the saturated vapor mixing ratio that depends on temperature $T$, determined by the Clausius-Clapeyron equation [5]. For convenience, we rescale supersaturation $s$ to the scalar field $c$ in the following way:

$$0 \le c(x, t) = \frac{s(x, t) + 1}{s_{\max} + 1} \le 1. \qquad (2)$$

The transport of $c$ is determined by

$$\frac{\partial c}{\partial t} + (\mathbf{u} \cdot \nabla)c = D\nabla^2 c, \qquad (3)$$

where $D$ is the scalar diffusivity and $\mathbf{u}$ the advecting (turbulent) velocity field. The dynamics of the latter field is given by the incompressible Navier-Stokes equations,

$$\nabla \cdot \mathbf{u} = 0 \quad \text{and} \quad \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho_0} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}. \qquad (4)$$

Here, $p(x, t)$ is the pressure field, $\rho_0$ the constant mass density, and $\nu$ the kinematic viscosity. Furthermore, $f(x, t)$ is a large-scale volume forcing which injects turbulent kinetic energy into the fluid flow at a fixed rate, see further below in this section for details. The dimensionless form of (3) and (4) is obtained by defining characteristic values of all physical quantities. This is the length of the cubic simulation domain $L$, the reference scalar amplitude $c_0$, and the root-mean-square value of the velocity field $U$, which leads to

$$\nabla \cdot \mathbf{u} = 0, \qquad (5)$$

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \frac{1}{\text{Re}}\nabla^2 \mathbf{u} + \mathbf{f}, \qquad (6)$$

$$\frac{\partial c}{\partial t} + (\mathbf{u} \cdot \nabla)c = \frac{1}{\text{Re Sc}}\nabla^2 c. \qquad (7)$$

Here, $\text{Re} = UL/\nu$ is the large-scale Reynolds number and $\text{Sc} = \nu/D$ is the Schmidt number. The Schmidt number $\text{Sc}$ relates consequently momentum to scalar diffusion.

In three-dimensional turbulent flows, the turbulent kinetic energy, which is transferred from large to small scales, is converted into heat starting at the Kolmogorov scale $\eta_K = (\nu^3/(\epsilon))^{1/4}$ with the mean rate of kinetic energy dissipation, $\langle \epsilon \rangle$ [34]. This cascade picture can be applied to the scalar variance cascade depending on $\text{Sc}$ [17]. For $\text{Sc} \le 1$, the scalar variance cascades down to smaller scales in the inertial-convective range, which is bounded from below by the Corrsin length $\eta_C = (D^3/(\epsilon))^{1/4} \gtrsim \eta_K$, as suggested by Obukhov [35] and Corrsin [36] independently. For $\text{Sc} \ge 1$, the scalar variance cascades down the inertial range to $\eta_K$. For even smaller scales $r$, the viscous-convective range is established up to the Batchelor scale $\eta_B = \eta_K/\sqrt{\text{Sc}} = \sqrt{D/\gamma}$ with the mean stretching rate $\gamma$ [20,21]. In this range between $\eta_K \gtrsim r \gtrsim \eta_B$, the scalar is stirred by a spatially smooth flow, which can be considered as a chaotic mixing regime [23,37]. It becomes clear now that
<table><caption>TABLE I. Parameters of the direct numerical simulations. Here, $N^3$ is the number of grid points, $\Delta x/\eta_K$ is the ratio of the uniform grid spacing to Kolmogorov length. The energy injection rate is $\epsilon_{inj} \approx 0.004 \text{ m}^2\text{s}^{-3}$, which results in Kolmogorov microscale $\eta_K = 1$ mm. The root mean square velocity is $U = 0.097 \text{ m/s}$. Thus, the large-scale eddy turnover time follows to $\tau_L = L/U = 2.64 \text{ s}$ for both runs. The number of diffuselets is $N_{diff} = 256^2$.</caption><thead><tr><th>Run</th><th>L [m]</th><th>N<sup>3</sup></th><th>Δx/η<sub>K</sub></th><th>δA<sub>0</sub> [mm<sup>2</sup>]</th><th>ℓ<sub>0</sub> [m]</th><th>Re</th><th>Re<sub>λ</sub></th><th>Sc</th><th>τ<sub>η</sub> [s]</th><th>D [m<sup>2</sup>s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>T1</td><td>0.256</td><td>512<sup>3</sup></td><td>0.5</td><td>1(=η<sub>K</sub><sup>2</sup>)</td><td>0.005</td><td>1634</td><td>52</td><td>0.7</td><td>0.07</td><td>2.52 × 10<sup>-5</sup></td></tr><tr><td>T2</td><td>0.256</td><td>2048<sup>3</sup></td><td>0.125</td><td>1(=η<sub>K</sub><sup>2</sup>)</td><td>0.005</td><td>1634</td><td>52</td><td>70</td><td>0.07</td><td>2.52 × 10<sup>-6</sup></td></tr></tbody></table>

the DNS at $Sc = 0.7$ is assigned to Corrsin-Obukhov regime of passive scalar turbulence, the one at $Sc = 70$ to the Batchelor-Kraichnan regime.

The DNS applies the pseudospectral method. Turbulent velocity and passive scalar field $u(x, t)$ and $c(x, t)$ are expanded in Fourier modes with respect to each of the three space directions relying on volumetric fast Fourier transformations of the software package P3DFFT [38]. Time stepping is done with a second-order predictor-corrector scheme [39], see also Ref. [31]. Further details on the required spectral resolution and the resulting Reynolds numbers will be listed in Tables I and II. The computational mesh is uniform in all three space directions with a mesh size $\Delta x$.

The volume forcing $f$ is applied in Fourier space for a selected set of low-wavenumber Fourier modes, following the scheme introduced in Ref. [32]. More specifically, the forcing acts on a small subset $\mathcal{K}$ of wavevectors:

$$ \mathcal{K} = \left\{ k_j^f = \frac{2\pi}{L}(\pm 1, \pm 1, \pm 2) \text{ plus permutations} \right\}, \quad (8) $$

ensuring an isotropic kinetic energy injection at the largest scales. The forcing is given by

$$ f_i(k_j, t) = \epsilon_{\text{inj}} \frac{\hat{u}_i(k_j, t)}{\sum_{k_j^f \in \mathcal{K}} |\hat{u}_i(k_j^f, t)|^2} \delta_{k_j k_j^f}. \quad (9) $$

In the statistically steady state, this balances the mean energy dissipation rate, such that $\langle \epsilon \rangle = \epsilon_{inj}$, resulting in a turbulent velocity field with a Kolmogorov length scale which can be prescribed. The kinetic energy injection rate, $\epsilon_{inj} = \langle u \cdot f \rangle$, is chosen such that $\eta_K = 1$ mm, which is typical of cloud turbulence. So far, we have discussed the treatment of Eulerian fields (velocity, passive scalar). The Lagrangian part follows now in Sec. II B.

<table><caption>TABLE II. Parameters of the simulations for turbulent scalar transport with phase changes. Here $n_d$ refers to the number density of liquid droplets, $K$ is the evaporation constant, $r_0$ is the initial droplet radius. Here, the initial thickness of cloud profile is always $\ell_0 = 0.01 \text{ m}$. The number of diffuselets is $N_{diff} = 256^2$ for runs D1, D2, D3, D7, and D8 and $N_{diff} = 512^2$ for runs D4, D5, and D6.</caption><thead><tr><th>Run</th><th>L [m]</th><th>N<sup>3</sup></th><th>n<sub>d</sub> [cm<sup>-3</sup>]</th><th>K [m<sup>2</sup>/s]</th><th>r<sub>0</sub> [µm]</th></tr></thead><tbody><tr><td>D1</td><td>0.256</td><td>512<sup>3</sup></td><td>100</td><td>1.2 × 10<sup>-11</sup></td><td>20</td></tr><tr><td>D2</td><td>0.256</td><td>512<sup>3</sup></td><td>500</td><td>1.2 × 10<sup>-11</sup></td><td>20</td></tr><tr><td>D3</td><td>0.256</td><td>512<sup>3</sup></td><td>1000</td><td>1.2 × 10<sup>-11</sup></td><td>20</td></tr><tr><td>D4</td><td>0.512</td><td>1024<sup>3</sup></td><td>100</td><td>1.2 × 10<sup>-11</sup></td><td>20</td></tr><tr><td>D5</td><td>0.512</td><td>1024<sup>3</sup></td><td>500</td><td>1.2 × 10<sup>-11</sup></td><td>20</td></tr><tr><td>D6</td><td>0.512</td><td>1024<sup>3</sup></td><td>1000</td><td>1.2 × 10<sup>-11</sup></td><td>20</td></tr><tr><td>D7</td><td>0.256</td><td>512<sup>3</sup></td><td>100</td><td>1.2 × 10<sup>-10</sup></td><td>20</td></tr><tr><td>D8</td><td>0.256</td><td>512<sup>3</sup></td><td>100</td><td>1.2 × 10<sup>-11</sup></td><td>10</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>FIG. 1. The sheet configuration of the initial concentration field $c(x, 0)$. (a) The sheet is placed in the triply periodic simulation domain $V = L^3$ parallel to the y-z plane and centered at $x = L/2$. (b) The cross-sheet concentration profile is a Gaussian, see Eq. (10), characterized by a thickness $\ell_0$ in the x direction. The thicknesses in panels (a) and (b) are magnified for better visualization. (c) Decomposition of the sheet into diffuselets. Center points $X_k(t = 0)$ are indicated by black dots as well as the normal vector $n_0$. Index $k$ ranges from 1 to $N_{\text{diff}}$.</figcaption></figure>

## B. Lagrangian diffuselet framework of turbulent mixing

In accordance with Ref. [25], a scalar sheet is divided into $N_{\text{diff}}$ finite-size elements, as illustrated in Fig. 1. The sheet is very thin, such that the well-mixed equilibrium will correspond to a mean vapor content below the saturation level. This implies that all droplets are evaporated. For simplicity, we select a Gaussian transversal profile for the initial scalar field,

$$c(x, t = 0) = \exp\left(-\frac{(x - L/2)^2}{\ell_0^2}\right), \qquad (10)$$

with a sheet thickness $\ell_0$. This profile translates to the same Gaussian shape for the supersaturation field $s$. The latter field is positive, i.e., supersaturated in a cross section of approximately $\ell_0/2$ with $s_{\text{max}} = 0.05$ in the sheet center, see Eq. (2). Each finite-size element is associated with a Lagrangian tracer, whose coordinates correspond to the coordinates of the geometrical center of the element, see Fig. 1(c). Lagrangian tracers are seeded in the center plane of the volume with initial positions $\{X_i^0\}_{i=1,...,N_{\text{diff}}}$ Their trajectories $\{X_i(t)\}_{i=1,...,N_{\text{diff}}}$ are given by

$$\frac{dX_i(t)}{dt} = u(X_i, t), \quad i = 1, \dots, N_{\text{diff}} \quad \text{with} \quad X_i(t = 0) = X_i^0. \qquad (11)$$

These elements are referred as *diffuselets*, and the transport equation (3) can be solved for each diffuselet individually. The resulting data will then be employed in the reconstruction of the global statistical properties of the scalar field. Time integration of the particle (or diffuselet) trajectories is performed with the same second-order predictor-corrector method as the Eulerian fields. The interpolation of the velocity field to each diffuselet position is trilinear.

The concentration field of each diffuselet can be calculated using the Ranz coordinate transformation [40]. The maximal concentration $c_i$ of the scalar field inside the diffuselet $i$ with the surface $\delta A_i$ and thickness $\ell_i$ is simply deduced from a new dimensionless time $\tau_i$ by $c_i = 1/\sqrt{\tau_i}$. This dimensionless time is determined numerically by integrating $(D$ and $\ell_i$ carry their physical dimension)

$$\frac{d\tau_i}{dt} = \frac{4D}{\ell_i^2}, \quad \text{with} \quad \tau_i(t = 0) = 1. \qquad (12)$$

The advection-diffusion equation (3) is simplified to pure diffusion equation in the new coordinates $(\xi, \tau)$, which is given by

$$\frac{\partial c_i}{\partial \tau} = \frac{1}{4} \frac{\partial^2 c_i}{\partial \xi^2}. \qquad (13)$$
Here, $\xi = n_i/\ell_i$ is a dimensionless coordinate normal to the diffuselet surface. Each diffuselet is characterized by two key parameters: (1) the initial orientation, designated as $\mathbf{n}_0$, and (2) the initial surface vector, $\delta\mathbf{A}_0 = \mathbf{n}_0\delta A_0$. All elements have an identical initial surface $\delta A_0 \sim \eta_K^2$. Under the influence of the local velocity gradients, the surface elements are subject to stretching and compression. A local stretching factor is defined as $\rho_i = \delta A_i/\delta A_0$. Incompressibility of the velocity field requires $\delta A_i(t)\ell_i(t) = \delta A_0\ell_0$ for the stretching history with $\delta A_i(t)$ and $\ell_i(t)$ being the surface and thickness of $i$th diffuselet at time $t$ initiating from $\delta A_0$ and $\ell_0$. More specifically, the kinematic stretching dynamics is determined by the velocity gradient tensor along the Lagrangian tracer track $\mathbf{X}_i(t)$ that was initiated in the center of diffuselet $i$,

$$ \frac{d\delta\mathbf{A}_i}{dt} = -(\nabla u)^T [\mathbf{X}_i(t), t] \delta\mathbf{A}_i. \qquad (14) $$

It is actually more powerful to define a surface dispersion tensor $\mathcal{L}_i$ for each tracer, which is a solution of

$$ \frac{d\mathcal{L}_i(t)}{dt} = -(\nabla u)^T [\mathbf{X}_i(t), t] \mathcal{L}_i(t) \quad \text{with} \quad \mathcal{L}_i(0) = \mathbf{1}. \qquad (15) $$

Multiplying this equation on the right by $\delta\mathbf{A}_0$ indicates that $\mathcal{L}_i(t)\delta\mathbf{A}_0$ is a solution of (14). The surface vector is thus given by

$$ \delta\mathbf{A}_i(t) = \mathcal{L}_i(t)\delta\mathbf{A}_0, \qquad (16) $$

whatever the initial surface vector $\delta\mathbf{A}_0$. The stretching factor of the surface element $\rho_i(t) = (\delta A_i^T \delta A_i)^{1/2}/\delta A_0$ simply writes

$$ \rho_i(t) = [\mathbf{n}_{0,i}^T \mathcal{L}_i(t)^T \mathcal{L}_i(t) \mathbf{n}_{0,i}]^{1/2}. \qquad (17) $$

It depends only on the initial orientation $\mathbf{n}_{0,i}$ of the diffuselet. It is possible to average $\rho_i$ over all initial orientations in order to accelerate the convergence of the statistics [25]. However, in this paper, we chose to fix $\mathbf{n}_{0,i}$ along the $x$ direction to be closer to the DNS results. The difference between the two results was found to be very small.

Using the incompressibility, the evolution of the diffuselet thickness $\ell_i(t)$ can be written as

$$ \ell_i(t) = \frac{\ell_0 \delta A_0}{\delta A_i(t)} = \frac{\ell_0}{[\mathbf{n}_{0,i}^T \mathcal{L}_i(t)^T \mathcal{L}_i(t) \mathbf{n}_{0,i}]^{1/2}}. \qquad (18) $$

It is then possible to generalize the dimensionless time to a dimensionless tensor, see again Eq. (12),

$$ \tau_i(t) = 1 + 4D \int_0^t \frac{1}{\ell_i^2(t')} dt' \rightarrow \mathcal{T}_i = \mathbf{1} + \frac{4D}{\ell_0^2} \int_0^t (\mathcal{L}_i^T \mathcal{L}_i)[\mathbf{X}_i(t'), t'] dt', \qquad (19) $$

from which the dimensionless time is obtained as

$$ \tau_i(t) = \mathbf{n}_{0,i}^T \mathcal{T}_i(t) \mathbf{n}_{0,i} \qquad (20) $$

The maximal concentration of each diffuselet is then equal to $c_i = 1/\sqrt{\tau_i}$. The mean square of each diffuselet is, to a good approximation, the product of the squared maximal concentration $c_i^2$ and the volume of the diffuselet at time $t$, given by $\delta A_0\ell_0/\ell_i \times \ell_i\sqrt{\tau_i} = \delta A_0\ell_0\sqrt{\tau_i}$ and multiplied by a factor of $\sqrt{\pi}/2$ due to the Gaussian profile. This leads to

$$ \langle c_i^2(t) \rangle \approx \frac{\sqrt{\pi} \delta A_0 \ell_0}{\sqrt{2 \tau_i(t)}}. \qquad (21) $$
Note that $\ell_i\sqrt{\tau_i(t)}$ is the diffusive thickness used in Eq. (21). The total mean square of all diffuselets is consequently given by sum of (21) over the entire diffuselet ensemble,

$$ \langle c^2(t) \rangle = \sum_{i=1}^{N_{\text{diff}}} \langle c_i^2(t) \rangle \approx \sum_{i=1}^{N_{\text{diff}}} \frac{\sqrt{\pi} \delta A_0 \ell_0}{\sqrt{2 \tau_i(t)}}. \quad (22) $$

Following a similar approach as for the mean square, each diffuselet contributes to the total concentration probability density function (PDF) at the level $c = c_i$ with a weight given by the diffuselet volume $\delta A_0 \ell_0/c_i$. The PDF of maximal concentration is then approximately obtained by $\sqrt{\tau_i/\pi}$ for a top-hat profile. A subsequent convolution with the Gaussian profile (which is used here) yields an additional factor of $1/(c\sqrt{-\ln(c\sqrt{\tau_i})})$. Thus, the PDF amplitude at concentration $c$, denoted as $P(c)$, is calculated as the sum over all diffuselets.

$$ P(c) \approx \sum_{i=1}^{N_{\text{diff}}} \frac{\sqrt{\tau_i/\pi}}{c\sqrt{-\ln(c\sqrt{\tau_i})}}. \quad (23) $$

Equations (22) and (23) will be numerically evaluated with respect to time and compared with the DNS results. More details on the mathematical formulation of the mathematical framework (including the diagonalization of $\mathcal{L}_i^T \mathcal{L}_i$) can be found in Ref. [25]. In the Appendix B, we also report resolution studies, where differently fine decompositions of the filament into diffuselets have been compared to each other. This was done to confirm that mean-squared scalar and scalar PDF remain unaffected for different resolutions. Also in Appendix A, we have validated the diffuselet model with an analytical flow example.

### III. RESULTS

#### A. Passive scalar mixing in three-dimensional homogeneous isotropic flow

##### 1. Stretching ratio statistics in turbulent flow

The turbulent mixing in a homogeneous isotropic three-dimensional velocity field adds significant complexity to the problem at hand, caused by the continuum of turbulent eddy scales, which are subject to vortex stretching and other vortex-vortex interactions, see, e.g., Refs. [41–43]. The dynamics of this turbulent flow requires DNS that solve numerically the incompressible Navier-Stokes equations together with the passive scalar advection-diffusion equation and the Lagrangian diffuselets.

The statistics of the stretching factors $\rho_i$, see Eq. (17), is shown in Fig. 2. We plot the PDF of $\tilde{\rho}_i = (\ln \rho_i - \langle \ln \rho_i \rangle)/\sigma (\ln \rho_i)$ in Fig. 2(a). The evolution is similar to that in the seminal work by Girimaji and Pope [42]. Their numerical analysis was conducted under similar but not exactly the same conditions in a homogeneous isotropic flow. The PDFs are very close to Gaussian, except for a small positive skewness that decreases with time.

In Figs. 2(b)–2(d) of the same figure, we list the first three moments of $\ln \rho_i$ versus time. The first moment corresponds to the mean Lyapunov exponent

$$ \lambda = \frac{\langle \ln \rho_i \rangle}{t}. \quad (24) $$

It is initially large and then decreases to a value close to $\lambda = 0.159/\tau_\eta$, where $\tau_\eta = 0.07$ s is the Kolmogorov timescale. This is in fair agreement with Girimaji and Pope [42] who found a value $\lambda = 0.165/\tau_\eta$ and who explained the initial decay of the stretching rate by the alignment of the surfaces with the shear. The second central moment is

$$ V = \langle (\ln \rho_i - \lambda t)^2 \rangle. \quad (25) $$

Since $\ln \rho_i$ follows a random process, the mean square increases linearly in time. The rescaled mean square $V/t$ reaches a saturation value equal to $V/t = 0.18/\tau_\eta$. This is larger than the value
<figure><img src="image_4.png" /><figcaption>FIG. 2. Statistics of the surface stretching ratios $\rho_i = \delta A_i(t)/\delta A_0$. (a) PDF of the normalized stretching $\tilde{\rho}_i = (\ln \rho_i - \langle \ln \rho_i \rangle)/\sigma(\ln \rho_i)$ for different moments of time. The black line represents the normal distribution $N(0, 1)$. (b) Mean finite-time Lyapunov exponent $\lambda = \langle \ln \rho_i \rangle/t$ versus time. (c) Temporal evolution of the time-rescaled mean square $V/t$ of $\ln \rho_i$, as defined by Eq. (25). (d) Temporal evolution of the skewness of the PDF of $\ln \rho_i$, as defined by Eq. (26).</figcaption></figure>

$V/t = 0.1/\tau_\eta$ found in Girimaji and Pope [42]. This difference can be caused by the slightly different Reynolds numbers, differences in the spectral resolutions, and the volume forcing methods $f$ in the Navier-Stokes equations (4). Note that we inject kinetic energy at a constant rate into the flow, which determines $\langle \epsilon \rangle$ [32]. This differs from the stochastic forcing used in Ref. [42]. The third central moment of $\ln \rho_i$ rescaled by $V^{3/2}$ corresponds to the skewness of the PDF of $\ln \rho_i$:

$$S = \frac{\langle (\ln \rho_i - \lambda t)^3 \rangle}{V^{3/2}}. \qquad (26)$$

In Fig. 2(d), the skewness of $\ln \rho_i$ is initially weakly positive and decreases with time. It is curious to see that the PDFs of the surface stretching are skewed, whereas the PDFs of line stretching (for pair dispersion) are not skewed. This effect is visible in the results of Girimaji and Pope [42]. It should be noted that the PDFs of surface stretching are not skewed for a sine flow [25].

The Taylor-microscale Reynolds number $Re_\lambda$ of the turbulent flow is given by

$$Re_\lambda = \sqrt{\frac{5}{3\nu\langle\epsilon\rangle}} U^2, \qquad (27)$$

and has a moderate value of 52, as provided in Table I. The spectral resolution is $k_{max}\eta_K \approx 6$ for $Sc = 0.7$ and $k_{max}\eta_B \approx 3$ for $Sc = 70$ with the maximally resolved wavenumber $k_{max} = 2\pi\sqrt{2N/(3L)}$. This guarantees a sufficient resolution of the smallest scales of the turbulent flow. The DNSs are performed for two Schmidt numbers of 0.7 and 70.

## 2. Schmidt number $Sc = 0.7$

The comparison of DNS run T1 with the diffuselet model is presented in Fig. 3. The central part of the figure displays the mean squared scalar $\langle c^2(t) \rangle$, which corresponds to the supersaturation field in the cloud slab. The top row shows contour plots of the concentration field at selected times, illustrating the dynamics of the mixing process. The bottom row compares the PDFs of the passive
<figure><img src="image_5.png" /><figcaption>FIG. 3. Comparison of the diffuselet method with DNS for a turbulent flow at $Sc = 0.7$ (Run T1). The central panel presents the evolution of mean-squared scalar $\langle c^2 \rangle$ over time, with selected time points I–IV marked by red dots. The top row shows contour plots of the scalar field (logarithm of concentration) at these times, illustrating the effect of high diffusivity on scalar structures. The bottom row compares PDF $P(c)$ of scalar concentration $c$ for both methods, showing their agreement across different stages of scalar mixing. The green dashed-dotted line represents a gamma distribution, see Eq. (29), around the mean value.</figcaption></figure>

scalar, obtained by both methods, at the same time instants. It is observed how the filament is increasingly deformed and convoluted while simultaneously being subject to molecular diffusion. Both the DNS and diffuselet methods initially show a very good agreement, as demonstrated by the mean square and PDF. Deviations arise for $t \gtrsim 1$ s, which corresponds to times larger than $\approx 0.3\tau_L$ and thus to instant III in the figure. The DNS data converge to the perfectly mixed state $\langle c^2 \rangle \rightarrow \langle c \rangle^2$ for times $t \gg 1$ s. This is a consequence of the conserved scalar dynamics, which implies that the volume mean $\langle c \rangle = \text{const.}$ with respect to time $t$ in the present setup with periodic boundary conditions in all three spatial directions. Such an offset is absent in the diffuselet model (which is formulated without boundary conditions), such that the mean square continues here to decay exponentially with respect to time. The mixing time definition (A2) is adapted to the turbulent case where surfaces grow exponentially in time:

$$t_s = \frac{1}{2\gamma_{\text{turb}}} \ln \left( \frac{\gamma_{\text{turb}} \ell_0^2}{D} \right) \quad \text{with} \quad \gamma_{\text{turb}} = \sqrt{\frac{\epsilon_{\text{in}}}{2\nu}} = \sqrt{\frac{\langle \epsilon \rangle}{2\nu}}. \qquad (28)$$

Here, the stretching rate $\gamma_{\text{turb}}$ is obtained for a statistically steady regime of the fluid turbulence with $\epsilon_{\text{inj}} = \langle \epsilon \rangle$. For run T1, this results in $\gamma_{\text{turb}} = 10.6 \text{ s}^{-1}$ and thus a short mixing time of $t_s \approx 0.11 \text{ s}$ since $D = 2.5 \times 10^{-5} \text{ m/s}^2$. The mixing time corresponds approximately to time instant II in Fig. 3. It explains why at time instant II the mean square has started to decay. It also explains why the peak at $c = 1$ has disappeared in the concentration PDF. There is a very good agreement between the diffuselet and the DNS at this time instant.
However, there is a discrepancy at later times for the PDFs between the DNS and the diffuselet method due to periodic boundary conditions and the aggregation process. The scalar concentration PDF at time instant IV can be fitted fairly well to a gamma distribution (except for the far tails), which suggests that a scalar sheet aggregation process is at work, despite the absence of a viscous-convective cascade range for Sc < 1. The gamma distribution has the following form:

$$P(c) = \frac{c^{\alpha-1}\theta^{\alpha}\exp(-\theta c)}{\Gamma(\alpha)}, \qquad (29)$$

where $\alpha = \langle c \rangle^2/\langle c^2 \rangle$ and $\theta = \langle c \rangle/\langle c^2 \rangle = \alpha/\langle c \rangle$. By application of the substitution, $P(c)dc = P(\tilde{c})d\tilde{c}$, we can convert the PDF in (29) into

$$P(\tilde{c} = \frac{c}{\langle c \rangle}) = \frac{\tilde{c}^{\alpha-1}\alpha^{\alpha}\exp(-\alpha\tilde{c})}{\Gamma(\alpha)}, \qquad (30)$$

with only one time-dependent fit parameter $\alpha$ left, which is interpreted as the average number of sheet aggregations [22]. This is at time instant IV, roughly just $\alpha \gtrsim 1$ for Sc = 0.7.

Finally, we note that the convergence into the well-mixed final state can be postponed by a reduction of $\langle c \rangle$, i.e., by a further reduction of the initial filament sheet thickness $\ell_0$.

### 3. Schmidt number Sc = 70

To proceed with a comparison of the DNS and diffuselet method, we now consider a Schmidt number of Sc = 70, which is two orders of magnitude larger than in Sec. III A 2. It is clear that this case is not relevant for cloud turbulence. However, it allows us to further explore the capabilities of the diffuselet framework. The high Schmidt number causes a Batchelor length $\eta_B$ of the passive scalar, which is by approximately an order of magnitude smaller than $\eta_K$. As shown in Fig. 4, structures of the scalar field are less blurred and show much finer striations in comparison to those of Fig. 3 at Sc = 0.7. Consequently, the PDF and mean square of the scalar field require a longer time interval to converge to the well-mixed equilibrium state. As is also evident from Fig. 4, both quantities agree for a longer time period. However, eventually this DNS case also starts to deviate from the predictions of the diffuselet model.

The DNSs for the present case were expensive since the full viscous-convective range had to be resolved for the high-Schmidt-number case. While leaving the outer box dimension L at the same value as for Sc = 0.7, the uniform mesh was by a factor of $4^3 = 64$ larger, as seen in Table I. Therefore, with a further increase of the Schmidt number to values Sc $\approx 10^3$, the diffuselet method can serve as a potentially efficient way to model the turbulent mixing without the necessity to resolve all features of the scalar all the way down to the Batchelor scale. The model would then require information on the local flow stretching only, i.e., a resolution of the Kolmogorov length is sufficient. From this perspective, the diffuselet approach can be considered as an efficient Lagrangian subgrid-scale parametrization of the chaotic mixing in the viscous-convective scale range. Subgrid scale models for Sc $\gg 1$ are rare. We mention the nonlinear large-eddy simulation approach of Burton and Dahm [44] that applies the multifractal model of kinetic energy and scalar dissipation rates to close the subgrid scales stresses. This closure model was extended in Ref. [45] to high Schmidt numbers.

## B. Cloudy air filament mixing in three-dimensional homogeneous isotropic flow

### 1. Cloud water droplet growth model

We extend the passive scalar mixing study to the case with phase changes in the following. As mentioned in the introduction, the interaction between cloud water droplets and small-scale turbulence plays a crucial role in the dynamics and lifetime of the (warm) turbulent cloud as a whole. To this end, we extend our physical model by cloud water droplets. Droplets are initially seeded randomly inside the supersaturated region, where $s > 0$. Recall that the liquid water mixing
<figure><img src="image_6.png" /><figcaption>FIG. 4. Comparison of the diffuselet method with DNS for a turbulent flow at Sc = 70 (Run T2). The central panel presents the evolution of the mean-squared scalar $\langle c^2 \rangle$ over time, with selected time points (I–IV) marked by red dots. The top row shows contour plots of the scalar field (logarithm of concentration) at these times, illustrating the effect of high diffusivity on scalar structures. The bottom row compares the PDF of scalar concentration c for both methods, showing their agreement across different stages of scalar mixing.</figcaption></figure>

ratio $q_l$ is not included as an Eulerian field, but as an ensemble of $N_l$ point-like Lagrangian tracers that carry a radius coordinate $r_j(t)$ with $j = 1, \dots, N_l$. The droplet evaporation and growth add a sink and source term to the advection-diffusion equation (3), respectively [9,11,30,36]. This leads to

$$ \frac{\partial s}{\partial t} + (\mathbf{u} \cdot \nabla)s = D_s \nabla^2 s - A \frac{4\pi \rho_l K}{V_{\text{cell}}} \sum_{m=1}^{\hat{N}} r_m(t)s[\mathbf{X}_m(t), t], \quad (31) $$

where A (measured in m³/kg) is a parameter which depends on thermodynamic and physical parameters, such as saturation water vapor pressure and latent heat [31]. Furthermore, $\rho_l$ is the liquid water density, K the evaporation constant (measured in m²/s), $V_{\text{cell}} = (\Delta x)^3$ the mesh cell volume, $r_j$ the radius of the jth cloud water droplet with spatial coordinates $\mathbf{X}_j$. The summation is over $\hat{N}$ droplets inside $V_{\text{cell}}$ at position $\mathbf{x}$. Furthermore, $D_s \approx D$. The constant K carries the physical dimension of a diffusion constant; it is a function of the latent heat, the saturation pressure, the environmental temperature in the vicinity of the droplet, see Refs. [5] or [27] for a detailed derivation. The radius changes of the jth droplet depend on the fluctuations of supersaturation in the vicinity of the droplet. This leads to [5]

$$ \frac{dr_j^2}{dt} = 2Ks[\mathbf{X}_j(t), t] \quad \text{for} \quad j = 1, \dots, N_t. \quad (32) $$

The equation describes a diffusional growth in the Lagrangian frame. The initial size of the droplets in the filament for the present case is uniformly set to $r_j(0) = r_0 = 20$ µm (except run D8 with $r_0 = 10$ µm). The droplets are thus too large to incorporate curvature and solutal effects in (32) on
the one hand. On the other hand, they are still too small for collision-coalescence. For simplicity, they will thus be considered as Lagrangian tracers. The droplet radius evolves due to the evaporation process and serves as a source of water vapor. This source term is the condensation date $C_d$, the last term on the right-hand size of Eq. (31), which is given by

$$C_d(\mathbf{x}, t) = \frac{4\pi \rho_l K}{V_{\text{cell}}} \sum_{m=1}^{\hat{N}} r_m(t) s[\mathbf{X}_m(t), t]. \quad (33)$$

The condensation rate enters as an Eulerian field in (31) and is obtained from the Lagrangian values along each discrete droplet trajectory $\mathbf{X}_m(t)$ by reshuffling to the eight neighboring sites $\mathbf{x}$ of the uniform computational grid, loosely speaking by an inverse interpolation. In this way, we couple Eulerian and Lagrangian model parts. The Schmidt number is again $Sc \approx 0.7$. Table II summarizes all cases that are studied in this subsection. Characteristic cloud water droplet number densities are $n \approx 10^2 \text{ cm}^{-3}$; we will, however, enhance this density to artificially high values to investigate the impact of a dense spray-like configuration on the mixing process. Following Ref. [26], we will also enhance the growth diffusion constant $K$ by an order of magnitude to quantify the impact of a stronger droplet response to subsaturation. Together with size effects, this results in eight DNS runs to compare DNS and diffuselet approaches.

The monodisperse droplets are seeded randomly in the center of the cloudy air filament for time $t = 0$, at places where the vapor mixing ratio $q_v > q_{vs}$, i.e., where the supersaturation $s > 0$. Note also that the supersaturation field $s(\mathbf{x}, t)$ is rescaled in the same way as $c(\mathbf{x}, t)$ by Eq. (2) as for the droplet-free mixing case.

## 2. Comparison to diffuselet model for different cloud water droplet number densities

The numerical simulations demonstrate that the additional source term due to cloud water droplets influences the statistical properties of the scalar field, as shown in Fig. 5. For a number density $n_d = 100 \text{ cm}^{-3}$ in run D1, which is typical for a warm cloud, the difference to the previous droplet-free case in run T2 remains small for the whole time evolution that is considered. Deviations from the droplet-free case become larger for higher number densities in runs D2 and D3. More droplets are subject to evaporation, which enhances the magnitude of the condensation rate; $C_d < 0$ holds for evaporation, which is a source of vapor content. Furthermore, the droplet size distribution is broadest for the highest number density, which is visible in the bottom row of Fig. 5. More droplets remain in the supersaturated fractions of the filaments such that their evaporation is delayed. The top row of the figure shows vertical cuts through the simulation volume displaying contours of the scalar field together with droplets in this cross section plane which are colored with respect to their individual size. Three time instants of $t = 0.2$, 1, and 3 s are shown. Since the droplets are Lagrangian tracers, they do not decouple from the scalar filaments in this study. This is different from Ref. [46] where gravitational settling and particle inertia were included. Blue colored satellite droplets are thus observable at the subsaturated outer scalar filaments that are subject to the strongest entrainment of clear subsaturated air. To conclude, the comparison with the diffuselet method is basically as for the droplet-free run T1 for the realistic value of the cloud droplet number density. Deviations increase with increasing $n_d$.

To verify the influence of initial droplet size on the mixing process, an additional simulation with a reduced initial droplet radius of $r_0 = 10 \text{ µm}$ was performed (run D8 in Table II) for $n_d = 100 \text{ cm}^{-3}$. The mean square of the scalar and the PDF showed negligible differences between the two cases, indicating that scalar mixing remains largely unaffected by droplet size at this range. However, the droplet radius distributions increasingly deviate from each other with progressing time. A consequence of the droplet evaporation law (32) is that the initially smaller droplets shrink faster and thus lead to a broader size distribution at later times. This result is not shown but is consistent with the trends discussed above.
<figure><img src="image_7.png" /></figure>

FIG. 5. Impact of additional cloud water droplets on the turbulent mixing process. (first row) Vertical cuts through simulation domain snapshots at different time instants. Droplets and scalar contours are shown for run D1 from Table II. The rescaled scalar supersaturation field of the filament and the seeded droplets are shown. The scalar is plotted in logarithmic contour levels. The droplets are colored with respect to their radius from red (largest) to blue (smallest). (second row) Mean squared scalar versus time $\langle c^2(t) \rangle$. (third row) Scalar PDFs at the three time instants from the top row. The dashed black line represents the critical value $c_c = 1/(1+s_{max})$, which distinguishes regions where droplets grow from those where they shrink. (fourth row) Corresponding droplet size distributions. The number densities and the diffuselet case are indicated by the legend in the second row. These are runs D1, D2, and D3 from Table II.

### 3. Comparison with diffuselet model for enhanced droplet evaporation

The main relative humidity (or saturation) for the specified configuration, even in a well-mixed equilibrium, is considerably less than 100%, thereby enabling droplets to evaporate further. Nevertheless, the rate of evaporation for realistic physical parameters does not permit the evaporation of more than 7% of the initial size for the time interval that we considered. The evaporation constant $K$ (which has the physical dimension of a diffusion constant) determines the rate at which droplets exchange mass with the surrounding air, see again Eq. (32).

To further enhance evaporation of cloud water droplets, $K$ can be increased beyond realistic values which was done for run D7 from Table II. An increase of $K$ by one order of magnitude is sufficient to produce notable changes in the evolution of droplet size distribution, as shown in Ref. [26]. Figure 6 summarizes the comparison of runs D2 and D7. Starting with $t = 2$ s, the number of droplets begins to decrease due to the complete droplet evaporation in D7. It can be observed in the bottom row of the figure that the PDF of the droplet size distribution evolves from an initial delta-type distribution to a broader one, similar to what was shown in Ref. [26]. Droplets outside the core of the cloudy-air slab shrink more rapidly. It is also seen that the deviations to the diffuselet model in D7 increase. The effect of the additional condensation rate term is again clearly visible
<figure><img src="image_8.png" /><figcaption>FIG. 6. Comparison of results between realistic (run D1) and artificially enhanced (run D7) constant K in the droplet growth equation (32). (top row) Mean-squared scalar versus time, $\langle c^2(t) \rangle$. (second row) Scalar PDFs at the three time instants from the top row. The dashed black line represents again the critical value $c_c = 1/(1+s_{max})$, which distinguishes regions where droplets grow from those where they shrink. (bottom row) Corresponding droplet size distributions. The values of K and the diffuselet case are indicated by the legend in the top row.</figcaption></figure>

in the scalar PDFs by a saturated tail at the smallest scalar concentrations, similar to the artificially enhanced droplet number densities in the third row of Fig. 5.

### 4. Effect of bigger simulation domain size on comparison to diffuselet model

From the previous discussions, it becomes evident that, for the given scalar diffusivity, the ratio of filament to simulation volume (which is effectively $\ell_0/L$) is a factor that has to be included in the comparison of DNS and diffuselet methods. This ratio also determines the magnitude of the scalar mean value to which the mean-squared scalar relaxes in the final state of the turbulent mixing process. Runs D4 to D6 thus repeat runs D1 to D3, but in a simulation box $V = L^3$ that is twice as big in each space direction. Figure 7 confirms an improvement of the agreement of DNS results to those of the diffuselet method for this series when comparing to Fig. 5.

### 5. Time evolution for individual cloud water droplets

The evaporation and size changes of the droplets can be approximated by the well-known $d^2$ law, which suggests that the square of the droplet diameter changes linearly with time. This result can be obtained from Eq. (32) as $d^2(t) \approx d_0^2 - 8Kt \langle s \rangle_V$ when the supersaturation $s[X_i(t)] \rightarrow \langle s \rangle_V$, i.e., a constant. Here $d$ is the droplet diameter, $d_0$ is the initial value [$r = d/2$ in Eq. (32)]. The lifetime of a single droplet in such a uniform field is

$$t_1 = \frac{d_0^2}{8K \langle s \rangle_V}, \qquad (34)$$

where $\langle s \rangle_V$ is the volume mean of the supersaturation field. This law assumes an idealized system, where evaporation is specified only by water vapor diffusion in a homogeneous background,
<figure><img src="image_9.png" /><figcaption>FIG. 7. Same as in Fig. 5, but for runs D4, D5, and D6 of Table II with box size $L = 0.512$ m.</figcaption></figure>

given here by the volume average of the supersaturation field, $\langle s \rangle_V$; fluctuations of humidity are disregarded. The $d^2$ law is likely to be violated when a droplet experiences large fluctuations of temperature and water vapor content in the field in which it evaporates, like those resulting from turbulent mixing in clouds, as first anticipated by Srivastava [47].

Indeed, selected tracks of individual cloud water droplets $r(t)$ from the DNS data display a substantial evaporation delay as compared to the $d^2$ expectation, as can be seen in Fig. 8. The reason is easy to understand given our starting strip topology: close-packed droplets do not evaporate when surrounded by their own saturating vapor but do so promptly (on a timescale given by the $d^2$ law).

<figure><img src="image_10.png" /><figcaption>FIG. 8. Evaporation process along individual Lagrangian cloud water droplet tracks. The figure shows 1310 individual particle tracks colored with respect to supersaturation (see legend). (a) Run D1. (b) Run D7. See Table II. The black dashed line displays the $d^2$ law for the evaporation in a homogeneous supersaturation field. The red dashed line represents the mixing time $t_s$.</figcaption></figure>
when their interstitial vapor concentration falls below saturation. The lifetime of a whole droplet strip of width $l_0$ is thus given by the time $t_v$ it takes for the vapor to evacuate the droplets assembly, and thus coincides with the mixing time $t_s$ of the vapor constructed on the initial size of the strip [33,48] which reads

$$ t_v = \frac{1}{\gamma_{\text{turb}}} \ln(1 + \phi^{-1}), \\ \text{with } \phi = \frac{\rho_s - \rho_{\infty}}{\varphi \rho_l + (1 - \varphi) \rho_s} \frac{2}{\sqrt{\text{Pe}}}, \qquad (35) $$

$$ \text{and } \text{Pe} = \frac{\gamma_{\text{turb}} l_0^2}{D}, $$

where $\rho_s$ is the saturated vapor density, $\rho_l$ the liquid water density, $\rho_\infty$ the vapor density outside of the spray, and $\varphi$ the liquid water volume fraction, $\text{Pe} = 4208$ is the Péclet number. Recall that the stretching rate is found $\gamma_{\text{turb}} = 10.6 \text{ s}^{-1}$, see directly below Eq. (28). The volume fraction $\varphi$ is defined as the ratio of liquid water volume $V_l$ to vapor volume $V_v$:

$$ \varphi = \frac{V_l}{V_v} = \frac{\sum_{i=1}^{N_i} 4\pi r_i^3 / 3}{V_v}, \qquad (36) $$

where $N_i$ is the total number of droplets. For case D1 of the present simulations, we have $\varphi \approx 10^{-6}$, leading to $\varphi \approx 4 \times 10^{-4}$, giving a spray lifetime calculated from Eq. (35) of is $t_v = 0.73 \text{ s}$ in very fair agreement with the observations reported in Fig. 8, showing that it takes roughly somewhat less than a second for most of the droplets to start to decrease rapidly. Evaporation of dense sprays is a mixing problem.

## 6. Evolution of individual cloud water droplets from the diffuselet model

In Sec. II B, we explained that the maximum concentration of the diffuselet $i$ is connected to the dimensionless time $\tau_i(t)$, see Eq. (12), by the following relation:

$$ \tilde{c}(X_i, t) \approx \frac{1}{\sqrt{\tau_i}}. \qquad (37) $$

This suggests investigating how this direct evaluation of the supersaturation at the droplet position $X_i(t)$ compares to the one which follows from the numerically simulated scalar field. We conduct this comparison for simulation run D1 at a Schmidt number of $Sc = 0.7$ in the following. Thus, we determine $\tilde{c}(X_i, t)$ from $\tilde{c}(X_i, t)$ and insert this into Eq. (32) for each droplet. We note at this point that the original diffuselet model does not include the phase changes which are studied here. Thus, we can expect deviations in the course of time.

Figure 9 shows the comparison between both methods. As can be seen for all three examples, up to the mixing time $t_s$, both methods agree very well, even for individual droplets. However, differences become visible at later times which either lead to a slightly faster or slightly slower decrease of the droplet radius with time. Similarly, the scalar concentration starts to differ, most probably caused by additional deformations of the diffuselet sheet in the two space directions normal to the cross-sheet coordinate. These effects are not fully captured by the model. When these data for individual droplets are averaged over the full droplet ensemble, we obtain again the results that have been reported already in Fig. 5. They are replotted here in the last column of Fig. 9.

Finally, Fig. 10 shows the PDFs of the droplet radii at different times, comparing both, diffuselet method and DNS. At initial times for $t = 0.25 \text{ s}$ and $0.5 \text{ s}$, the PDFs agree fairly well. The distributions are still localized near initial radius $r_0 = 20 \mu\text{m}$. As time progresses, for $t = 0.75 \text{ s}$ and $1 \text{ s}$, differences between both methods start to occur, particularly in the left tail of the distribution for the droplets that evaporate fastest. For the last two instants at $t = 2 \text{ s}$ and $3 \text{ s}$, a pronounced right tail
<figure><img src="image_11.png" /><figcaption>FIG. 9. Individual droplet evaporation histories directly obtained from the diffuselet model by Eq. (37) in comparison with the simulation results. This comparison is done for individual droplets. Three example cases are shown together with the average, which is taken over the whole droplet ensemble. (upper row) Time evolution of droplet radius. (bottom row) Evolution of scalar concentration. Quantities with a tilde correspond to results from the diffuselet model. Data are for run D1.</figcaption></figure>

is observed for the diffuselet method. The distribution obtained from the DNS model is narrower on both sides of the tails, which is a result of the saturation of the mean square to the well-mixed limit.

<figure><img src="image_12.png" /><figcaption>FIG. 10. Comparison of probability density function of the cloud water droplet radius obtained by diffuselet method (blue lines) and DNS (orange line) for six different time instants, which are indicated in each panel. Data are for run D1.</figcaption></figure>
## IV. CONCLUSION

The primary focus of this work was a comparison of the turbulent mixing of a scalar filament in a three-dimensional incompressible turbulent Navier-Stokes flow by two methods: the Lagrangian kinematic diffuselet model [25] and direct numerical simulations of scalar turbulence. The configuration of interest is a sheet-like thin scalar filament at the beginning of the mixing process. The turbulence is in a statistically stationary regime. This filament (or sheet) is covered by little area elements, the diffuselets; their center is attached to a Lagrangian tracer which is advected in the turbulent flow. Along the individual diffuselet element trajectories, we then collect the individual stretching history, which allows us to predict the time evolution of the mean squared scalar and the probability density function PDF of the scalar concentration $c$ (and correspondingly $s$) of the mixing process in the control volume $V = L^3$. The results of the diffuselet framework and the DNS agree as long as the limitations of approximating the process in a box with periodic boundary conditions in the numerical simulations become important. The scalar mixing process is compared for two Schmidt numbers, $Sc = 0.7$ and $70$. We consider single-phase and two-phase turbulent mixing processes for this study. The latter case is a simple model for the entrainment and subsequent turbulent mixing of clear air with cloudy air (which is represented as the simple sheet-like filament) in a small control volume $V$ at the edge of a turbulent atmospheric cloud. In the two-phase case, cloud water droplets are seeded to the filament subject to slow shrinkage and eventual evaporation. The droplets are treated here simply as Lagrangian tracers; gravitational settling and droplet inertia are neglected in the present study. For a discussion of this simplification, we refer to Ref. [31].

The major findings can be summarized as follows. The results of the diffuselet method and the DNS for the single-phase case agree well up to the mixing timescale $t_s$. Afterwards, deviations start to grow systematically. They are caused by the relaxation of the concentration fluctuations to the finite volumetric mean value of the passive scalar concentration field—the well-mixed final state, which is not included in the diffuselet method. Here, the finite mean is a consequence of the conservation of concentration in the triply periodic control volume. The agreement improves when the Schmidt number is enhanced and when the initial filament thickness is smaller (not shown).

The comparison of the diffuselet model to DNS for the cloudy air filament showed a level of agreement similar to that of a single-phase scalar field. This holds particularly for values of droplet number density and diffusion constant $K$ in the droplet growth equation that correspond to realistic values inside the cloud. The scalar field in this case is the supersaturation field $s$, which determines the local growth conditions of the cloud water droplets. For enhanced number density and droplet response to subsaturation, the deviations become larger since the magnitude of the source (or sink) term – the condensation rate $C_d$ – in the scalar transport Eq. (31) grows. Again, we can conclude that a smaller initial cloudy air filament thickness relative to the box size improves the agreement to the diffuselet model. Turbulent fluctuations of the supersaturation field due to entrainment cause also deviations from the $d^2$ law for droplet evaporation in sprays. Our results thus agree with a recent mixing experiment in a laboratory flow [33]. To summarize this part, we demonstrated a first extension of the diffuselet method to an evaporating phase of a collection of closely packed droplets and anticipated the overall evaporation time of a strip of a relatively dense spray, see Eqs. (35).

To conclude, we have thus demonstrated that the diffuselet method is applicable even for the case of low Schmidt numbers $Sc \approx 1$ including phase changes. Our work shows clearly that the diffuselet method has its limitations, in particular when longer evolution times $t > t_s$ have to be considered. By construction, the agreement with DNS is best for high-Schmidt-number mixing processes. Here, the model can be even used as a Lagrangian subgrid scale approach that models the mixing in the viscous-convective scale range. A further extension of the diffuselet method would consist in a direct incorporation of the phase changes, i.e., the integration of loss of liquid water for evaporation or the loss of vapor content for condensation into the kinematic Lagrangian framework. These processes would then introduce an active scalar component into the model with latent heat release and evaporative cooling. These efforts can be expected to be reported elsewhere.
ACKNOWLEDGMENTS

The work of V.P. and S.S. was supported by ITN CoPerMix. The project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 956457. Supercomputing time has been provided at the University Computer Center (UniRZ) of the TU Ilmenau. The authors also gratefully acknowledge the Gauss Centre for Supercomputing e.V. [49] for funding this project by providing computing time on the GCS Supercomputer Super MUC-NG at Leibniz Supercomputing Center [50].

DATA AVAILABILITY

The data supporting this publication are available on Zenodo at [51].

APPENDIX A: PASSIVE SCALAR MIXING IN STEADY ANALYTICAL FLOW

Here, we validate the diffuselet method for a simple analytical steady velocity field in a three-dimensional cubic domain with periodic boundary conditions:

$$u_x(\mathbf{x}, t) = A_1 \cos(\omega_1 y), \quad u_y(\mathbf{x}, t) = A_2 \sin(\omega_2 x), \quad u_z(\mathbf{x}, t) = A_3, \qquad (A1)$$

where $A_1$, $A_2$, $A_3$, as well as $\omega_1$, $\omega_2$ are constants. For simplicity, we consider $A_1 = A_2 = A = 0.1 \text{ m/s}$, $A_3 = 0$, together with $\omega_1 = \omega_2 = \omega = 2\pi/L$. For convenience, we also provide scales and velocities often in physical units, such that the reader gets a better impression about the dimensions of the considered mixing volumes. The diffusivity coefficient $D$ was chosen to be five times smaller than that of water vapor, so $D = 5 \times 10^{-6} \text{ m/s}^2$. As mentioned above, the initial scalar field takes the form of a sheet with a transversal Gaussian profile and an initial thickness $\ell_0 = 0.005 \text{ m}$. This is about 2% of the box size $L = 0.512 \text{ m}$. The sheet is divided into $512 \times 512$ surface elements with surface $\delta A_0$ and sheet thickness $\ell_0$; these are the diffuselets (see again Fig. 1). This implies an advection of 262 144 Lagrangian tracers simultaneously with the analytical flow.

The comparison of the results for mean square and PDF of the scalar field $c$, obtained by the diffuselet method and the DNS, is presented in Fig. 11. The mean squared scalar in the DNS, which is calculated as a volume average at time $t$, coincides perfectly with the one from the diffuselet. The mean square starts to decay after $t_s \approx 1 \text{ s}$. At this mixing time, scalar diffusion is enhanced because the thickness of the sheet becomes smaller than the Batchelor scale. Here, the flow is two-dimensional (2D) and steady, which leads to a linear growth in time of the sheet's surface as in a pure shear flow. The mixing time is thus given by [19]

$$t_s = \frac{1}{\gamma}(3Pe)^{1/3}, \qquad (A2)$$

where the Péclet number $Pe = \gamma \ell_0^2/D$ is based on the initial thickness $\ell_0$ of the scalar sheet. Here, the shear $\gamma = 2\pi A_1/L \approx 1.2 \text{ s}^{-1}$, see again Eqs. (A1), which results in a mixing time $t_s \approx 2 \text{ s}$, in good agreement with the numerical results.

Figure 11 displays the concentration PDFs, calculated using the DNS and using the diffuselet model with Eq. (23). There is an excellent agreement between the two methods up to times $t \lesssim 10 \text{ s}$. Initially, the PDF exhibits a peak at $c = 1$, corresponding to the maximum concentration inside the sheet of scalar. After the mixing time, this peak disappears and the PDFs become monotonic decreasing functions of $c$. After $t_{II} = 5 \text{ s}$, the maximal concentration starts to decrease as a function of time, meaning that all parts of the sheet have started to diffuse. The maximum concentration is equal to $c_{max} = 0.5$ at $t_{III} = 5 \text{ s}$. At a later time ($t_{IV} = 10 \text{ s}$), the maximal concentration predicted by the diffuselet model is slightly larger than in the DNS. Indeed, this maximal concentration corresponds to a weakly stretched scalar sheet for which the diffuselet method is inappropriate, since it models the diffusion as a one-dimensional (1D) transverse diffusive process, whereas these points diffuse in two directions in the DNS.
<figure><img src="image_13.png" /><figcaption>FIG. 11. Comparison of the diffuselet method with DNS for scalar transport in a steady analytical flow given by Eq. (A1). The central panel shows the time evolution of the mean squared scalar ($c^2(t)$) in a doubly logarithmic plot. The red markers indicate selected time instants I–IV. The top row shows two-dimensional contour plots of the scalar field (logarithm of concentration) at these times. The bottom row of panels provides the corresponding PDFs of scalar concentration $c$ for both methods at the four times.</figcaption></figure>

## APPENDIX B: TESTS OF NUMERICAL RESOLUTION

The diffuselet method might be sensitive to the number of diffuselets that cover the filament sheet. In this appendix, we present additional runs (see Table III) to investigate the performance for differently fine diffuselet decompositions. To this end, we compare three DNS runs R1, R2, and R3, with different spectral resolution and correspondingly different diffuselet numbers, all advecting the scalar in the same turbulent flow. We can conclude that the different decompositions lead to the same statistics, as seen in Fig. 12, for both the mean-squared scalar and scalar distributions.

<table><caption>TABLE III. Parameters of three DNS runs to compare differently fine decompositions of the initial filament into diffuselets together with differently fine numerical grids. The energy injection rate is $\epsilon_{inj} \approx 0.004 \text{ m}^2/\text{s}^3$, which results in Kolmogorov microscale $\eta_K = 1 \text{ mm}$.</caption><thead><tr><th>Run</th><th>L [m]</th><th>N<sup>3</sup></th><th>Δx/η<sub>K</sub></th><th>δA<sub>0</sub></th><th>ℓ<sub>0</sub> [m]</th><th>Re</th><th>Re<sub>λ</sub></th><th>Sc</th></tr></thead><tbody><tr><td>R1</td><td>0.256</td><td>256<sup>3</sup></td><td>1.0</td><td>Δx × Δx</td><td>0.01</td><td>1634</td><td>52</td><td>0.7</td></tr><tr><td>R2</td><td>0.256</td><td>512<sup>3</sup></td><td>0.5</td><td>2Δx × 2Δx</td><td>0.01</td><td>1634</td><td>52</td><td>0.7</td></tr><tr><td>R3</td><td>0.256</td><td>1024<sup>3</sup></td><td>0.25</td><td>4Δx × 4Δx</td><td>0.01</td><td>1634</td><td>52</td><td>0.7</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>FIG. 12. Comparing results of the diffuselet method with differently coarse decompositions of the initial filament sheet. (top row) Mean-squared scalar versus time. (bottom row) PDFs of the scalar concentration for four different instants. The initial diffuselet surface for runs R1 to R3 is indicated in Table III.</figcaption></figure>

[1] R. A. Shaw, Particle-turbulence interactions in atmospheric clouds, Annu. Rev. Fluid Mech. 35, 183 (2003).

[2] E. Bodenschatz, S. P. Malinowski, R. A. Shaw, and F. Stratmann, Can we understand clouds without turbulence? Science 327, 970 (2010).

[3] B. Stevens and S. Bony, Water in the atmosphere, Phys. Today 66(6), 29 (2013).

[4] R. Vogel, A. L. Albright, J. Vial, G. George, B. Stevens, and S. Bony, Strong cloud-circulation coupling explains weak trade cumulus feedback, Nature (London) 612, 696 (2022).

[5] R. R. Rogers and M. K. Yau, A Short Course in Cloud Physics (Butterworth-Heinemann, Oxford, UK, 1989).

[6] M. B. Baker, R. E. Breidenthal, T. W. Choularton, and J. Latham, The effects of turbulent mixing in clouds, J. Atmos. Sci. 41, 299 (1984).

[7] M. Andrejczuk, W. W. Grabowski, S. P. Malinowski, and P. K. Smolarkiewicz, Numerical simulation of cloud-clear air interfacial mixing: Effects on cloud microphysics, J. Atmos. Sci. 63, 3204 (2006).

[8] A. Celani, A. Mazzino, and M. Tizzi, The equivalent size of cloud condensation nuclei, New J. Phys. 10, 075021 (2008).

[9] G. Sardina, F. Picano, L. Brandt, and R. Caballero, Continuous growth of droplet size variance due to condensation in turbulent clouds, Phys. Rev. Lett. 115, 184501 (2015).

[10] B. Kumar, P. Götzfried, N. Suresh, J. Schumacher, and R. A. Shaw, Scale dependence of cloud microphysical response to turbulent entrainment and mixing, J. Adv. Model. Earth Syst. 10, 2777 (2018).

[11] J. Fries, G. Sardina, G. Svensson, and B. Mehlig, Key parameters for droplet evaporation and mixing at the cloud edge, Q. J. R. Meteorol. Soc. 147, 2160 (2021).

[12] B. Kumar, R. Ranjan, M. K. Yau, S. Bera, and S. A. Rao, Impact of high and low vorticity turbulence on cloud environment mixing and cloud microphysics processes, Atmos. Chem. Phys. 21, 12317 (2021).

[13] W. W. Grabowski, L. Thomas, and B. Kumar, Impact of cloud-base turbulence on CCN activation: single-size CCN, J. Atmos. Sci. 79, 551 (2022).

[14] J. Fries, G. Sardina, G. Svensson, A. Pumir, and B. Mehlig, Lagrangian supersaturation fluctuations at the cloud edge, Phys. Rev. Lett. 131, 254201 (2023).

[15] L. Magaritz-Ronen, M. Pinsky, and A. Khain, Drizzle formation in stratocumulus clouds: Effects of turbulent mixing, Atmos. Chem. Phys. 16, 1849 (2016).
[16] P. E. Dimotakis, Turbulent mixing, Annu. Rev. Fluid Mech. 37, 329 (2005).

[17] K. R. Sreenivasan and J. Schumacher, Lagrangian views on turbulent mixing of passive scalars, Phil. Trans. R. Soc. A 368, 1561 (2010).

[18] K. R. Sreenivasan, Turbulent mixing: A perspective, Proc. Natl. Acad. Sci. USA 116, 18175 (2019).

[19] E. Villermaux, Mixing versus stirring, Annu. Rev. Fluid Mech. 51, 245 (2019).

[20] G. K. Batchelor, Small-scale variation of convected quantities like temperature in turbulent fluid Part 1. General discussion and the case of small conductivity, J. Fluid Mech. 5, 113 (1959).

[21] R. H. Kraichnan, Small-scale structure of a scalar field convected by turbulence, Phys. Fluids 11, 945 (1968).

[22] E. Villermaux and J. Duplat, Mixing as an aggregation process, Phys. Rev. Lett. 91, 184501 (2003).

[23] P. Meunier and E. Villermaux, How vortices mix, J. Fluid Mech. 476, 213 (2003).

[24] P. Meunier and E. Villermaux, The diffusive strip method for scalar mixing in two dimensions, J. Fluid Mech. 662, 134 (2010).

[25] P. Meunier and E. Villermaux, The diffuselet concept for scalar mixing, J. Fluid Mech. 951, A33 (2022).

[26] B. Kumar, F. Janetzko, J. Schumacher, and R. A. Shaw, Extreme responses of a coupled scalar-particle system during turbulent mixing, New J. Phys. 14, 115020 (2012).

[27] B. Kumar, J. Schumacher, and R. A. Shaw, Cloud microphysical effects of turbulent mixing and entrainment, Theor. Comput. Fluid Dyn. 27, 361 (2013).

[28] B. Kumar, J. Schumacher, and R. A. Shaw, Lagrangian mixing dynamics at the cloudy-clear air interface, J. Atmos. Sci. 71, 2564 (2014).

[29] P. Götzfried, B. Kumar, R. A. Shaw, and J. Schumacher, Droplet dynamics and fine-scale structure in a shearless turbulent mixing layer with phase changes, J. Fluid Mech. 814, 452 (2017).

[30] A. Celani, G. Falkovich, A. Mazzino, and A. Seminara, Droplet condensation in turbulent flows, Europhys. Lett. 70, 775 (2005).

[31] V. Pushenko and J. Schumacher, Connecting finite-time Lyapunov exponents with supersaturation and droplet dynamics in a turbulent bulk flow, Phys. Rev. E 109, 045101 (2024).

[32] J. Schumacher, Sub-Kolmogorov-scale fluctuations in fluid turbulence, Europhys. Lett. 80, 54001 (2007).

[33] A. de Rivas and E. Villermaux, Dense spray evaporation as a mixing process, Phys. Rev. Fluids 1, 014201 (2016).

[34] A. N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Dokl. Akad. Nauk SSSR 30, 299 (1941).

[35] A. M. Obukhov, Structure of the temperature field in a turbulent flow, Izv. Akad. Nauk SSSR, Ser. Geogr. and Geofiz. 13, 281 (1949).

[36] S. Corrsin, On the spectrum of isotropic temperature fluctuations in an isotropic turbulence, J. Appl. Phys. 22, 469 (1951).

[37] J. M. Ottino, Mixing, chaotic advection, and turbulence, Annu. Rev. Fluid Mech. 22, 207 (1990).

[38] D. Pekurovsky, P3DFFT: A framework for parallel computations of Fourier transforms in three dimensions, SIAM J. Sci. Comput. 34, C192 (2012).

[39] P. K. Yeung and S. B. Pope, An algorithm for tracking fluid particles in numerical simulations of homogeneous turbulence, J. Comput. Phys. 79, 373 (1988).

[40] W. E. Ranz, Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows, AIChE J. 25, 41 (1979).

[41] W. T. Ashurst, A. R. Kerstein, R. M. Kerr, and C. H. Gibson, Alignment of vorticity and scalar gradient with strain rate in simulated Navier-Stokes turbulence, Phys. Fluids 30, 2343 (1987).

[42] S. S. Girimaji and S. B. Pope, Material-element deformation in isotropic turbulence, J. Fluid Mech. 220, 427 (1990).

[43] G. Zinchenko, V. Pushenko, and J. Schumacher, Local precursors to anomalous dissipation in Navier-Stokes turbulence: Burgers vortex-type models and simulation analysis, Phys. Rev. Fluids 9, 114608 (2024).

[44] G. C. Burton and W. J. A. Dahm, Multifractal subgrid-scale modeling for large-eddy simulation. II. Backscatter limiting and a *posteriori* evaluation, Phys. Fluids 17, 075112 (2005).
[45] G. C. Burton, The nonlinear large-eddy simulation method applied to Sc ≈ 1 and Sc ≫ 1 passive-scalar mixing, *Phys. Fluids* **20**, 035103 (2008).

[46] P. Götzfried, M. S. Emran, E. Villermaux, and J. Schumacher, Comparison of Lagrangian and Eulerian frames of passive scalar turbulent mixing, *Phys. Rev. Fluids* **4**, 044607 (2019).

[47] R. C. Srivastava, Growth of cloud drops by condensation: A criticism of currently accepted theory and a new approach, *J. Atmos. Sci.* **46**, 869 (1989).

[48] E. Villermaux, A. Moutte, M. Amielh, and P. Meunier, Fine structure of the vapor field in evaporating dense sprays, *Phys. Rev. Fluids* **2**, 074501 (2017).

[49] https://www.gauss-centre.eu.

[50] https://www.lrz.de.

[51] 10.5281/zenodo.15853362.